xref: /haiku/src/system/libroot/posix/glibc/stdio-common/printf_fp.c (revision c237c4ce593ee823d9867fd997e51e4c447f5623)
1 /* Floating point output for `printf'.
2    Copyright (C) 1995-1999, 2000, 2001 Free Software Foundation, Inc.
3    This file is part of the GNU C Library.
4    Written by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
5 
6    The GNU C Library is free software; you can redistribute it and/or
7    modify it under the terms of the GNU Lesser General Public
8    License as published by the Free Software Foundation; either
9    version 2.1 of the License, or (at your option) any later version.
10 
11    The GNU C Library is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14    Lesser General Public License for more details.
15 
16    You should have received a copy of the GNU Lesser General Public
17    License along with the GNU C Library; if not, write to the Free
18    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
19    02111-1307 USA.  */
20 
21 /* The gmp headers need some configuration frobs.  */
22 #define HAVE_ALLOCA 1
23 
24 #ifdef USE_IN_LIBIO
25 #  include <libioP.h>
26 #else
27 #  include <stdio.h>
28 #endif
29 #include <alloca.h>
30 #include <ctype.h>
31 #include <float.h>
32 #include <gmp-mparam.h>
33 #include <gmp.h>
34 #include <stdlib/gmp-impl.h>
35 #include <stdlib/longlong.h>
36 #include <stdlib/fpioconst.h>
37 #include <locale/localeinfo.h>
38 #include <limits.h>
39 #include <math.h>
40 #include <printf.h>
41 #include <string.h>
42 #include <unistd.h>
43 #include <stdlib.h>
44 #include <wchar.h>
45 
46 #ifndef NDEBUG
47 # define NDEBUG			/* Undefine this for debugging assertions.  */
48 #endif
49 #include <assert.h>
50 
51 /* This defines make it possible to use the same code for GNU C library and
52    the GNU I/O library.	 */
53 #ifdef USE_IN_LIBIO
54 # define PUT(f, s, n) _IO_sputn (f, s, n)
55 # define PAD(f, c, n) (wide ? _IO_wpadn (f, c, n) : _IO_padn (f, c, n))
56 /* We use this file GNU C library and GNU I/O library.	So make
57    names equal.	 */
58 # undef putc
59 # define putc(c, f) (wide \
60 		      ? (int)_IO_putwc_unlocked (c, f) : _IO_putc_unlocked (c, f))
61 # define size_t     _IO_size_t
62 # define FILE	     _IO_FILE
63 #else	/* ! USE_IN_LIBIO */
64 # define PUT(f, s, n) fwrite (s, 1, n, f)
65 # define PAD(f, c, n) __printf_pad (f, c, n)
66 ssize_t __printf_pad __P ((FILE *, char pad, int n)); /* In vfprintf.c.  */
67 #endif	/* USE_IN_LIBIO */
68 
69 /* Macros for doing the actual output.  */
70 
71 #define outchar(ch)							      \
72   do									      \
73     {									      \
74       register const int outc = (ch);					      \
75       if (putc (outc, fp) == EOF)					      \
76 	return -1;							      \
77       ++done;								      \
78     } while (0)
79 
80 #define PRINT(ptr, wptr, len)						      \
81   do									      \
82     {									      \
83       register size_t outlen = (len);					      \
84       if (len > 20)							      \
85 	{								      \
86 	  if (PUT (fp, wide ? (const char *) wptr : ptr, outlen) != outlen)   \
87 	    return -1;							      \
88 	  ptr += outlen;						      \
89 	  done += outlen;						      \
90 	}								      \
91       else								      \
92 	{								      \
93 	  if (wide)							      \
94 	    while (outlen-- > 0)					      \
95 	      outchar (*wptr++);					      \
96 	  else								      \
97 	    while (outlen-- > 0)					      \
98 	      outchar (*ptr++);						      \
99 	}								      \
100     } while (0)
101 
102 #define PADN(ch, len)							      \
103   do									      \
104     {									      \
105       if (PAD (fp, ch, len) != len)					      \
106 	return -1;							      \
107       done += len;							      \
108     }									      \
109   while (0)
110 
111 /* We use the GNU MP library to handle large numbers.
112 
113    An MP variable occupies a varying number of entries in its array.  We keep
114    track of this number for efficiency reasons.  Otherwise we would always
115    have to process the whole array.  */
116 #define MPN_VAR(name) mp_limb_t *name; mp_size_t name##size
117 
118 #define MPN_ASSIGN(dst,src)						      \
119   memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
120 #define MPN_GE(u,v) \
121   (u##size > v##size || (u##size == v##size && __mpn_cmp (u, v, u##size) >= 0))
122 
123 extern mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
124 				       int *expt, int *is_neg,
125 				       double value);
126 extern mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
127 					    int *expt, int *is_neg,
128 					    long double value);
129 extern unsigned int __guess_grouping (unsigned int intdig_max,
130 				      const char *grouping);
131 
132 
133 static wchar_t *group_number (wchar_t *buf, wchar_t *bufend,
134 			      unsigned int intdig_no, const char *grouping,
135 			      wchar_t thousands_sep, int ngroups)
136      internal_function;
137 
138 
139 int
140 __printf_fp (FILE *fp,
141 	     const struct printf_info *info,
142 	     const void *const *args)
143 {
144   /* The floating-point value to output.  */
145   union
146     {
147       double dbl;
148       __long_double_t ldbl;
149     }
150   fpnum;
151 
152   /* Locale-dependent representation of decimal point.	*/
153   const char *decimal;
154   wchar_t decimalwc;
155 
156   /* Locale-dependent thousands separator and grouping specification.  */
157   const char *thousands_sep = NULL;
158   wchar_t thousands_sepwc = 0;
159   const char *grouping;
160 
161   /* "NaN" or "Inf" for the special cases.  */
162   const char *special = NULL;
163   const wchar_t *wspecial = NULL;
164 
165   /* We need just a few limbs for the input before shifting to the right
166      position.	*/
167   mp_limb_t fp_input[(LDBL_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB];
168   /* We need to shift the contents of fp_input by this amount of bits.	*/
169   int to_shift = 0;
170 
171   /* The fraction of the floting-point value in question  */
172   MPN_VAR(frac);
173   /* and the exponent.	*/
174   int exponent;
175   /* Sign of the exponent.  */
176   int expsign = 0;
177   /* Sign of float number.  */
178   int is_neg = 0;
179 
180   /* Scaling factor.  */
181   MPN_VAR(scale);
182 
183   /* Temporary bignum value.  */
184   MPN_VAR(tmp);
185 
186   /* Digit which is result of last hack_digit() call.  */
187   wchar_t digit;
188 
189   /* The type of output format that will be used: 'e'/'E' or 'f'/'F'.  */
190   int type;
191 
192   /* Counter for number of written characters.	*/
193   int done = 0;
194 
195   /* General helper (carry limb).  */
196   mp_limb_t cy;
197 
198   /* Nonzero if this is output on a wide character stream.  */
199   int wide = info->wide;
200 
201   wchar_t hack_digit_ret;
202   int hack_digit_callee;
203 
204   while (0)
205     {
206       mp_limb_t hi;
207 
208 hack_digit:
209       if (expsign != 0 && _tolower (type) == 'f' && exponent-- > 0)
210 	hi = 0;
211       else if (scalesize == 0)
212 	{
213 	  hi = frac[fracsize - 1];
214 	  cy = __mpn_mul_1 (frac, frac, fracsize - 1, 10);
215 	  frac[fracsize - 1] = cy;
216 	}
217       else
218 	{
219 	  if (fracsize < scalesize)
220 	    hi = 0;
221 	  else
222 	    {
223 	      hi = mpn_divmod (tmp, frac, fracsize, scale, scalesize);
224 	      tmp[fracsize - scalesize] = hi;
225 	      hi = tmp[0];
226 
227 	      fracsize = scalesize;
228 	      while (fracsize != 0 && frac[fracsize - 1] == 0)
229 		--fracsize;
230 	      if (fracsize == 0)
231 		{
232 		  /* We're not prepared for an mpn variable with zero
233 		     limbs.  */
234 		  fracsize = 1;
235 		  hack_digit_ret = L'0' + hi;
236 		  goto hack_digit_end;
237 		}
238 	    }
239 
240 	  cy = __mpn_mul_1 (frac, frac, fracsize, 10);
241 	  if (cy != 0)
242 	    frac[fracsize++] = cy;
243 	}
244 
245       hack_digit_ret = L'0' + hi;
246 hack_digit_end:
247       switch (hack_digit_callee)
248         {
249 	  case 1: goto hack_digit_callee1;
250 	  case 2: goto hack_digit_callee2;
251 	  case 3: goto hack_digit_callee3;
252 	  default: abort();
253 	}
254     }
255 
256 
257   /* Figure out the decimal point character.  */
258   if (info->extra == 0)
259     {
260       decimal = _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
261       decimalwc = _NL_CURRENT_WORD (LC_NUMERIC, _NL_NUMERIC_DECIMAL_POINT_WC);
262     }
263   else
264     {
265       decimal = _NL_CURRENT (LC_MONETARY, MON_DECIMAL_POINT);
266       if (*decimal == '\0')
267 	decimal = _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
268       decimalwc = _NL_CURRENT_WORD (LC_MONETARY,
269 				    _NL_MONETARY_DECIMAL_POINT_WC);
270       if (decimalwc == L'\0')
271 	decimalwc = _NL_CURRENT_WORD (LC_NUMERIC,
272 				      _NL_NUMERIC_DECIMAL_POINT_WC);
273     }
274   /* The decimal point character must not be zero.  */
275   assert (*decimal != '\0');
276   assert (decimalwc != L'\0');
277 
278   if (info->group)
279     {
280       if (info->extra == 0)
281 	grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
282       else
283 	grouping = _NL_CURRENT (LC_MONETARY, MON_GROUPING);
284 
285       if (*grouping <= 0 || *grouping == CHAR_MAX)
286 	grouping = NULL;
287       else
288 	{
289 	  /* Figure out the thousands separator character.  */
290 	  if (wide)
291 	    {
292 	      if (info->extra == 0)
293 		thousands_sepwc =
294 		  _NL_CURRENT_WORD (LC_NUMERIC, _NL_NUMERIC_THOUSANDS_SEP_WC);
295 	      else
296 		thousands_sepwc =
297 		  _NL_CURRENT_WORD (LC_MONETARY,
298 				    _NL_MONETARY_THOUSANDS_SEP_WC);
299 	    }
300 	  else
301 	    {
302 	      if (info->extra == 0)
303 		thousands_sep = _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP);
304 	      else
305 		thousands_sep = _NL_CURRENT (LC_MONETARY, MON_THOUSANDS_SEP);
306 	    }
307 
308 	  if ((wide && thousands_sepwc == L'\0')
309 	      || (! wide && *thousands_sep == '\0'))
310 	    grouping = NULL;
311 	  else if (thousands_sepwc == L'\0')
312 	    /* If we are printing multibyte characters and there is a
313 	       multibyte representation for the thousands separator,
314 	       we must ensure the wide character thousands separator
315 	       is available, even if it is fake.  */
316 	    thousands_sepwc = 0xfffffffe;
317 	}
318     }
319   else
320     grouping = NULL;
321 
322   /* Fetch the argument value.	*/
323 #ifndef __NO_LONG_DOUBLE_MATH
324   if (info->is_long_double && sizeof (long double) > sizeof (double))
325     {
326       fpnum.ldbl = *(const long double *) args[0];
327 
328       /* Check for special values: not a number or infinity.  */
329       if (isnan (fpnum.ldbl))
330 	{
331 	  if (isupper (info->spec))
332 	    {
333 	      special = "NAN";
334 	      wspecial = L"NAN";
335 	    }
336 	    else
337 	      {
338 		special = "nan";
339 		wspecial = L"nan";
340 	      }
341 	  is_neg = 0;
342 	}
343       else if (isinf (fpnum.ldbl))
344 	{
345 	  if (isupper (info->spec))
346 	    {
347 	      special = "INF";
348 	      wspecial = L"INF";
349 	    }
350 	  else
351 	    {
352 	      special = "inf";
353 	      wspecial = L"inf";
354 	    }
355 	  is_neg = fpnum.ldbl < 0;
356 	}
357       else
358 	{
359 	  fracsize = __mpn_extract_long_double (fp_input,
360 						(sizeof (fp_input) /
361 						 sizeof (fp_input[0])),
362 						&exponent, &is_neg,
363 						fpnum.ldbl);
364 	  to_shift = 1 + fracsize * BITS_PER_MP_LIMB - LDBL_MANT_DIG;
365 	}
366     }
367   else
368 #endif	/* no long double */
369     {
370       fpnum.dbl = *(const double *) args[0];
371 
372       /* Check for special values: not a number or infinity.  */
373       if (isnan (fpnum.dbl))
374 	{
375 	  if (isupper (info->spec))
376 	    {
377 	      special = "NAN";
378 	      wspecial = L"NAN";
379 	    }
380 	  else
381 	    {
382 	      special = "nan";
383 	      wspecial = L"nan";
384 	    }
385 	  is_neg = 0;
386 	}
387       else if (isinf (fpnum.dbl))
388 	{
389 	  if (isupper (info->spec))
390 	    {
391 	      special = "INF";
392 	      wspecial = L"INF";
393 	    }
394 	  else
395 	    {
396 	      special = "inf";
397 	      wspecial = L"inf";
398 	    }
399 	  is_neg = fpnum.dbl < 0;
400 	}
401       else
402 	{
403 	  fracsize = __mpn_extract_double (fp_input,
404 					   (sizeof (fp_input)
405 					    / sizeof (fp_input[0])),
406 					   &exponent, &is_neg, fpnum.dbl);
407 	  to_shift = 1 + fracsize * BITS_PER_MP_LIMB - DBL_MANT_DIG;
408 	}
409     }
410 
411   if (special)
412     {
413       int width = info->width;
414 
415       if (is_neg || info->showsign || info->space)
416 	--width;
417       width -= 3;
418 
419       if (!info->left && width > 0)
420 	PADN (' ', width);
421 
422       if (is_neg)
423 	outchar ('-');
424       else if (info->showsign)
425 	outchar ('+');
426       else if (info->space)
427 	outchar (' ');
428 
429       PRINT (special, wspecial, 3);
430 
431       if (info->left && width > 0)
432 	PADN (' ', width);
433 
434       return done;
435     }
436 
437 
438   /* We need three multiprecision variables.  Now that we have the exponent
439      of the number we can allocate the needed memory.  It would be more
440      efficient to use variables of the fixed maximum size but because this
441      would be really big it could lead to memory problems.  */
442   {
443     mp_size_t bignum_size = ((ABS (exponent) + BITS_PER_MP_LIMB - 1)
444 			     / BITS_PER_MP_LIMB + 4) * sizeof (mp_limb_t);
445     frac = (mp_limb_t *) alloca (bignum_size);
446     tmp = (mp_limb_t *) alloca (bignum_size);
447     scale = (mp_limb_t *) alloca (bignum_size);
448   }
449 
450   /* We now have to distinguish between numbers with positive and negative
451      exponents because the method used for the one is not applicable/efficient
452      for the other.  */
453   scalesize = 0;
454   if (exponent > 2)
455     {
456       /* |FP| >= 8.0.  */
457       int scaleexpo = 0;
458       int explog = LDBL_MAX_10_EXP_LOG;
459       int exp10 = 0;
460       const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
461       int cnt_h, cnt_l, i;
462 
463       if ((exponent + to_shift) % BITS_PER_MP_LIMB == 0)
464 	{
465 	  MPN_COPY_DECR (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
466 			 fp_input, fracsize);
467 	  fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
468 	}
469       else
470 	{
471 	  cy = __mpn_lshift (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
472 			     fp_input, fracsize,
473 			     (exponent + to_shift) % BITS_PER_MP_LIMB);
474 	  fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
475 	  if (cy)
476 	    frac[fracsize++] = cy;
477 	}
478       MPN_ZERO (frac, (exponent + to_shift) / BITS_PER_MP_LIMB);
479 
480       assert (powers > &_fpioconst_pow10[0]);
481       do
482 	{
483 	  --powers;
484 
485 	  /* The number of the product of two binary numbers with n and m
486 	     bits respectively has m+n or m+n-1 bits.	*/
487 	  if (exponent >= scaleexpo + powers->p_expo - 1)
488 	    {
489 	      if (scalesize == 0)
490 		{
491 #ifndef __NO_LONG_DOUBLE_MATH
492 		  if (LDBL_MANT_DIG > _FPIO_CONST_OFFSET * BITS_PER_MP_LIMB
493 		      && info->is_long_double)
494 		    {
495 #define _FPIO_CONST_SHIFT \
496   (((LDBL_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) \
497    - _FPIO_CONST_OFFSET)
498 		      /* 64bit const offset is not enough for
499 			 IEEE quad long double.  */
500 		      tmpsize = powers->arraysize + _FPIO_CONST_SHIFT;
501 		      memcpy (tmp + _FPIO_CONST_SHIFT,
502 			      &__tens[powers->arrayoff],
503 			      tmpsize * sizeof (mp_limb_t));
504 		      MPN_ZERO (tmp, _FPIO_CONST_SHIFT);
505 		    }
506 		  else
507 #endif
508 		    {
509 		      tmpsize = powers->arraysize;
510 		      memcpy (tmp, &__tens[powers->arrayoff],
511 			      tmpsize * sizeof (mp_limb_t));
512 		    }
513 		}
514 	      else
515 		{
516 		  cy = __mpn_mul (tmp, scale, scalesize,
517 				  &__tens[powers->arrayoff
518 					 + _FPIO_CONST_OFFSET],
519 				  powers->arraysize - _FPIO_CONST_OFFSET);
520 		  tmpsize = scalesize + powers->arraysize - _FPIO_CONST_OFFSET;
521 		  if (cy == 0)
522 		    --tmpsize;
523 		}
524 
525 	      if (MPN_GE (frac, tmp))
526 		{
527 		  int cnt;
528 		  MPN_ASSIGN (scale, tmp);
529 		  count_leading_zeros (cnt, scale[scalesize - 1]);
530 		  scaleexpo = (scalesize - 2) * BITS_PER_MP_LIMB - cnt - 1;
531 		  exp10 |= 1 << explog;
532 		}
533 	    }
534 	  --explog;
535 	}
536       while (powers > &_fpioconst_pow10[0]);
537       exponent = exp10;
538 
539       /* Optimize number representations.  We want to represent the numbers
540 	 with the lowest number of bytes possible without losing any
541 	 bytes. Also the highest bit in the scaling factor has to be set
542 	 (this is a requirement of the MPN division routines).  */
543       if (scalesize > 0)
544 	{
545 	  /* Determine minimum number of zero bits at the end of
546 	     both numbers.  */
547 	  for (i = 0; scale[i] == 0 && frac[i] == 0; i++)
548 	    ;
549 
550 	  /* Determine number of bits the scaling factor is misplaced.	*/
551 	  count_leading_zeros (cnt_h, scale[scalesize - 1]);
552 
553 	  if (cnt_h == 0)
554 	    {
555 	      /* The highest bit of the scaling factor is already set.	So
556 		 we only have to remove the trailing empty limbs.  */
557 	      if (i > 0)
558 		{
559 		  MPN_COPY_INCR (scale, scale + i, scalesize - i);
560 		  scalesize -= i;
561 		  MPN_COPY_INCR (frac, frac + i, fracsize - i);
562 		  fracsize -= i;
563 		}
564 	    }
565 	  else
566 	    {
567 	      if (scale[i] != 0)
568 		{
569 		  count_trailing_zeros (cnt_l, scale[i]);
570 		  if (frac[i] != 0)
571 		    {
572 		      int cnt_l2;
573 		      count_trailing_zeros (cnt_l2, frac[i]);
574 		      if (cnt_l2 < cnt_l)
575 			cnt_l = cnt_l2;
576 		    }
577 		}
578 	      else
579 		count_trailing_zeros (cnt_l, frac[i]);
580 
581 	      /* Now shift the numbers to their optimal position.  */
582 	      if (i == 0 && BITS_PER_MP_LIMB - cnt_h > cnt_l)
583 		{
584 		  /* We cannot save any memory.	 So just roll both numbers
585 		     so that the scaling factor has its highest bit set.  */
586 
587 		  (void) __mpn_lshift (scale, scale, scalesize, cnt_h);
588 		  cy = __mpn_lshift (frac, frac, fracsize, cnt_h);
589 		  if (cy != 0)
590 		    frac[fracsize++] = cy;
591 		}
592 	      else if (BITS_PER_MP_LIMB - cnt_h <= cnt_l)
593 		{
594 		  /* We can save memory by removing the trailing zero limbs
595 		     and by packing the non-zero limbs which gain another
596 		     free one. */
597 
598 		  (void) __mpn_rshift (scale, scale + i, scalesize - i,
599 				       BITS_PER_MP_LIMB - cnt_h);
600 		  scalesize -= i + 1;
601 		  (void) __mpn_rshift (frac, frac + i, fracsize - i,
602 				       BITS_PER_MP_LIMB - cnt_h);
603 		  fracsize -= frac[fracsize - i - 1] == 0 ? i + 1 : i;
604 		}
605 	      else
606 		{
607 		  /* We can only save the memory of the limbs which are zero.
608 		     The non-zero parts occupy the same number of limbs.  */
609 
610 		  (void) __mpn_rshift (scale, scale + (i - 1),
611 				       scalesize - (i - 1),
612 				       BITS_PER_MP_LIMB - cnt_h);
613 		  scalesize -= i;
614 		  (void) __mpn_rshift (frac, frac + (i - 1),
615 				       fracsize - (i - 1),
616 				       BITS_PER_MP_LIMB - cnt_h);
617 		  fracsize -= frac[fracsize - (i - 1) - 1] == 0 ? i : i - 1;
618 		}
619 	    }
620 	}
621     }
622   else if (exponent < 0)
623     {
624       /* |FP| < 1.0.  */
625       int exp10 = 0;
626       int explog = LDBL_MAX_10_EXP_LOG;
627       const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
628       mp_size_t used_limbs = fracsize - 1;
629 
630       /* Now shift the input value to its right place.	*/
631       cy = __mpn_lshift (frac, fp_input, fracsize, to_shift);
632       frac[fracsize++] = cy;
633       assert (cy == 1 || (frac[fracsize - 2] == 0 && frac[0] == 0));
634 
635       expsign = 1;
636       exponent = -exponent;
637 
638       assert (powers != &_fpioconst_pow10[0]);
639       do
640 	{
641 	  --powers;
642 
643 	  if (exponent >= powers->m_expo)
644 	    {
645 	      int i, incr, cnt_h, cnt_l;
646 	      mp_limb_t topval[2];
647 
648 	      /* The __mpn_mul function expects the first argument to be
649 		 bigger than the second.  */
650 	      if (fracsize < powers->arraysize - _FPIO_CONST_OFFSET)
651 		cy = __mpn_mul (tmp, &__tens[powers->arrayoff
652 					    + _FPIO_CONST_OFFSET],
653 				powers->arraysize - _FPIO_CONST_OFFSET,
654 				frac, fracsize);
655 	      else
656 		cy = __mpn_mul (tmp, frac, fracsize,
657 				&__tens[powers->arrayoff + _FPIO_CONST_OFFSET],
658 				powers->arraysize - _FPIO_CONST_OFFSET);
659 	      tmpsize = fracsize + powers->arraysize - _FPIO_CONST_OFFSET;
660 	      if (cy == 0)
661 		--tmpsize;
662 
663 	      count_leading_zeros (cnt_h, tmp[tmpsize - 1]);
664 	      incr = (tmpsize - fracsize) * BITS_PER_MP_LIMB
665 		     + BITS_PER_MP_LIMB - 1 - cnt_h;
666 
667 	      assert (incr <= powers->p_expo);
668 
669 	      /* If we increased the exponent by exactly 3 we have to test
670 		 for overflow.	This is done by comparing with 10 shifted
671 		 to the right position.	 */
672 	      if (incr == exponent + 3)
673 		{
674 		  if (cnt_h <= BITS_PER_MP_LIMB - 4)
675 		    {
676 		      topval[0] = 0;
677 		      topval[1]
678 			= ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4 - cnt_h);
679 		    }
680 		  else
681 		    {
682 		      topval[0] = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4);
683 		      topval[1] = 0;
684 		      (void) __mpn_lshift (topval, topval, 2,
685 					   BITS_PER_MP_LIMB - cnt_h);
686 		    }
687 		}
688 
689 	      /* We have to be careful when multiplying the last factor.
690 		 If the result is greater than 1.0 be have to test it
691 		 against 10.0.  If it is greater or equal to 10.0 the
692 		 multiplication was not valid.  This is because we cannot
693 		 determine the number of bits in the result in advance.  */
694 	      if (incr < exponent + 3
695 		  || (incr == exponent + 3 &&
696 		      (tmp[tmpsize - 1] < topval[1]
697 		       || (tmp[tmpsize - 1] == topval[1]
698 			   && tmp[tmpsize - 2] < topval[0]))))
699 		{
700 		  /* The factor is right.  Adapt binary and decimal
701 		     exponents.	 */
702 		  exponent -= incr;
703 		  exp10 |= 1 << explog;
704 
705 		  /* If this factor yields a number greater or equal to
706 		     1.0, we must not shift the non-fractional digits down. */
707 		  if (exponent < 0)
708 		    cnt_h += -exponent;
709 
710 		  /* Now we optimize the number representation.	 */
711 		  for (i = 0; tmp[i] == 0; ++i);
712 		  if (cnt_h == BITS_PER_MP_LIMB - 1)
713 		    {
714 		      MPN_COPY (frac, tmp + i, tmpsize - i);
715 		      fracsize = tmpsize - i;
716 		    }
717 		  else
718 		    {
719 		      count_trailing_zeros (cnt_l, tmp[i]);
720 
721 		      /* Now shift the numbers to their optimal position.  */
722 		      if (i == 0 && BITS_PER_MP_LIMB - 1 - cnt_h > cnt_l)
723 			{
724 			  /* We cannot save any memory.	 Just roll the
725 			     number so that the leading digit is in a
726 			     separate limb.  */
727 
728 			  cy = __mpn_lshift (frac, tmp, tmpsize, cnt_h + 1);
729 			  fracsize = tmpsize + 1;
730 			  frac[fracsize - 1] = cy;
731 			}
732 		      else if (BITS_PER_MP_LIMB - 1 - cnt_h <= cnt_l)
733 			{
734 			  (void) __mpn_rshift (frac, tmp + i, tmpsize - i,
735 					       BITS_PER_MP_LIMB - 1 - cnt_h);
736 			  fracsize = tmpsize - i;
737 			}
738 		      else
739 			{
740 			  /* We can only save the memory of the limbs which
741 			     are zero.	The non-zero parts occupy the same
742 			     number of limbs.  */
743 
744 			  (void) __mpn_rshift (frac, tmp + (i - 1),
745 					       tmpsize - (i - 1),
746 					       BITS_PER_MP_LIMB - 1 - cnt_h);
747 			  fracsize = tmpsize - (i - 1);
748 			}
749 		    }
750 		  used_limbs = fracsize - 1;
751 		}
752 	    }
753 	  --explog;
754 	}
755       while (powers != &_fpioconst_pow10[1] && exponent > 0);
756       /* All factors but 10^-1 are tested now.	*/
757       if (exponent > 0)
758 	{
759 	  int cnt_l;
760 
761 	  cy = __mpn_mul_1 (tmp, frac, fracsize, 10);
762 	  tmpsize = fracsize;
763 	  assert (cy == 0 || tmp[tmpsize - 1] < 20);
764 
765 	  count_trailing_zeros (cnt_l, tmp[0]);
766 	  if (cnt_l < MIN (4, exponent))
767 	    {
768 	      cy = __mpn_lshift (frac, tmp, tmpsize,
769 				 BITS_PER_MP_LIMB - MIN (4, exponent));
770 	      if (cy != 0)
771 		frac[tmpsize++] = cy;
772 	    }
773 	  else
774 	    (void) __mpn_rshift (frac, tmp, tmpsize, MIN (4, exponent));
775 	  fracsize = tmpsize;
776 	  exp10 |= 1;
777 	  assert (frac[fracsize - 1] < 10);
778 	}
779       exponent = exp10;
780     }
781   else
782     {
783       /* This is a special case.  We don't need a factor because the
784 	 numbers are in the range of 0.0 <= fp < 8.0.  We simply
785 	 shift it to the right place and divide it by 1.0 to get the
786 	 leading digit.	 (Of course this division is not really made.)	*/
787       assert (0 <= exponent && exponent < 3 &&
788 	      exponent + to_shift < BITS_PER_MP_LIMB);
789 
790       /* Now shift the input value to its right place.	*/
791       cy = __mpn_lshift (frac, fp_input, fracsize, (exponent + to_shift));
792       frac[fracsize++] = cy;
793       exponent = 0;
794     }
795 
796   {
797     int width = info->width;
798     wchar_t *wbuffer, *wstartp, *wcp;
799     int buffer_malloced;
800     int chars_needed;
801     int expscale;
802     int intdig_max, intdig_no = 0;
803     int fracdig_min, fracdig_max, fracdig_no = 0;
804     int dig_max;
805     int significant;
806     int ngroups = 0;
807 
808     if (_tolower (info->spec) == 'e')
809       {
810 	type = info->spec;
811 	intdig_max = 1;
812 	fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
813 	chars_needed = 1 + 1 + fracdig_max + 1 + 1 + 4;
814 	/*	       d   .	 ddd	     e	 +-  ddd  */
815 	dig_max = INT_MAX;		/* Unlimited.  */
816 	significant = 1;		/* Does not matter here.  */
817       }
818     else if (_tolower (info->spec) == 'f')
819       {
820 	type = info->spec;
821 	fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
822 	if (expsign == 0)
823 	  {
824 	    intdig_max = exponent + 1;
825 	    /* This can be really big!	*/  /* XXX Maybe malloc if too big? */
826 	    chars_needed = exponent + 1 + 1 + fracdig_max;
827 	  }
828 	else
829 	  {
830 	    intdig_max = 1;
831 	    chars_needed = 1 + 1 + fracdig_max;
832 	  }
833 	dig_max = INT_MAX;		/* Unlimited.  */
834 	significant = 1;		/* Does not matter here.  */
835       }
836     else
837       {
838 	dig_max = info->prec < 0 ? 6 : (info->prec == 0 ? 1 : info->prec);
839 	if ((expsign == 0 && exponent >= dig_max)
840 	    || (expsign != 0 && exponent > 4))
841 	  {
842 	    if ('g' - 'G' == 'e' - 'E')
843 	      type = 'E' + (info->spec - 'G');
844 	    else
845 	      type = isupper (info->spec) ? 'E' : 'e';
846 	    fracdig_max = dig_max - 1;
847 	    intdig_max = 1;
848 	    chars_needed = 1 + 1 + fracdig_max + 1 + 1 + 4;
849 	  }
850 	else
851 	  {
852 	    type = 'f';
853 	    intdig_max = expsign == 0 ? exponent + 1 : 0;
854 	    fracdig_max = dig_max - intdig_max;
855 	    /* We need space for the significant digits and perhaps
856 	       for leading zeros when < 1.0.  The number of leading
857 	       zeros can be as many as would be required for
858 	       exponential notation with a negative two-digit
859 	       exponent, which is 4.  */
860 	    chars_needed = dig_max + 1 + 4;
861 	  }
862 	fracdig_min = info->alt ? fracdig_max : 0;
863 	significant = 0;		/* We count significant digits.	 */
864       }
865 
866     if (grouping)
867       {
868 	/* Guess the number of groups we will make, and thus how
869 	   many spaces we need for separator characters.  */
870 	ngroups = __guess_grouping (intdig_max, grouping);
871 	chars_needed += ngroups;
872       }
873 
874     /* Allocate buffer for output.  We need two more because while rounding
875        it is possible that we need two more characters in front of all the
876        other output.  If the amount of memory we have to allocate is too
877        large use `malloc' instead of `alloca'.  */
878     buffer_malloced = chars_needed > 5000;
879     if (buffer_malloced)
880       {
881 	wbuffer = (wchar_t *) malloc ((2 + chars_needed) * sizeof (wchar_t));
882 	if (wbuffer == NULL)
883 	  /* Signal an error to the caller.  */
884 	  return -1;
885       }
886     else
887       wbuffer = (wchar_t *) alloca ((2 + chars_needed) * sizeof (wchar_t));
888     wcp = wstartp = wbuffer + 2;	/* Let room for rounding.  */
889 
890     /* Do the real work: put digits in allocated buffer.  */
891     if (expsign == 0 || _tolower (type) != 'f')
892       {
893 	assert (expsign == 0 || intdig_max == 1);
894 	while (intdig_no < intdig_max)
895 	  {
896 	    ++intdig_no;
897 	    hack_digit_callee = 1;
898 	    goto hack_digit;
899 hack_digit_callee1:
900 	    *wcp++ = hack_digit_ret;
901 	  }
902 	significant = 1;
903 	if (info->alt
904 	    || fracdig_min > 0
905 	    || (fracdig_max > 0 && (fracsize > 1 || frac[0] != 0)))
906 	  *wcp++ = decimalwc;
907       }
908     else
909       {
910 	/* |fp| < 1.0 and the selected type is 'f', so put "0."
911 	   in the buffer.  */
912 	*wcp++ = L'0';
913 	--exponent;
914 	*wcp++ = decimalwc;
915       }
916 
917     /* Generate the needed number of fractional digits.	 */
918     while (fracdig_no < fracdig_min
919 	   || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0)))
920       {
921 	++fracdig_no;
922 	hack_digit_callee = 2;
923 	goto hack_digit;
924 hack_digit_callee2:
925 	*wcp = hack_digit_ret;
926 	if (*wcp != L'0')
927 	  significant = 1;
928 	else if (significant == 0)
929 	  {
930 	    ++fracdig_max;
931 	    if (fracdig_min > 0)
932 	      ++fracdig_min;
933 	  }
934 	++wcp;
935       }
936 
937     /* Do rounding.  */
938     hack_digit_callee = 3;
939     goto hack_digit;
940 hack_digit_callee3:
941     digit = hack_digit_ret;
942     if (digit > L'4')
943       {
944 	wchar_t *wtp = wcp;
945 
946 	if (digit == L'5'
947 	    && ((*(wcp - 1) != decimalwc && (*(wcp - 1) & 1) == 0)
948 		|| ((*(wcp - 1) == decimalwc && (*(wcp - 2) & 1) == 0))))
949 	  {
950 	    /* This is the critical case.	 */
951 	    if (fracsize == 1 && frac[0] == 0)
952 	      /* Rest of the number is zero -> round to even.
953 		 (IEEE 754-1985 4.1 says this is the default rounding.)  */
954 	      goto do_expo;
955 	    else if (scalesize == 0)
956 	      {
957 		/* Here we have to see whether all limbs are zero since no
958 		   normalization happened.  */
959 		size_t lcnt = fracsize;
960 		while (lcnt >= 1 && frac[lcnt - 1] == 0)
961 		  --lcnt;
962 		if (lcnt == 0)
963 		  /* Rest of the number is zero -> round to even.
964 		     (IEEE 754-1985 4.1 says this is the default rounding.)  */
965 		  goto do_expo;
966 	      }
967 	  }
968 
969 	if (fracdig_no > 0)
970 	  {
971 	    /* Process fractional digits.  Terminate if not rounded or
972 	       radix character is reached.  */
973 	    while (*--wtp != decimalwc && *wtp == L'9')
974 	      *wtp = '0';
975 	    if (*wtp != decimalwc)
976 	      /* Round up.  */
977 	      (*wtp)++;
978 	  }
979 
980 	if (fracdig_no == 0 || *wtp == decimalwc)
981 	  {
982 	    /* Round the integer digits.  */
983 	    if (*(wtp - 1) == decimalwc)
984 	      --wtp;
985 
986 	    while (--wtp >= wstartp && *wtp == L'9')
987 	      *wtp = L'0';
988 
989 	    if (wtp >= wstartp)
990 	      /* Round up.  */
991 	      (*wtp)++;
992 	    else
993 	      /* It is more critical.  All digits were 9's.  */
994 	      {
995 		if (_tolower (type) != 'f')
996 		  {
997 		    *wstartp = '1';
998 		    exponent += expsign == 0 ? 1 : -1;
999 		  }
1000 		else if (intdig_no == dig_max)
1001 		  {
1002 		    /* This is the case where for type %g the number fits
1003 		       really in the range for %f output but after rounding
1004 		       the number of digits is too big.	 */
1005 		    *--wstartp = decimalwc;
1006 		    *--wstartp = L'1';
1007 
1008 		    if (info->alt || fracdig_no > 0)
1009 		      {
1010 			/* Overwrite the old radix character.  */
1011 			wstartp[intdig_no + 2] = L'0';
1012 			++fracdig_no;
1013 		      }
1014 
1015 		    fracdig_no += intdig_no;
1016 		    intdig_no = 1;
1017 		    fracdig_max = intdig_max - intdig_no;
1018 		    ++exponent;
1019 		    /* Now we must print the exponent.	*/
1020 		    type = isupper (info->spec) ? 'E' : 'e';
1021 		  }
1022 		else
1023 		  {
1024 		    /* We can simply add another another digit before the
1025 		       radix.  */
1026 		    *--wstartp = L'1';
1027 		    ++intdig_no;
1028 		  }
1029 
1030 		/* While rounding the number of digits can change.
1031 		   If the number now exceeds the limits remove some
1032 		   fractional digits.  */
1033 		if (intdig_no + fracdig_no > dig_max)
1034 		  {
1035 		    wcp -= intdig_no + fracdig_no - dig_max;
1036 		    fracdig_no -= intdig_no + fracdig_no - dig_max;
1037 		  }
1038 	      }
1039 	  }
1040       }
1041 
1042   do_expo:
1043     /* Now remove unnecessary '0' at the end of the string.  */
1044     while (fracdig_no > fracdig_min && *(wcp - 1) == L'0')
1045       {
1046 	--wcp;
1047 	--fracdig_no;
1048       }
1049     /* If we eliminate all fractional digits we perhaps also can remove
1050        the radix character.  */
1051     if (fracdig_no == 0 && !info->alt && *(wcp - 1) == decimalwc)
1052       --wcp;
1053 
1054     if (grouping)
1055       /* Add in separator characters, overwriting the same buffer.  */
1056       wcp = group_number (wstartp, wcp, intdig_no, grouping, thousands_sepwc,
1057 			  ngroups);
1058 
1059     /* Write the exponent if it is needed.  */
1060     if (_tolower (type) != 'f')
1061       {
1062 	*wcp++ = (wchar_t) type;
1063 	*wcp++ = expsign ? L'-' : L'+';
1064 
1065 	/* Find the magnitude of the exponent.	*/
1066 	expscale = 10;
1067 	while (expscale <= exponent)
1068 	  expscale *= 10;
1069 
1070 	if (exponent < 10)
1071 	  /* Exponent always has at least two digits.  */
1072 	  *wcp++ = L'0';
1073 	else
1074 	  do
1075 	    {
1076 	      expscale /= 10;
1077 	      *wcp++ = L'0' + (exponent / expscale);
1078 	      exponent %= expscale;
1079 	    }
1080 	  while (expscale > 10);
1081 	*wcp++ = L'0' + exponent;
1082       }
1083 
1084     /* Compute number of characters which must be filled with the padding
1085        character.  */
1086     if (is_neg || info->showsign || info->space)
1087       --width;
1088     width -= wcp - wstartp;
1089 
1090     if (!info->left && info->pad != '0' && width > 0)
1091       PADN (info->pad, width);
1092 
1093     if (is_neg)
1094       outchar ('-');
1095     else if (info->showsign)
1096       outchar ('+');
1097     else if (info->space)
1098       outchar (' ');
1099 
1100     if (!info->left && info->pad == '0' && width > 0)
1101       PADN ('0', width);
1102 
1103     {
1104       char *buffer = NULL;
1105       char *cp = NULL;
1106       char *tmpptr;
1107 
1108       if (! wide)
1109 	{
1110 	  /* Create the single byte string.  */
1111 	  size_t decimal_len;
1112 	  size_t thousands_sep_len;
1113 	  wchar_t *copywc;
1114 
1115 	  decimal_len = strlen (decimal);
1116 
1117 	  if (thousands_sep == NULL)
1118 	    thousands_sep_len = 0;
1119 	  else
1120 	    thousands_sep_len = strlen (thousands_sep);
1121 
1122 	  if (buffer_malloced)
1123 	    {
1124 	      buffer = (char *) malloc (2 + chars_needed + decimal_len
1125 					+ ngroups * thousands_sep_len);
1126 	      if (buffer == NULL)
1127 		/* Signal an error to the caller.  */
1128 		return -1;
1129 	    }
1130 	  else
1131 	    buffer = (char *) alloca (2 + chars_needed + decimal_len
1132 				      + ngroups * thousands_sep_len);
1133 
1134 	  /* Now copy the wide character string.  Since the character
1135 	     (except for the decimal point and thousands separator) must
1136 	     be coming from the ASCII range we can esily convert the
1137 	     string without mapping tables.  */
1138 	  for (cp = buffer, copywc = wstartp; copywc < wcp; ++copywc)
1139 	    if (*copywc == decimalwc)
1140 	      cp = (char *) __mempcpy (cp, decimal, decimal_len);
1141 	    else if (*copywc == thousands_sepwc)
1142 	      cp = (char *) __mempcpy (cp, thousands_sep, thousands_sep_len);
1143 	    else
1144 	      *cp++ = (char) *copywc;
1145 	}
1146 
1147       tmpptr = buffer;
1148       PRINT (tmpptr, wstartp, wide ? wcp - wstartp : cp - tmpptr);
1149 
1150       /* Free the memory if necessary.  */
1151       if (buffer_malloced)
1152 	{
1153 	  free (buffer);
1154 	  free (wbuffer);
1155 	}
1156     }
1157 
1158     if (info->left && width > 0)
1159       PADN (info->pad, width);
1160   }
1161   return done;
1162 }
1163 
1164 /* Return the number of extra grouping characters that will be inserted
1165    into a number with INTDIG_MAX integer digits.  */
1166 
1167 unsigned int
1168 __guess_grouping (unsigned int intdig_max, const char *grouping)
1169 {
1170   unsigned int groups;
1171 
1172   /* We treat all negative values like CHAR_MAX.  */
1173 
1174   if (*grouping == CHAR_MAX || *grouping <= 0)
1175     /* No grouping should be done.  */
1176     return 0;
1177 
1178   groups = 0;
1179   while (intdig_max > (unsigned int) *grouping)
1180     {
1181       ++groups;
1182       intdig_max -= *grouping++;
1183 
1184       if (*grouping == CHAR_MAX
1185 #if CHAR_MIN < 0
1186 	  || *grouping < 0
1187 #endif
1188 	  )
1189 	/* No more grouping should be done.  */
1190 	break;
1191       else if (*grouping == 0)
1192 	{
1193 	  /* Same grouping repeats.  */
1194 	  groups += (intdig_max - 1) / grouping[-1];
1195 	  break;
1196 	}
1197     }
1198 
1199   return groups;
1200 }
1201 
1202 /* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND).
1203    There is guaranteed enough space past BUFEND to extend it.
1204    Return the new end of buffer.  */
1205 
1206 static wchar_t *
1207 internal_function
1208 group_number (wchar_t *buf, wchar_t *bufend, unsigned int intdig_no,
1209 	      const char *grouping, wchar_t thousands_sep, int ngroups)
1210 {
1211   wchar_t *p;
1212 
1213   if (ngroups == 0)
1214     return bufend;
1215 
1216   /* Move the fractional part down.  */
1217   __wmemmove (buf + intdig_no + ngroups, buf + intdig_no,
1218 	      bufend - (buf + intdig_no));
1219 
1220   p = buf + intdig_no + ngroups - 1;
1221   do
1222     {
1223       unsigned int len = *grouping++;
1224       do
1225 	*p-- = buf[--intdig_no];
1226       while (--len > 0);
1227       *p-- = thousands_sep;
1228 
1229       if (*grouping == CHAR_MAX
1230 #if CHAR_MIN < 0
1231 	  || *grouping < 0
1232 #endif
1233 	  )
1234 	/* No more grouping should be done.  */
1235 	break;
1236       else if (*grouping == 0)
1237 	/* Same grouping repeats.  */
1238 	--grouping;
1239     } while (intdig_no > (unsigned int) *grouping);
1240 
1241   /* Copy the remaining ungrouped digits.  */
1242   do
1243     *p-- = buf[--intdig_no];
1244   while (p > buf);
1245 
1246   return bufend + ngroups;
1247 }
1248