1/* Copyright (c) 2007, 2012, Oracle and/or its affiliates. All rights reserved.
2 2016,2018 MariaDB Corporation AB
3
4 This library is free software; you can redistribute it and/or
5 modify it under the terms of the GNU Library General Public
6 License as published by the Free Software Foundation; version 2
7 of the License.
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 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */
17
18/****************************************************************
19
20 This file incorporates work covered by the following copyright and
21 permission notice:
22
23 The author of this software is David M. Gay.
24
25 Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
26
27 Permission to use, copy, modify, and distribute this software for any
28 purpose without fee is hereby granted, provided that this entire notice
29 is included in all copies of any software which is or includes a copy
30 or modification of this software and in all copies of the supporting
31 documentation for such software.
32
33 THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
34 WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
35 REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
36 OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
37
38 ***************************************************************/
39
40//#include "strings_def.h"
41//#include <my_base.h> /* for EOVERFLOW on Windows */
42#include <ma_global.h>
43#include <memory.h>
44#include "ma_string.h"
45
46/**
47 Appears to suffice to not call malloc() in most cases.
48 @todo
49 see if it is possible to get rid of malloc().
50 this constant is sufficient to avoid malloc() on all inputs I have tried.
51*/
52#define DTOA_BUFF_SIZE (460 * sizeof(void *))
53
54/* Magic value returned by dtoa() to indicate overflow */
55#define DTOA_OVERFLOW 9999
56
57static char *dtoa(double, int, int, int *, int *, char **, char *, size_t);
58static void dtoa_free(char *, char *, size_t);
59
60/**
61 @brief
62 Converts a given floating point number to a zero-terminated string
63 representation using the 'f' format.
64
65 @details
66 This function is a wrapper around dtoa() to do the same as
67 sprintf(to, "%-.*f", precision, x), though the conversion is usually more
68 precise. The only difference is in handling [-,+]infinity and nan values,
69 in which case we print '0\0' to the output string and indicate an overflow.
70
71 @param x the input floating point number.
72 @param precision the number of digits after the decimal point.
73 All properties of sprintf() apply:
74 - if the number of significant digits after the decimal
75 point is less than precision, the resulting string is
76 right-padded with zeros
77 - if the precision is 0, no decimal point appears
78 - if a decimal point appears, at least one digit appears
79 before it
80 @param to pointer to the output buffer. The longest string which
81 my_fcvt() can return is FLOATING_POINT_BUFFER bytes
82 (including the terminating '\0').
83 @param error if not NULL, points to a location where the status of
84 conversion is stored upon return.
85 FALSE successful conversion
86 TRUE the input number is [-,+]infinity or nan.
87 The output string in this case is always '0'.
88 @return number of written characters (excluding terminating '\0')
89*/
90
91size_t ma_fcvt(double x, int precision, char *to, my_bool *error)
92{
93 int decpt, sign, len, i;
94 char *res, *src, *end, *dst= to;
95 char buf[DTOA_BUFF_SIZE];
96 DBUG_ASSERT(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
97
98 res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
99
100 if (decpt == DTOA_OVERFLOW)
101 {
102 dtoa_free(res, buf, sizeof(buf));
103 *to++= '0';
104 *to= '\0';
105 if (error != NULL)
106 *error= TRUE;
107 return 1;
108 }
109
110 src= res;
111 len= (int)(end - src);
112
113 if (sign)
114 *dst++= '-';
115
116 if (decpt <= 0)
117 {
118 *dst++= '0';
119 *dst++= '.';
120 for (i= decpt; i < 0; i++)
121 *dst++= '0';
122 }
123
124 for (i= 1; i <= len; i++)
125 {
126 *dst++= *src++;
127 if (i == decpt && i < len)
128 *dst++= '.';
129 }
130 while (i++ <= decpt)
131 *dst++= '0';
132
133 if (precision > 0)
134 {
135 if (len <= decpt)
136 *dst++= '.';
137
138 for (i= precision - MAX(0, (len - decpt)); i > 0; i--)
139 *dst++= '0';
140 }
141
142 *dst= '\0';
143 if (error != NULL)
144 *error= FALSE;
145
146 dtoa_free(res, buf, sizeof(buf));
147
148 return dst - to;
149}
150
151/**
152 @brief
153 Converts a given floating point number to a zero-terminated string
154 representation with a given field width using the 'e' format
155 (aka scientific notation) or the 'f' one.
156
157 @details
158 The format is chosen automatically to provide the most number of significant
159 digits (and thus, precision) with a given field width. In many cases, the
160 result is similar to that of sprintf(to, "%g", x) with a few notable
161 differences:
162 - the conversion is usually more precise than C library functions.
163 - there is no 'precision' argument. instead, we specify the number of
164 characters available for conversion (i.e. a field width).
165 - the result never exceeds the specified field width. If the field is too
166 short to contain even a rounded decimal representation, ma_gcvt()
167 indicates overflow and truncates the output string to the specified width.
168 - float-type arguments are handled differently than double ones. For a
169 float input number (i.e. when the 'type' argument is MY_GCVT_ARG_FLOAT)
170 we deliberately limit the precision of conversion by FLT_DIG digits to
171 avoid garbage past the significant digits.
172 - unlike sprintf(), in cases where the 'e' format is preferred, we don't
173 zero-pad the exponent to save space for significant digits. The '+' sign
174 for a positive exponent does not appear for the same reason.
175
176 @param x the input floating point number.
177 @param type is either MY_GCVT_ARG_FLOAT or MY_GCVT_ARG_DOUBLE.
178 Specifies the type of the input number (see notes above).
179 @param width field width in characters. The minimal field width to
180 hold any number representation (albeit rounded) is 7
181 characters ("-Ne-NNN").
182 @param to pointer to the output buffer. The result is always
183 zero-terminated, and the longest returned string is thus
184 'width + 1' bytes.
185 @param error if not NULL, points to a location where the status of
186 conversion is stored upon return.
187 FALSE successful conversion
188 TRUE the input number is [-,+]infinity or nan.
189 The output string in this case is always '0'.
190 @return number of written characters (excluding terminating '\0')
191
192 @todo
193 Check if it is possible and makes sense to do our own rounding on top of
194 dtoa() instead of calling dtoa() twice in (rare) cases when the resulting
195 string representation does not fit in the specified field width and we want
196 to re-round the input number with fewer significant digits. Examples:
197
198 ma_gcvt(-9e-3, ..., 4, ...);
199 ma_gcvt(-9e-3, ..., 2, ...);
200 ma_gcvt(1.87e-3, ..., 4, ...);
201 ma_gcvt(55, ..., 1, ...);
202
203 We do our best to minimize such cases by:
204
205 - passing to dtoa() the field width as the number of significant digits
206
207 - removing the sign of the number early (and decreasing the width before
208 passing it to dtoa())
209
210 - choosing the proper format to preserve the most number of significant
211 digits.
212*/
213
214size_t ma_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
215 my_bool *error)
216{
217 int decpt, sign, len, exp_len;
218 char *res, *src, *end, *dst= to, *dend= dst + width;
219 char buf[DTOA_BUFF_SIZE];
220 my_bool have_space, force_e_format;
221 DBUG_ASSERT(width > 0 && to != NULL);
222
223 /* We want to remove '-' from equations early */
224 if (x < 0.)
225 width--;
226
227 res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : MIN(width, FLT_DIG),
228 &decpt, &sign, &end, buf, sizeof(buf));
229 if (decpt == DTOA_OVERFLOW)
230 {
231 dtoa_free(res, buf, sizeof(buf));
232 *to++= '0';
233 *to= '\0';
234 if (error != NULL)
235 *error= TRUE;
236 return 1;
237 }
238
239 if (error != NULL)
240 *error= FALSE;
241
242 src= res;
243 len= (int)(end - res);
244
245 /*
246 Number of digits in the exponent from the 'e' conversion.
247 The sign of the exponent is taken into account separetely, we don't need
248 to count it here.
249 */
250 exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
251
252 /*
253 Do we have enough space for all digits in the 'f' format?
254 Let 'len' be the number of significant digits returned by dtoa,
255 and F be the length of the resulting decimal representation.
256 Consider the following cases:
257 1. decpt <= 0, i.e. we have "0.NNN" => F = len - decpt + 2
258 2. 0 < decpt < len, i.e. we have "NNN.NNN" => F = len + 1
259 3. len <= decpt, i.e. we have "NNN00" => F = decpt
260 */
261 have_space= (decpt <= 0 ? len - decpt + 2 :
262 decpt > 0 && decpt < len ? len + 1 :
263 decpt) <= width;
264 /*
265 The following is true when no significant digits can be placed with the
266 specified field width using the 'f' format, and the 'e' format
267 will not be truncated.
268 */
269 force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len);
270 /*
271 Assume that we don't have enough space to place all significant digits in
272 the 'f' format. We have to choose between the 'e' format and the 'f' one
273 to keep as many significant digits as possible.
274 Let E and F be the lengths of decimal representation in the 'e' and 'f'
275 formats, respectively. We want to use the 'f' format if, and only if F <= E.
276 Consider the following cases:
277 1. decpt <= 0.
278 F = len - decpt + 2 (see above)
279 E = len + (len > 1) + 1 + 1 (decpt <= -99) + (decpt <= -9) + 1
280 ("N.NNe-MMM")
281 (F <= E) <=> (len == 1 && decpt >= -1) || (len > 1 && decpt >= -2)
282 We also need to ensure that if the 'f' format is chosen,
283 the field width allows us to place at least one significant digit
284 (i.e. width > 2 - decpt). If not, we prefer the 'e' format.
285 2. 0 < decpt < len
286 F = len + 1 (see above)
287 E = len + 1 + 1 + ... ("N.NNeMMM")
288 F is always less than E.
289 3. len <= decpt <= width
290 In this case we have enough space to represent the number in the 'f'
291 format, so we prefer it with some exceptions.
292 4. width < decpt
293 The number cannot be represented in the 'f' format at all, always use
294 the 'e' 'one.
295 */
296 if ((have_space ||
297 /*
298 Not enough space, let's see if the 'f' format provides the most number
299 of significant digits.
300 */
301 ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
302 (len > 1 || !force_e_format)))) &&
303 !force_e_format)) &&
304
305 /*
306 Use the 'e' format in some cases even if we have enough space for the
307 'f' one. See comment for DBL_DIG.
308 */
309 (!have_space || (decpt >= -DBL_DIG + 1 &&
310 (decpt <= DBL_DIG || len > decpt))))
311 {
312 /* 'f' format */
313 int i;
314
315 width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0);
316
317 /* Do we have to truncate any digits? */
318 if (width < len)
319 {
320 if (width < decpt)
321 {
322 if (error != NULL)
323 *error= TRUE;
324 width= decpt;
325 }
326
327 /*
328 We want to truncate (len - width) least significant digits after the
329 decimal point. For this we are calling dtoa with mode=5, passing the
330 number of significant digits = (len-decpt) - (len-width) = width-decpt
331 */
332 dtoa_free(res, buf, sizeof(buf));
333 res= dtoa(x, 5, width - decpt, &decpt, &sign, &end, buf, sizeof(buf));
334 src= res;
335 len= (int)(end - res);
336 }
337
338 if (len == 0)
339 {
340 /* Underflow. Just print '0' and exit */
341 *dst++= '0';
342 goto end;
343 }
344
345 /*
346 At this point we are sure we have enough space to put all digits
347 returned by dtoa
348 */
349 if (sign && dst < dend)
350 *dst++= '-';
351 if (decpt <= 0)
352 {
353 if (dst < dend)
354 *dst++= '0';
355 if (len > 0 && dst < dend)
356 *dst++= '.';
357 for (; decpt < 0 && dst < dend; decpt++)
358 *dst++= '0';
359 }
360
361 for (i= 1; i <= len && dst < dend; i++)
362 {
363 *dst++= *src++;
364 if (i == decpt && i < len && dst < dend)
365 *dst++= '.';
366 }
367 while (i++ <= decpt && dst < dend)
368 *dst++= '0';
369 }
370 else
371 {
372 /* 'e' format */
373 int decpt_sign= 0;
374
375 if (--decpt < 0)
376 {
377 decpt= -decpt;
378 width--;
379 decpt_sign= 1;
380 }
381 width-= 1 + exp_len; /* eNNN */
382
383 if (len > 1)
384 width--;
385
386 if (width <= 0)
387 {
388 /* Overflow */
389 if (error != NULL)
390 *error= TRUE;
391 width= 0;
392 }
393
394 /* Do we have to truncate any digits? */
395 if (width < len)
396 {
397 /* Yes, re-convert with a smaller width */
398 dtoa_free(res, buf, sizeof(buf));
399 res= dtoa(x, 4, width, &decpt, &sign, &end, buf, sizeof(buf));
400 src= res;
401 len= (int)(end - res);
402 if (--decpt < 0)
403 decpt= -decpt;
404 }
405 /*
406 At this point we are sure we have enough space to put all digits
407 returned by dtoa
408 */
409 if (sign && dst < dend)
410 *dst++= '-';
411 if (dst < dend)
412 *dst++= *src++;
413 if (len > 1 && dst < dend)
414 {
415 *dst++= '.';
416 while (src < end && dst < dend)
417 *dst++= *src++;
418 }
419 if (dst < dend)
420 *dst++= 'e';
421 if (decpt_sign && dst < dend)
422 *dst++= '-';
423
424 if (decpt >= 100 && dst < dend)
425 {
426 *dst++= decpt / 100 + '0';
427 decpt%= 100;
428 if (dst < dend)
429 *dst++= decpt / 10 + '0';
430 }
431 else if (decpt >= 10 && dst < dend)
432 *dst++= decpt / 10 + '0';
433 if (dst < dend)
434 *dst++= decpt % 10 + '0';
435
436 }
437
438end:
439 dtoa_free(res, buf, sizeof(buf));
440 *dst= '\0';
441
442 return dst - to;
443}
444
445/****************************************************************
446 *
447 * The author of this software is David M. Gay.
448 *
449 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
450 *
451 * Permission to use, copy, modify, and distribute this software for any
452 * purpose without fee is hereby granted, provided that this entire notice
453 * is included in all copies of any software which is or includes a copy
454 * or modification of this software and in all copies of the supporting
455 * documentation for such software.
456 *
457 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
458 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
459 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
460 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
461 *
462 ***************************************************************/
463/* Please send bug reports to David M. Gay (dmg at acm dot org,
464 * with " at " changed at "@" and " dot " changed to "."). */
465
466/*
467 Original copy of the software is located at http://www.netlib.org/fp/dtoa.c
468 It was adjusted to serve MySQL server needs:
469 * strtod() was modified to not expect a zero-terminated string.
470 It now honors 'se' (end of string) argument as the input parameter,
471 not just as the output one.
472 * in dtoa(), in case of overflow/underflow/NaN result string now contains "0";
473 decpt is set to DTOA_OVERFLOW to indicate overflow.
474 * support for VAX, IBM mainframe and 16-bit hardware removed
475 * we always assume that 64-bit integer type is available
476 * support for Kernigan-Ritchie style headers (pre-ANSI compilers)
477 removed
478 * all gcc warnings ironed out
479 * we always assume multithreaded environment, so we had to change
480 memory allocation procedures to use stack in most cases;
481 malloc is used as the last resort.
482 * pow5mult rewritten to use pre-calculated pow5 list instead of
483 the one generated on the fly.
484*/
485
486
487/*
488 On a machine with IEEE extended-precision registers, it is
489 necessary to specify double-precision (53-bit) rounding precision
490 before invoking strtod or dtoa. If the machine uses (the equivalent
491 of) Intel 80x87 arithmetic, the call
492 _control87(PC_53, MCW_PC);
493 does this with many compilers. Whether this or another call is
494 appropriate depends on the compiler; for this to work, it may be
495 necessary to #include "float.h" or another system-dependent header
496 file.
497*/
498
499/*
500 #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
501 and dtoa should round accordingly.
502 #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
503 and Honor_FLT_ROUNDS is not #defined.
504
505 TODO: check if we can get rid of the above two
506*/
507
508typedef int32 Long;
509typedef uint32 ULong;
510typedef int64 LLong;
511typedef uint64 ULLong;
512
513typedef union { double d; ULong L[2]; } U;
514
515#if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) && \
516 (__FLOAT_WORD_ORDER == __BIG_ENDIAN))
517#define word0(x) (x)->L[0]
518#define word1(x) (x)->L[1]
519#else
520#define word0(x) (x)->L[1]
521#define word1(x) (x)->L[0]
522#endif
523
524#define dval(x) (x)->d
525
526/* #define P DBL_MANT_DIG */
527/* Ten_pmax= floor(P*log(2)/log(5)) */
528/* Bletch= (highest power of 2 < DBL_MAX_10_EXP) / 16 */
529/* Quick_max= floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
530/* Int_max= floor(P*log(FLT_RADIX)/log(10) - 1) */
531
532#define Exp_shift 20
533#define Exp_shift1 20
534#define Exp_msk1 0x100000
535#define Exp_mask 0x7ff00000
536#define P 53
537#define Bias 1023
538#define Emin (-1022)
539#define Exp_1 0x3ff00000
540#define Exp_11 0x3ff00000
541#define Ebits 11
542#define Frac_mask 0xfffff
543#define Frac_mask1 0xfffff
544#define Ten_pmax 22
545#define Bletch 0x10
546#define Bndry_mask 0xfffff
547#define Bndry_mask1 0xfffff
548#define LSB 1
549#define Sign_bit 0x80000000
550#define Log2P 1
551#define Tiny1 1
552#define Quick_max 14
553#define Int_max 14
554
555#ifndef Flt_Rounds
556#ifdef FLT_ROUNDS
557#define Flt_Rounds FLT_ROUNDS
558#else
559#define Flt_Rounds 1
560#endif
561#endif /*Flt_Rounds*/
562
563#ifdef Honor_FLT_ROUNDS
564#define Rounding rounding
565#undef Check_FLT_ROUNDS
566#define Check_FLT_ROUNDS
567#else
568#define Rounding Flt_Rounds
569#endif
570
571#define rounded_product(a,b) a*= b
572#define rounded_quotient(a,b) a/= b
573
574#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
575#define Big1 0xffffffff
576#define FFFFFFFF 0xffffffffUL
577
578/* This is tested to be enough for dtoa */
579
580#define Kmax 15
581
582#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
583 2*sizeof(int) + y->wds*sizeof(ULong))
584
585/* Arbitrary-length integer */
586
587typedef struct Bigint
588{
589 union {
590 ULong *x; /* points right after this Bigint object */
591 struct Bigint *next; /* to maintain free lists */
592 } p;
593 int k; /* 2^k = maxwds */
594 int maxwds; /* maximum length in 32-bit words */
595 int sign; /* not zero if number is negative */
596 int wds; /* current length in 32-bit words */
597} Bigint;
598
599
600/* A simple stack-memory based allocator for Bigints */
601
602typedef struct Stack_alloc
603{
604 char *begin;
605 char *free;
606 char *end;
607 /*
608 Having list of free blocks lets us reduce maximum required amount
609 of memory from ~4000 bytes to < 1680 (tested on x86).
610 */
611 Bigint *freelist[Kmax+1];
612} Stack_alloc;
613
614
615/*
616 Try to allocate object on stack, and resort to malloc if all
617 stack memory is used. Ensure allocated objects to be aligned by the pointer
618 size in order to not break the alignment rules when storing a pointer to a
619 Bigint.
620*/
621
622static Bigint *Balloc(int k, Stack_alloc *alloc)
623{
624 Bigint *rv;
625 DBUG_ASSERT(k <= Kmax);
626 if (k <= Kmax && alloc->freelist[k])
627 {
628 rv= alloc->freelist[k];
629 alloc->freelist[k]= rv->p.next;
630 }
631 else
632 {
633 int x, len;
634
635 x= 1 << k;
636 len= MY_ALIGN(sizeof(Bigint) + x * sizeof(ULong), SIZEOF_CHARP);
637
638 if (alloc->free + len <= alloc->end)
639 {
640 rv= (Bigint*) alloc->free;
641 alloc->free+= len;
642 }
643 else
644 rv= (Bigint*) malloc(len);
645
646 rv->k= k;
647 rv->maxwds= x;
648 }
649 rv->sign= rv->wds= 0;
650 rv->p.x= (ULong*) (rv + 1);
651 return rv;
652}
653
654
655/*
656 If object was allocated on stack, try putting it to the free
657 list. Otherwise call free().
658*/
659
660static void Bfree(Bigint *v, Stack_alloc *alloc)
661{
662 char *gptr= (char*) v; /* generic pointer */
663 if (gptr < alloc->begin || gptr >= alloc->end)
664 free(gptr);
665 else if (v->k <= Kmax)
666 {
667 /*
668 Maintain free lists only for stack objects: this way we don't
669 have to bother with freeing lists in the end of dtoa;
670 heap should not be used normally anyway.
671 */
672 v->p.next= alloc->freelist[v->k];
673 alloc->freelist[v->k]= v;
674 }
675}
676
677
678/*
679 This is to place return value of dtoa in: tries to use stack
680 as well, but passes by free lists management and just aligns len by
681 the pointer size in order to not break the alignment rules when storing a
682 pointer to a Bigint.
683*/
684
685static char *dtoa_alloc(int i, Stack_alloc *alloc)
686{
687 char *rv;
688 int aligned_size= MY_ALIGN(i, SIZEOF_CHARP);
689 if (alloc->free + aligned_size <= alloc->end)
690 {
691 rv= alloc->free;
692 alloc->free+= aligned_size;
693 }
694 else
695 rv= malloc(i);
696 return rv;
697}
698
699
700/*
701 dtoa_free() must be used to free values s returned by dtoa()
702 This is the counterpart of dtoa_alloc()
703*/
704
705static void dtoa_free(char *gptr, char *buf, size_t buf_size)
706{
707 if (gptr < buf || gptr >= buf + buf_size)
708 free(gptr);
709}
710
711
712/* Bigint arithmetic functions */
713
714/* Multiply by m and add a */
715
716static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc)
717{
718 int i, wds;
719 ULong *x;
720 ULLong carry, y;
721 Bigint *b1;
722
723 wds= b->wds;
724 x= b->p.x;
725 i= 0;
726 carry= a;
727 do
728 {
729 y= *x * (ULLong)m + carry;
730 carry= y >> 32;
731 *x++= (ULong)(y & FFFFFFFF);
732 }
733 while (++i < wds);
734 if (carry)
735 {
736 if (wds >= b->maxwds)
737 {
738 b1= Balloc(b->k+1, alloc);
739 Bcopy(b1, b);
740 Bfree(b, alloc);
741 b= b1;
742 }
743 b->p.x[wds++]= (ULong) carry;
744 b->wds= wds;
745 }
746 return b;
747}
748
749
750static int hi0bits(register ULong x)
751{
752 register int k= 0;
753
754 if (!(x & 0xffff0000))
755 {
756 k= 16;
757 x<<= 16;
758 }
759 if (!(x & 0xff000000))
760 {
761 k+= 8;
762 x<<= 8;
763 }
764 if (!(x & 0xf0000000))
765 {
766 k+= 4;
767 x<<= 4;
768 }
769 if (!(x & 0xc0000000))
770 {
771 k+= 2;
772 x<<= 2;
773 }
774 if (!(x & 0x80000000))
775 {
776 k++;
777 if (!(x & 0x40000000))
778 return 32;
779 }
780 return k;
781}
782
783
784static int lo0bits(ULong *y)
785{
786 register int k;
787 register ULong x= *y;
788
789 if (x & 7)
790 {
791 if (x & 1)
792 return 0;
793 if (x & 2)
794 {
795 *y= x >> 1;
796 return 1;
797 }
798 *y= x >> 2;
799 return 2;
800 }
801 k= 0;
802 if (!(x & 0xffff))
803 {
804 k= 16;
805 x>>= 16;
806 }
807 if (!(x & 0xff))
808 {
809 k+= 8;
810 x>>= 8;
811 }
812 if (!(x & 0xf))
813 {
814 k+= 4;
815 x>>= 4;
816 }
817 if (!(x & 0x3))
818 {
819 k+= 2;
820 x>>= 2;
821 }
822 if (!(x & 1))
823 {
824 k++;
825 x>>= 1;
826 if (!x)
827 return 32;
828 }
829 *y= x;
830 return k;
831}
832
833
834/* Convert integer to Bigint number */
835
836static Bigint *i2b(int i, Stack_alloc *alloc)
837{
838 Bigint *b;
839
840 b= Balloc(1, alloc);
841 b->p.x[0]= i;
842 b->wds= 1;
843 return b;
844}
845
846
847/* Multiply two Bigint numbers */
848
849static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc)
850{
851 Bigint *c;
852 int k, wa, wb, wc;
853 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
854 ULong y;
855 ULLong carry, z;
856
857 if (a->wds < b->wds)
858 {
859 c= a;
860 a= b;
861 b= c;
862 }
863 k= a->k;
864 wa= a->wds;
865 wb= b->wds;
866 wc= wa + wb;
867 if (wc > a->maxwds)
868 k++;
869 c= Balloc(k, alloc);
870 for (x= c->p.x, xa= x + wc; x < xa; x++)
871 *x= 0;
872 xa= a->p.x;
873 xae= xa + wa;
874 xb= b->p.x;
875 xbe= xb + wb;
876 xc0= c->p.x;
877 for (; xb < xbe; xc0++)
878 {
879 if ((y= *xb++))
880 {
881 x= xa;
882 xc= xc0;
883 carry= 0;
884 do
885 {
886 z= *x++ * (ULLong)y + *xc + carry;
887 carry= z >> 32;
888 *xc++= (ULong) (z & FFFFFFFF);
889 }
890 while (x < xae);
891 *xc= (ULong) carry;
892 }
893 }
894 for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ;
895 c->wds= wc;
896 return c;
897}
898
899
900/*
901 Precalculated array of powers of 5: tested to be enough for
902 vasting majority of dtoa_r cases.
903*/
904
905static ULong powers5[]=
906{
907 625UL,
908
909 390625UL,
910
911 2264035265UL, 35UL,
912
913 2242703233UL, 762134875UL, 1262UL,
914
915 3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL,
916
917 781532673UL, 64985353UL, 253049085UL, 594863151UL, 3553621484UL,
918 3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL,
919
920 2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL,
921 3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL,
922 1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL,
923 2161952759UL, 4100910556UL, 1608314830UL, 349175UL
924};
925
926
927static Bigint p5_a[]=
928{
929 /* { x } - k - maxwds - sign - wds */
930 { { powers5 }, 1, 1, 0, 1 },
931 { { powers5 + 1 }, 1, 1, 0, 1 },
932 { { powers5 + 2 }, 1, 2, 0, 2 },
933 { { powers5 + 4 }, 2, 3, 0, 3 },
934 { { powers5 + 7 }, 3, 5, 0, 5 },
935 { { powers5 + 12 }, 4, 10, 0, 10 },
936 { { powers5 + 22 }, 5, 19, 0, 19 }
937};
938
939#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
940
941static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc)
942{
943 Bigint *b1, *p5, *p51=NULL;
944 int i;
945 static int p05[3]= { 5, 25, 125 };
946 my_bool overflow= FALSE;
947
948 if ((i= k & 3))
949 b= multadd(b, p05[i-1], 0, alloc);
950
951 if (!(k>>= 2))
952 return b;
953 p5= p5_a;
954 for (;;)
955 {
956 if (k & 1)
957 {
958 b1= mult(b, p5, alloc);
959 Bfree(b, alloc);
960 b= b1;
961 }
962 if (!(k>>= 1))
963 break;
964 /* Calculate next power of 5 */
965 if (overflow)
966 {
967 p51= mult(p5, p5, alloc);
968 Bfree(p5, alloc);
969 p5= p51;
970 }
971 else if (p5 < p5_a + P5A_MAX)
972 ++p5;
973 else if (p5 == p5_a + P5A_MAX)
974 {
975 p5= mult(p5, p5, alloc);
976 overflow= TRUE;
977 }
978 }
979 if (p51)
980 Bfree(p51, alloc);
981 return b;
982}
983
984
985static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc)
986{
987 int i, k1, n, n1;
988 Bigint *b1;
989 ULong *x, *x1, *xe, z;
990
991 n= k >> 5;
992 k1= b->k;
993 n1= n + b->wds + 1;
994 for (i= b->maxwds; n1 > i; i<<= 1)
995 k1++;
996 b1= Balloc(k1, alloc);
997 x1= b1->p.x;
998 for (i= 0; i < n; i++)
999 *x1++= 0;
1000 x= b->p.x;
1001 xe= x + b->wds;
1002 if (k&= 0x1f)
1003 {
1004 k1= 32 - k;
1005 z= 0;
1006 do
1007 {
1008 *x1++= *x << k | z;
1009 z= *x++ >> k1;
1010 }
1011 while (x < xe);
1012 if ((*x1= z))
1013 ++n1;
1014 }
1015 else
1016 do
1017 *x1++= *x++;
1018 while (x < xe);
1019 b1->wds= n1 - 1;
1020 Bfree(b, alloc);
1021 return b1;
1022}
1023
1024
1025static int cmp(Bigint *a, Bigint *b)
1026{
1027 ULong *xa, *xa0, *xb, *xb0;
1028 int i, j;
1029
1030 i= a->wds;
1031 j= b->wds;
1032 if (i-= j)
1033 return i;
1034 xa0= a->p.x;
1035 xa= xa0 + j;
1036 xb0= b->p.x;
1037 xb= xb0 + j;
1038 for (;;)
1039 {
1040 if (*--xa != *--xb)
1041 return *xa < *xb ? -1 : 1;
1042 if (xa <= xa0)
1043 break;
1044 }
1045 return 0;
1046}
1047
1048
1049static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc)
1050{
1051 Bigint *c;
1052 int i, wa, wb;
1053 ULong *xa, *xae, *xb, *xbe, *xc;
1054 ULLong borrow, y;
1055
1056 i= cmp(a,b);
1057 if (!i)
1058 {
1059 c= Balloc(0, alloc);
1060 c->wds= 1;
1061 c->p.x[0]= 0;
1062 return c;
1063 }
1064 if (i < 0)
1065 {
1066 c= a;
1067 a= b;
1068 b= c;
1069 i= 1;
1070 }
1071 else
1072 i= 0;
1073 c= Balloc(a->k, alloc);
1074 c->sign= i;
1075 wa= a->wds;
1076 xa= a->p.x;
1077 xae= xa + wa;
1078 wb= b->wds;
1079 xb= b->p.x;
1080 xbe= xb + wb;
1081 xc= c->p.x;
1082 borrow= 0;
1083 do
1084 {
1085 y= (ULLong)*xa++ - *xb++ - borrow;
1086 borrow= y >> 32 & (ULong)1;
1087 *xc++= (ULong) (y & FFFFFFFF);
1088 }
1089 while (xb < xbe);
1090 while (xa < xae)
1091 {
1092 y= *xa++ - borrow;
1093 borrow= y >> 32 & (ULong)1;
1094 *xc++= (ULong) (y & FFFFFFFF);
1095 }
1096 while (!*--xc)
1097 wa--;
1098 c->wds= wa;
1099 return c;
1100}
1101
1102
1103static Bigint *d2b(U *d, int *e, int *bits, Stack_alloc *alloc)
1104{
1105 Bigint *b;
1106 int de, k;
1107 ULong *x, y, z;
1108 int i;
1109#define d0 word0(d)
1110#define d1 word1(d)
1111
1112 b= Balloc(1, alloc);
1113 x= b->p.x;
1114
1115 z= d0 & Frac_mask;
1116 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1117 if ((de= (int)(d0 >> Exp_shift)))
1118 z|= Exp_msk1;
1119 if ((y= d1))
1120 {
1121 if ((k= lo0bits(&y)))
1122 {
1123 x[0]= y | z << (32 - k);
1124 z>>= k;
1125 }
1126 else
1127 x[0]= y;
1128 i= b->wds= (x[1]= z) ? 2 : 1;
1129 }
1130 else
1131 {
1132 k= lo0bits(&z);
1133 x[0]= z;
1134 i= b->wds= 1;
1135 k+= 32;
1136 }
1137 if (de)
1138 {
1139 *e= de - Bias - (P-1) + k;
1140 *bits= P - k;
1141 }
1142 else
1143 {
1144 *e= de - Bias - (P-1) + 1 + k;
1145 *bits= 32*i - hi0bits(x[i-1]);
1146 }
1147 return b;
1148#undef d0
1149#undef d1
1150}
1151
1152
1153static const double tens[] =
1154{
1155 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1156 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1157 1e20, 1e21, 1e22
1158};
1159
1160static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 };
1161/*
1162 The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
1163 flag unnecessarily. It leads to a song and dance at the end of strtod.
1164*/
1165#define Scale_Bit 0x10
1166#define n_bigtens 5
1167
1168
1169static int quorem(Bigint *b, Bigint *S)
1170{
1171 int n;
1172 ULong *bx, *bxe, q, *sx, *sxe;
1173 ULLong borrow, carry, y, ys;
1174
1175 n= S->wds;
1176 if (b->wds < n)
1177 return 0;
1178 sx= S->p.x;
1179 sxe= sx + --n;
1180 bx= b->p.x;
1181 bxe= bx + n;
1182 q= *bxe / (*sxe + 1); /* ensure q <= true quotient */
1183 if (q)
1184 {
1185 borrow= 0;
1186 carry= 0;
1187 do
1188 {
1189 ys= *sx++ * (ULLong)q + carry;
1190 carry= ys >> 32;
1191 y= *bx - (ys & FFFFFFFF) - borrow;
1192 borrow= y >> 32 & (ULong)1;
1193 *bx++= (ULong) (y & FFFFFFFF);
1194 }
1195 while (sx <= sxe);
1196 if (!*bxe)
1197 {
1198 bx= b->p.x;
1199 while (--bxe > bx && !*bxe)
1200 --n;
1201 b->wds= n;
1202 }
1203 }
1204 if (cmp(b, S) >= 0)
1205 {
1206 q++;
1207 borrow= 0;
1208 carry= 0;
1209 bx= b->p.x;
1210 sx= S->p.x;
1211 do
1212 {
1213 ys= *sx++ + carry;
1214 carry= ys >> 32;
1215 y= *bx - (ys & FFFFFFFF) - borrow;
1216 borrow= y >> 32 & (ULong)1;
1217 *bx++= (ULong) (y & FFFFFFFF);
1218 }
1219 while (sx <= sxe);
1220 bx= b->p.x;
1221 bxe= bx + n;
1222 if (!*bxe)
1223 {
1224 while (--bxe > bx && !*bxe)
1225 --n;
1226 b->wds= n;
1227 }
1228 }
1229 return q;
1230}
1231
1232
1233/*
1234 dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1235
1236 Inspired by "How to Print Floating-Point Numbers Accurately" by
1237 Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
1238
1239 Modifications:
1240 1. Rather than iterating, we use a simple numeric overestimate
1241 to determine k= floor(log10(d)). We scale relevant
1242 quantities using O(log2(k)) rather than O(k) multiplications.
1243 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1244 try to generate digits strictly left to right. Instead, we
1245 compute with fewer bits and propagate the carry if necessary
1246 when rounding the final digit up. This is often faster.
1247 3. Under the assumption that input will be rounded nearest,
1248 mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1249 That is, we allow equality in stopping tests when the
1250 round-nearest rule will give the same floating-point value
1251 as would satisfaction of the stopping test with strict
1252 inequality.
1253 4. We remove common factors of powers of 2 from relevant
1254 quantities.
1255 5. When converting floating-point integers less than 1e16,
1256 we use floating-point arithmetic rather than resorting
1257 to multiple-precision integers.
1258 6. When asked to produce fewer than 15 digits, we first try
1259 to get by with floating-point arithmetic; we resort to
1260 multiple-precision integer arithmetic only if we cannot
1261 guarantee that the floating-point calculation has given
1262 the correctly rounded result. For k requested digits and
1263 "uniformly" distributed input, the probability is
1264 something like 10^(k-15) that we must resort to the Long
1265 calculation.
1266 */
1267
1268static char *dtoa(double dd, int mode, int ndigits, int *decpt, int *sign,
1269 char **rve, char *buf, size_t buf_size)
1270{
1271 /*
1272 Arguments ndigits, decpt, sign are similar to those
1273 of ecvt and fcvt; trailing zeros are suppressed from
1274 the returned string. If not null, *rve is set to point
1275 to the end of the return value. If d is +-Infinity or NaN,
1276 then *decpt is set to DTOA_OVERFLOW.
1277
1278 mode:
1279 0 ==> shortest string that yields d when read in
1280 and rounded to nearest.
1281 1 ==> like 0, but with Steele & White stopping rule;
1282 e.g. with IEEE P754 arithmetic , mode 0 gives
1283 1e23 whereas mode 1 gives 9.999999999999999e22.
1284 2 ==> MAX(1,ndigits) significant digits. This gives a
1285 return value similar to that of ecvt, except
1286 that trailing zeros are suppressed.
1287 3 ==> through ndigits past the decimal point. This
1288 gives a return value similar to that from fcvt,
1289 except that trailing zeros are suppressed, and
1290 ndigits can be negative.
1291 4,5 ==> similar to 2 and 3, respectively, but (in
1292 round-nearest mode) with the tests of mode 0 to
1293 possibly return a shorter string that rounds to d.
1294 With IEEE arithmetic and compilation with
1295 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
1296 as modes 2 and 3 when FLT_ROUNDS != 1.
1297 6-9 ==> Debugging modes similar to mode - 4: don't try
1298 fast floating-point estimate (if applicable).
1299
1300 Values of mode other than 0-9 are treated as mode 0.
1301
1302 Sufficient space is allocated to the return value
1303 to hold the suppressed trailing zeros.
1304 */
1305
1306 int bbits, b2, b5, be, dig, i, ieps, UNINIT_VAR(ilim), ilim0,
1307 UNINIT_VAR(ilim1), j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1308 spec_case, try_quick;
1309 Long L;
1310 int denorm;
1311 ULong x;
1312 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1313 U d2, eps, u;
1314 double ds;
1315 char *s, *s0;
1316#ifdef Honor_FLT_ROUNDS
1317 int rounding;
1318#endif
1319 Stack_alloc alloc;
1320
1321 alloc.begin= alloc.free= buf;
1322 alloc.end= buf + buf_size;
1323 memset(alloc.freelist, 0, sizeof(alloc.freelist));
1324
1325 u.d= dd;
1326 if (word0(&u) & Sign_bit)
1327 {
1328 /* set sign for everything, including 0's and NaNs */
1329 *sign= 1;
1330 word0(&u) &= ~Sign_bit; /* clear sign bit */
1331 }
1332 else
1333 *sign= 0;
1334
1335 /* If infinity, set decpt to DTOA_OVERFLOW, if 0 set it to 1 */
1336 if (((word0(&u) & Exp_mask) == Exp_mask && (*decpt= DTOA_OVERFLOW)) ||
1337 (!dval(&u) && (*decpt= 1)))
1338 {
1339 /* Infinity, NaN, 0 */
1340 char *res= (char*) dtoa_alloc(2, &alloc);
1341 res[0]= '0';
1342 res[1]= '\0';
1343 if (rve)
1344 *rve= res + 1;
1345 return res;
1346 }
1347
1348#ifdef Honor_FLT_ROUNDS
1349 if ((rounding= Flt_Rounds) >= 2)
1350 {
1351 if (*sign)
1352 rounding= rounding == 2 ? 0 : 2;
1353 else
1354 if (rounding != 2)
1355 rounding= 0;
1356 }
1357#endif
1358
1359 b= d2b(&u, &be, &bbits, &alloc);
1360 if ((i= (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1361 {
1362 dval(&d2)= dval(&u);
1363 word0(&d2) &= Frac_mask1;
1364 word0(&d2) |= Exp_11;
1365
1366 /*
1367 log(x) ~=~ log(1.5) + (x-1.5)/1.5
1368 log10(x) = log(x) / log(10)
1369 ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1370 log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
1371
1372 This suggests computing an approximation k to log10(d) by
1373
1374 k= (i - Bias)*0.301029995663981
1375 + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1376
1377 We want k to be too large rather than too small.
1378 The error in the first-order Taylor series approximation
1379 is in our favor, so we just round up the constant enough
1380 to compensate for any error in the multiplication of
1381 (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1382 and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1383 adding 1e-13 to the constant term more than suffices.
1384 Hence we adjust the constant term to 0.1760912590558.
1385 (We could get a more accurate k by invoking log10,
1386 but this is probably not worthwhile.)
1387 */
1388
1389 i-= Bias;
1390 denorm= 0;
1391 }
1392 else
1393 {
1394 /* d is denormalized */
1395
1396 i= bbits + be + (Bias + (P-1) - 1);
1397 x= i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
1398 : word1(&u) << (32 - i);
1399 dval(&d2)= x;
1400 word0(&d2)-= 31*Exp_msk1; /* adjust exponent */
1401 i-= (Bias + (P-1) - 1) + 1;
1402 denorm= 1;
1403 }
1404 ds= (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1405 k= (int)ds;
1406 if (ds < 0. && ds != k)
1407 k--; /* want k= floor(ds) */
1408 k_check= 1;
1409 if (k >= 0 && k <= Ten_pmax)
1410 {
1411 if (dval(&u) < tens[k])
1412 k--;
1413 k_check= 0;
1414 }
1415 j= bbits - i - 1;
1416 if (j >= 0)
1417 {
1418 b2= 0;
1419 s2= j;
1420 }
1421 else
1422 {
1423 b2= -j;
1424 s2= 0;
1425 }
1426 if (k >= 0)
1427 {
1428 b5= 0;
1429 s5= k;
1430 s2+= k;
1431 }
1432 else
1433 {
1434 b2-= k;
1435 b5= -k;
1436 s5= 0;
1437 }
1438 if (mode < 0 || mode > 9)
1439 mode= 0;
1440
1441#ifdef Check_FLT_ROUNDS
1442 try_quick= Rounding == 1;
1443#else
1444 try_quick= 1;
1445#endif
1446
1447 if (mode > 5)
1448 {
1449 mode-= 4;
1450 try_quick= 0;
1451 }
1452 leftright= 1;
1453 switch (mode) {
1454 case 0:
1455 case 1:
1456 ilim= ilim1= -1;
1457 i= 18;
1458 ndigits= 0;
1459 break;
1460 case 2:
1461 leftright= 0;
1462 /* fall through */
1463 case 4:
1464 if (ndigits <= 0)
1465 ndigits= 1;
1466 ilim= ilim1= i= ndigits;
1467 break;
1468 case 3:
1469 leftright= 0;
1470 /* fall through */
1471 case 5:
1472 i= ndigits + k + 1;
1473 ilim= i;
1474 ilim1= i - 1;
1475 if (i <= 0)
1476 i= 1;
1477 }
1478 s= s0= dtoa_alloc(i, &alloc);
1479
1480#ifdef Honor_FLT_ROUNDS
1481 if (mode > 1 && rounding != 1)
1482 leftright= 0;
1483#endif
1484
1485 if (ilim >= 0 && ilim <= Quick_max && try_quick)
1486 {
1487 /* Try to get by with floating-point arithmetic. */
1488 i= 0;
1489 dval(&d2)= dval(&u);
1490 k0= k;
1491 ilim0= ilim;
1492 ieps= 2; /* conservative */
1493 if (k > 0)
1494 {
1495 ds= tens[k&0xf];
1496 j= k >> 4;
1497 if (j & Bletch)
1498 {
1499 /* prevent overflows */
1500 j&= Bletch - 1;
1501 dval(&u)/= bigtens[n_bigtens-1];
1502 ieps++;
1503 }
1504 for (; j; j>>= 1, i++)
1505 {
1506 if (j & 1)
1507 {
1508 ieps++;
1509 ds*= bigtens[i];
1510 }
1511 }
1512 dval(&u)/= ds;
1513 }
1514 else if ((j1= -k))
1515 {
1516 dval(&u)*= tens[j1 & 0xf];
1517 for (j= j1 >> 4; j; j>>= 1, i++)
1518 {
1519 if (j & 1)
1520 {
1521 ieps++;
1522 dval(&u)*= bigtens[i];
1523 }
1524 }
1525 }
1526 if (k_check && dval(&u) < 1. && ilim > 0)
1527 {
1528 if (ilim1 <= 0)
1529 goto fast_failed;
1530 ilim= ilim1;
1531 k--;
1532 dval(&u)*= 10.;
1533 ieps++;
1534 }
1535 dval(&eps)= ieps*dval(&u) + 7.;
1536 word0(&eps)-= (P-1)*Exp_msk1;
1537 if (ilim == 0)
1538 {
1539 S= mhi= 0;
1540 dval(&u)-= 5.;
1541 if (dval(&u) > dval(&eps))
1542 goto one_digit;
1543 if (dval(&u) < -dval(&eps))
1544 goto no_digits;
1545 goto fast_failed;
1546 }
1547 if (leftright)
1548 {
1549 /* Use Steele & White method of only generating digits needed. */
1550 dval(&eps)= 0.5/tens[ilim-1] - dval(&eps);
1551 for (i= 0;;)
1552 {
1553 L= (Long) dval(&u);
1554 dval(&u)-= L;
1555 *s++= '0' + (int)L;
1556 if (dval(&u) < dval(&eps))
1557 goto ret1;
1558 if (1. - dval(&u) < dval(&eps))
1559 goto bump_up;
1560 if (++i >= ilim)
1561 break;
1562 dval(&eps)*= 10.;
1563 dval(&u)*= 10.;
1564 }
1565 }
1566 else
1567 {
1568 /* Generate ilim digits, then fix them up. */
1569 dval(&eps)*= tens[ilim-1];
1570 for (i= 1;; i++, dval(&u)*= 10.)
1571 {
1572 L= (Long)(dval(&u));
1573 if (!(dval(&u)-= L))
1574 ilim= i;
1575 *s++= '0' + (int)L;
1576 if (i == ilim)
1577 {
1578 if (dval(&u) > 0.5 + dval(&eps))
1579 goto bump_up;
1580 else if (dval(&u) < 0.5 - dval(&eps))
1581 {
1582 while (*--s == '0');
1583 s++;
1584 goto ret1;
1585 }
1586 break;
1587 }
1588 }
1589 }
1590 fast_failed:
1591 s= s0;
1592 dval(&u)= dval(&d2);
1593 k= k0;
1594 ilim= ilim0;
1595 }
1596
1597 /* Do we have a "small" integer? */
1598
1599 if (be >= 0 && k <= Int_max)
1600 {
1601 /* Yes. */
1602 ds= tens[k];
1603 if (ndigits < 0 && ilim <= 0)
1604 {
1605 S= mhi= 0;
1606 if (ilim < 0 || dval(&u) <= 5*ds)
1607 goto no_digits;
1608 goto one_digit;
1609 }
1610 for (i= 1;; i++, dval(&u)*= 10.)
1611 {
1612 L= (Long)(dval(&u) / ds);
1613 dval(&u)-= L*ds;
1614#ifdef Check_FLT_ROUNDS
1615 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1616 if (dval(&u) < 0)
1617 {
1618 L--;
1619 dval(&u)+= ds;
1620 }
1621#endif
1622 *s++= '0' + (int)L;
1623 if (!dval(&u))
1624 {
1625 break;
1626 }
1627 if (i == ilim)
1628 {
1629#ifdef Honor_FLT_ROUNDS
1630 if (mode > 1)
1631 {
1632 switch (rounding) {
1633 case 0: goto ret1;
1634 case 2: goto bump_up;
1635 }
1636 }
1637#endif
1638 dval(&u)+= dval(&u);
1639 if (dval(&u) > ds || (dval(&u) == ds && L & 1))
1640 {
1641bump_up:
1642 while (*--s == '9')
1643 if (s == s0)
1644 {
1645 k++;
1646 *s= '0';
1647 break;
1648 }
1649 ++*s++;
1650 }
1651 break;
1652 }
1653 }
1654 goto ret1;
1655 }
1656
1657 m2= b2;
1658 m5= b5;
1659 mhi= mlo= 0;
1660 if (leftright)
1661 {
1662 i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
1663 b2+= i;
1664 s2+= i;
1665 mhi= i2b(1, &alloc);
1666 }
1667 if (m2 > 0 && s2 > 0)
1668 {
1669 i= m2 < s2 ? m2 : s2;
1670 b2-= i;
1671 m2-= i;
1672 s2-= i;
1673 }
1674 if (b5 > 0)
1675 {
1676 if (leftright)
1677 {
1678 if (m5 > 0)
1679 {
1680 mhi= pow5mult(mhi, m5, &alloc);
1681 b1= mult(mhi, b, &alloc);
1682 Bfree(b, &alloc);
1683 b= b1;
1684 }
1685 if ((j= b5 - m5))
1686 b= pow5mult(b, j, &alloc);
1687 }
1688 else
1689 b= pow5mult(b, b5, &alloc);
1690 }
1691 S= i2b(1, &alloc);
1692 if (s5 > 0)
1693 S= pow5mult(S, s5, &alloc);
1694
1695 /* Check for special case that d is a normalized power of 2. */
1696
1697 spec_case= 0;
1698 if ((mode < 2 || leftright)
1699#ifdef Honor_FLT_ROUNDS
1700 && rounding == 1
1701#endif
1702 )
1703 {
1704 if (!word1(&u) && !(word0(&u) & Bndry_mask) &&
1705 word0(&u) & (Exp_mask & ~Exp_msk1)
1706 )
1707 {
1708 /* The special case */
1709 b2+= Log2P;
1710 s2+= Log2P;
1711 spec_case= 1;
1712 }
1713 }
1714
1715 /*
1716 Arrange for convenient computation of quotients:
1717 shift left if necessary so divisor has 4 leading 0 bits.
1718
1719 Perhaps we should just compute leading 28 bits of S once
1720 a nd for all and pass them and a shift to quorem, so it
1721 can do shifts and ors to compute the numerator for q.
1722 */
1723 if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
1724 i= 32 - i;
1725 if (i > 4)
1726 {
1727 i-= 4;
1728 b2+= i;
1729 m2+= i;
1730 s2+= i;
1731 }
1732 else if (i < 4)
1733 {
1734 i+= 28;
1735 b2+= i;
1736 m2+= i;
1737 s2+= i;
1738 }
1739 if (b2 > 0)
1740 b= lshift(b, b2, &alloc);
1741 if (s2 > 0)
1742 S= lshift(S, s2, &alloc);
1743 if (k_check)
1744 {
1745 if (cmp(b,S) < 0)
1746 {
1747 k--;
1748 /* we botched the k estimate */
1749 b= multadd(b, 10, 0, &alloc);
1750 if (leftright)
1751 mhi= multadd(mhi, 10, 0, &alloc);
1752 ilim= ilim1;
1753 }
1754 }
1755 if (ilim <= 0 && (mode == 3 || mode == 5))
1756 {
1757 if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
1758 {
1759 /* no digits, fcvt style */
1760no_digits:
1761 k= -1 - ndigits;
1762 goto ret;
1763 }
1764one_digit:
1765 *s++= '1';
1766 k++;
1767 goto ret;
1768 }
1769 if (leftright)
1770 {
1771 if (m2 > 0)
1772 mhi= lshift(mhi, m2, &alloc);
1773
1774 /*
1775 Compute mlo -- check for special case that d is a normalized power of 2.
1776 */
1777
1778 mlo= mhi;
1779 if (spec_case)
1780 {
1781 mhi= Balloc(mhi->k, &alloc);
1782 Bcopy(mhi, mlo);
1783 mhi= lshift(mhi, Log2P, &alloc);
1784 }
1785
1786 for (i= 1;;i++)
1787 {
1788 dig= quorem(b,S) + '0';
1789 /* Do we yet have the shortest decimal string that will round to d? */
1790 j= cmp(b, mlo);
1791 delta= diff(S, mhi, &alloc);
1792 j1= delta->sign ? 1 : cmp(b, delta);
1793 Bfree(delta, &alloc);
1794 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
1795#ifdef Honor_FLT_ROUNDS
1796 && rounding >= 1
1797#endif
1798 )
1799 {
1800 if (dig == '9')
1801 goto round_9_up;
1802 if (j > 0)
1803 dig++;
1804 *s++= dig;
1805 goto ret;
1806 }
1807 if (j < 0 || (j == 0 && mode != 1 && !(word1(&u) & 1)))
1808 {
1809 if (!b->p.x[0] && b->wds <= 1)
1810 {
1811 goto accept_dig;
1812 }
1813#ifdef Honor_FLT_ROUNDS
1814 if (mode > 1)
1815 switch (rounding) {
1816 case 0: goto accept_dig;
1817 case 2: goto keep_dig;
1818 }
1819#endif /*Honor_FLT_ROUNDS*/
1820 if (j1 > 0)
1821 {
1822 b= lshift(b, 1, &alloc);
1823 j1= cmp(b, S);
1824 if ((j1 > 0 || (j1 == 0 && dig & 1))
1825 && dig++ == '9')
1826 goto round_9_up;
1827 }
1828accept_dig:
1829 *s++= dig;
1830 goto ret;
1831 }
1832 if (j1 > 0)
1833 {
1834#ifdef Honor_FLT_ROUNDS
1835 if (!rounding)
1836 goto accept_dig;
1837#endif
1838 if (dig == '9')
1839 { /* possible if i == 1 */
1840round_9_up:
1841 *s++= '9';
1842 goto roundoff;
1843 }
1844 *s++= dig + 1;
1845 goto ret;
1846 }
1847#ifdef Honor_FLT_ROUNDS
1848keep_dig:
1849#endif
1850 *s++= dig;
1851 if (i == ilim)
1852 break;
1853 b= multadd(b, 10, 0, &alloc);
1854 if (mlo == mhi)
1855 mlo= mhi= multadd(mhi, 10, 0, &alloc);
1856 else
1857 {
1858 mlo= multadd(mlo, 10, 0, &alloc);
1859 mhi= multadd(mhi, 10, 0, &alloc);
1860 }
1861 }
1862 }
1863 else
1864 for (i= 1;; i++)
1865 {
1866 *s++= dig= quorem(b,S) + '0';
1867 if (!b->p.x[0] && b->wds <= 1)
1868 {
1869 goto ret;
1870 }
1871 if (i >= ilim)
1872 break;
1873 b= multadd(b, 10, 0, &alloc);
1874 }
1875
1876 /* Round off last digit */
1877
1878#ifdef Honor_FLT_ROUNDS
1879 switch (rounding) {
1880 case 0: goto trimzeros;
1881 case 2: goto roundoff;
1882 }
1883#endif
1884 b= lshift(b, 1, &alloc);
1885 j= cmp(b, S);
1886 if (j > 0 || (j == 0 && dig & 1))
1887 {
1888roundoff:
1889 while (*--s == '9')
1890 if (s == s0)
1891 {
1892 k++;
1893 *s++= '1';
1894 goto ret;
1895 }
1896 ++*s++;
1897 }
1898 else
1899 {
1900#ifdef Honor_FLT_ROUNDS
1901trimzeros:
1902#endif
1903 while (*--s == '0');
1904 s++;
1905 }
1906ret:
1907 Bfree(S, &alloc);
1908 if (mhi)
1909 {
1910 if (mlo && mlo != mhi)
1911 Bfree(mlo, &alloc);
1912 Bfree(mhi, &alloc);
1913 }
1914ret1:
1915 Bfree(b, &alloc);
1916 *s= 0;
1917 *decpt= k + 1;
1918 if (rve)
1919 *rve= s;
1920 return s0;
1921}
1922