Annotation of OpenXM_contrib/gmp/mpf/get_str.c, Revision 1.1.1.2
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:
1.1.1.2 ! maekawa 7: Copyright (C) 1993, 1994, 1995, 1996, 1997, 2000 Free Software Foundation,
! 8: Inc.
1.1 maekawa 9:
10: This file is part of the GNU MP Library.
11:
12: The GNU MP Library is free software; you can redistribute it and/or modify
1.1.1.2 ! maekawa 13: it under the terms of the GNU Lesser General Public License as published by
! 14: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 15: option) any later version.
16:
17: The GNU MP Library is distributed in the hope that it will be useful, but
18: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! maekawa 19: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 20: License for more details.
21:
1.1.1.2 ! maekawa 22: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 23: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
24: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25: MA 02111-1307, USA. */
26:
27: #include "gmp.h"
28: #include "gmp-impl.h"
29: #include "longlong.h"
30:
31: /*
1.1.1.2 ! maekawa 32: The conversion routine works like this:
! 33:
! 34: 1. If U >= 1, compute U' = U / base**n, where n is chosen such that U' is
! 35: the largest number smaller than 1.
! 36: 2. Else, if U < 1, compute U' = U * base**n, where n is chosen such that U'
! 37: is the largest number smaller than 1.
! 38: 3. Convert U' (by repeatedly multiplying it by base). This process can
! 39: easily be interrupted when the needed number of digits are generated.
1.1 maekawa 40: */
41:
42: char *
43: #if __STDC__
44: mpf_get_str (char *digit_ptr, mp_exp_t *exp, int base, size_t n_digits, mpf_srcptr u)
45: #else
46: mpf_get_str (digit_ptr, exp, base, n_digits, u)
47: char *digit_ptr;
48: mp_exp_t *exp;
49: int base;
50: size_t n_digits;
51: mpf_srcptr u;
52: #endif
53: {
1.1.1.2 ! maekawa 54: mp_ptr up;
1.1 maekawa 55: mp_size_t usize;
56: mp_exp_t uexp;
1.1.1.2 ! maekawa 57: mp_size_t prec;
1.1 maekawa 58: unsigned char *str;
59: char *num_to_text;
60: mp_ptr rp;
1.1.1.2 ! maekawa 61: mp_size_t rsize;
1.1 maekawa 62: mp_limb_t big_base;
63: size_t digits_computed_so_far;
64: int dig_per_u;
65: unsigned char *tstr;
66: mp_exp_t exp_in_base;
1.1.1.2 ! maekawa 67: int cnt;
1.1 maekawa 68: TMP_DECL (marker);
69:
70: TMP_MARK (marker);
71: usize = u->_mp_size;
72: uexp = u->_mp_exp;
1.1.1.2 ! maekawa 73: prec = u->_mp_prec + 1;
1.1 maekawa 74:
75: if (base >= 0)
76: {
77: if (base == 0)
78: base = 10;
79: num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz";
80: }
81: else
82: {
83: base = -base;
84: num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
85: }
86:
87: /* Don't compute more digits than U can accurately represent.
88: Also, if 0 digits were requested, give *exactly* as many digits
89: as can be accurately represented. */
90: {
1.1.1.2 ! maekawa 91: size_t max_digits = 2 + (size_t) (((prec - 2) * BITS_PER_MP_LIMB)
! 92: * __mp_bases[base].chars_per_bit_exactly);
1.1 maekawa 93: if (n_digits == 0 || n_digits > max_digits)
94: n_digits = max_digits;
1.1.1.2 ! maekawa 95: #if 0
! 96: /* This seems to work, but check it better before enabling it. */
! 97: else
! 98: /* Limit the computation precision if only a limited digits are
! 99: desired. We could probably decrease both this, and avoid the +1
! 100: for setting prec above. */
! 101: prec = 2 + (mp_size_t)
! 102: (n_digits / (BITS_PER_MP_LIMB * __mp_bases[base].chars_per_bit_exactly));
! 103: #endif
1.1 maekawa 104: }
105:
106: if (digit_ptr == 0)
107: {
108: /* We didn't get a string from the user. Allocate one (and return
109: a pointer to it) with space for `-' and terminating null. */
110: digit_ptr = (char *) (*_mp_allocate_func) (n_digits + 2);
111: }
112:
113: if (usize == 0)
114: {
115: *exp = 0;
116: *digit_ptr = 0;
117: return digit_ptr;
118: }
119:
120: str = (unsigned char *) digit_ptr;
121:
122: if (usize < 0)
123: {
124: *digit_ptr = '-';
125: str++;
126: usize = -usize;
127: }
128:
1.1.1.2 ! maekawa 129: up = PTR (u);
1.1 maekawa 130:
1.1.1.2 ! maekawa 131: if (uexp > 0)
1.1 maekawa 132: {
1.1.1.2 ! maekawa 133: /* U >= 1. Compute U' = U / base**n, where n is chosen such that U' < 1. */
! 134: mp_size_t ralloc;
! 135: mp_ptr tp;
! 136: int i;
1.1 maekawa 137:
1.1.1.2 ! maekawa 138: /* Limit the number of digits to develop for small integers. */
! 139: #if 0
! 140: if (exp_in_base < n_digits)
! 141: n_digits = exp_in_base;
! 142: #endif
1.1 maekawa 143:
1.1.1.2 ! maekawa 144: count_leading_zeros (cnt, up[usize - 1]);
! 145: exp_in_base = (((double) uexp * BITS_PER_MP_LIMB - cnt)
! 146: * __mp_bases[base].chars_per_bit_exactly);
! 147: exp_in_base += 1;
! 148:
! 149: ralloc = (prec + 1) * 2;
! 150: rp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB);
! 151: tp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB);
! 152:
! 153: rp[0] = base;
! 154: rsize = 1;
! 155: count_leading_zeros (cnt, exp_in_base);
! 156: for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
! 157: {
! 158: mpn_mul_n (tp, rp, rp, rsize);
! 159: rsize = 2 * rsize;
! 160: rsize -= tp[rsize - 1] == 0;
1.1 maekawa 161:
1.1.1.2 ! maekawa 162: if (rsize > prec)
! 163: {
! 164: MPN_COPY (rp, tp + rsize - prec, prec + 1);
! 165: rsize = prec;
! 166: }
! 167: else
! 168: MP_PTR_SWAP (rp, tp);
1.1 maekawa 169:
1.1.1.2 ! maekawa 170: if (((exp_in_base >> i) & 1) != 0)
! 171: {
! 172: mp_limb_t cy;
! 173: cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
! 174: rp[rsize] = cy;
! 175: rsize += cy != 0;
! 176: }
! 177: }
! 178:
! 179: count_leading_zeros (cnt, rp[rsize - 1]);
! 180: if (cnt != 0)
1.1 maekawa 181: {
1.1.1.2 ! maekawa 182: mpn_lshift (rp, rp, rsize, cnt);
! 183:
! 184: if (usize < rsize)
! 185: {
! 186: /* Pad out U to the size of R while shifting it.
! 187: (Reuse temporary space at tp.) */
! 188: mp_limb_t cy;
! 189:
! 190: MPN_ZERO (tp, rsize - usize);
! 191: cy = mpn_lshift (tp + rsize - usize, up, usize, cnt);
! 192: up = tp;
! 193: usize = rsize;
! 194: if (cy)
! 195: up[usize++] = cy;
! 196: ASSERT_ALWAYS (usize <= ralloc); /* sufficient space? */
! 197: }
! 198: else
! 199: {
! 200: /* Copy U to temporary space. */
! 201: /* FIXME: Allocate more space for tp above, and reuse it here. */
! 202: mp_limb_t cy;
! 203: mp_ptr tup = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB);
! 204:
! 205: cy = mpn_lshift (tup, up, usize, cnt);
! 206: up = tup;
! 207: if (cy)
! 208: up[usize++] = cy;
! 209: }
1.1 maekawa 210: }
211: else
212: {
1.1.1.2 ! maekawa 213: if (usize < rsize)
! 214: {
! 215: /* Pad out U to the size of R. (Reuse temporary space at tp.) */
! 216: MPN_ZERO (tp, rsize - usize);
! 217: MPN_COPY (tp + rsize - usize, up, usize);
! 218: up = tp;
! 219: usize = rsize;
! 220: }
! 221: else
! 222: {
! 223: /* Copy U to temporary space. */
! 224: mp_ptr tmp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
! 225: MPN_COPY (tmp, up, usize);
! 226: up = tmp;
! 227: }
! 228: }
! 229:
! 230: {
! 231: mp_ptr qp;
! 232: qp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
! 233: mpn_divrem (qp, prec - (usize - rsize), up, usize, rp, rsize);
! 234: rsize = prec;
! 235: rp = qp;
! 236: }
! 237: }
! 238: else
! 239: {
! 240: /* U < 1. Compute U' = U * base**n, where n is chosen such that U' is
! 241: the greatest number that still satisfies U' < 1. */
! 242: mp_size_t ralloc;
! 243: mp_ptr tp;
! 244: int i;
! 245:
! 246: uexp = -uexp;
! 247: count_leading_zeros (cnt, up[usize - 1]);
! 248: exp_in_base = (((double) uexp * BITS_PER_MP_LIMB + cnt - 1)
! 249: * __mp_bases[base].chars_per_bit_exactly);
! 250: if (exp_in_base < 0)
! 251: exp_in_base = 0;
! 252:
! 253: if (exp_in_base != 0)
! 254: {
! 255: ralloc = (prec + 1) * 2;
! 256: rp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB);
! 257: tp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB);
1.1 maekawa 258:
259: rp[0] = base;
260: rsize = 1;
261: count_leading_zeros (cnt, exp_in_base);
262: for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
263: {
264: mpn_mul_n (tp, rp, rp, rsize);
265: rsize = 2 * rsize;
266: rsize -= tp[rsize - 1] == 0;
1.1.1.2 ! maekawa 267: if (rsize > prec)
! 268: {
! 269: MPN_COPY (rp, tp + rsize - prec, prec + 1);
! 270: rsize = prec;
! 271: }
! 272: else
! 273: MP_PTR_SWAP (rp, tp);
1.1 maekawa 274:
275: if (((exp_in_base >> i) & 1) != 0)
276: {
1.1.1.2 ! maekawa 277: mp_limb_t cy;
1.1 maekawa 278: cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
279: rp[rsize] = cy;
280: rsize += cy != 0;
281: }
282: }
283:
284: {
1.1.1.2 ! maekawa 285: mp_limb_t cy;
! 286: tp = (mp_ptr) TMP_ALLOC ((rsize + usize) * BYTES_PER_MP_LIMB);
! 287: if (rsize > usize)
! 288: cy = mpn_mul (tp, rp, rsize, up, usize);
1.1 maekawa 289: else
1.1.1.2 ! maekawa 290: cy = mpn_mul (tp, up, usize, rp, rsize);
! 291: rsize += usize;
! 292: rsize -= cy == 0;
! 293: rp = tp;
1.1 maekawa 294: }
1.1.1.2 ! maekawa 295: exp_in_base = -exp_in_base;
1.1 maekawa 296: }
1.1.1.2 ! maekawa 297: else
1.1 maekawa 298: {
1.1.1.2 ! maekawa 299: rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
! 300: MPN_COPY (rp, up, usize);
! 301: rsize = usize;
1.1 maekawa 302: }
303: }
304:
1.1.1.2 ! maekawa 305: big_base = __mp_bases[base].big_base;
! 306: dig_per_u = __mp_bases[base].chars_per_limb;
1.1 maekawa 307:
1.1.1.2 ! maekawa 308: /* Hack for correctly (although not optimally) converting to bases that are
! 309: powers of 2. If we deem it important, we could handle powers of 2 by
! 310: shifting and masking (just like mpn_get_str). */
! 311: if (big_base < 10) /* logarithm of base when power of two */
1.1 maekawa 312: {
1.1.1.2 ! maekawa 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: }
1.1 maekawa 319:
1.1.1.2 ! maekawa 320: /* Now that we have normalized the number, develop the digits, essentially by
! 321: multiplying it by BASE. We initially develop at least 3 extra digits,
! 322: since the two leading digits might become zero, and we need one extra for
! 323: rounding the output properly. */
1.1 maekawa 324:
1.1.1.2 ! maekawa 325: /* Allocate temporary digit space. We can't put digits directly in the user
! 326: area, since we generate more digits than requested. (We allocate
! 327: BITS_PER_MP_LIMB extra bytes because of the digit block nature of the
! 328: conversion.) */
! 329: tstr = (unsigned char *) TMP_ALLOC (n_digits + BITS_PER_MP_LIMB + 3);
1.1 maekawa 330:
1.1.1.2 ! maekawa 331: for (digits_computed_so_far = 0; digits_computed_so_far < n_digits + 3;
! 332: digits_computed_so_far += dig_per_u)
! 333: {
! 334: mp_limb_t cy;
! 335: /* For speed: skip trailing zeroes. */
! 336: if (rp[0] == 0)
1.1 maekawa 337: {
1.1.1.2 ! maekawa 338: rp++;
! 339: rsize--;
! 340: if (rsize == 0)
1.1 maekawa 341: {
1.1.1.2 ! maekawa 342: n_digits = digits_computed_so_far;
! 343: break;
1.1 maekawa 344: }
345: }
346:
1.1.1.2 ! maekawa 347: cy = mpn_mul_1 (rp, rp, rsize, big_base);
1.1 maekawa 348:
1.1.1.2 ! maekawa 349: ASSERT_ALWAYS (! (digits_computed_so_far == 0 && cy == 0));
1.1 maekawa 350:
1.1.1.2 ! maekawa 351: /* Convert N1 from BIG_BASE to a string of digits in BASE
! 352: using single precision operations. */
1.1 maekawa 353: {
1.1.1.2 ! maekawa 354: int i;
! 355: unsigned char *s = tstr + digits_computed_so_far + dig_per_u;
! 356: for (i = dig_per_u - 1; i >= 0; i--)
1.1 maekawa 357: {
1.1.1.2 ! maekawa 358: *--s = cy % base;
! 359: cy /= base;
1.1 maekawa 360: }
361: }
1.1.1.2 ! maekawa 362: }
1.1 maekawa 363:
1.1.1.2 ! maekawa 364: /* We can have at most two leading 0. Remove them. */
1.1 maekawa 365: if (tstr[0] == 0)
366: {
367: tstr++;
368: digits_computed_so_far--;
369: exp_in_base--;
370:
1.1.1.2 ! maekawa 371: if (tstr[0] == 0)
1.1 maekawa 372: {
1.1.1.2 ! maekawa 373: tstr++;
! 374: digits_computed_so_far--;
! 375: exp_in_base--;
! 376:
! 377: if (tstr[0] == 0)
! 378: abort ();
1.1 maekawa 379: }
380: }
381:
1.1.1.2 ! maekawa 382: {
! 383: size_t i;
! 384:
! 385: /* We should normally have computed too many digits. Round the result
! 386: at the point indicated by n_digits. */
! 387: if (digits_computed_so_far > n_digits)
! 388: {
! 389: /* Round the result. */
! 390: if (tstr[n_digits] * 2 >= base)
! 391: {
! 392: digits_computed_so_far = n_digits;
! 393: for (i = n_digits - 1;; i--)
! 394: {
! 395: unsigned int x;
! 396: x = ++(tstr[i]);
! 397: if (x != base)
! 398: break;
! 399: digits_computed_so_far--;
! 400: if (i == 0)
! 401: {
! 402: /* We had something like `9999999...9d', where 2*d >= base.
! 403: This rounds up to `1', increasing the exponent. */
! 404: tstr[0] = 1;
! 405: digits_computed_so_far = 1;
! 406: exp_in_base++;
! 407: break;
! 408: }
! 409: }
! 410: }
! 411: }
! 412:
! 413: /* We might have fewer digits than requested as a result of rounding above,
! 414: (i.e. 0.999999 => 1.0) or because we have a number that simply doesn't
! 415: need many digits in this base (i.e., 0.125 in base 10). */
! 416: if (n_digits > digits_computed_so_far)
! 417: n_digits = digits_computed_so_far;
! 418:
! 419: /* Remove trailing 0. There can be many zeros. */
! 420: while (n_digits != 0 && tstr[n_digits - 1] == 0)
! 421: n_digits--;
! 422:
! 423: /* Translate to ascii and null-terminate. */
! 424: for (i = 0; i < n_digits; i++)
! 425: *str++ = num_to_text[tstr[i]];
! 426: }
1.1 maekawa 427: *str = 0;
428: *exp = exp_in_base;
429: TMP_FREE (marker);
430: return digit_ptr;
431: }
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