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