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)
66ssize_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
123extern mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
124				       int *expt, int *is_neg,
125				       double value);
126extern 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);
129extern unsigned int __guess_grouping (unsigned int intdig_max,
130				      const char *grouping);
131
132
133static 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
139int
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
208hack_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;
246hack_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;
899hack_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;
924hack_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;
940hack_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
1167unsigned 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
1206static wchar_t *
1207internal_function
1208group_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