[BACK]Return to get_str.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2-ssh-2 / mpf

Annotation of OpenXM/src/kan96xx/gmp-2.0.2-ssh-2/mpf/get_str.c, Revision 1.1.1.1

1.1       takayama    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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>