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