Annotation of OpenXM_contrib/gmp/mpf/get_str.c, Revision 1.1.1.3
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;
1.1.1.3 ! maekawa 335: /* For speed: skip trailing zero limbs. */
1.1.1.2 maekawa 336: if (rp[0] == 0)
1.1 maekawa 337: {
1.1.1.2 maekawa 338: rp++;
339: rsize--;
340: if (rsize == 0)
1.1.1.3 ! maekawa 341: break;
1.1 maekawa 342: }
343:
1.1.1.2 maekawa 344: cy = mpn_mul_1 (rp, rp, rsize, big_base);
1.1 maekawa 345:
1.1.1.2 maekawa 346: ASSERT_ALWAYS (! (digits_computed_so_far == 0 && cy == 0));
1.1 maekawa 347:
1.1.1.2 maekawa 348: /* Convert N1 from BIG_BASE to a string of digits in BASE
349: using single precision operations. */
1.1 maekawa 350: {
1.1.1.2 maekawa 351: int i;
352: unsigned char *s = tstr + digits_computed_so_far + dig_per_u;
353: for (i = dig_per_u - 1; i >= 0; i--)
1.1 maekawa 354: {
1.1.1.2 maekawa 355: *--s = cy % base;
356: cy /= base;
1.1 maekawa 357: }
358: }
1.1.1.2 maekawa 359: }
1.1 maekawa 360:
1.1.1.2 maekawa 361: /* We can have at most two leading 0. Remove them. */
1.1 maekawa 362: if (tstr[0] == 0)
363: {
364: tstr++;
365: digits_computed_so_far--;
366: exp_in_base--;
367:
1.1.1.2 maekawa 368: if (tstr[0] == 0)
1.1 maekawa 369: {
1.1.1.2 maekawa 370: tstr++;
371: digits_computed_so_far--;
372: exp_in_base--;
373:
1.1.1.3 ! maekawa 374: ASSERT_ALWAYS (tstr[0] != 0);
1.1 maekawa 375: }
376: }
377:
1.1.1.3 ! maekawa 378: /* We should normally have computed too many digits. Round the result
! 379: at the point indicated by n_digits. */
! 380: if (digits_computed_so_far > n_digits)
! 381: {
! 382: size_t i;
! 383: /* Round the result. */
! 384: if (tstr[n_digits] * 2 >= base)
! 385: {
! 386: digits_computed_so_far = n_digits;
! 387: for (i = n_digits - 1;; i--)
! 388: {
! 389: unsigned int x;
! 390: x = ++(tstr[i]);
! 391: if (x != base)
! 392: break;
! 393: digits_computed_so_far--;
! 394: if (i == 0)
! 395: {
! 396: /* We had something like `bbbbbbb...bd', where 2*d >= base
! 397: and `b' denotes digit with significance base - 1.
! 398: This rounds up to `1', increasing the exponent. */
! 399: tstr[0] = 1;
! 400: digits_computed_so_far = 1;
! 401: exp_in_base++;
1.1.1.2 maekawa 402: break;
1.1.1.3 ! maekawa 403: }
! 404: }
! 405: }
! 406: }
1.1.1.2 maekawa 407:
1.1.1.3 ! maekawa 408: /* We might have fewer digits than requested as a result of rounding above,
! 409: (i.e. 0.999999 => 1.0) or because we have a number that simply doesn't
! 410: need many digits in this base (e.g., 0.125 in base 10). */
! 411: if (n_digits > digits_computed_so_far)
! 412: n_digits = digits_computed_so_far;
! 413:
! 414: /* Remove trailing 0. There can be many zeros. */
! 415: while (n_digits != 0 && tstr[n_digits - 1] == 0)
! 416: n_digits--;
1.1.1.2 maekawa 417:
1.1.1.3 ! maekawa 418: /* Translate to ASCII and copy to result string. */
! 419: while (n_digits != 0)
! 420: {
! 421: *str++ = num_to_text[*tstr++];
1.1.1.2 maekawa 422: n_digits--;
1.1.1.3 ! maekawa 423: }
1.1.1.2 maekawa 424:
1.1 maekawa 425: *str = 0;
426: *exp = exp_in_base;
427: TMP_FREE (marker);
428: return digit_ptr;
429: }
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