checkin of 0.3.0
[pspp-builds.git] / lib / gmp / mpf / get_str.c
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 <config.h>
27 #include "gmp.h"
28 #include "gmp-impl.h"
29 #include "longlong.h"
30
31 /*
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.
38    2. Compute B^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.
43
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.
50    2. Compute B^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.
56 */
57
58 char *
59 #if __STDC__
60 mpf_get_str (char *digit_ptr, mp_exp_t *exp, int base, size_t n_digits, mpf_srcptr u)
61 #else
62 mpf_get_str (digit_ptr, exp, base, n_digits, u)
63      char *digit_ptr;
64      mp_exp_t *exp;
65      int base;
66      size_t n_digits;
67      mpf_srcptr u;
68 #endif
69 {
70   mp_size_t usize;
71   mp_exp_t uexp;
72   unsigned char *str;
73   size_t str_size;
74   char *num_to_text;
75   long i;                       /* should be size_t */
76   mp_ptr rp;
77   mp_limb_t big_base;
78   size_t digits_computed_so_far;
79   int dig_per_u;
80   mp_srcptr up;
81   unsigned char *tstr;
82   mp_exp_t exp_in_base;
83   TMP_DECL (marker);
84
85   TMP_MARK (marker);
86   usize = u->_mp_size;
87   uexp = u->_mp_exp;
88
89   if (base >= 0)
90     {
91       if (base == 0)
92         base = 10;
93       num_to_text = (char *) "0123456789abcdefghijklmnopqrstuvwxyz";
94     }
95   else
96     {
97       base = -base;
98       num_to_text = (char *) "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
99     }
100
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.  */
104   {
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;
109   }
110
111   if (digit_ptr == 0)
112     {
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);
116     }
117
118   if (usize == 0)
119     {
120       *exp = 0;
121       *digit_ptr = 0;
122       return digit_ptr;
123     }
124
125   str = (unsigned char *) digit_ptr;
126
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);
130
131   if (usize < 0)
132     {
133       *digit_ptr = '-';
134       str++;
135       usize = -usize;
136     }
137
138   digits_computed_so_far = 0;
139
140   if (uexp > usize)
141     {
142       /* The number has just an integral part.  */
143       mp_size_t rsize;
144       mp_size_t exp_in_limbs;
145       mp_size_t msize;
146       mp_ptr tp, xp, mp;
147       int cnt;
148       mp_limb_t cy;
149       mp_size_t start_str;
150       mp_size_t n_limbs;
151
152       n_limbs = 2 + ((mp_size_t) (n_digits / __mp_bases[base].chars_per_bit_exactly)
153                      / BITS_PER_MP_LIMB);
154
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.)  */
158
159       exp_in_limbs = uexp;
160
161       if (n_limbs >= exp_in_limbs)
162         {
163           /* The number is so small that we convert the entire number.  */
164           exp_in_base = 0;
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;
169         }
170       else
171         {
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);
175
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);
179
180           rp[0] = base;
181           rsize = 1;
182
183           count_leading_zeros (cnt, exp_in_base);
184           for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
185             {
186               mpn_mul_n (tp, rp, rp, rsize);
187               rsize = 2 * rsize;
188               rsize -= tp[rsize - 1] == 0;
189               xp = tp; tp = rp; rp = xp;
190
191               if (((exp_in_base >> i) & 1) != 0)
192                 {
193                   cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
194                   rp[rsize] = cy;
195                   rsize += cy != 0;
196                 }
197             }
198
199           mp = u->_mp_d;
200           msize = usize;
201
202           {
203             mp_ptr qp;
204             mp_limb_t qflag;
205             mp_size_t xtra;
206             if (msize < rsize)
207               {
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);
211                 mp = tmp;
212                 msize = rsize;
213               }
214             else
215               {
216                 mp_ptr tmp = (mp_ptr) TMP_ALLOC ((msize+1)* BYTES_PER_MP_LIMB);
217                 MPN_COPY (tmp, mp, msize);
218                 mp = tmp;
219               }
220             count_leading_zeros (cnt, rp[rsize - 1]);
221             cy = 0;
222             if (cnt != 0)
223               {
224                 mpn_lshift (rp, rp, rsize, cnt);
225                 cy = mpn_lshift (mp, mp, msize, cnt);
226                 if (cy)
227                   mp[msize++] = cy;
228               }
229
230             {
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);
235               qp[qsize] = qflag;
236               rsize = qsize + qflag;
237               rp = qp;
238             }
239           }
240         }
241
242       str_size = mpn_get_str (tstr, base, rp, rsize);
243
244       if (str_size > n_digits + 3 * BITS_PER_MP_LIMB)
245         abort ();
246
247       start_str = 0;
248       while (tstr[start_str] == 0)
249         start_str++;
250
251       for (i = start_str; i < (int) str_size; i++)
252         {
253           tstr[digits_computed_so_far++] = tstr[i];
254           if (digits_computed_so_far > n_digits)
255             break;
256         }
257       exp_in_base = exp_in_base + str_size - start_str;
258       goto finish_up;
259     }
260
261   exp_in_base = 0;
262
263   if (uexp > 0)
264     {
265       /* The number has an integral part, convert that first.
266          If there is a fractional part too, it will be handled later.  */
267       mp_size_t start_str;
268
269       rp = (mp_ptr) TMP_ALLOC (uexp * BYTES_PER_MP_LIMB);
270       up = u->_mp_d + usize - uexp;
271       MPN_COPY (rp, up, uexp);
272
273       str_size = mpn_get_str (tstr, base, rp, uexp);
274
275       start_str = 0;
276       while (tstr[start_str] == 0)
277         start_str++;
278
279       for (i = start_str; i < (int) str_size; i++)
280         {
281           tstr[digits_computed_so_far++] = tstr[i];
282           if (digits_computed_so_far > n_digits)
283             {
284               exp_in_base = str_size - start_str;
285               goto finish_up;
286             }
287         }
288
289       exp_in_base = str_size - start_str;
290       /* Modify somewhat and fall out to convert fraction... */
291       usize -= uexp;
292       uexp = 0;
293     }
294
295   if (usize <= 0)
296     goto finish_up;
297
298   /* Convert the fraction.  */
299   {
300     mp_size_t rsize, msize;
301     mp_ptr rp, tp, xp, mp;
302     int cnt;
303     mp_limb_t cy;
304     mp_exp_t nexp;
305
306     big_base = __mp_bases[base].big_base;
307     dig_per_u = __mp_bases[base].chars_per_limb;
308
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 */
313       {
314         int logbase = big_base;
315         if (dig_per_u * logbase == BITS_PER_MP_LIMB)
316           dig_per_u--;
317         big_base = (mp_limb_t) 1 << (dig_per_u * logbase);
318         /* fall out to general code... */
319       }
320
321 #if 0
322     if (0 && uexp == 0)
323       {
324         rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
325         up = u->_mp_d;
326         MPN_COPY (rp, up, usize);
327         rsize = usize;
328         nexp = 0;
329       }
330     else
331       {}
332 #endif
333     uexp = -uexp;
334     if (u->_mp_d[usize - 1] == 0)
335       cnt = 0;
336     else
337       count_leading_zeros (cnt, u->_mp_d[usize - 1]);
338
339     nexp = ((uexp * BITS_PER_MP_LIMB) + cnt)
340       * __mp_bases[base].chars_per_bit_exactly;
341
342     if (nexp == 0)
343       {
344         rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
345         up = u->_mp_d;
346         MPN_COPY (rp, up, usize);
347         rsize = usize;
348       }
349     else
350       {
351         rsize = uexp + 2;
352         rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
353         tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
354
355         rp[0] = base;
356         rsize = 1;
357
358         count_leading_zeros (cnt, nexp);
359         for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
360           {
361             mpn_mul_n (tp, rp, rp, rsize);
362             rsize = 2 * rsize;
363             rsize -= tp[rsize - 1] == 0;
364             xp = tp; tp = rp; rp = xp;
365
366             if (((nexp >> i) & 1) != 0)
367               {
368                 cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base);
369                 rp[rsize] = cy;
370                 rsize += cy != 0;
371               }
372           }
373
374         /* Did our multiplier (base^nexp) cancel with uexp?  */
375 #if 0
376         if (uexp != rsize)
377           {
378             do
379               {
380                 cy = mpn_mul_1 (rp, rp, rsize, big_base);
381                 nexp += dig_per_u;
382               }
383             while (cy == 0);
384             rp[rsize++] = cy;
385           }
386 #endif
387         mp = u->_mp_d;
388         msize = usize;
389
390         tp = (mp_ptr) TMP_ALLOC ((rsize + msize) * BYTES_PER_MP_LIMB);
391         if (rsize > msize)
392           cy = mpn_mul (tp, rp, rsize, mp, msize);
393         else
394           cy = mpn_mul (tp, mp, msize, rp, rsize);
395         rsize += msize;
396         rsize -= cy == 0;
397         rp = tp;
398
399         /* If we already output digits (for an integral part) pad
400            leading zeros.  */
401         if (digits_computed_so_far != 0)
402           for (i = 0; i < nexp; i++)
403             tstr[digits_computed_so_far++] = 0;
404       }
405
406     while (digits_computed_so_far <= n_digits)
407       {
408         /* For speed: skip trailing zeroes.  */
409         if (rp[0] == 0)
410           {
411             rp++;
412             rsize--;
413             if (rsize == 0)
414               {
415                 n_digits = digits_computed_so_far;
416                 break;
417               }
418           }
419
420         cy = mpn_mul_1 (rp, rp, rsize, big_base);
421         if (digits_computed_so_far == 0 && cy == 0)
422           {
423             abort ();
424             nexp += dig_per_u;
425             continue;
426           }
427         /* Convert N1 from BIG_BASE to a string of digits in BASE
428            using single precision operations.  */
429         {
430           unsigned char *s = tstr + digits_computed_so_far + dig_per_u;
431           for (i = dig_per_u - 1; i >= 0; i--)
432             {
433               *--s = cy % base;
434               cy /= base;
435             }
436         }
437         digits_computed_so_far += dig_per_u;
438       }
439     if (exp_in_base == 0)
440       exp_in_base = -nexp;
441   }
442
443  finish_up:
444
445   /* We can have at most one leading 0.  Remove it.  */
446   if (tstr[0] == 0)
447     {
448       tstr++;
449       digits_computed_so_far--;
450       exp_in_base--;
451     }
452
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)
456     {
457       /* Round the result.  */
458       if (tstr[n_digits] * 2 >= base)
459         {
460           digits_computed_so_far = n_digits;
461           for (i = n_digits - 1; i >= 0; i--)
462             {
463               unsigned int x;
464               x = ++(tstr[i]);
465               if ((int) x < base)
466                 goto rounded_ok;
467               digits_computed_so_far--;
468             }
469           tstr[0] = 1;
470           digits_computed_so_far = 1;
471           exp_in_base++;
472         rounded_ok:;
473         }
474     }
475
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;
481
482   /* Remove trailing 0.  There can be many zeros. */
483   while (n_digits != 0 && tstr[n_digits - 1] == 0)
484     n_digits--;
485
486   /* Translate to ascii and null-terminate.  */
487   for (i = 0; i < (int) n_digits; i++)
488     *str++ = num_to_text[tstr[i]];
489   *str = 0;
490   *exp = exp_in_base;
491   TMP_FREE (marker);
492   return digit_ptr;
493 }
494
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;
501 #endif