Annotation of OpenXM_contrib/gmp/mpf/get_str.c, Revision 1.1.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