contrib/libgmp/mpf/get_str.c

   1 /* mpf_get_str (digit_ptr, exp, base, n_digits, a) -- Convert the floating
   2   point number A to a base BASE number and store N_DIGITS raw digits at
   3   DIGIT_PTR, and the base BASE exponent in the word pointed to by EXP.  For
   4   example, the number 3.1416 would be returned as "31416" in DIGIT_PTR and
   5   1 in EXP.
   6
   7 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
   8
   9 This file is part of the GNU MP Library.
  10
  11 The GNU MP Library is free software; you can redistribute it and/or modify
  12 it under the terms of the GNU Library General Public License as published by
  13 the Free Software Foundation; either version 2 of the License, or (at your
  14 option) any later version.
  15
  16 The GNU MP Library is distributed in the hope that it will be useful, but
  17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
  19 License for more details.
  20
  21 You should have received a copy of the GNU Library General Public License
  22 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  23 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  24 MA 02111-1307, USA. */
  25
  26 #include "gmp.h"
  27 #include "gmp-impl.h"
  28 #include "longlong.h"
  29
  30 /*
  31    New algorithm for converting fractions (951019):
  32    0. Call the fraction to convert F.
  33    1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e.,
  34       [exp * BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly].  Exp is
  35       the number of limbs between the limb point and the most significant
  36       non-zero limb.  Call this result n.
  37    2. Compute B^n.
  38    3. F*B^n will now be just below 1, which can be converted easily.  (Just
  39       multiply by B repeatedly, and see the digits fall out as integers.)
  40    We should interrupt the conversion process of F*B^n as soon as the number
  41    of digits requested have been generated.
  42
  43    New algorithm for converting integers (951019):
  44    0. Call the integer to convert I.
  45    1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e.,
  46       [exp BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly].  Exp is
  47       the number of limbs between the limb point and the least significant
  48       non-zero limb.  Call this result n.
  49    2. Compute B^n.
  50    3. I/B^n can be converted easily.  (Just divide by B repeatedly.  In GMP,
  51       this is best done by calling mpn_get_str.)
  52    Note that converting I/B^n could yield more digits than requested.  For
  53    efficiency, the variable n above should be set larger in such cases, to
  54    kill all undesired digits in the division in step 3.
  55 */
  56
  57 char *
  58 #if __STDC__
  59 mpf_get_str (char *digit_ptr, mp_exp_t *exp, int base, size_t n_digits, mpf_srcptr u)
  60 #else
  61 mpf_get_str (digit_ptr, exp, base, n_digits, u)
  62      char *digit_ptr;
  63      mp_exp_t *exp;
  64      int base;
  65      size_t n_digits;
  66      mpf_srcptr u;
  67 #endif
  68 {
  69   mp_size_t usize;
  70   mp_exp_t uexp;
  71   unsigned char *str;
  72   size_t str_size;
  73   char *num_to_text;
  74   long i;                       /* should be size_t */
  75   mp_ptr rp;
  76   mp_limb_t big_base;
  77   size_t digits_computed_so_far;
  78   int dig_per_u;
  79   mp_srcptr up;
  80   unsigned char *tstr;
  81   mp_exp_t exp_in_base;
  82   TMP_DECL (marker);
  83
  84   TMP_MARK (marker);
  85   usize = u->_mp_size;
  86   uexp = u->_mp_exp;
  87
  88   if (base >= 0)
  89     {
  90       if (base == 0)
  91         base = 10;
  92       num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz";
  93     }
  94   else
  95     {
  96       base = -base;
  97       num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
  98     }
  99
 100   /* Don't compute more digits than U can accurately represent.
 101      Also, if 0 digits were requested, give *exactly* as many digits
 102      as can be accurately represented.  */
 103   {
 104     size_t max_digits = (((u->_mp_prec - 1) * BITS_PER_MP_LIMB)
 105                          * __mp_bases[base].chars_per_bit_exactly);
 106     if (n_digits == 0 || n_digits > max_digits)
 107       n_digits = max_digits;
 108   }
 109
 110   if (digit_ptr == 0)
 111     {
 112       /* We didn't get a string from the user.  Allocate one (and return
 113          a pointer to it) with space for `-' and terminating null.  */
 114       digit_ptr = (char *) (*_mp_allocate_func) (n_digits + 2);
 115     }
 116
 117   if (usize == 0)
 118     {
 119       *exp = 0;
 120       *digit_ptr = 0;
 121       return digit_ptr;
 122     }
 123
 124   str = (unsigned char *) digit_ptr;
 125
 126   /* Allocate temporary digit space.  We can't put digits directly in the user
 127      area, since we almost always generate more digits than requested.  */
 128   tstr = (unsigned char *) TMP_ALLOC (n_digits + 3 * BITS_PER_MP_LIMB);
 129
 130   if (usize < 0)
 131     {
 132       *digit_ptr = '-';
 133       str++;
 134       usize = -usize;
 135     }
 136
 137   digits_computed_so_far = 0;
 138
 139   if (uexp > usize)
 140     {
 141       /* The number has just an integral part.  */
 142       mp_size_t rsize;
 143       mp_size_t exp_in_limbs;
 144       mp_size_t msize;
 145       mp_ptr tp, xp, mp;
 146       int cnt;
 147       mp_limb_t cy;
 148       mp_size_t start_str;
 149       mp_size_t n_limbs;
 150
 151       n_limbs = 2 + ((mp_size_t) (n_digits / __mp_bases[base].chars_per_bit_exactly)
 152                      / BITS_PER_MP_LIMB);
 153
 154       /* Compute n such that [u/B^n] contains (somewhat) more than n_digits
 155          digits.  (We compute less than that only if that is an exact number,
 156          i.e., exp is small enough.)  */
 157
 158       exp_in_limbs = uexp;
 159
 160       if (n_limbs >= exp_in_limbs)
 161         {
 162           /* The number is so small that we convert the entire number.  */
 163           exp_in_base = 0;
 164           rp = (mp_ptr) TMP_ALLOC (exp_in_limbs * BYTES_PER_MP_LIMB);
 165           MPN_ZERO (rp, exp_in_limbs - usize);
 166           MPN_COPY (rp + (exp_in_limbs - usize), u->_mp_d, usize);
 167           rsize = exp_in_limbs;
 168         }
 169       else
 170         {
 171           exp_in_limbs -= n_limbs;
 172           exp_in_base = (((exp_in_limbs * BITS_PER_MP_LIMB - 1))
 173                          * __mp_bases[base].chars_per_bit_exactly);
 174
 175           rsize = exp_in_limbs + 1;
 176           rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
 177           tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
 178
 179           rp[0] = base;
 180           rsize = 1;
 181
 182           count_leading_zeros (cnt, exp_in_base);
 183           for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
 184             {
 185               mpn_mul_n (tp, rp, rp, rsize);
 186               rsize = 2 * rsize;
 187               rsize -= tp[rsize - 1] == 0;
 188               xp = tp; tp = rp; rp = xp;
 189
 190               if (((exp_in_base >> i) & 1) != 0)
 191                 {
 192                   cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
 193                   rp[rsize] = cy;
 194                   rsize += cy != 0;
 195                 }
 196             }
 197
 198           mp = u->_mp_d;
 199           msize = usize;
 200
 201           {
 202             mp_ptr qp;
 203             mp_limb_t qflag;
 204             mp_size_t xtra;
 205             if (msize < rsize)
 206               {
 207                 mp_ptr tmp = (mp_ptr) TMP_ALLOC ((rsize+1)* BYTES_PER_MP_LIMB);
 208                 MPN_ZERO (tmp, rsize - msize);
 209                 MPN_COPY (tmp + rsize - msize, mp, msize);
 210                 mp = tmp;
 211                 msize = rsize;
 212               }
 213             else
 214               {
 215                 mp_ptr tmp = (mp_ptr) TMP_ALLOC ((msize+1)* BYTES_PER_MP_LIMB);
 216                 MPN_COPY (tmp, mp, msize);
 217                 mp = tmp;
 218               }
 219             count_leading_zeros (cnt, rp[rsize - 1]);
 220             cy = 0;
 221             if (cnt != 0)
 222               {
 223                 mpn_lshift (rp, rp, rsize, cnt);
 224                 cy = mpn_lshift (mp, mp, msize, cnt);
 225                 if (cy)
 226                   mp[msize++] = cy;
 227               }
 228
 229             {
 230               mp_size_t qsize = n_limbs + (cy != 0);
 231               qp = (mp_ptr) TMP_ALLOC ((qsize + 1) * BYTES_PER_MP_LIMB);
 232               xtra = qsize - (msize - rsize);
 233               qflag = mpn_divrem (qp, xtra, mp, msize, rp, rsize);
 234               qp[qsize] = qflag;
 235               rsize = qsize + qflag;
 236               rp = qp;
 237             }
 238           }
 239         }
 240
 241       str_size = mpn_get_str (tstr, base, rp, rsize);
 242
 243       if (str_size > n_digits + 3 * BITS_PER_MP_LIMB)
 244         abort ();
 245
 246       start_str = 0;
 247       while (tstr[start_str] == 0)
 248         start_str++;
 249
 250       for (i = start_str; i < str_size; i++)
 251         {
 252           tstr[digits_computed_so_far++] = tstr[i];
 253           if (digits_computed_so_far > n_digits)
 254             break;
 255         }
 256       exp_in_base = exp_in_base + str_size - start_str;
 257       goto finish_up;
 258     }
 259
 260   exp_in_base = 0;
 261
 262   if (uexp > 0)
 263     {
 264       /* The number has an integral part, convert that first.
 265          If there is a fractional part too, it will be handled later.  */
 266       mp_size_t start_str;
 267
 268       rp = (mp_ptr) TMP_ALLOC (uexp * BYTES_PER_MP_LIMB);
 269       up = u->_mp_d + usize - uexp;
 270       MPN_COPY (rp, up, uexp);
 271
 272       str_size = mpn_get_str (tstr, base, rp, uexp);
 273
 274       start_str = 0;
 275       while (tstr[start_str] == 0)
 276         start_str++;
 277
 278       for (i = start_str; i < str_size; i++)
 279         {
 280           tstr[digits_computed_so_far++] = tstr[i];
 281           if (digits_computed_so_far > n_digits)
 282             {
 283               exp_in_base = str_size - start_str;
 284               goto finish_up;
 285             }
 286         }
 287
 288       exp_in_base = str_size - start_str;
 289       /* Modify somewhat and fall out to convert fraction... */
 290       usize -= uexp;
 291       uexp = 0;
 292     }
 293
 294   if (usize <= 0)
 295     goto finish_up;
 296
 297   /* Convert the fraction.  */
 298   {
 299     mp_size_t rsize, msize;
 300     mp_ptr rp, tp, xp, mp;
 301     int cnt;
 302     mp_limb_t cy;
 303     mp_exp_t nexp;
 304
 305     big_base = __mp_bases[base].big_base;
 306     dig_per_u = __mp_bases[base].chars_per_limb;
 307
 308     /* Hack for correctly (although not efficiently) converting to bases that
 309        are powers of 2.  If we deem it important, we could handle powers of 2
 310        by shifting and masking (just like mpn_get_str).  */
 311     if (big_base < 10)          /* logarithm of base when power of two */
 312       {
 313         int logbase = big_base;
 314         if (dig_per_u * logbase == BITS_PER_MP_LIMB)
 315           dig_per_u--;
 316         big_base = (mp_limb_t) 1 << (dig_per_u * logbase);
 317         /* fall out to general code... */
 318       }
 319
 320 #if 0
 321     if (0 && uexp == 0)
 322       {
 323         rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
 324         up = u->_mp_d;
 325         MPN_COPY (rp, up, usize);
 326         rsize = usize;
 327         nexp = 0;
 328       }
 329     else
 330       {}
 331 #endif
 332     uexp = -uexp;
 333     if (u->_mp_d[usize - 1] == 0)
 334       cnt = 0;
 335     else
 336       count_leading_zeros (cnt, u->_mp_d[usize - 1]);
 337
 338     nexp = ((uexp * BITS_PER_MP_LIMB) + cnt)
 339       * __mp_bases[base].chars_per_bit_exactly;
 340
 341     if (nexp == 0)
 342       {
 343         rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
 344         up = u->_mp_d;
 345         MPN_COPY (rp, up, usize);
 346         rsize = usize;
 347       }
 348     else
 349       {
 350         rsize = uexp + 2;
 351         rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
 352         tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
 353
 354         rp[0] = base;
 355         rsize = 1;
 356
 357         count_leading_zeros (cnt, nexp);
 358         for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
 359           {
 360             mpn_mul_n (tp, rp, rp, rsize);
 361             rsize = 2 * rsize;
 362             rsize -= tp[rsize - 1] == 0;
 363             xp = tp; tp = rp; rp = xp;
 364
 365             if (((nexp >> i) & 1) != 0)
 366               {
 367                 cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
 368                 rp[rsize] = cy;
 369                 rsize += cy != 0;
 370               }
 371           }
 372
 373         /* Did our multiplier (base^nexp) cancel with uexp?  */
 374 #if 0
 375         if (uexp != rsize)
 376           {
 377             do
 378               {
 379                 cy = mpn_mul_1 (rp, rp, rsize, big_base);
 380                 nexp += dig_per_u;
 381               }
 382             while (cy == 0);
 383             rp[rsize++] = cy;
 384           }
 385 #endif
 386         mp = u->_mp_d;
 387         msize = usize;
 388
 389         tp = (mp_ptr) TMP_ALLOC ((rsize + msize) * BYTES_PER_MP_LIMB);
 390         if (rsize > msize)
 391           cy = mpn_mul (tp, rp, rsize, mp, msize);
 392         else
 393           cy = mpn_mul (tp, mp, msize, rp, rsize);
 394         rsize += msize;
 395         rsize -= cy == 0;
 396         rp = tp;
 397
 398         /* If we already output digits (for an integral part) pad
 399            leading zeros.  */
 400         if (digits_computed_so_far != 0)
 401           for (i = 0; i < nexp; i++)
 402             tstr[digits_computed_so_far++] = 0;
 403       }
 404
 405     while (digits_computed_so_far <= n_digits)
 406       {
 407         /* For speed: skip trailing zeroes.  */
 408         if (rp[0] == 0)
 409           {
 410             rp++;
 411             rsize--;
 412             if (rsize == 0)
 413               {
 414                 n_digits = digits_computed_so_far;
 415                 break;
 416               }
 417           }
 418
 419         cy = mpn_mul_1 (rp, rp, rsize, big_base);
 420         if (digits_computed_so_far == 0 && cy == 0)
 421           {
 422             abort ();
 423             nexp += dig_per_u;
 424             continue;
 425           }
 426         /* Convert N1 from BIG_BASE to a string of digits in BASE
 427            using single precision operations.  */
 428         {
 429           unsigned char *s = tstr + digits_computed_so_far + dig_per_u;
 430           for (i = dig_per_u - 1; i >= 0; i--)
 431             {
 432               *--s = cy % base;
 433               cy /= base;
 434             }
 435         }
 436         digits_computed_so_far += dig_per_u;
 437       }
 438     if (exp_in_base == 0)
 439       exp_in_base = -nexp;
 440   }
 441
 442  finish_up:
 443
 444   /* We can have at most one leading 0.  Remove it.  */
 445   if (tstr[0] == 0)
 446     {
 447       tstr++;
 448       digits_computed_so_far--;
 449       exp_in_base--;
 450     }
 451
 452   /* We should normally have computed too many digits.  Round the result
 453      at the point indicated by n_digits.  */
 454   if (digits_computed_so_far > n_digits)
 455     {
 456       /* Round the result.  */
 457       if (tstr[n_digits] * 2 >= base)
 458         {
 459           digits_computed_so_far = n_digits;
 460           for (i = n_digits - 1; i >= 0; i--)
 461             {
 462               unsigned int x;
 463               x = ++(tstr[i]);
 464               if (x < base)
 465                 goto rounded_ok;
 466               digits_computed_so_far--;
 467             }
 468           tstr[0] = 1;
 469           digits_computed_so_far = 1;
 470           exp_in_base++;
 471         rounded_ok:;
 472         }
 473     }
 474
 475   /* We might have fewer digits than requested as a result of rounding above,
 476      (i.e. 0.999999 => 1.0) or because we have a number that simply doesn't
 477      need many digits in this base (i.e., 0.125 in base 10).  */
 478   if (n_digits > digits_computed_so_far)
 479     n_digits = digits_computed_so_far;
 480
 481   /* Remove trailing 0.  There can be many zeros. */
 482   while (n_digits != 0 && tstr[n_digits - 1] == 0)
 483     n_digits--;
 484
 485   /* Translate to ascii and null-terminate.  */
 486   for (i = 0; i < n_digits; i++)
 487     *str++ = num_to_text[tstr[i]];
 488   *str = 0;
 489   *exp = exp_in_base;
 490   TMP_FREE (marker);
 491   return digit_ptr;
 492 }
 493
 494 #if COPY_THIS_TO_OTHER_PLACES
 495       /* Use this expression in lots of places in the library instead of the
 496          count_leading_zeros+expression that is used currently.  This expression
 497          is much more accurate and will save odles of memory.  */
 498       rsize = ((mp_size_t) (exp_in_base / __mp_bases[base].chars_per_bit_exactly)
 499                + BITS_PER_MP_LIMB) / BITS_PER_MP_LIMB;
 500 #endif