Annotation of OpenXM_contrib/gmp/mpf/get_str.c, Revision 1.1
1.1 ! maekawa 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>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to get_str.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpf