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
7 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
9 This file is part of the GNU MP Library.
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.
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.
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. */
32 New algorithm for converting fractions (951019):
33 0. Call the fraction to convert F.
34 1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e.,
35 [exp * BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly]. Exp is
36 the number of limbs between the limb point and the most significant
37 non-zero limb. Call this result n.
39 3. F*B^n will now be just below 1, which can be converted easily. (Just
40 multiply by B repeatedly, and see the digits fall out as integers.)
41 We should interrupt the conversion process of F*B^n as soon as the number
42 of digits requested have been generated.
44 New algorithm for converting integers (951019):
45 0. Call the integer to convert I.
46 1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e.,
47 [exp BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly]. Exp is
48 the number of limbs between the limb point and the least significant
49 non-zero limb. Call this result n.
51 3. I/B^n can be converted easily. (Just divide by B repeatedly. In GMP,
52 this is best done by calling mpn_get_str.)
53 Note that converting I/B^n could yield more digits than requested. For
54 efficiency, the variable n above should be set larger in such cases, to
55 kill all undesired digits in the division in step 3.
60 mpf_get_str (char *digit_ptr, mp_exp_t *exp, int base, size_t n_digits, mpf_srcptr u)
62 mpf_get_str (digit_ptr, exp, base, n_digits, u)
75 long i; /* should be size_t */
78 size_t digits_computed_so_far;
93 num_to_text = (char *) "0123456789abcdefghijklmnopqrstuvwxyz";
98 num_to_text = (char *) "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
101 /* Don't compute more digits than U can accurately represent.
102 Also, if 0 digits were requested, give *exactly* as many digits
103 as can be accurately represented. */
105 size_t max_digits = (((u->_mp_prec - 1) * BITS_PER_MP_LIMB)
106 * __mp_bases[base].chars_per_bit_exactly);
107 if (n_digits == 0 || n_digits > max_digits)
108 n_digits = max_digits;
113 /* We didn't get a string from the user. Allocate one (and return
114 a pointer to it) with space for `-' and terminating null. */
115 digit_ptr = (char *) (*_mp_allocate_func) (n_digits + 2);
125 str = (unsigned char *) digit_ptr;
127 /* Allocate temporary digit space. We can't put digits directly in the user
128 area, since we almost always generate more digits than requested. */
129 tstr = (unsigned char *) TMP_ALLOC (n_digits + 3 * BITS_PER_MP_LIMB);
138 digits_computed_so_far = 0;
142 /* The number has just an integral part. */
144 mp_size_t exp_in_limbs;
152 n_limbs = 2 + ((mp_size_t) (n_digits / __mp_bases[base].chars_per_bit_exactly)
155 /* Compute n such that [u/B^n] contains (somewhat) more than n_digits
156 digits. (We compute less than that only if that is an exact number,
157 i.e., exp is small enough.) */
161 if (n_limbs >= exp_in_limbs)
163 /* The number is so small that we convert the entire number. */
165 rp = (mp_ptr) TMP_ALLOC (exp_in_limbs * BYTES_PER_MP_LIMB);
166 MPN_ZERO (rp, exp_in_limbs - usize);
167 MPN_COPY (rp + (exp_in_limbs - usize), u->_mp_d, usize);
168 rsize = exp_in_limbs;
172 exp_in_limbs -= n_limbs;
173 exp_in_base = (((exp_in_limbs * BITS_PER_MP_LIMB - 1))
174 * __mp_bases[base].chars_per_bit_exactly);
176 rsize = exp_in_limbs + 1;
177 rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
178 tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
183 count_leading_zeros (cnt, exp_in_base);
184 for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
186 mpn_mul_n (tp, rp, rp, rsize);
188 rsize -= tp[rsize - 1] == 0;
189 xp = tp; tp = rp; rp = xp;
191 if (((exp_in_base >> i) & 1) != 0)
193 cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
208 mp_ptr tmp = (mp_ptr) TMP_ALLOC ((rsize+1)* BYTES_PER_MP_LIMB);
209 MPN_ZERO (tmp, rsize - msize);
210 MPN_COPY (tmp + rsize - msize, mp, msize);
216 mp_ptr tmp = (mp_ptr) TMP_ALLOC ((msize+1)* BYTES_PER_MP_LIMB);
217 MPN_COPY (tmp, mp, msize);
220 count_leading_zeros (cnt, rp[rsize - 1]);
224 mpn_lshift (rp, rp, rsize, cnt);
225 cy = mpn_lshift (mp, mp, msize, cnt);
231 mp_size_t qsize = n_limbs + (cy != 0);
232 qp = (mp_ptr) TMP_ALLOC ((qsize + 1) * BYTES_PER_MP_LIMB);
233 xtra = qsize - (msize - rsize);
234 qflag = mpn_divrem (qp, xtra, mp, msize, rp, rsize);
236 rsize = qsize + qflag;
242 str_size = mpn_get_str (tstr, base, rp, rsize);
244 if (str_size > n_digits + 3 * BITS_PER_MP_LIMB)
248 while (tstr[start_str] == 0)
251 for (i = start_str; i < (int) str_size; i++)
253 tstr[digits_computed_so_far++] = tstr[i];
254 if (digits_computed_so_far > n_digits)
257 exp_in_base = exp_in_base + str_size - start_str;
265 /* The number has an integral part, convert that first.
266 If there is a fractional part too, it will be handled later. */
269 rp = (mp_ptr) TMP_ALLOC (uexp * BYTES_PER_MP_LIMB);
270 up = u->_mp_d + usize - uexp;
271 MPN_COPY (rp, up, uexp);
273 str_size = mpn_get_str (tstr, base, rp, uexp);
276 while (tstr[start_str] == 0)
279 for (i = start_str; i < (int) str_size; i++)
281 tstr[digits_computed_so_far++] = tstr[i];
282 if (digits_computed_so_far > n_digits)
284 exp_in_base = str_size - start_str;
289 exp_in_base = str_size - start_str;
290 /* Modify somewhat and fall out to convert fraction... */
298 /* Convert the fraction. */
300 mp_size_t rsize, msize;
301 mp_ptr rp, tp, xp, mp;
306 big_base = __mp_bases[base].big_base;
307 dig_per_u = __mp_bases[base].chars_per_limb;
309 /* Hack for correctly (although not efficiently) converting to bases that
310 are powers of 2. If we deem it important, we could handle powers of 2
311 by shifting and masking (just like mpn_get_str). */
312 if (big_base < 10) /* logarithm of base when power of two */
314 int logbase = big_base;
315 if (dig_per_u * logbase == BITS_PER_MP_LIMB)
317 big_base = (mp_limb_t) 1 << (dig_per_u * logbase);
318 /* fall out to general code... */
324 rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
326 MPN_COPY (rp, up, usize);
334 if (u->_mp_d[usize - 1] == 0)
337 count_leading_zeros (cnt, u->_mp_d[usize - 1]);
339 nexp = ((uexp * BITS_PER_MP_LIMB) + cnt)
340 * __mp_bases[base].chars_per_bit_exactly;
344 rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
346 MPN_COPY (rp, up, usize);
352 rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
353 tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
358 count_leading_zeros (cnt, nexp);
359 for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
361 mpn_mul_n (tp, rp, rp, rsize);
363 rsize -= tp[rsize - 1] == 0;
364 xp = tp; tp = rp; rp = xp;
366 if (((nexp >> i) & 1) != 0)
368 cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
374 /* Did our multiplier (base^nexp) cancel with uexp? */
380 cy = mpn_mul_1 (rp, rp, rsize, big_base);
390 tp = (mp_ptr) TMP_ALLOC ((rsize + msize) * BYTES_PER_MP_LIMB);
392 cy = mpn_mul (tp, rp, rsize, mp, msize);
394 cy = mpn_mul (tp, mp, msize, rp, rsize);
399 /* If we already output digits (for an integral part) pad
401 if (digits_computed_so_far != 0)
402 for (i = 0; i < nexp; i++)
403 tstr[digits_computed_so_far++] = 0;
406 while (digits_computed_so_far <= n_digits)
408 /* For speed: skip trailing zeroes. */
415 n_digits = digits_computed_so_far;
420 cy = mpn_mul_1 (rp, rp, rsize, big_base);
421 if (digits_computed_so_far == 0 && cy == 0)
427 /* Convert N1 from BIG_BASE to a string of digits in BASE
428 using single precision operations. */
430 unsigned char *s = tstr + digits_computed_so_far + dig_per_u;
431 for (i = dig_per_u - 1; i >= 0; i--)
437 digits_computed_so_far += dig_per_u;
439 if (exp_in_base == 0)
445 /* We can have at most one leading 0. Remove it. */
449 digits_computed_so_far--;
453 /* We should normally have computed too many digits. Round the result
454 at the point indicated by n_digits. */
455 if (digits_computed_so_far > n_digits)
457 /* Round the result. */
458 if (tstr[n_digits] * 2 >= base)
460 digits_computed_so_far = n_digits;
461 for (i = n_digits - 1; i >= 0; i--)
467 digits_computed_so_far--;
470 digits_computed_so_far = 1;
476 /* We might have fewer digits than requested as a result of rounding above,
477 (i.e. 0.999999 => 1.0) or because we have a number that simply doesn't
478 need many digits in this base (i.e., 0.125 in base 10). */
479 if (n_digits > digits_computed_so_far)
480 n_digits = digits_computed_so_far;
482 /* Remove trailing 0. There can be many zeros. */
483 while (n_digits != 0 && tstr[n_digits - 1] == 0)
486 /* Translate to ascii and null-terminate. */
487 for (i = 0; i < (int) n_digits; i++)
488 *str++ = num_to_text[tstr[i]];
495 #if COPY_THIS_TO_OTHER_PLACES
496 /* Use this expression in lots of places in the library instead of the
497 count_leading_zeros+expression that is used currently. This expression
498 is much more accurate and will save odles of memory. */
499 rsize = ((mp_size_t) (exp_in_base / __mp_bases[base].chars_per_bit_exactly)
500 + BITS_PER_MP_LIMB) / BITS_PER_MP_LIMB;