Update gnulib files.
[gnutls:gnutls.git] / lgl / vasnprintf.c
1 /* vsprintf with automatic memory allocation.
2    Copyright (C) 1999, 2002-2008 Free Software Foundation, Inc.
3
4    This program is free software; you can redistribute it and/or modify
5    it under the terms of the GNU Lesser General Public License as published by
6    the Free Software Foundation; either version 2.1, or (at your option)
7    any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU Lesser General Public License for more details.
13
14    You should have received a copy of the GNU Lesser General Public License along
15    with this program; if not, write to the Free Software Foundation,
16    Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
17
18 /* This file can be parametrized with the following macros:
19      VASNPRINTF         The name of the function being defined.
20      FCHAR_T            The element type of the format string.
21      DCHAR_T            The element type of the destination (result) string.
22      FCHAR_T_ONLY_ASCII Set to 1 to enable verification that all characters
23                         in the format string are ASCII. MUST be set if
24                         FCHAR_T and DCHAR_T are not the same type.
25      DIRECTIVE          Structure denoting a format directive.
26                         Depends on FCHAR_T.
27      DIRECTIVES         Structure denoting the set of format directives of a
28                         format string.  Depends on FCHAR_T.
29      PRINTF_PARSE       Function that parses a format string.
30                         Depends on FCHAR_T.
31      DCHAR_CPY          memcpy like function for DCHAR_T[] arrays.
32      DCHAR_SET          memset like function for DCHAR_T[] arrays.
33      DCHAR_MBSNLEN      mbsnlen like function for DCHAR_T[] arrays.
34      SNPRINTF           The system's snprintf (or similar) function.
35                         This may be either snprintf or swprintf.
36      TCHAR_T            The element type of the argument and result string
37                         of the said SNPRINTF function.  This may be either
38                         char or wchar_t.  The code exploits that
39                         sizeof (TCHAR_T) | sizeof (DCHAR_T) and
40                         alignof (TCHAR_T) <= alignof (DCHAR_T).
41      DCHAR_IS_TCHAR     Set to 1 if DCHAR_T and TCHAR_T are the same type.
42      DCHAR_CONV_FROM_ENCODING A function to convert from char[] to DCHAR[].
43      DCHAR_IS_UINT8_T   Set to 1 if DCHAR_T is uint8_t.
44      DCHAR_IS_UINT16_T  Set to 1 if DCHAR_T is uint16_t.
45      DCHAR_IS_UINT32_T  Set to 1 if DCHAR_T is uint32_t.  */
46
47 /* Tell glibc's <stdio.h> to provide a prototype for snprintf().
48    This must come before <config.h> because <config.h> may include
49    <features.h>, and once <features.h> has been included, it's too late.  */
50 #ifndef _GNU_SOURCE
51 # define _GNU_SOURCE    1
52 #endif
53
54 #ifndef VASNPRINTF
55 # include <config.h>
56 #endif
57 #ifndef IN_LIBINTL
58 # include <alloca.h>
59 #endif
60
61 /* Specification.  */
62 #ifndef VASNPRINTF
63 # if WIDE_CHAR_VERSION
64 #  include "vasnwprintf.h"
65 # else
66 #  include "vasnprintf.h"
67 # endif
68 #endif
69
70 #include <locale.h>     /* localeconv() */
71 #include <stdio.h>      /* snprintf(), sprintf() */
72 #include <stdlib.h>     /* abort(), malloc(), realloc(), free() */
73 #include <string.h>     /* memcpy(), strlen() */
74 #include <errno.h>      /* errno */
75 #include <limits.h>     /* CHAR_BIT */
76 #include <float.h>      /* DBL_MAX_EXP, LDBL_MAX_EXP */
77 #if HAVE_NL_LANGINFO
78 # include <langinfo.h>
79 #endif
80 #ifndef VASNPRINTF
81 # if WIDE_CHAR_VERSION
82 #  include "wprintf-parse.h"
83 # else
84 #  include "printf-parse.h"
85 # endif
86 #endif
87
88 /* Checked size_t computations.  */
89 #include "xsize.h"
90
91 #if (NEED_PRINTF_DOUBLE || NEED_PRINTF_LONG_DOUBLE) && !defined IN_LIBINTL
92 # include <math.h>
93 # include "float+.h"
94 #endif
95
96 #if (NEED_PRINTF_DOUBLE || NEED_PRINTF_INFINITE_DOUBLE) && !defined IN_LIBINTL
97 # include <math.h>
98 # include "isnand.h"
99 #endif
100
101 #if (NEED_PRINTF_LONG_DOUBLE || NEED_PRINTF_INFINITE_LONG_DOUBLE) && !defined IN_LIBINTL
102 # include <math.h>
103 # include "isnanl-nolibm.h"
104 # include "fpucw.h"
105 #endif
106
107 #if (NEED_PRINTF_DIRECTIVE_A || NEED_PRINTF_DOUBLE) && !defined IN_LIBINTL
108 # include <math.h>
109 # include "isnand.h"
110 # include "printf-frexp.h"
111 #endif
112
113 #if (NEED_PRINTF_DIRECTIVE_A || NEED_PRINTF_LONG_DOUBLE) && !defined IN_LIBINTL
114 # include <math.h>
115 # include "isnanl-nolibm.h"
116 # include "printf-frexpl.h"
117 # include "fpucw.h"
118 #endif
119
120 #if HAVE_WCHAR_T
121 # if HAVE_WCSLEN
122 #  define local_wcslen wcslen
123 # else
124    /* Solaris 2.5.1 has wcslen() in a separate library libw.so. To avoid
125       a dependency towards this library, here is a local substitute.
126       Define this substitute only once, even if this file is included
127       twice in the same compilation unit.  */
128 #  ifndef local_wcslen_defined
129 #   define local_wcslen_defined 1
130 static size_t
131 local_wcslen (const wchar_t *s)
132 {
133   const wchar_t *ptr;
134
135   for (ptr = s; *ptr != (wchar_t) 0; ptr++)
136     ;
137   return ptr - s;
138 }
139 #  endif
140 # endif
141 #endif
142
143 /* Default parameters.  */
144 #ifndef VASNPRINTF
145 # if WIDE_CHAR_VERSION
146 #  define VASNPRINTF vasnwprintf
147 #  define FCHAR_T wchar_t
148 #  define DCHAR_T wchar_t
149 #  define TCHAR_T wchar_t
150 #  define DCHAR_IS_TCHAR 1
151 #  define DIRECTIVE wchar_t_directive
152 #  define DIRECTIVES wchar_t_directives
153 #  define PRINTF_PARSE wprintf_parse
154 #  define DCHAR_CPY wmemcpy
155 # else
156 #  define VASNPRINTF vasnprintf
157 #  define FCHAR_T char
158 #  define DCHAR_T char
159 #  define TCHAR_T char
160 #  define DCHAR_IS_TCHAR 1
161 #  define DIRECTIVE char_directive
162 #  define DIRECTIVES char_directives
163 #  define PRINTF_PARSE printf_parse
164 #  define DCHAR_CPY memcpy
165 # endif
166 #endif
167 #if WIDE_CHAR_VERSION
168   /* TCHAR_T is wchar_t.  */
169 # define USE_SNPRINTF 1
170 # if HAVE_DECL__SNWPRINTF
171    /* On Windows, the function swprintf() has a different signature than
172       on Unix; we use the _snwprintf() function instead.  */
173 #  define SNPRINTF _snwprintf
174 # else
175    /* Unix.  */
176 #  define SNPRINTF swprintf
177 # endif
178 #else
179   /* TCHAR_T is char.  */
180 # /* Use snprintf if it exists under the name 'snprintf' or '_snprintf'.
181      But don't use it on BeOS, since BeOS snprintf produces no output if the
182      size argument is >= 0x3000000.  */
183 # if (HAVE_DECL__SNPRINTF || HAVE_SNPRINTF) && !defined __BEOS__
184 #  define USE_SNPRINTF 1
185 # else
186 #  define USE_SNPRINTF 0
187 # endif
188 # if HAVE_DECL__SNPRINTF
189    /* Windows.  */
190 #  define SNPRINTF _snprintf
191 # else
192    /* Unix.  */
193 #  define SNPRINTF snprintf
194    /* Here we need to call the native snprintf, not rpl_snprintf.  */
195 #  undef snprintf
196 # endif
197 #endif
198 /* Here we need to call the native sprintf, not rpl_sprintf.  */
199 #undef sprintf
200
201 #if (NEED_PRINTF_DIRECTIVE_A || NEED_PRINTF_LONG_DOUBLE || NEED_PRINTF_DOUBLE || NEED_PRINTF_INFINITE_DOUBLE) && !defined IN_LIBINTL
202 /* Determine the decimal-point character according to the current locale.  */
203 # ifndef decimal_point_char_defined
204 #  define decimal_point_char_defined 1
205 static char
206 decimal_point_char ()
207 {
208   const char *point;
209   /* Determine it in a multithread-safe way.  We know nl_langinfo is
210      multithread-safe on glibc systems, but is not required to be multithread-
211      safe by POSIX.  sprintf(), however, is multithread-safe.  localeconv()
212      is rarely multithread-safe.  */
213 #  if HAVE_NL_LANGINFO && __GLIBC__
214   point = nl_langinfo (RADIXCHAR);
215 #  elif 1
216   char pointbuf[5];
217   sprintf (pointbuf, "%#.0f", 1.0);
218   point = &pointbuf[1];
219 #  else
220   point = localeconv () -> decimal_point;
221 #  endif
222   /* The decimal point is always a single byte: either '.' or ','.  */
223   return (point[0] != '\0' ? point[0] : '.');
224 }
225 # endif
226 #endif
227
228 #if NEED_PRINTF_INFINITE_DOUBLE && !NEED_PRINTF_DOUBLE && !defined IN_LIBINTL
229
230 /* Equivalent to !isfinite(x) || x == 0, but does not require libm.  */
231 static int
232 is_infinite_or_zero (double x)
233 {
234   return isnand (x) || x + x == x;
235 }
236
237 #endif
238
239 #if NEED_PRINTF_INFINITE_LONG_DOUBLE && !NEED_PRINTF_LONG_DOUBLE && !defined IN_LIBINTL
240
241 /* Equivalent to !isfinite(x), but does not require libm.  */
242 static int
243 is_infinitel (long double x)
244 {
245   return isnanl (x) || (x + x == x && x != 0.0L);
246 }
247
248 #endif
249
250 #if (NEED_PRINTF_LONG_DOUBLE || NEED_PRINTF_DOUBLE) && !defined IN_LIBINTL
251
252 /* Converting 'long double' to decimal without rare rounding bugs requires
253    real bignums.  We use the naming conventions of GNU gmp, but vastly simpler
254    (and slower) algorithms.  */
255
256 typedef unsigned int mp_limb_t;
257 # define GMP_LIMB_BITS 32
258 typedef int mp_limb_verify[2 * (sizeof (mp_limb_t) * CHAR_BIT == GMP_LIMB_BITS) - 1];
259
260 typedef unsigned long long mp_twolimb_t;
261 # define GMP_TWOLIMB_BITS 64
262 typedef int mp_twolimb_verify[2 * (sizeof (mp_twolimb_t) * CHAR_BIT == GMP_TWOLIMB_BITS) - 1];
263
264 /* Representation of a bignum >= 0.  */
265 typedef struct
266 {
267   size_t nlimbs;
268   mp_limb_t *limbs; /* Bits in little-endian order, allocated with malloc().  */
269 } mpn_t;
270
271 /* Compute the product of two bignums >= 0.
272    Return the allocated memory in case of success, NULL in case of memory
273    allocation failure.  */
274 static void *
275 multiply (mpn_t src1, mpn_t src2, mpn_t *dest)
276 {
277   const mp_limb_t *p1;
278   const mp_limb_t *p2;
279   size_t len1;
280   size_t len2;
281
282   if (src1.nlimbs <= src2.nlimbs)
283     {
284       len1 = src1.nlimbs;
285       p1 = src1.limbs;
286       len2 = src2.nlimbs;
287       p2 = src2.limbs;
288     }
289   else
290     {
291       len1 = src2.nlimbs;
292       p1 = src2.limbs;
293       len2 = src1.nlimbs;
294       p2 = src1.limbs;
295     }
296   /* Now 0 <= len1 <= len2.  */
297   if (len1 == 0)
298     {
299       /* src1 or src2 is zero.  */
300       dest->nlimbs = 0;
301       dest->limbs = (mp_limb_t *) malloc (1);
302     }
303   else
304     {
305       /* Here 1 <= len1 <= len2.  */
306       size_t dlen;
307       mp_limb_t *dp;
308       size_t k, i, j;
309
310       dlen = len1 + len2;
311       dp = (mp_limb_t *) malloc (dlen * sizeof (mp_limb_t));
312       if (dp == NULL)
313         return NULL;
314       for (k = len2; k > 0; )
315         dp[--k] = 0;
316       for (i = 0; i < len1; i++)
317         {
318           mp_limb_t digit1 = p1[i];
319           mp_twolimb_t carry = 0;
320           for (j = 0; j < len2; j++)
321             {
322               mp_limb_t digit2 = p2[j];
323               carry += (mp_twolimb_t) digit1 * (mp_twolimb_t) digit2;
324               carry += dp[i + j];
325               dp[i + j] = (mp_limb_t) carry;
326               carry = carry >> GMP_LIMB_BITS;
327             }
328           dp[i + len2] = (mp_limb_t) carry;
329         }
330       /* Normalise.  */
331       while (dlen > 0 && dp[dlen - 1] == 0)
332         dlen--;
333       dest->nlimbs = dlen;
334       dest->limbs = dp;
335     }
336   return dest->limbs;
337 }
338
339 /* Compute the quotient of a bignum a >= 0 and a bignum b > 0.
340    a is written as  a = q * b + r  with 0 <= r < b.  q is the quotient, r
341    the remainder.
342    Finally, round-to-even is performed: If r > b/2 or if r = b/2 and q is odd,
343    q is incremented.
344    Return the allocated memory in case of success, NULL in case of memory
345    allocation failure.  */
346 static void *
347 divide (mpn_t a, mpn_t b, mpn_t *q)
348 {
349   /* Algorithm:
350      First normalise a and b: a=[a[m-1],...,a[0]], b=[b[n-1],...,b[0]]
351      with m>=0 and n>0 (in base beta = 2^GMP_LIMB_BITS).
352      If m<n, then q:=0 and r:=a.
353      If m>=n=1, perform a single-precision division:
354        r:=0, j:=m,
355        while j>0 do
356          {Here (q[m-1]*beta^(m-1)+...+q[j]*beta^j) * b[0] + r*beta^j =
357                = a[m-1]*beta^(m-1)+...+a[j]*beta^j und 0<=r<b[0]<beta}
358          j:=j-1, r:=r*beta+a[j], q[j]:=floor(r/b[0]), r:=r-b[0]*q[j].
359        Normalise [q[m-1],...,q[0]], yields q.
360      If m>=n>1, perform a multiple-precision division:
361        We have a/b < beta^(m-n+1).
362        s:=intDsize-1-(hightest bit in b[n-1]), 0<=s<intDsize.
363        Shift a and b left by s bits, copying them. r:=a.
364        r=[r[m],...,r[0]], b=[b[n-1],...,b[0]] with b[n-1]>=beta/2.
365        For j=m-n,...,0: {Here 0 <= r < b*beta^(j+1).}
366          Compute q* :
367            q* := floor((r[j+n]*beta+r[j+n-1])/b[n-1]).
368            In case of overflow (q* >= beta) set q* := beta-1.
369            Compute c2 := ((r[j+n]*beta+r[j+n-1]) - q* * b[n-1])*beta + r[j+n-2]
370            and c3 := b[n-2] * q*.
371            {We have 0 <= c2 < 2*beta^2, even 0 <= c2 < beta^2 if no overflow
372             occurred.  Furthermore 0 <= c3 < beta^2.
373             If there was overflow and
374             r[j+n]*beta+r[j+n-1] - q* * b[n-1] >= beta, i.e. c2 >= beta^2,
375             the next test can be skipped.}
376            While c3 > c2, {Here 0 <= c2 < c3 < beta^2}
377              Put q* := q* - 1, c2 := c2 + b[n-1]*beta, c3 := c3 - b[n-2].
378            If q* > 0:
379              Put r := r - b * q* * beta^j. In detail:
380                [r[n+j],...,r[j]] := [r[n+j],...,r[j]] - q* * [b[n-1],...,b[0]].
381                hence: u:=0, for i:=0 to n-1 do
382                               u := u + q* * b[i],
383                               r[j+i]:=r[j+i]-(u mod beta) (+ beta, if carry),
384                               u:=u div beta (+ 1, if carry in subtraction)
385                       r[n+j]:=r[n+j]-u.
386                {Since always u = (q* * [b[i-1],...,b[0]] div beta^i) + 1
387                                < q* + 1 <= beta,
388                 the carry u does not overflow.}
389              If a negative carry occurs, put q* := q* - 1
390                and [r[n+j],...,r[j]] := [r[n+j],...,r[j]] + [0,b[n-1],...,b[0]].
391          Set q[j] := q*.
392        Normalise [q[m-n],..,q[0]]; this yields the quotient q.
393        Shift [r[n-1],...,r[0]] right by s bits and normalise; this yields the
394        rest r.
395        The room for q[j] can be allocated at the memory location of r[n+j].
396      Finally, round-to-even:
397        Shift r left by 1 bit.
398        If r > b or if r = b and q[0] is odd, q := q+1.
399    */
400   const mp_limb_t *a_ptr = a.limbs;
401   size_t a_len = a.nlimbs;
402   const mp_limb_t *b_ptr = b.limbs;
403   size_t b_len = b.nlimbs;
404   mp_limb_t *roomptr;
405   mp_limb_t *tmp_roomptr = NULL;
406   mp_limb_t *q_ptr;
407   size_t q_len;
408   mp_limb_t *r_ptr;
409   size_t r_len;
410
411   /* Allocate room for a_len+2 digits.
412      (Need a_len+1 digits for the real division and 1 more digit for the
413      final rounding of q.)  */
414   roomptr = (mp_limb_t *) malloc ((a_len + 2) * sizeof (mp_limb_t));
415   if (roomptr == NULL)
416     return NULL;
417
418   /* Normalise a.  */
419   while (a_len > 0 && a_ptr[a_len - 1] == 0)
420     a_len--;
421
422   /* Normalise b.  */
423   for (;;)
424     {
425       if (b_len == 0)
426         /* Division by zero.  */
427         abort ();
428       if (b_ptr[b_len - 1] == 0)
429         b_len--;
430       else
431         break;
432     }
433
434   /* Here m = a_len >= 0 and n = b_len > 0.  */
435
436   if (a_len < b_len)
437     {
438       /* m<n: trivial case.  q=0, r := copy of a.  */
439       r_ptr = roomptr;
440       r_len = a_len;
441       memcpy (r_ptr, a_ptr, a_len * sizeof (mp_limb_t));
442       q_ptr = roomptr + a_len;
443       q_len = 0;
444     }
445   else if (b_len == 1)
446     {
447       /* n=1: single precision division.
448          beta^(m-1) <= a < beta^m  ==>  beta^(m-2) <= a/b < beta^m  */
449       r_ptr = roomptr;
450       q_ptr = roomptr + 1;
451       {
452         mp_limb_t den = b_ptr[0];
453         mp_limb_t remainder = 0;
454         const mp_limb_t *sourceptr = a_ptr + a_len;
455         mp_limb_t *destptr = q_ptr + a_len;
456         size_t count;
457         for (count = a_len; count > 0; count--)
458           {
459             mp_twolimb_t num =
460               ((mp_twolimb_t) remainder << GMP_LIMB_BITS) | *--sourceptr;
461             *--destptr = num / den;
462             remainder = num % den;
463           }
464         /* Normalise and store r.  */
465         if (remainder > 0)
466           {
467             r_ptr[0] = remainder;
468             r_len = 1;
469           }
470         else
471           r_len = 0;
472         /* Normalise q.  */
473         q_len = a_len;
474         if (q_ptr[q_len - 1] == 0)
475           q_len--;
476       }
477     }
478   else
479     {
480       /* n>1: multiple precision division.
481          beta^(m-1) <= a < beta^m, beta^(n-1) <= b < beta^n  ==>
482          beta^(m-n-1) <= a/b < beta^(m-n+1).  */
483       /* Determine s.  */
484       size_t s;
485       {
486         mp_limb_t msd = b_ptr[b_len - 1]; /* = b[n-1], > 0 */
487         s = 31;
488         if (msd >= 0x10000)
489           {
490             msd = msd >> 16;
491             s -= 16;
492           }
493         if (msd >= 0x100)
494           {
495             msd = msd >> 8;
496             s -= 8;
497           }
498         if (msd >= 0x10)
499           {
500             msd = msd >> 4;
501             s -= 4;
502           }
503         if (msd >= 0x4)
504           {
505             msd = msd >> 2;
506             s -= 2;
507           }
508         if (msd >= 0x2)
509           {
510             msd = msd >> 1;
511             s -= 1;
512           }
513       }
514       /* 0 <= s < GMP_LIMB_BITS.
515          Copy b, shifting it left by s bits.  */
516       if (s > 0)
517         {
518           tmp_roomptr = (mp_limb_t *) malloc (b_len * sizeof (mp_limb_t));
519           if (tmp_roomptr == NULL)
520             {
521               free (roomptr);
522               return NULL;
523             }
524           {
525             const mp_limb_t *sourceptr = b_ptr;
526             mp_limb_t *destptr = tmp_roomptr;
527             mp_twolimb_t accu = 0;
528             size_t count;
529             for (count = b_len; count > 0; count--)
530               {
531                 accu += (mp_twolimb_t) *sourceptr++ << s;
532                 *destptr++ = (mp_limb_t) accu;
533                 accu = accu >> GMP_LIMB_BITS;
534               }
535             /* accu must be zero, since that was how s was determined.  */
536             if (accu != 0)
537               abort ();
538           }
539           b_ptr = tmp_roomptr;
540         }
541       /* Copy a, shifting it left by s bits, yields r.
542          Memory layout:
543          At the beginning: r = roomptr[0..a_len],
544          at the end: r = roomptr[0..b_len-1], q = roomptr[b_len..a_len]  */
545       r_ptr = roomptr;
546       if (s == 0)
547         {
548           memcpy (r_ptr, a_ptr, a_len * sizeof (mp_limb_t));
549           r_ptr[a_len] = 0;
550         }
551       else
552         {
553           const mp_limb_t *sourceptr = a_ptr;
554           mp_limb_t *destptr = r_ptr;
555           mp_twolimb_t accu = 0;
556           size_t count;
557           for (count = a_len; count > 0; count--)
558             {
559               accu += (mp_twolimb_t) *sourceptr++ << s;
560               *destptr++ = (mp_limb_t) accu;
561               accu = accu >> GMP_LIMB_BITS;
562             }
563           *destptr++ = (mp_limb_t) accu;
564         }
565       q_ptr = roomptr + b_len;
566       q_len = a_len - b_len + 1; /* q will have m-n+1 limbs */
567       {
568         size_t j = a_len - b_len; /* m-n */
569         mp_limb_t b_msd = b_ptr[b_len - 1]; /* b[n-1] */
570         mp_limb_t b_2msd = b_ptr[b_len - 2]; /* b[n-2] */
571         mp_twolimb_t b_msdd = /* b[n-1]*beta+b[n-2] */
572           ((mp_twolimb_t) b_msd << GMP_LIMB_BITS) | b_2msd;
573         /* Division loop, traversed m-n+1 times.
574            j counts down, b is unchanged, beta/2 <= b[n-1] < beta.  */
575         for (;;)
576           {
577             mp_limb_t q_star;
578             mp_limb_t c1;
579             if (r_ptr[j + b_len] < b_msd) /* r[j+n] < b[n-1] ? */
580               {
581                 /* Divide r[j+n]*beta+r[j+n-1] by b[n-1], no overflow.  */
582                 mp_twolimb_t num =
583                   ((mp_twolimb_t) r_ptr[j + b_len] << GMP_LIMB_BITS)
584                   | r_ptr[j + b_len - 1];
585                 q_star = num / b_msd;
586                 c1 = num % b_msd;
587               }
588             else
589               {
590                 /* Overflow, hence r[j+n]*beta+r[j+n-1] >= beta*b[n-1].  */
591                 q_star = (mp_limb_t)~(mp_limb_t)0; /* q* = beta-1 */
592                 /* Test whether r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] >= beta
593                    <==> r[j+n]*beta+r[j+n-1] + b[n-1] >= beta*b[n-1]+beta
594                    <==> b[n-1] < floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta)
595                         {<= beta !}.
596                    If yes, jump directly to the subtraction loop.
597                    (Otherwise, r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] < beta
598                     <==> floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta) = b[n-1] ) */
599                 if (r_ptr[j + b_len] > b_msd
600                     || (c1 = r_ptr[j + b_len - 1] + b_msd) < b_msd)
601                   /* r[j+n] >= b[n-1]+1 or
602                      r[j+n] = b[n-1] and the addition r[j+n-1]+b[n-1] gives a
603                      carry.  */
604                   goto subtract;
605               }
606             /* q_star = q*,
607                c1 = (r[j+n]*beta+r[j+n-1]) - q* * b[n-1] (>=0, <beta).  */
608             {
609               mp_twolimb_t c2 = /* c1*beta+r[j+n-2] */
610                 ((mp_twolimb_t) c1 << GMP_LIMB_BITS) | r_ptr[j + b_len - 2];
611               mp_twolimb_t c3 = /* b[n-2] * q* */
612                 (mp_twolimb_t) b_2msd * (mp_twolimb_t) q_star;
613               /* While c2 < c3, increase c2 and decrease c3.
614                  Consider c3-c2.  While it is > 0, decrease it by
615                  b[n-1]*beta+b[n-2].  Because of b[n-1]*beta+b[n-2] >= beta^2/2
616                  this can happen only twice.  */
617               if (c3 > c2)
618                 {
619                   q_star = q_star - 1; /* q* := q* - 1 */
620                   if (c3 - c2 > b_msdd)
621                     q_star = q_star - 1; /* q* := q* - 1 */
622                 }
623             }
624             if (q_star > 0)
625               subtract:
626               {
627                 /* Subtract r := r - b * q* * beta^j.  */
628                 mp_limb_t cr;
629                 {
630                   const mp_limb_t *sourceptr = b_ptr;
631                   mp_limb_t *destptr = r_ptr + j;
632                   mp_twolimb_t carry = 0;
633                   size_t count;
634                   for (count = b_len; count > 0; count--)
635                     {
636                       /* Here 0 <= carry <= q*.  */
637                       carry =
638                         carry
639                         + (mp_twolimb_t) q_star * (mp_twolimb_t) *sourceptr++
640                         + (mp_limb_t) ~(*destptr);
641                       /* Here 0 <= carry <= beta*q* + beta-1.  */
642                       *destptr++ = ~(mp_limb_t) carry;
643                       carry = carry >> GMP_LIMB_BITS; /* <= q* */
644                     }
645                   cr = (mp_limb_t) carry;
646                 }
647                 /* Subtract cr from r_ptr[j + b_len], then forget about
648                    r_ptr[j + b_len].  */
649                 if (cr > r_ptr[j + b_len])
650                   {
651                     /* Subtraction gave a carry.  */
652                     q_star = q_star - 1; /* q* := q* - 1 */
653                     /* Add b back.  */
654                     {
655                       const mp_limb_t *sourceptr = b_ptr;
656                       mp_limb_t *destptr = r_ptr + j;
657                       mp_limb_t carry = 0;
658                       size_t count;
659                       for (count = b_len; count > 0; count--)
660                         {
661                           mp_limb_t source1 = *sourceptr++;
662                           mp_limb_t source2 = *destptr;
663                           *destptr++ = source1 + source2 + carry;
664                           carry =
665                             (carry
666                              ? source1 >= (mp_limb_t) ~source2
667                              : source1 > (mp_limb_t) ~source2);
668                         }
669                     }
670                     /* Forget about the carry and about r[j+n].  */
671                   }
672               }
673             /* q* is determined.  Store it as q[j].  */
674             q_ptr[j] = q_star;
675             if (j == 0)
676               break;
677             j--;
678           }
679       }
680       r_len = b_len;
681       /* Normalise q.  */
682       if (q_ptr[q_len - 1] == 0)
683         q_len--;
684 # if 0 /* Not needed here, since we need r only to compare it with b/2, and
685           b is shifted left by s bits.  */
686       /* Shift r right by s bits.  */
687       if (s > 0)
688         {
689           mp_limb_t ptr = r_ptr + r_len;
690           mp_twolimb_t accu = 0;
691           size_t count;
692           for (count = r_len; count > 0; count--)
693             {
694               accu = (mp_twolimb_t) (mp_limb_t) accu << GMP_LIMB_BITS;
695               accu += (mp_twolimb_t) *--ptr << (GMP_LIMB_BITS - s);
696               *ptr = (mp_limb_t) (accu >> GMP_LIMB_BITS);
697             }
698         }
699 # endif
700       /* Normalise r.  */
701       while (r_len > 0 && r_ptr[r_len - 1] == 0)
702         r_len--;
703     }
704   /* Compare r << 1 with b.  */
705   if (r_len > b_len)
706     goto increment_q;
707   {
708     size_t i;
709     for (i = b_len;;)
710       {
711         mp_limb_t r_i =
712           (i <= r_len && i > 0 ? r_ptr[i - 1] >> (GMP_LIMB_BITS - 1) : 0)
713           | (i < r_len ? r_ptr[i] << 1 : 0);
714         mp_limb_t b_i = (i < b_len ? b_ptr[i] : 0);
715         if (r_i > b_i)
716           goto increment_q;
717         if (r_i < b_i)
718           goto keep_q;
719         if (i == 0)
720           break;
721         i--;
722       }
723   }
724   if (q_len > 0 && ((q_ptr[0] & 1) != 0))
725     /* q is odd.  */
726     increment_q:
727     {
728       size_t i;
729       for (i = 0; i < q_len; i++)
730         if (++(q_ptr[i]) != 0)
731           goto keep_q;
732       q_ptr[q_len++] = 1;
733     }
734   keep_q:
735   if (tmp_roomptr != NULL)
736     free (tmp_roomptr);
737   q->limbs = q_ptr;
738   q->nlimbs = q_len;
739   return roomptr;
740 }
741
742 /* Convert a bignum a >= 0, multiplied with 10^extra_zeroes, to decimal
743    representation.
744    Destroys the contents of a.
745    Return the allocated memory - containing the decimal digits in low-to-high
746    order, terminated with a NUL character - in case of success, NULL in case
747    of memory allocation failure.  */
748 static char *
749 convert_to_decimal (mpn_t a, size_t extra_zeroes)
750 {
751   mp_limb_t *a_ptr = a.limbs;
752   size_t a_len = a.nlimbs;
753   /* 0.03345 is slightly larger than log(2)/(9*log(10)).  */
754   size_t c_len = 9 * ((size_t)(a_len * (GMP_LIMB_BITS * 0.03345f)) + 1);
755   char *c_ptr = (char *) malloc (xsum (c_len, extra_zeroes));
756   if (c_ptr != NULL)
757     {
758       char *d_ptr = c_ptr;
759       for (; extra_zeroes > 0; extra_zeroes--)
760         *d_ptr++ = '0';
761       while (a_len > 0)
762         {
763           /* Divide a by 10^9, in-place.  */
764           mp_limb_t remainder = 0;
765           mp_limb_t *ptr = a_ptr + a_len;
766           size_t count;
767           for (count = a_len; count > 0; count--)
768             {
769               mp_twolimb_t num =
770                 ((mp_twolimb_t) remainder << GMP_LIMB_BITS) | *--ptr;
771               *ptr = num / 1000000000;
772               remainder = num % 1000000000;
773             }
774           /* Store the remainder as 9 decimal digits.  */
775           for (count = 9; count > 0; count--)
776             {
777               *d_ptr++ = '0' + (remainder % 10);
778               remainder = remainder / 10;
779             }
780           /* Normalize a.  */
781           if (a_ptr[a_len - 1] == 0)
782             a_len--;
783         }
784       /* Remove leading zeroes.  */
785       while (d_ptr > c_ptr && d_ptr[-1] == '0')
786         d_ptr--;
787       /* But keep at least one zero.  */
788       if (d_ptr == c_ptr)
789         *d_ptr++ = '0';
790       /* Terminate the string.  */
791       *d_ptr = '\0';
792     }
793   return c_ptr;
794 }
795
796 # if NEED_PRINTF_LONG_DOUBLE
797
798 /* Assuming x is finite and >= 0:
799    write x as x = 2^e * m, where m is a bignum.
800    Return the allocated memory in case of success, NULL in case of memory
801    allocation failure.  */
802 static void *
803 decode_long_double (long double x, int *ep, mpn_t *mp)
804 {
805   mpn_t m;
806   int exp;
807   long double y;
808   size_t i;
809
810   /* Allocate memory for result.  */
811   m.nlimbs = (LDBL_MANT_BIT + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
812   m.limbs = (mp_limb_t *) malloc (m.nlimbs * sizeof (mp_limb_t));
813   if (m.limbs == NULL)
814     return NULL;
815   /* Split into exponential part and mantissa.  */
816   y = frexpl (x, &exp);
817   if (!(y >= 0.0L && y < 1.0L))
818     abort ();
819   /* x = 2^exp * y = 2^(exp - LDBL_MANT_BIT) * (y * LDBL_MANT_BIT), and the
820      latter is an integer.  */
821   /* Convert the mantissa (y * LDBL_MANT_BIT) to a sequence of limbs.
822      I'm not sure whether it's safe to cast a 'long double' value between
823      2^31 and 2^32 to 'unsigned int', therefore play safe and cast only
824      'long double' values between 0 and 2^16 (to 'unsigned int' or 'int',
825      doesn't matter).  */
826 #  if (LDBL_MANT_BIT % GMP_LIMB_BITS) != 0
827 #   if (LDBL_MANT_BIT % GMP_LIMB_BITS) > GMP_LIMB_BITS / 2
828     {
829       mp_limb_t hi, lo;
830       y *= (mp_limb_t) 1 << (LDBL_MANT_BIT % (GMP_LIMB_BITS / 2));
831       hi = (int) y;
832       y -= hi;
833       if (!(y >= 0.0L && y < 1.0L))
834         abort ();
835       y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
836       lo = (int) y;
837       y -= lo;
838       if (!(y >= 0.0L && y < 1.0L))
839         abort ();
840       m.limbs[LDBL_MANT_BIT / GMP_LIMB_BITS] = (hi << (GMP_LIMB_BITS / 2)) | lo;
841     }
842 #   else
843     {
844       mp_limb_t d;
845       y *= (mp_limb_t) 1 << (LDBL_MANT_BIT % GMP_LIMB_BITS);
846       d = (int) y;
847       y -= d;
848       if (!(y >= 0.0L && y < 1.0L))
849         abort ();
850       m.limbs[LDBL_MANT_BIT / GMP_LIMB_BITS] = d;
851     }
852 #   endif
853 #  endif
854   for (i = LDBL_MANT_BIT / GMP_LIMB_BITS; i > 0; )
855     {
856       mp_limb_t hi, lo;
857       y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
858       hi = (int) y;
859       y -= hi;
860       if (!(y >= 0.0L && y < 1.0L))
861         abort ();
862       y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
863       lo = (int) y;
864       y -= lo;
865       if (!(y >= 0.0L && y < 1.0L))
866         abort ();
867       m.limbs[--i] = (hi << (GMP_LIMB_BITS / 2)) | lo;
868     }
869 #if 0 /* On FreeBSD 6.1/x86, 'long double' numbers sometimes have excess
870          precision.  */
871   if (!(y == 0.0L))
872     abort ();
873 #endif
874   /* Normalise.  */
875   while (m.nlimbs > 0 && m.limbs[m.nlimbs - 1] == 0)
876     m.nlimbs--;
877   *mp = m;
878   *ep = exp - LDBL_MANT_BIT;
879   return m.limbs;
880 }
881
882 # endif
883
884 # if NEED_PRINTF_DOUBLE
885
886 /* Assuming x is finite and >= 0:
887    write x as x = 2^e * m, where m is a bignum.
888    Return the allocated memory in case of success, NULL in case of memory
889    allocation failure.  */
890 static void *
891 decode_double (double x, int *ep, mpn_t *mp)
892 {
893   mpn_t m;
894   int exp;
895   double y;
896   size_t i;
897
898   /* Allocate memory for result.  */
899   m.nlimbs = (DBL_MANT_BIT + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
900   m.limbs = (mp_limb_t *) malloc (m.nlimbs * sizeof (mp_limb_t));
901   if (m.limbs == NULL)
902     return NULL;
903   /* Split into exponential part and mantissa.  */
904   y = frexp (x, &exp);
905   if (!(y >= 0.0 && y < 1.0))
906     abort ();
907   /* x = 2^exp * y = 2^(exp - DBL_MANT_BIT) * (y * DBL_MANT_BIT), and the
908      latter is an integer.  */
909   /* Convert the mantissa (y * DBL_MANT_BIT) to a sequence of limbs.
910      I'm not sure whether it's safe to cast a 'double' value between
911      2^31 and 2^32 to 'unsigned int', therefore play safe and cast only
912      'double' values between 0 and 2^16 (to 'unsigned int' or 'int',
913      doesn't matter).  */
914 #  if (DBL_MANT_BIT % GMP_LIMB_BITS) != 0
915 #   if (DBL_MANT_BIT % GMP_LIMB_BITS) > GMP_LIMB_BITS / 2
916     {
917       mp_limb_t hi, lo;
918       y *= (mp_limb_t) 1 << (DBL_MANT_BIT % (GMP_LIMB_BITS / 2));
919       hi = (int) y;
920       y -= hi;
921       if (!(y >= 0.0 && y < 1.0))
922         abort ();
923       y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
924       lo = (int) y;
925       y -= lo;
926       if (!(y >= 0.0 && y < 1.0))
927         abort ();
928       m.limbs[DBL_MANT_BIT / GMP_LIMB_BITS] = (hi << (GMP_LIMB_BITS / 2)) | lo;
929     }
930 #   else
931     {
932       mp_limb_t d;
933       y *= (mp_limb_t) 1 << (DBL_MANT_BIT % GMP_LIMB_BITS);
934       d = (int) y;
935       y -= d;
936       if (!(y >= 0.0 && y < 1.0))
937         abort ();
938       m.limbs[DBL_MANT_BIT / GMP_LIMB_BITS] = d;
939     }
940 #   endif
941 #  endif
942   for (i = DBL_MANT_BIT / GMP_LIMB_BITS; i > 0; )
943     {
944       mp_limb_t hi, lo;
945       y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
946       hi = (int) y;
947       y -= hi;
948       if (!(y >= 0.0 && y < 1.0))
949         abort ();
950       y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
951       lo = (int) y;
952       y -= lo;
953       if (!(y >= 0.0 && y < 1.0))
954         abort ();
955       m.limbs[--i] = (hi << (GMP_LIMB_BITS / 2)) | lo;
956     }
957   if (!(y == 0.0))
958     abort ();
959   /* Normalise.  */
960   while (m.nlimbs > 0 && m.limbs[m.nlimbs - 1] == 0)
961     m.nlimbs--;
962   *mp = m;
963   *ep = exp - DBL_MANT_BIT;
964   return m.limbs;
965 }
966
967 # endif
968
969 /* Assuming x = 2^e * m is finite and >= 0, and n is an integer:
970    Returns the decimal representation of round (x * 10^n).
971    Return the allocated memory - containing the decimal digits in low-to-high
972    order, terminated with a NUL character - in case of success, NULL in case
973    of memory allocation failure.  */
974 static char *
975 scale10_round_decimal_decoded (int e, mpn_t m, void *memory, int n)
976 {
977   int s;
978   size_t extra_zeroes;
979   unsigned int abs_n;
980   unsigned int abs_s;
981   mp_limb_t *pow5_ptr;
982   size_t pow5_len;
983   unsigned int s_limbs;
984   unsigned int s_bits;
985   mpn_t pow5;
986   mpn_t z;
987   void *z_memory;
988   char *digits;
989
990   if (memory == NULL)
991     return NULL;
992   /* x = 2^e * m, hence
993      y = round (2^e * 10^n * m) = round (2^(e+n) * 5^n * m)
994        = round (2^s * 5^n * m).  */
995   s = e + n;
996   extra_zeroes = 0;
997   /* Factor out a common power of 10 if possible.  */
998   if (s > 0 && n > 0)
999     {
1000       extra_zeroes = (s < n ? s : n);
1001       s -= extra_zeroes;
1002       n -= extra_zeroes;
1003     }
1004   /* Here y = round (2^s * 5^n * m) * 10^extra_zeroes.
1005      Before converting to decimal, we need to compute
1006      z = round (2^s * 5^n * m).  */
1007   /* Compute 5^|n|, possibly shifted by |s| bits if n and s have the same
1008      sign.  2.322 is slightly larger than log(5)/log(2).  */
1009   abs_n = (n >= 0 ? n : -n);
1010   abs_s = (s >= 0 ? s : -s);
1011   pow5_ptr = (mp_limb_t *) malloc (((int)(abs_n * (2.322f / GMP_LIMB_BITS)) + 1
1012                                     + abs_s / GMP_LIMB_BITS + 1)
1013                                    * sizeof (mp_limb_t));
1014   if (pow5_ptr == NULL)
1015     {
1016       free (memory);
1017       return NULL;
1018     }
1019   /* Initialize with 1.  */
1020   pow5_ptr[0] = 1;
1021   pow5_len = 1;
1022   /* Multiply with 5^|n|.  */
1023   if (abs_n > 0)
1024     {
1025       static mp_limb_t const small_pow5[13 + 1] =
1026         {
1027           1, 5, 25, 125, 625, 3125, 15625, 78125, 390625, 1953125, 9765625,
1028           48828125, 244140625, 1220703125
1029         };
1030       unsigned int n13;
1031       for (n13 = 0; n13 <= abs_n; n13 += 13)
1032         {
1033           mp_limb_t digit1 = small_pow5[n13 + 13 <= abs_n ? 13 : abs_n - n13];
1034           size_t j;
1035           mp_twolimb_t carry = 0;
1036           for (j = 0; j < pow5_len; j++)
1037             {
1038               mp_limb_t digit2 = pow5_ptr[j];
1039               carry += (mp_twolimb_t) digit1 * (mp_twolimb_t) digit2;
1040               pow5_ptr[j] = (mp_limb_t) carry;
1041               carry = carry >> GMP_LIMB_BITS;
1042             }
1043           if (carry > 0)
1044             pow5_ptr[pow5_len++] = (mp_limb_t) carry;
1045         }
1046     }
1047   s_limbs = abs_s / GMP_LIMB_BITS;
1048   s_bits = abs_s % GMP_LIMB_BITS;
1049   if (n >= 0 ? s >= 0 : s <= 0)
1050     {
1051       /* Multiply with 2^|s|.  */
1052       if (s_bits > 0)
1053         {
1054           mp_limb_t *ptr = pow5_ptr;
1055           mp_twolimb_t accu = 0;
1056           size_t count;
1057           for (count = pow5_len; count > 0; count--)
1058             {
1059               accu += (mp_twolimb_t) *ptr << s_bits;
1060               *ptr++ = (mp_limb_t) accu;
1061               accu = accu >> GMP_LIMB_BITS;
1062             }
1063           if (accu > 0)
1064             {
1065               *ptr = (mp_limb_t) accu;
1066               pow5_len++;
1067             }
1068         }
1069       if (s_limbs > 0)
1070         {
1071           size_t count;
1072           for (count = pow5_len; count > 0;)
1073             {
1074               count--;
1075               pow5_ptr[s_limbs + count] = pow5_ptr[count];
1076             }
1077           for (count = s_limbs; count > 0;)
1078             {
1079               count--;
1080               pow5_ptr[count] = 0;
1081             }
1082           pow5_len += s_limbs;
1083         }
1084       pow5.limbs = pow5_ptr;
1085       pow5.nlimbs = pow5_len;
1086       if (n >= 0)
1087         {
1088           /* Multiply m with pow5.  No division needed.  */
1089           z_memory = multiply (m, pow5, &z);
1090         }
1091       else
1092         {
1093           /* Divide m by pow5 and round.  */
1094           z_memory = divide (m, pow5, &z);
1095         }
1096     }
1097   else
1098     {
1099       pow5.limbs = pow5_ptr;
1100       pow5.nlimbs = pow5_len;
1101       if (n >= 0)
1102         {
1103           /* n >= 0, s < 0.
1104              Multiply m with pow5, then divide by 2^|s|.  */
1105           mpn_t numerator;
1106           mpn_t denominator;
1107           void *tmp_memory;
1108           tmp_memory = multiply (m, pow5, &numerator);
1109           if (tmp_memory == NULL)
1110             {
1111               free (pow5_ptr);
1112               free (memory);
1113               return NULL;
1114             }
1115           /* Construct 2^|s|.  */
1116           {
1117             mp_limb_t *ptr = pow5_ptr + pow5_len;
1118             size_t i;
1119             for (i = 0; i < s_limbs; i++)
1120               ptr[i] = 0;
1121             ptr[s_limbs] = (mp_limb_t) 1 << s_bits;
1122             denominator.limbs = ptr;
1123             denominator.nlimbs = s_limbs + 1;
1124           }
1125           z_memory = divide (numerator, denominator, &z);
1126           free (tmp_memory);
1127         }
1128       else
1129         {
1130           /* n < 0, s > 0.
1131              Multiply m with 2^s, then divide by pow5.  */
1132           mpn_t numerator;
1133           mp_limb_t *num_ptr;
1134           num_ptr = (mp_limb_t *) malloc ((m.nlimbs + s_limbs + 1)
1135                                           * sizeof (mp_limb_t));
1136           if (num_ptr == NULL)
1137             {
1138               free (pow5_ptr);
1139               free (memory);
1140               return NULL;
1141             }
1142           {
1143             mp_limb_t *destptr = num_ptr;
1144             {
1145               size_t i;
1146               for (i = 0; i < s_limbs; i++)
1147                 *destptr++ = 0;
1148             }
1149             if (s_bits > 0)
1150               {
1151                 const mp_limb_t *sourceptr = m.limbs;
1152                 mp_twolimb_t accu = 0;
1153                 size_t count;
1154                 for (count = m.nlimbs; count > 0; count--)
1155                   {
1156                     accu += (mp_twolimb_t) *sourceptr++ << s_bits;
1157                     *destptr++ = (mp_limb_t) accu;
1158                     accu = accu >> GMP_LIMB_BITS;
1159                   }
1160                 if (accu > 0)
1161                   *destptr++ = (mp_limb_t) accu;
1162               }
1163             else
1164               {
1165                 const mp_limb_t *sourceptr = m.limbs;
1166                 size_t count;
1167                 for (count = m.nlimbs; count > 0; count--)
1168                   *destptr++ = *sourceptr++;
1169               }
1170             numerator.limbs = num_ptr;
1171             numerator.nlimbs = destptr - num_ptr;
1172           }
1173           z_memory = divide (numerator, pow5, &z);
1174           free (num_ptr);
1175         }
1176     }
1177   free (pow5_ptr);
1178   free (memory);
1179
1180   /* Here y = round (x * 10^n) = z * 10^extra_zeroes.  */
1181
1182   if (z_memory == NULL)
1183     return NULL;
1184   digits = convert_to_decimal (z, extra_zeroes);
1185   free (z_memory);
1186   return digits;
1187 }
1188
1189 # if NEED_PRINTF_LONG_DOUBLE
1190
1191 /* Assuming x is finite and >= 0, and n is an integer:
1192    Returns the decimal representation of round (x * 10^n).
1193    Return the allocated memory - containing the decimal digits in low-to-high
1194    order, terminated with a NUL character - in case of success, NULL in case
1195    of memory allocation failure.  */
1196 static char *
1197 scale10_round_decimal_long_double (long double x, int n)
1198 {
1199   int e;
1200   mpn_t m;
1201   void *memory = decode_long_double (x, &e, &m);
1202   return scale10_round_decimal_decoded (e, m, memory, n);
1203 }
1204
1205 # endif
1206
1207 # if NEED_PRINTF_DOUBLE
1208
1209 /* Assuming x is finite and >= 0, and n is an integer:
1210    Returns the decimal representation of round (x * 10^n).
1211    Return the allocated memory - containing the decimal digits in low-to-high
1212    order, terminated with a NUL character - in case of success, NULL in case
1213    of memory allocation failure.  */
1214 static char *
1215 scale10_round_decimal_double (double x, int n)
1216 {
1217   int e;
1218   mpn_t m;
1219   void *memory = decode_double (x, &e, &m);
1220   return scale10_round_decimal_decoded (e, m, memory, n);
1221 }
1222
1223 # endif
1224
1225 # if NEED_PRINTF_LONG_DOUBLE
1226
1227 /* Assuming x is finite and > 0:
1228    Return an approximation for n with 10^n <= x < 10^(n+1).
1229    The approximation is usually the right n, but may be off by 1 sometimes.  */
1230 static int
1231 floorlog10l (long double x)
1232 {
1233   int exp;
1234   long double y;
1235   double z;
1236   double l;
1237
1238   /* Split into exponential part and mantissa.  */
1239   y = frexpl (x, &exp);
1240   if (!(y >= 0.0L && y < 1.0L))
1241     abort ();
1242   if (y == 0.0L)
1243     return INT_MIN;
1244   if (y < 0.5L)
1245     {
1246       while (y < (1.0L / (1 << (GMP_LIMB_BITS / 2)) / (1 << (GMP_LIMB_BITS / 2))))
1247         {
1248           y *= 1.0L * (1 << (GMP_LIMB_BITS / 2)) * (1 << (GMP_LIMB_BITS / 2));
1249           exp -= GMP_LIMB_BITS;
1250         }
1251       if (y < (1.0L / (1 << 16)))
1252         {
1253           y *= 1.0L * (1 << 16);
1254           exp -= 16;
1255         }
1256       if (y < (1.0L / (1 << 8)))
1257         {
1258           y *= 1.0L * (1 << 8);
1259           exp -= 8;
1260         }
1261       if (y < (1.0L / (1 << 4)))
1262         {
1263           y *= 1.0L * (1 << 4);
1264           exp -= 4;
1265         }
1266       if (y < (1.0L / (1 << 2)))
1267         {
1268           y *= 1.0L * (1 << 2);
1269           exp -= 2;
1270         }
1271       if (y < (1.0L / (1 << 1)))
1272         {
1273           y *= 1.0L * (1 << 1);
1274           exp -= 1;
1275         }
1276     }
1277   if (!(y >= 0.5L && y < 1.0L))
1278     abort ();
1279   /* Compute an approximation for l = log2(x) = exp + log2(y).  */
1280   l = exp;
1281   z = y;
1282   if (z < 0.70710678118654752444)
1283     {
1284       z *= 1.4142135623730950488;
1285       l -= 0.5;
1286     }
1287   if (z < 0.8408964152537145431)
1288     {
1289       z *= 1.1892071150027210667;
1290       l -= 0.25;
1291     }
1292   if (z < 0.91700404320467123175)
1293     {
1294       z *= 1.0905077326652576592;
1295       l -= 0.125;
1296     }
1297   if (z < 0.9576032806985736469)
1298     {
1299       z *= 1.0442737824274138403;
1300       l -= 0.0625;
1301     }
1302   /* Now 0.95 <= z <= 1.01.  */
1303   z = 1 - z;
1304   /* log(1-z) = - z - z^2/2 - z^3/3 - z^4/4 - ...
1305      Four terms are enough to get an approximation with error < 10^-7.  */
1306   l -= z * (1.0 + z * (0.5 + z * ((1.0 / 3) + z * 0.25)));
1307   /* Finally multiply with log(2)/log(10), yields an approximation for
1308      log10(x).  */
1309   l *= 0.30102999566398119523;
1310   /* Round down to the next integer.  */
1311   return (int) l + (l < 0 ? -1 : 0);
1312 }
1313
1314 # endif
1315
1316 # if NEED_PRINTF_DOUBLE
1317
1318 /* Assuming x is finite and > 0:
1319    Return an approximation for n with 10^n <= x < 10^(n+1).
1320    The approximation is usually the right n, but may be off by 1 sometimes.  */
1321 static int
1322 floorlog10 (double x)
1323 {
1324   int exp;
1325   double y;
1326   double z;
1327   double l;
1328
1329   /* Split into exponential part and mantissa.  */
1330   y = frexp (x, &exp);
1331   if (!(y >= 0.0 && y < 1.0))
1332     abort ();
1333   if (y == 0.0)
1334     return INT_MIN;
1335   if (y < 0.5)
1336     {
1337       while (y < (1.0 / (1 << (GMP_LIMB_BITS / 2)) / (1 << (GMP_LIMB_BITS / 2))))
1338         {
1339           y *= 1.0 * (1 << (GMP_LIMB_BITS / 2)) * (1 << (GMP_LIMB_BITS / 2));
1340           exp -= GMP_LIMB_BITS;
1341         }
1342       if (y < (1.0 / (1 << 16)))
1343         {
1344           y *= 1.0 * (1 << 16);
1345           exp -= 16;
1346         }
1347       if (y < (1.0 / (1 << 8)))
1348         {
1349           y *= 1.0 * (1 << 8);
1350           exp -= 8;
1351         }
1352       if (y < (1.0 / (1 << 4)))
1353         {
1354           y *= 1.0 * (1 << 4);
1355           exp -= 4;
1356         }
1357       if (y < (1.0 / (1 << 2)))
1358         {
1359           y *= 1.0 * (1 << 2);
1360           exp -= 2;
1361         }
1362       if (y < (1.0 / (1 << 1)))
1363         {
1364           y *= 1.0 * (1 << 1);
1365           exp -= 1;
1366         }
1367     }
1368   if (!(y >= 0.5 && y < 1.0))
1369     abort ();
1370   /* Compute an approximation for l = log2(x) = exp + log2(y).  */
1371   l = exp;
1372   z = y;
1373   if (z < 0.70710678118654752444)
1374     {
1375       z *= 1.4142135623730950488;
1376       l -= 0.5;
1377     }
1378   if (z < 0.8408964152537145431)
1379     {
1380       z *= 1.1892071150027210667;
1381       l -= 0.25;
1382     }
1383   if (z < 0.91700404320467123175)
1384     {
1385       z *= 1.0905077326652576592;
1386       l -= 0.125;
1387     }
1388   if (z < 0.9576032806985736469)
1389     {
1390       z *= 1.0442737824274138403;
1391       l -= 0.0625;
1392     }
1393   /* Now 0.95 <= z <= 1.01.  */
1394   z = 1 - z;
1395   /* log(1-z) = - z - z^2/2 - z^3/3 - z^4/4 - ...
1396      Four terms are enough to get an approximation with error < 10^-7.  */
1397   l -= z * (1.0 + z * (0.5 + z * ((1.0 / 3) + z * 0.25)));
1398   /* Finally multiply with log(2)/log(10), yields an approximation for
1399      log10(x).  */
1400   l *= 0.30102999566398119523;
1401   /* Round down to the next integer.  */
1402   return (int) l + (l < 0 ? -1 : 0);
1403 }
1404
1405 # endif
1406
1407 #endif
1408
1409 DCHAR_T *
1410 VASNPRINTF (DCHAR_T *resultbuf, size_t *lengthp,
1411             const FCHAR_T *format, va_list args)
1412 {
1413   DIRECTIVES d;
1414   arguments a;
1415
1416   if (PRINTF_PARSE (format, &d, &a) < 0)
1417     /* errno is already set.  */
1418     return NULL;
1419
1420 #define CLEANUP() \
1421   free (d.dir);                                                         \
1422   if (a.arg)                                                            \
1423     free (a.arg);
1424
1425   if (PRINTF_FETCHARGS (args, &a) < 0)
1426     {
1427       CLEANUP ();
1428       errno = EINVAL;
1429       return NULL;
1430     }
1431
1432   {
1433     size_t buf_neededlength;
1434     TCHAR_T *buf;
1435     TCHAR_T *buf_malloced;
1436     const FCHAR_T *cp;
1437     size_t i;
1438     DIRECTIVE *dp;
1439     /* Output string accumulator.  */
1440     DCHAR_T *result;
1441     size_t allocated;
1442     size_t length;
1443
1444     /* Allocate a small buffer that will hold a directive passed to
1445        sprintf or snprintf.  */
1446     buf_neededlength =
1447       xsum4 (7, d.max_width_length, d.max_precision_length, 6);
1448 #if HAVE_ALLOCA
1449     if (buf_neededlength < 4000 / sizeof (TCHAR_T))
1450       {
1451         buf = (TCHAR_T *) alloca (buf_neededlength * sizeof (TCHAR_T));
1452         buf_malloced = NULL;
1453       }
1454     else
1455 #endif
1456       {
1457         size_t buf_memsize = xtimes (buf_neededlength, sizeof (TCHAR_T));
1458         if (size_overflow_p (buf_memsize))
1459           goto out_of_memory_1;
1460         buf = (TCHAR_T *) malloc (buf_memsize);
1461         if (buf == NULL)
1462           goto out_of_memory_1;
1463         buf_malloced = buf;
1464       }
1465
1466     if (resultbuf != NULL)
1467       {
1468         result = resultbuf;
1469         allocated = *lengthp;
1470       }
1471     else
1472       {
1473         result = NULL;
1474         allocated = 0;
1475       }
1476     length = 0;
1477     /* Invariants:
1478        result is either == resultbuf or == NULL or malloc-allocated.
1479        If length > 0, then result != NULL.  */
1480
1481     /* Ensures that allocated >= needed.  Aborts through a jump to
1482        out_of_memory if needed is SIZE_MAX or otherwise too big.  */
1483 #define ENSURE_ALLOCATION(needed) \
1484     if ((needed) > allocated)                                                \
1485       {                                                                      \
1486         size_t memory_size;                                                  \
1487         DCHAR_T *memory;                                                     \
1488                                                                              \
1489         allocated = (allocated > 0 ? xtimes (allocated, 2) : 12);            \
1490         if ((needed) > allocated)                                            \
1491           allocated = (needed);                                              \
1492         memory_size = xtimes (allocated, sizeof (DCHAR_T));                  \
1493         if (size_overflow_p (memory_size))                                   \
1494           goto out_of_memory;                                                \
1495         if (result == resultbuf || result == NULL)                           \
1496           memory = (DCHAR_T *) malloc (memory_size);                         \
1497         else                                                                 \
1498           memory = (DCHAR_T *) realloc (result, memory_size);                \
1499         if (memory == NULL)                                                  \
1500           goto out_of_memory;                                                \
1501         if (result == resultbuf && length > 0)                               \
1502           DCHAR_CPY (memory, result, length);                                \
1503         result = memory;                                                     \
1504       }
1505
1506     for (cp = format, i = 0, dp = &d.dir[0]; ; cp = dp->dir_end, i++, dp++)
1507       {
1508         if (cp != dp->dir_start)
1509           {
1510             size_t n = dp->dir_start - cp;
1511             size_t augmented_length = xsum (length, n);
1512
1513             ENSURE_ALLOCATION (augmented_length);
1514             /* This copies a piece of FCHAR_T[] into a DCHAR_T[].  Here we
1515                need that the format string contains only ASCII characters
1516                if FCHAR_T and DCHAR_T are not the same type.  */
1517             if (sizeof (FCHAR_T) == sizeof (DCHAR_T))
1518               {
1519                 DCHAR_CPY (result + length, (const DCHAR_T *) cp, n);
1520                 length = augmented_length;
1521               }
1522             else
1523               {
1524                 do
1525                   result[length++] = (unsigned char) *cp++;
1526                 while (--n > 0);
1527               }
1528           }
1529         if (i == d.count)
1530           break;
1531
1532         /* Execute a single directive.  */
1533         if (dp->conversion == '%')
1534           {
1535             size_t augmented_length;
1536
1537             if (!(dp->arg_index == ARG_NONE))
1538               abort ();
1539             augmented_length = xsum (length, 1);
1540             ENSURE_ALLOCATION (augmented_length);
1541             result[length] = '%';
1542             length = augmented_length;
1543           }
1544         else
1545           {
1546             if (!(dp->arg_index != ARG_NONE))
1547               abort ();
1548
1549             if (dp->conversion == 'n')
1550               {
1551                 switch (a.arg[dp->arg_index].type)
1552                   {
1553                   case TYPE_COUNT_SCHAR_POINTER:
1554                     *a.arg[dp->arg_index].a.a_count_schar_pointer = length;
1555                     break;
1556                   case TYPE_COUNT_SHORT_POINTER:
1557                     *a.arg[dp->arg_index].a.a_count_short_pointer = length;
1558                     break;
1559                   case TYPE_COUNT_INT_POINTER:
1560                     *a.arg[dp->arg_index].a.a_count_int_pointer = length;
1561                     break;
1562                   case TYPE_COUNT_LONGINT_POINTER:
1563                     *a.arg[dp->arg_index].a.a_count_longint_pointer = length;
1564                     break;
1565 #if HAVE_LONG_LONG_INT
1566                   case TYPE_COUNT_LONGLONGINT_POINTER:
1567                     *a.arg[dp->arg_index].a.a_count_longlongint_pointer = length;
1568                     break;
1569 #endif
1570                   default:
1571                     abort ();
1572                   }
1573               }
1574 #if ENABLE_UNISTDIO
1575             /* The unistdio extensions.  */
1576             else if (dp->conversion == 'U')
1577               {
1578                 arg_type type = a.arg[dp->arg_index].type;
1579                 int flags = dp->flags;
1580                 int has_width;
1581                 size_t width;
1582                 int has_precision;
1583                 size_t precision;
1584
1585                 has_width = 0;
1586                 width = 0;
1587                 if (dp->width_start != dp->width_end)
1588                   {
1589                     if (dp->width_arg_index != ARG_NONE)
1590                       {
1591                         int arg;
1592
1593                         if (!(a.arg[dp->width_arg_index].type == TYPE_INT))
1594                           abort ();
1595                         arg = a.arg[dp->width_arg_index].a.a_int;
1596                         if (arg < 0)
1597                           {
1598                             /* "A negative field width is taken as a '-' flag
1599                                 followed by a positive field width."  */
1600                             flags |= FLAG_LEFT;
1601                             width = (unsigned int) (-arg);
1602                           }
1603                         else
1604                           width = arg;
1605                       }
1606                     else
1607                       {
1608                         const FCHAR_T *digitp = dp->width_start;
1609
1610                         do
1611                           width = xsum (xtimes (width, 10), *digitp++ - '0');
1612                         while (digitp != dp->width_end);
1613                       }
1614                     has_width = 1;
1615                   }
1616
1617                 has_precision = 0;
1618                 precision = 0;
1619                 if (dp->precision_start != dp->precision_end)
1620                   {
1621                     if (dp->precision_arg_index != ARG_NONE)
1622                       {
1623                         int arg;
1624
1625                         if (!(a.arg[dp->precision_arg_index].type == TYPE_INT))
1626                           abort ();
1627                         arg = a.arg[dp->precision_arg_index].a.a_int;
1628                         /* "A negative precision is taken as if the precision
1629                             were omitted."  */
1630                         if (arg >= 0)
1631                           {
1632                             precision = arg;
1633                             has_precision = 1;
1634                           }
1635                       }
1636                     else
1637                       {
1638                         const FCHAR_T *digitp = dp->precision_start + 1;
1639
1640                         precision = 0;
1641                         while (digitp != dp->precision_end)
1642                           precision = xsum (xtimes (precision, 10), *digitp++ - '0');
1643                         has_precision = 1;
1644                       }
1645                   }
1646
1647                 switch (type)
1648                   {
1649                   case TYPE_U8_STRING:
1650                     {
1651                       const uint8_t *arg = a.arg[dp->arg_index].a.a_u8_string;
1652                       const uint8_t *arg_end;
1653                       size_t characters;
1654
1655                       if (has_precision)
1656                         {
1657                           /* Use only PRECISION characters, from the left.  */
1658                           arg_end = arg;
1659                           characters = 0;
1660                           for (; precision > 0; precision--)
1661                             {
1662                               int count = u8_strmblen (arg_end);
1663                               if (count == 0)
1664                                 break;
1665                               if (count < 0)
1666                                 {
1667                                   if (!(result == resultbuf || result == NULL))
1668                                     free (result);
1669                                   if (buf_malloced != NULL)
1670                                     free (buf_malloced);
1671                                   CLEANUP ();
1672                                   errno = EILSEQ;
1673                                   return NULL;
1674                                 }
1675                               arg_end += count;
1676                               characters++;
1677                             }
1678                         }
1679                       else if (has_width)
1680                         {
1681                           /* Use the entire string, and count the number of
1682                              characters.  */
1683                           arg_end = arg;
1684                           characters = 0;
1685                           for (;;)
1686                             {
1687                               int count = u8_strmblen (arg_end);
1688                               if (count == 0)
1689                                 break;
1690                               if (count < 0)
1691                                 {
1692                                   if (!(result == resultbuf || result == NULL))
1693                                     free (result);
1694                                   if (buf_malloced != NULL)
1695                                     free (buf_malloced);
1696                                   CLEANUP ();
1697                                   errno = EILSEQ;
1698                                   return NULL;
1699                                 }
1700                               arg_end += count;
1701                               characters++;
1702                             }
1703                         }
1704                       else
1705                         {
1706                           /* Use the entire string.  */
1707                           arg_end = arg + u8_strlen (arg);
1708                           /* The number of characters doesn't matter.  */
1709                           characters = 0;
1710                         }
1711
1712                       if (has_width && width > characters
1713                           && !(dp->flags & FLAG_LEFT))
1714                         {
1715                           size_t n = width - characters;
1716                           ENSURE_ALLOCATION (xsum (length, n));
1717                           DCHAR_SET (result + length, ' ', n);
1718                           length += n;
1719                         }
1720
1721 # if DCHAR_IS_UINT8_T
1722                       {
1723                         size_t n = arg_end - arg;
1724                         ENSURE_ALLOCATION (xsum (length, n));
1725                         DCHAR_CPY (result + length, arg, n);
1726                         length += n;
1727                       }
1728 # else
1729                       { /* Convert.  */
1730                         DCHAR_T *converted = result + length;
1731                         size_t converted_len = allocated - length;
1732 #  if DCHAR_IS_TCHAR
1733                         /* Convert from UTF-8 to locale encoding.  */
1734                         if (u8_conv_to_encoding (locale_charset (),
1735                                                  iconveh_question_mark,
1736                                                  arg, arg_end - arg, NULL,
1737                                                  &converted, &converted_len)
1738                             < 0)
1739 #  else
1740                         /* Convert from UTF-8 to UTF-16/UTF-32.  */
1741                         converted =
1742                           U8_TO_DCHAR (arg, arg_end - arg,
1743                                        converted, &converted_len);
1744                         if (converted == NULL)
1745 #  endif
1746                           {
1747                             int saved_errno = errno;
1748                             if (!(result == resultbuf || result == NULL))
1749                               free (result);
1750                             if (buf_malloced != NULL)
1751                               free (buf_malloced);
1752                             CLEANUP ();
1753                             errno = saved_errno;
1754                             return NULL;
1755                           }
1756                         if (converted != result + length)
1757                           {
1758                             ENSURE_ALLOCATION (xsum (length, converted_len));
1759                             DCHAR_CPY (result + length, converted, converted_len);
1760                             free (converted);
1761                           }
1762                         length += converted_len;
1763                       }
1764 # endif
1765
1766                       if (has_width && width > characters
1767                           && (dp->flags & FLAG_LEFT))
1768                         {
1769                           size_t n = width - characters;
1770                           ENSURE_ALLOCATION (xsum (length, n));
1771                           DCHAR_SET (result + length, ' ', n);
1772                           length += n;
1773                         }
1774                     }
1775                     break;
1776
1777                   case TYPE_U16_STRING:
1778                     {
1779                       const uint16_t *arg = a.arg[dp->arg_index].a.a_u16_string;
1780                       const uint16_t *arg_end;
1781                       size_t characters;
1782
1783                       if (has_precision)
1784                         {
1785                           /* Use only PRECISION characters, from the left.  */
1786                           arg_end = arg;
1787                           characters = 0;
1788                           for (; precision > 0; precision--)
1789                             {
1790                               int count = u16_strmblen (arg_end);
1791                               if (count == 0)
1792                                 break;
1793                               if (count < 0)
1794                                 {
1795                                   if (!(result == resultbuf || result == NULL))
1796                                     free (result);
1797                                   if (buf_malloced != NULL)
1798                                     free (buf_malloced);
1799                                   CLEANUP ();
1800                                   errno = EILSEQ;
1801                                   return NULL;
1802                                 }
1803                               arg_end += count;
1804                               characters++;
1805                             }
1806                         }
1807                       else if (has_width)
1808                         {
1809                           /* Use the entire string, and count the number of
1810                              characters.  */
1811                           arg_end = arg;
1812                           characters = 0;
1813                           for (;;)
1814                             {
1815                               int count = u16_strmblen (arg_end);
1816                               if (count == 0)
1817                                 break;
1818                               if (count < 0)
1819                                 {
1820                                   if (!(result == resultbuf || result == NULL))
1821                                     free (result);
1822                                   if (buf_malloced != NULL)
1823                                     free (buf_malloced);
1824                                   CLEANUP ();
1825                                   errno = EILSEQ;
1826                                   return NULL;
1827                                 }
1828                               arg_end += count;
1829                               characters++;
1830                             }
1831                         }
1832                       else
1833                         {
1834                           /* Use the entire string.  */
1835                           arg_end = arg + u16_strlen (arg);
1836                           /* The number of characters doesn't matter.  */
1837                           characters = 0;
1838                         }
1839
1840                       if (has_width && width > characters
1841                           && !(dp->flags & FLAG_LEFT))
1842                         {
1843                           size_t n = width - characters;
1844                           ENSURE_ALLOCATION (xsum (length, n));
1845                           DCHAR_SET (result + length, ' ', n);
1846                           length += n;
1847                         }
1848
1849 # if DCHAR_IS_UINT16_T
1850                       {
1851                         size_t n = arg_end - arg;
1852                         ENSURE_ALLOCATION (xsum (length, n));
1853                         DCHAR_CPY (result + length, arg, n);
1854                         length += n;
1855                       }
1856 # else
1857                       { /* Convert.  */
1858                         DCHAR_T *converted = result + length;
1859                         size_t converted_len = allocated - length;
1860 #  if DCHAR_IS_TCHAR
1861                         /* Convert from UTF-16 to locale encoding.  */
1862                         if (u16_conv_to_encoding (locale_charset (),
1863                                                   iconveh_question_mark,
1864                                                   arg, arg_end - arg, NULL,
1865                                                   &converted, &converted_len)
1866                             < 0)
1867 #  else
1868                         /* Convert from UTF-16 to UTF-8/UTF-32.  */
1869                         converted =
1870                           U16_TO_DCHAR (arg, arg_end - arg,
1871                                         converted, &converted_len);
1872                         if (converted == NULL)
1873 #  endif
1874                           {
1875                             int saved_errno = errno;
1876                             if (!(result == resultbuf || result == NULL))
1877                               free (result);
1878                             if (buf_malloced != NULL)
1879                               free (buf_malloced);
1880                             CLEANUP ();
1881                             errno = saved_errno;
1882                             return NULL;
1883                           }
1884                         if (converted != result + length)
1885                           {
1886                             ENSURE_ALLOCATION (xsum (length, converted_len));
1887                             DCHAR_CPY (result + length, converted, converted_len);
1888                             free (converted);
1889                           }
1890                         length += converted_len;
1891                       }
1892 # endif
1893
1894                       if (has_width && width > characters
1895                           && (dp->flags & FLAG_LEFT))
1896                         {
1897                           size_t n = width - characters;
1898                           ENSURE_ALLOCATION (xsum (length, n));
1899                           DCHAR_SET (result + length, ' ', n);
1900                           length += n;
1901                         }
1902                     }
1903                     break;
1904
1905                   case TYPE_U32_STRING:
1906                     {
1907                       const uint32_t *arg = a.arg[dp->arg_index].a.a_u32_string;
1908                       const uint32_t *arg_end;
1909                       size_t characters;
1910
1911                       if (has_precision)
1912                         {
1913                           /* Use only PRECISION characters, from the left.  */
1914                           arg_end = arg;
1915                           characters = 0;
1916                           for (; precision > 0; precision--)
1917                             {
1918                               int count = u32_strmblen (arg_end);
1919                               if (count == 0)
1920                                 break;
1921                               if (count < 0)
1922                                 {
1923                                   if (!(result == resultbuf || result == NULL))
1924                                     free (result);
1925                                   if (buf_malloced != NULL)
1926                                     free (buf_malloced);
1927                                   CLEANUP ();
1928                                   errno = EILSEQ;
1929                                   return NULL;
1930                                 }
1931                               arg_end += count;
1932                               characters++;
1933                             }
1934                         }
1935                       else if (has_width)
1936                         {
1937                           /* Use the entire string, and count the number of
1938                              characters.  */
1939                           arg_end = arg;
1940                           characters = 0;
1941                           for (;;)
1942                             {
1943                               int count = u32_strmblen (arg_end);
1944                               if (count == 0)
1945                                 break;
1946                               if (count < 0)
1947                                 {
1948                                   if (!(result == resultbuf || result == NULL))
1949                                     free (result);
1950                                   if (buf_malloced != NULL)
1951                                     free (buf_malloced);
1952                                   CLEANUP ();
1953                                   errno = EILSEQ;
1954                                   return NULL;
1955                                 }
1956                               arg_end += count;
1957                               characters++;
1958                             }
1959                         }
1960                       else
1961                         {
1962                           /* Use the entire string.  */
1963                           arg_end = arg + u32_strlen (arg);
1964                           /* The number of characters doesn't matter.  */
1965                           characters = 0;
1966                         }
1967
1968                       if (has_width && width > characters
1969                           && !(dp->flags & FLAG_LEFT))
1970                         {
1971                           size_t n = width - characters;
1972                           ENSURE_ALLOCATION (xsum (length, n));
1973                           DCHAR_SET (result + length, ' ', n);
1974                           length += n;
1975                         }
1976
1977 # if DCHAR_IS_UINT32_T
1978                       {
1979                         size_t n = arg_end - arg;
1980                         ENSURE_ALLOCATION (xsum (length, n));
1981                         DCHAR_CPY (result + length, arg, n);
1982                         length += n;
1983                       }
1984 # else
1985                       { /* Convert.  */
1986                         DCHAR_T *converted = result + length;
1987                         size_t converted_len = allocated - length;
1988 #  if DCHAR_IS_TCHAR
1989                         /* Convert from UTF-32 to locale encoding.  */
1990                         if (u32_conv_to_encoding (locale_charset (),
1991                                                   iconveh_question_mark,
1992                                                   arg, arg_end - arg, NULL,
1993                                                   &converted, &converted_len)
1994                             < 0)
1995 #  else
1996                         /* Convert from UTF-32 to UTF-8/UTF-16.  */
1997                         converted =
1998                           U32_TO_DCHAR (arg, arg_end - arg,
1999                                         converted, &converted_len);
2000                         if (converted == NULL)
2001 #  endif
2002                           {
2003                             int saved_errno = errno;
2004                             if (!(result == resultbuf || result == NULL))
2005                               free (result);
2006                             if (buf_malloced != NULL)
2007                               free (buf_malloced);
2008                             CLEANUP ();
2009                             errno = saved_errno;
2010                             return NULL;
2011                           }
2012                         if (converted != result + length)
2013                           {
2014                             ENSURE_ALLOCATION (xsum (length, converted_len));
2015                             DCHAR_CPY (result + length, converted, converted_len);
2016                             free (converted);
2017                           }
2018                         length += converted_len;
2019                       }
2020 # endif
2021
2022                       if (has_width && width > characters
2023                           && (dp->flags & FLAG_LEFT))
2024                         {
2025                           size_t n = width - characters;
2026                           ENSURE_ALLOCATION (xsum (length, n));
2027                           DCHAR_SET (result + length, ' ', n);
2028                           length += n;
2029                         }
2030                     }
2031                     break;
2032
2033                   default:
2034                     abort ();
2035                   }
2036               }
2037 #endif
2038 #if (NEED_PRINTF_DIRECTIVE_A || NEED_PRINTF_LONG_DOUBLE || NEED_PRINTF_DOUBLE) && !defined IN_LIBINTL
2039             else if ((dp->conversion == 'a' || dp->conversion == 'A')
2040 # if !(NEED_PRINTF_DIRECTIVE_A || (NEED_PRINTF_LONG_DOUBLE && NEED_PRINTF_DOUBLE))
2041                      && (0
2042 #  if NEED_PRINTF_DOUBLE
2043                          || a.arg[dp->arg_index].type == TYPE_DOUBLE
2044 #  endif
2045 #  if NEED_PRINTF_LONG_DOUBLE
2046                          || a.arg[dp->arg_index].type == TYPE_LONGDOUBLE
2047 #  endif
2048                         )
2049 # endif
2050                     )
2051               {
2052                 arg_type type = a.arg[dp->arg_index].type;
2053                 int flags = dp->flags;
2054                 int has_width;
2055                 size_t width;
2056                 int has_precision;
2057                 size_t precision;
2058                 size_t tmp_length;
2059                 DCHAR_T tmpbuf[700];
2060                 DCHAR_T *tmp;
2061                 DCHAR_T *pad_ptr;
2062                 DCHAR_T *p;
2063
2064                 has_width = 0;
2065                 width = 0;
2066                 if (dp->width_start != dp->width_end)
2067                   {
2068                     if (dp->width_arg_index != ARG_NONE)
2069                       {
2070                         int arg;
2071
2072                         if (!(a.arg[dp->width_arg_index].type == TYPE_INT))
2073                           abort ();
2074                         arg = a.arg[dp->width_arg_index].a.a_int;
2075                         if (arg < 0)
2076                           {
2077                             /* "A negative field width is taken as a '-' flag
2078                                 followed by a positive field width."  */
2079                             flags |= FLAG_LEFT;
2080                             width = (unsigned int) (-arg);
2081                           }
2082                         else
2083                           width = arg;
2084                       }
2085                     else
2086                       {
2087                         const FCHAR_T *digitp = dp->width_start;
2088
2089                         do
2090                           width = xsum (xtimes (width, 10), *digitp++ - '0');
2091                         while (digitp != dp->width_end);
2092                       }
2093                     has_width = 1;
2094                   }
2095
2096                 has_precision = 0;
2097                 precision = 0;
2098                 if (dp->precision_start != dp->precision_end)
2099                   {
2100                     if (dp->precision_arg_index != ARG_NONE)
2101                       {
2102                         int arg;
2103
2104                         if (!(a.arg[dp->precision_arg_index].type == TYPE_INT))
2105                           abort ();
2106                         arg = a.arg[dp->precision_arg_index].a.a_int;
2107                         /* "A negative precision is taken as if the precision
2108                             were omitted."  */
2109                         if (arg >= 0)
2110                           {
2111                             precision = arg;
2112                             has_precision = 1;
2113                           }
2114                       }
2115                     else
2116                       {
2117                         const FCHAR_T *digitp = dp->precision_start + 1;
2118
2119                         precision = 0;
2120                         while (digitp != dp->precision_end)
2121                           precision = xsum (xtimes (precision, 10), *digitp++ - '0');
2122                         has_precision = 1;
2123                       }
2124                   }
2125
2126                 /* Allocate a temporary buffer of sufficient size.  */
2127                 if (type == TYPE_LONGDOUBLE)
2128                   tmp_length =
2129                     (unsigned int) ((LDBL_DIG + 1)
2130                                     * 0.831 /* decimal -> hexadecimal */
2131                                    )
2132                     + 1; /* turn floor into ceil */
2133                 else
2134                   tmp_length =
2135                     (unsigned int) ((DBL_DIG + 1)
2136                                     * 0.831 /* decimal -> hexadecimal */
2137                                    )
2138                     + 1; /* turn floor into ceil */
2139                 if (tmp_length < precision)
2140                   tmp_length = precision;
2141                 /* Account for sign, decimal point etc. */
2142                 tmp_length = xsum (tmp_length, 12);
2143
2144                 if (tmp_length < width)
2145                   tmp_length = width;
2146
2147                 tmp_length = xsum (tmp_length, 1); /* account for trailing NUL */
2148
2149                 if (tmp_length <= sizeof (tmpbuf) / sizeof (DCHAR_T))
2150                   tmp = tmpbuf;
2151                 else
2152                   {
2153                     size_t tmp_memsize = xtimes (tmp_length, sizeof (DCHAR_T));
2154
2155                     if (size_overflow_p (tmp_memsize))
2156                       /* Overflow, would lead to out of memory.  */
2157                       goto out_of_memory;
2158                     tmp = (DCHAR_T *) malloc (tmp_memsize);
2159                     if (tmp == NULL)
2160                       /* Out of memory.  */
2161                       goto out_of_memory;
2162                   }
2163
2164                 pad_ptr = NULL;
2165                 p = tmp;
2166                 if (type == TYPE_LONGDOUBLE)
2167                   {
2168 # if NEED_PRINTF_DIRECTIVE_A || NEED_PRINTF_LONG_DOUBLE
2169                     long double arg = a.arg[dp->arg_index].a.a_longdouble;
2170
2171                     if (isnanl (arg))
2172                       {
2173                         if (dp->conversion == 'A')
2174                           {
2175                             *p++ = 'N'; *p++ = 'A'; *p++ = 'N';
2176                           }
2177                         else
2178                           {
2179                             *p++ = 'n'; *p++ = 'a'; *p++ = 'n';
2180                           }
2181                       }
2182                     else
2183                       {
2184                         int sign = 0;
2185                         DECL_LONG_DOUBLE_ROUNDING
2186
2187                         BEGIN_LONG_DOUBLE_ROUNDING ();
2188
2189                         if (signbit (arg)) /* arg < 0.0L or negative zero */
2190                           {
2191                             sign = -1;
2192                             arg = -arg;
2193                           }
2194
2195                         if (sign < 0)
2196                           *p++ = '-';
2197                         else if (flags & FLAG_SHOWSIGN)
2198                           *p++ = '+';
2199                         else if (flags & FLAG_SPACE)
2200                           *p++ = ' ';
2201
2202                         if (arg > 0.0L && arg + arg == arg)
2203                           {
2204                             if (dp->conversion == 'A')
2205                               {
2206                                 *p++ = 'I'; *p++ = 'N'; *p++ = 'F';
2207                               }
2208                             else
2209                               {
2210                                 *p++ = 'i'; *p++ = 'n'; *p++ = 'f';
2211                               }
2212                           }
2213                         else
2214                           {
2215                             int exponent;
2216                             long double mantissa;
2217
2218                             if (arg > 0.0L)
2219                               mantissa = printf_frexpl (arg, &exponent);
2220                             else
2221                               {
2222                                 exponent = 0;
2223                                 mantissa = 0.0L;
2224                               }
2225
2226                             if (has_precision
2227                                 && precision < (unsigned int) ((LDBL_DIG + 1) * 0.831) + 1)
2228                               {
2229                                 /* Round the mantissa.  */
2230                                 long double tail = mantissa;
2231                                 size_t q;
2232
2233                                 for (q = precision; ; q--)
2234                                   {
2235                                     int digit = (int) tail;
2236                                     tail -= digit;
2237                                     if (q == 0)
2238                                       {
2239                                         if (digit & 1 ? tail >= 0.5L : tail > 0.5L)
2240                                           tail = 1 - tail;
2241                                         else
2242                                           tail = - tail;
2243                                         break;
2244                                       }
2245                                     tail *= 16.0L;
2246                                   }
2247                                 if (tail != 0.0L)
2248                                   for (q = precision; q > 0; q--)
2249                                     tail *= 0.0625L;
2250                                 mantissa += tail;
2251                               }
2252
2253                             *p++ = '0';
2254                             *p++ = dp->conversion - 'A' + 'X';
2255                             pad_ptr = p;
2256                             {
2257                               int digit;
2258
2259                               digit = (int) mantissa;
2260                               mantissa -= digit;
2261                               *p++ = '0' + digit;
2262                               if ((flags & FLAG_ALT)
2263                                   || mantissa > 0.0L || precision > 0)
2264                                 {
2265                                   *p++ = decimal_point_char ();
2266                                   /* This loop terminates because we assume
2267                                      that FLT_RADIX is a power of 2.  */
2268                                   while (mantissa > 0.0L)
2269                                     {
2270                                       mantissa *= 16.0L;
2271                                       digit = (int) mantissa;
2272                                       mantissa -= digit;
2273                                       *p++ = digit
2274                                              + (digit < 10
2275                                                 ? '0'
2276                                                 : dp->conversion - 10);
2277                                       if (precision > 0)
2278                                         precision--;
2279                                     }
2280                                   while (precision > 0)
2281                                     {
2282                                       *p++ = '0';
2283                                       precision--;
2284                                     }
2285                                 }
2286                               }
2287                               *p++ = dp->conversion - 'A' + 'P';
2288 #  if WIDE_CHAR_VERSION
2289                               {
2290                                 static const wchar_t decimal_format[] =
2291                                   { '%', '+', 'd', '\0' };
2292                                 SNPRINTF (p, 6 + 1, decimal_format, exponent);
2293                               }
2294                               while (*p != '\0')
2295                                 p++;
2296 #  else
2297                               if (sizeof (DCHAR_T) == 1)
2298                                 {
2299                                   sprintf ((char *) p, "%+d", exponent);
2300                                   while (*p != '\0')
2301                                     p++;
2302                                 }
2303                               else
2304                                 {
2305                                   char expbuf[6 + 1];
2306                                   const char *ep;
2307                                   sprintf (expbuf, "%+d", exponent);
2308                                   for (ep = expbuf; (*p = *ep) != '\0'; ep++)
2309                                     p++;
2310                                 }
2311 #  endif
2312                           }
2313
2314                         END_LONG_DOUBLE_ROUNDING ();
2315                       }
2316 # else
2317                     abort ();
2318 # endif
2319                   }
2320                 else
2321                   {
2322 # if NEED_PRINTF_DIRECTIVE_A || NEED_PRINTF_DOUBLE
2323                     double arg = a.arg[dp->arg_index].a.a_double;
2324
2325                     if (isnand (arg))
2326                       {
2327                         if (dp->conversion == 'A')
2328                           {
2329                             *p++ = 'N'; *p++ = 'A'; *p++ = 'N';
2330                           }
2331                         else
2332                           {
2333                             *p++ = 'n'; *p++ = 'a'; *p++ = 'n';
2334                           }
2335                       }
2336                     else
2337                       {
2338                         int sign = 0;
2339
2340                         if (signbit (arg)) /* arg < 0.0 or negative zero */
2341                           {
2342                             sign = -1;
2343                             arg = -arg;
2344                           }
2345
2346                         if (sign < 0)
2347                           *p++ = '-';
2348                         else if (flags & FLAG_SHOWSIGN)
2349                           *p++ = '+';
2350                         else if (flags & FLAG_SPACE)
2351                           *p++ = ' ';
2352
2353                         if (arg > 0.0 && arg + arg == arg)
2354                           {
2355                             if (dp->conversion == 'A')
2356                               {
2357                                 *p++ = 'I'; *p++ = 'N'; *p++ = 'F';
2358                               }
2359                             else
2360                               {
2361                                 *p++ = 'i'; *p++ = 'n'; *p++ = 'f';
2362                               }
2363                           }
2364                         else
2365                           {
2366                             int exponent;
2367                             double mantissa;
2368
2369                             if (arg > 0.0)
2370                               mantissa = printf_frexp (arg, &exponent);
2371                             else
2372                               {
2373                                 exponent = 0;
2374                                 mantissa = 0.0;
2375                               }
2376
2377                             if (has_precision
2378                                 && precision < (unsigned int) ((DBL_DIG + 1) * 0.831) + 1)
2379                               {
2380                                 /* Round the mantissa.  */
2381                                 double tail = mantissa;
2382                                 size_t q;
2383
2384                                 for (q = precision; ; q--)
2385                                   {
2386                                     int digit = (int) tail;
2387                                     tail -= digit;
2388                                     if (q == 0)
2389                                       {
2390                                         if (digit & 1 ? tail >= 0.5 : tail > 0.5)
2391                                           tail = 1 - tail;
2392                                         else
2393                                           tail = - tail;
2394                                         break;
2395                                       }
2396                                     tail *= 16.0;
2397                                   }
2398                                 if (tail != 0.0)
2399                                   for (q = precision; q > 0; q--)
2400                                     tail *= 0.0625;
2401                                 mantissa += tail;
2402                               }
2403
2404                             *p++ = '0';
2405                             *p++ = dp->conversion - 'A' + 'X';
2406                             pad_ptr = p;
2407                             {
2408                               int digit;
2409
2410                               digit = (int) mantissa;
2411                               mantissa -= digit;
2412                               *p++ = '0' + digit;
2413                               if ((flags & FLAG_ALT)
2414                                   || mantissa > 0.0 || precision > 0)
2415                                 {
2416                                   *p++ = decimal_point_char ();
2417                                   /* This loop terminates because we assume
2418                                      that FLT_RADIX is a power of 2.  */
2419                                   while (mantissa > 0.0)
2420                                     {
2421                                       mantissa *= 16.0;
2422                                       digit = (int) mantissa;
2423                                       mantissa -= digit;
2424                                       *p++ = digit
2425                                              + (digit < 10
2426                                                 ? '0'
2427                                                 : dp->conversion - 10);
2428                                       if (precision > 0)
2429                                         precision--;
2430                                     }
2431                                   while (precision > 0)
2432                                     {
2433                                       *p++ = '0';
2434                                       precision--;
2435                                     }
2436                                 }
2437                               }
2438                               *p++ = dp->conversion - 'A' + 'P';
2439 #  if WIDE_CHAR_VERSION
2440                               {
2441                                 static const wchar_t decimal_format[] =
2442                                   { '%', '+', 'd', '\0' };
2443                                 SNPRINTF (p, 6 + 1, decimal_format, exponent);
2444                               }
2445                               while (*p != '\0')
2446                                 p++;
2447 #  else
2448                               if (sizeof (DCHAR_T) == 1)
2449                                 {
2450                                   sprintf ((char *) p, "%+d", exponent);
2451                                   while (*p != '\0')
2452                                     p++;
2453                                 }
2454                               else
2455                                 {
2456                                   char expbuf[6 + 1];
2457                                   const char *ep;
2458                                   sprintf (expbuf, "%+d", exponent);
2459                                   for (ep = expbuf; (*p = *ep) != '\0'; ep++)
2460                                     p++;
2461                                 }
2462 #  endif
2463                           }
2464                       }
2465 # else
2466                     abort ();
2467 # endif
2468                   }
2469                 /* The generated string now extends from tmp to p, with the
2470                    zero padding insertion point being at pad_ptr.  */
2471                 if (has_width && p - tmp < width)
2472                   {
2473                     size_t pad = width - (p - tmp);
2474                     DCHAR_T *end = p + pad;
2475
2476                     if (flags & FLAG_LEFT)
2477                       {
2478                         /* Pad with spaces on the right.  */
2479                         for (; pad > 0; pad--)
2480                           *p++ = ' ';
2481                       }
2482                     else if ((flags & FLAG_ZERO) && pad_ptr != NULL)
2483                       {
2484                         /* Pad with zeroes.  */
2485                         DCHAR_T *q = end;
2486
2487                         while (p > pad_ptr)
2488                           *--q = *--p;
2489                         for (; pad > 0; pad--)
2490                           *p++ = '0';
2491                       }
2492                     else
2493                       {
2494                         /* Pad with spaces on the left.  */
2495                         DCHAR_T *q = end;
2496
2497                         while (p > tmp)
2498                           *--q = *--p;
2499                         for (; pad > 0; pad--)
2500                           *p++ = ' ';
2501                       }
2502
2503                     p = end;
2504                   }
2505
2506                 {
2507                   size_t count = p - tmp;
2508
2509                   if (count >= tmp_length)
2510                     /* tmp_length was incorrectly calculated - fix the
2511                        code above!  */
2512                     abort ();
2513
2514                   /* Make room for the result.  */
2515                   if (count >= allocated - length)
2516                     {
2517                       size_t n = xsum (length, count);
2518
2519                       ENSURE_ALLOCATION (n);
2520                     }
2521
2522                   /* Append the result.  */
2523                   memcpy (result + length, tmp, count * sizeof (DCHAR_T));
2524                   if (tmp != tmpbuf)
2525                     free (tmp);
2526                   length += count;
2527                 }
2528               }
2529 #endif
2530 #if (NEED_PRINTF_INFINITE_DOUBLE || NEED_PRINTF_DOUBLE || NEED_PRINTF_INFINITE_LONG_DOUBLE || NEED_PRINTF_LONG_DOUBLE) && !defined IN_LIBINTL
2531             else if ((dp->conversion == 'f' || dp->conversion == 'F'
2532                       || dp->conversion == 'e' || dp->conversion == 'E'
2533                       || dp->conversion == 'g' || dp->conversion == 'G'
2534                       || dp->conversion == 'a' || dp->conversion == 'A')
2535                      && (0
2536 # if NEED_PRINTF_DOUBLE
2537                          || a.arg[dp->arg_index].type == TYPE_DOUBLE
2538 # elif NEED_PRINTF_INFINITE_DOUBLE
2539                          || (a.arg[dp->arg_index].type == TYPE_DOUBLE
2540                              /* The systems (mingw) which produce wrong output
2541                                 for Inf, -Inf, and NaN also do so for -0.0.
2542                                 Therefore we treat this case here as well.  */
2543                              && is_infinite_or_zero (a.arg[dp->arg_index].a.a_double))
2544 # endif
2545 # if NEED_PRINTF_LONG_DOUBLE
2546                          || a.arg[dp->arg_index].type == TYPE_LONGDOUBLE
2547 # elif NEED_PRINTF_INFINITE_LONG_DOUBLE
2548                          || (a.arg[dp->arg_index].type == TYPE_LONGDOUBLE
2549                              /* Some systems produce wrong output for Inf,
2550                                 -Inf, and NaN.  */
2551                              && is_infinitel (a.arg[dp->arg_index].a.a_longdouble))
2552 # endif
2553                         ))
2554               {
2555 # if (NEED_PRINTF_DOUBLE || NEED_PRINTF_INFINITE_DOUBLE) && (NEED_PRINTF_LONG_DOUBLE || NEED_PRINTF_INFINITE_LONG_DOUBLE)
2556                 arg_type type = a.arg[dp->arg_index].type;
2557 # endif
2558                 int flags = dp->flags;
2559                 int has_width;
2560                 size_t width;
2561                 int has_precision;
2562                 size_t precision;
2563                 size_t tmp_length;
2564                 DCHAR_T tmpbuf[700];
2565                 DCHAR_T *tmp;
2566                 DCHAR_T *pad_ptr;
2567                 DCHAR_T *p;
2568
2569                 has_width = 0;
2570                 width = 0;
2571                 if (dp->width_start != dp->width_end)
2572                   {
2573                     if (dp->width_arg_index != ARG_NONE)
2574                       {
2575                         int arg;
2576
2577                         if (!(a.arg[dp->width_arg_index].type == TYPE_INT))
2578                           abort ();
2579                         arg = a.arg[dp->width_arg_index].a.a_int;
2580                         if (arg < 0)
2581                           {
2582                             /* "A negative field width is taken as a '-' flag
2583                                 followed by a positive field width."  */
2584                             flags |= FLAG_LEFT;
2585                             width = (unsigned int) (-arg);
2586                           }
2587                         else
2588                           width = arg;
2589                       }
2590                     else
2591                       {
2592                         const FCHAR_T *digitp = dp->width_start;
2593
2594                         do
2595                           width = xsum (xtimes (width, 10), *digitp++ - '0');
2596                         while (digitp != dp->width_end);
2597                       }
2598                     has_width = 1;
2599                   }
2600
2601                 has_precision = 0;
2602                 precision = 0;
2603                 if (dp->precision_start != dp->precision_end)
2604                   {
2605                     if (dp->precision_arg_index != ARG_NONE)
2606                       {
2607                         int arg;
2608
2609                         if (!(a.arg[dp->precision_arg_index].type == TYPE_INT))
2610                           abort ();
2611                         arg = a.arg[dp->precision_arg_index].a.a_int;
2612                         /* "A negative precision is taken as if the precision
2613                             were omitted."  */
2614                         if (arg >= 0)
2615                           {
2616                             precision = arg;
2617                             has_precision = 1;
2618                           }
2619                       }
2620                     else
2621                       {
2622                         const FCHAR_T *digitp = dp->precision_start + 1;
2623
2624                         precision = 0;
2625                         while (digitp != dp->precision_end)
2626                           precision = xsum (xtimes (precision, 10), *digitp++ - '0');
2627                         has_precision = 1;
2628                       }
2629                   }
2630
2631                 /* POSIX specifies the default precision to be 6 for %f, %F,
2632                    %e, %E, but not for %g, %G.  Implementations appear to use
2633                    the same default precision also for %g, %G.  */
2634                 if (!has_precision)
2635                   precision = 6;
2636
2637                 /* Allocate a temporary buffer of sufficient size.  */
2638 # if NEED_PRINTF_DOUBLE && NEED_PRINTF_LONG_DOUBLE
2639                 tmp_length = (type == TYPE_LONGDOUBLE ? LDBL_DIG + 1 : DBL_DIG + 1);
2640 # elif NEED_PRINTF_INFINITE_DOUBLE && NEED_PRINTF_LONG_DOUBLE
2641                 tmp_length = (type == TYPE_LONGDOUBLE ? LDBL_DIG + 1 : 0);
2642 # elif NEED_PRINTF_LONG_DOUBLE
2643                 tmp_length = LDBL_DIG + 1;
2644 # elif NEED_PRINTF_DOUBLE
2645                 tmp_length = DBL_DIG + 1;
2646 # else
2647                 tmp_length = 0;
2648 # endif
2649                 if (tmp_length < precision)
2650                   tmp_length = precision;
2651 # if NEED_PRINTF_LONG_DOUBLE
2652 #  if NEED_PRINTF_DOUBLE || NEED_PRINTF_INFINITE_DOUBLE
2653                 if (type == TYPE_LONGDOUBLE)
2654 #  endif
2655                   if (dp->conversion == 'f' || dp->conversion == 'F')
2656                     {
2657                       long double arg = a.arg[dp->arg_index].a.a_longdouble;
2658                       if (!(isnanl (arg) || arg + arg == arg))
2659                         {
2660                           /* arg is finite and nonzero.  */
2661                           int exponent = floorlog10l (arg < 0 ? -arg : arg);
2662                           if (exponent >= 0 && tmp_length < exponent + precision)
2663                             tmp_length = exponent + precision;
2664                         }
2665                     }
2666 # endif
2667 # if NEED_PRINTF_DOUBLE
2668 #  if NEED_PRINTF_LONG_DOUBLE || NEED_PRINTF_INFINITE_LONG_DOUBLE
2669                 if (type == TYPE_DOUBLE)
2670 #  endif
2671                   if (dp->conversion == 'f' || dp->conversion == 'F')
2672                     {
2673                       double arg = a.arg[dp->arg_index].a.a_double;
2674                       if (!(isnand (arg) || arg + arg == arg))
2675                         {
2676                           /* arg is finite and nonzero.  */