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(HAVE_BIGENDIAN) || defined(WORDS_BIGENDIAN) || \
516 (defined(__FLOAT_WORD_ORDER) && (__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 /* coverity[assign_where_compare_meant] */
1337 if (((word0(&u) & Exp_mask) == Exp_mask && (*decpt= DTOA_OVERFLOW)) ||
1338 /* coverity[assign_where_compare_meant] */
1339 (!dval(&u) && (*decpt= 1)))
1340 {
1341 /* Infinity, NaN, 0 */
1342 char *res= (char*) dtoa_alloc(2, &alloc);
1343 res[0]= '0';
1344 res[1]= '\0';
1345 if (rve)
1346 *rve= res + 1;
1347 return res;
1348 }
1349
1350#ifdef Honor_FLT_ROUNDS
1351 if ((rounding= Flt_Rounds) >= 2)
1352 {
1353 if (*sign)
1354 rounding= rounding == 2 ? 0 : 2;
1355 else
1356 if (rounding != 2)
1357 rounding= 0;
1358 }
1359#endif
1360
1361 b= d2b(&u, &be, &bbits, &alloc);
1362 if ((i= (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1363 {
1364 dval(&d2)= dval(&u);
1365 word0(&d2) &= Frac_mask1;
1366 word0(&d2) |= Exp_11;
1367
1368 /*
1369 log(x) ~=~ log(1.5) + (x-1.5)/1.5
1370 log10(x) = log(x) / log(10)
1371 ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1372 log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
1373
1374 This suggests computing an approximation k to log10(d) by
1375
1376 k= (i - Bias)*0.301029995663981
1377 + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1378
1379 We want k to be too large rather than too small.
1380 The error in the first-order Taylor series approximation
1381 is in our favor, so we just round up the constant enough
1382 to compensate for any error in the multiplication of
1383 (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1384 and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1385 adding 1e-13 to the constant term more than suffices.
1386 Hence we adjust the constant term to 0.1760912590558.
1387 (We could get a more accurate k by invoking log10,
1388 but this is probably not worthwhile.)
1389 */
1390
1391 i-= Bias;
1392 denorm= 0;
1393 }
1394 else
1395 {
1396 /* d is denormalized */
1397
1398 i= bbits + be + (Bias + (P-1) - 1);
1399 x= i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
1400 : word1(&u) << (32 - i);
1401 dval(&d2)= x;
1402 word0(&d2)-= 31*Exp_msk1; /* adjust exponent */
1403 i-= (Bias + (P-1) - 1) + 1;
1404 denorm= 1;
1405 }
1406 ds= (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1407 k= (int)ds;
1408 if (ds < 0. && ds != k)
1409 k--; /* want k= floor(ds) */
1410 k_check= 1;
1411 if (k >= 0 && k <= Ten_pmax)
1412 {
1413 if (dval(&u) < tens[k])
1414 k--;
1415 k_check= 0;
1416 }
1417 j= bbits - i - 1;
1418 if (j >= 0)
1419 {
1420 b2= 0;
1421 s2= j;
1422 }
1423 else
1424 {
1425 b2= -j;
1426 s2= 0;
1427 }
1428 if (k >= 0)
1429 {
1430 b5= 0;
1431 s5= k;
1432 s2+= k;
1433 }
1434 else
1435 {
1436 b2-= k;
1437 b5= -k;
1438 s5= 0;
1439 }
1440 if (mode < 0 || mode > 9)
1441 mode= 0;
1442
1443#ifdef Check_FLT_ROUNDS
1444 try_quick= Rounding == 1;
1445#else
1446 try_quick= 1;
1447#endif
1448
1449 if (mode > 5)
1450 {
1451 mode-= 4;
1452 try_quick= 0;
1453 }
1454 leftright= 1;
1455 switch (mode) {
1456 case 0:
1457 case 1:
1458 ilim= ilim1= -1;
1459 i= 18;
1460 ndigits= 0;
1461 break;
1462 case 2:
1463 leftright= 0;
1464 /* fall through */
1465 case 4:
1466 if (ndigits <= 0)
1467 ndigits= 1;
1468 ilim= ilim1= i= ndigits;
1469 break;
1470 case 3:
1471 leftright= 0;
1472 /* fall through */
1473 case 5:
1474 i= ndigits + k + 1;
1475 ilim= i;
1476 ilim1= i - 1;
1477 if (i <= 0)
1478 i= 1;
1479 }
1480 s= s0= dtoa_alloc(i, &alloc);
1481
1482#ifdef Honor_FLT_ROUNDS
1483 if (mode > 1 && rounding != 1)
1484 leftright= 0;
1485#endif
1486
1487 if (ilim >= 0 && ilim <= Quick_max && try_quick)
1488 {
1489 /* Try to get by with floating-point arithmetic. */
1490 i= 0;
1491 dval(&d2)= dval(&u);
1492 k0= k;
1493 ilim0= ilim;
1494 ieps= 2; /* conservative */
1495 if (k > 0)
1496 {
1497 ds= tens[k&0xf];
1498 j= k >> 4;
1499 if (j & Bletch)
1500 {
1501 /* prevent overflows */
1502 j&= Bletch - 1;
1503 dval(&u)/= bigtens[n_bigtens-1];
1504 ieps++;
1505 }
1506 for (; j; j>>= 1, i++)
1507 {
1508 if (j & 1)
1509 {
1510 ieps++;
1511 ds*= bigtens[i];
1512 }
1513 }
1514 dval(&u)/= ds;
1515 }
1516 else if ((j1= -k))
1517 {
1518 dval(&u)*= tens[j1 & 0xf];
1519 for (j= j1 >> 4; j; j>>= 1, i++)
1520 {
1521 if (j & 1)
1522 {
1523 ieps++;
1524 dval(&u)*= bigtens[i];
1525 }
1526 }
1527 }
1528 if (k_check && dval(&u) < 1. && ilim > 0)
1529 {
1530 if (ilim1 <= 0)
1531 goto fast_failed;
1532 ilim= ilim1;
1533 k--;
1534 dval(&u)*= 10.;
1535 ieps++;
1536 }
1537 dval(&eps)= ieps*dval(&u) + 7.;
1538 word0(&eps)-= (P-1)*Exp_msk1;
1539 if (ilim == 0)
1540 {
1541 S= mhi= 0;
1542 dval(&u)-= 5.;
1543 if (dval(&u) > dval(&eps))
1544 goto one_digit;
1545 if (dval(&u) < -dval(&eps))
1546 goto no_digits;
1547 goto fast_failed;
1548 }
1549 if (leftright)
1550 {
1551 /* Use Steele & White method of only generating digits needed. */
1552 dval(&eps)= 0.5/tens[ilim-1] - dval(&eps);
1553 for (i= 0;;)
1554 {
1555 L= (Long) dval(&u);
1556 dval(&u)-= L;
1557 *s++= '0' + (int)L;
1558 if (dval(&u) < dval(&eps))
1559 goto ret1;
1560 if (1. - dval(&u) < dval(&eps))
1561 goto bump_up;
1562 if (++i >= ilim)
1563 break;
1564 dval(&eps)*= 10.;
1565 dval(&u)*= 10.;
1566 }
1567 }
1568 else
1569 {
1570 /* Generate ilim digits, then fix them up. */
1571 dval(&eps)*= tens[ilim-1];
1572 for (i= 1;; i++, dval(&u)*= 10.)
1573 {
1574 L= (Long)(dval(&u));
1575 if (!(dval(&u)-= L))
1576 ilim= i;
1577 *s++= '0' + (int)L;
1578 if (i == ilim)
1579 {
1580 if (dval(&u) > 0.5 + dval(&eps))
1581 goto bump_up;
1582 else if (dval(&u) < 0.5 - dval(&eps))
1583 {
1584 while (*--s == '0');
1585 s++;
1586 goto ret1;
1587 }
1588 break;
1589 }
1590 }
1591 }
1592 fast_failed:
1593 s= s0;
1594 dval(&u)= dval(&d2);
1595 k= k0;
1596 ilim= ilim0;
1597 }
1598
1599 /* Do we have a "small" integer? */
1600
1601 if (be >= 0 && k <= Int_max)
1602 {
1603 /* Yes. */
1604 ds= tens[k];
1605 if (ndigits < 0 && ilim <= 0)
1606 {
1607 S= mhi= 0;
1608 if (ilim < 0 || dval(&u) <= 5*ds)
1609 goto no_digits;
1610 goto one_digit;
1611 }
1612 for (i= 1;; i++, dval(&u)*= 10.)
1613 {
1614 L= (Long)(dval(&u) / ds);
1615 dval(&u)-= L*ds;
1616#ifdef Check_FLT_ROUNDS
1617 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1618 if (dval(&u) < 0)
1619 {
1620 L--;
1621 dval(&u)+= ds;
1622 }
1623#endif
1624 *s++= '0' + (int)L;
1625 if (!dval(&u))
1626 {
1627 break;
1628 }
1629 if (i == ilim)
1630 {
1631#ifdef Honor_FLT_ROUNDS
1632 if (mode > 1)
1633 {
1634 switch (rounding) {
1635 case 0: goto ret1;
1636 case 2: goto bump_up;
1637 }
1638 }
1639#endif
1640 dval(&u)+= dval(&u);
1641 if (dval(&u) > ds || (dval(&u) == ds && L & 1))
1642 {
1643bump_up:
1644 while (*--s == '9')
1645 if (s == s0)
1646 {
1647 k++;
1648 *s= '0';
1649 break;
1650 }
1651 ++*s++;
1652 }
1653 break;
1654 }
1655 }
1656 goto ret1;
1657 }
1658
1659 m2= b2;
1660 m5= b5;
1661 mhi= mlo= 0;
1662 if (leftright)
1663 {
1664 i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
1665 b2+= i;
1666 s2+= i;
1667 mhi= i2b(1, &alloc);
1668 }
1669 if (m2 > 0 && s2 > 0)
1670 {
1671 i= m2 < s2 ? m2 : s2;
1672 b2-= i;
1673 m2-= i;
1674 s2-= i;
1675 }
1676 if (b5 > 0)
1677 {
1678 if (leftright)
1679 {
1680 if (m5 > 0)
1681 {
1682 mhi= pow5mult(mhi, m5, &alloc);
1683 b1= mult(mhi, b, &alloc);
1684 Bfree(b, &alloc);
1685 b= b1;
1686 }
1687 if ((j= b5 - m5))
1688 b= pow5mult(b, j, &alloc);
1689 }
1690 else
1691 b= pow5mult(b, b5, &alloc);
1692 }
1693 S= i2b(1, &alloc);
1694 if (s5 > 0)
1695 S= pow5mult(S, s5, &alloc);
1696
1697 /* Check for special case that d is a normalized power of 2. */
1698
1699 spec_case= 0;
1700 if ((mode < 2 || leftright)
1701#ifdef Honor_FLT_ROUNDS
1702 && rounding == 1
1703#endif
1704 )
1705 {
1706 if (!word1(&u) && !(word0(&u) & Bndry_mask) &&
1707 word0(&u) & (Exp_mask & ~Exp_msk1)
1708 )
1709 {
1710 /* The special case */
1711 b2+= Log2P;
1712 s2+= Log2P;
1713 spec_case= 1;
1714 }
1715 }
1716
1717 /*
1718 Arrange for convenient computation of quotients:
1719 shift left if necessary so divisor has 4 leading 0 bits.
1720
1721 Perhaps we should just compute leading 28 bits of S once
1722 a nd for all and pass them and a shift to quorem, so it
1723 can do shifts and ors to compute the numerator for q.
1724 */
1725 if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
1726 i= 32 - i;
1727 if (i > 4)
1728 {
1729 i-= 4;
1730 b2+= i;
1731 m2+= i;
1732 s2+= i;
1733 }
1734 else if (i < 4)
1735 {
1736 i+= 28;
1737 b2+= i;
1738 m2+= i;
1739 s2+= i;
1740 }
1741 if (b2 > 0)
1742 b= lshift(b, b2, &alloc);
1743 if (s2 > 0)
1744 S= lshift(S, s2, &alloc);
1745 if (k_check)
1746 {
1747 if (cmp(b,S) < 0)
1748 {
1749 k--;
1750 /* we botched the k estimate */
1751 b= multadd(b, 10, 0, &alloc);
1752 if (leftright)
1753 mhi= multadd(mhi, 10, 0, &alloc);
1754 ilim= ilim1;
1755 }
1756 }
1757 if (ilim <= 0 && (mode == 3 || mode == 5))
1758 {
1759 if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
1760 {
1761 /* no digits, fcvt style */
1762no_digits:
1763 k= -1 - ndigits;
1764 goto ret;
1765 }
1766one_digit:
1767 *s++= '1';
1768 k++;
1769 goto ret;
1770 }
1771 if (leftright)
1772 {
1773 if (m2 > 0)
1774 mhi= lshift(mhi, m2, &alloc);
1775
1776 /*
1777 Compute mlo -- check for special case that d is a normalized power of 2.
1778 */
1779
1780 mlo= mhi;
1781 if (spec_case)
1782 {
1783 mhi= Balloc(mhi->k, &alloc);
1784 Bcopy(mhi, mlo);
1785 mhi= lshift(mhi, Log2P, &alloc);
1786 }
1787
1788 for (i= 1;;i++)
1789 {
1790 dig= quorem(b,S) + '0';
1791 /* Do we yet have the shortest decimal string that will round to d? */
1792 j= cmp(b, mlo);
1793 delta= diff(S, mhi, &alloc);
1794 j1= delta->sign ? 1 : cmp(b, delta);
1795 Bfree(delta, &alloc);
1796 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
1797#ifdef Honor_FLT_ROUNDS
1798 && rounding >= 1
1799#endif
1800 )
1801 {
1802 if (dig == '9')
1803 goto round_9_up;
1804 if (j > 0)
1805 dig++;
1806 *s++= dig;
1807 goto ret;
1808 }
1809 if (j < 0 || (j == 0 && mode != 1 && !(word1(&u) & 1)))
1810 {
1811 if (!b->p.x[0] && b->wds <= 1)
1812 {
1813 goto accept_dig;
1814 }
1815#ifdef Honor_FLT_ROUNDS
1816 if (mode > 1)
1817 switch (rounding) {
1818 case 0: goto accept_dig;
1819 case 2: goto keep_dig;
1820 }
1821#endif /*Honor_FLT_ROUNDS*/
1822 if (j1 > 0)
1823 {
1824 b= lshift(b, 1, &alloc);
1825 j1= cmp(b, S);
1826 if ((j1 > 0 || (j1 == 0 && dig & 1))
1827 && dig++ == '9')
1828 goto round_9_up;
1829 }
1830accept_dig:
1831 *s++= dig;
1832 goto ret;
1833 }
1834 if (j1 > 0)
1835 {
1836#ifdef Honor_FLT_ROUNDS
1837 if (!rounding)
1838 goto accept_dig;
1839#endif
1840 if (dig == '9')
1841 { /* possible if i == 1 */
1842round_9_up:
1843 *s++= '9';
1844 goto roundoff;
1845 }
1846 *s++= dig + 1;
1847 goto ret;
1848 }
1849#ifdef Honor_FLT_ROUNDS
1850keep_dig:
1851#endif
1852 *s++= dig;
1853 if (i == ilim)
1854 break;
1855 b= multadd(b, 10, 0, &alloc);
1856 if (mlo == mhi)
1857 mlo= mhi= multadd(mhi, 10, 0, &alloc);
1858 else
1859 {
1860 mlo= multadd(mlo, 10, 0, &alloc);
1861 mhi= multadd(mhi, 10, 0, &alloc);
1862 }
1863 }
1864 }
1865 else
1866 for (i= 1;; i++)
1867 {
1868 *s++= dig= quorem(b,S) + '0';
1869 if (!b->p.x[0] && b->wds <= 1)
1870 {
1871 goto ret;
1872 }
1873 if (i >= ilim)
1874 break;
1875 b= multadd(b, 10, 0, &alloc);
1876 }
1877
1878 /* Round off last digit */
1879
1880#ifdef Honor_FLT_ROUNDS
1881 switch (rounding) {
1882 case 0: goto trimzeros;
1883 case 2: goto roundoff;
1884 }
1885#endif
1886 b= lshift(b, 1, &alloc);
1887 j= cmp(b, S);
1888 if (j > 0 || (j == 0 && dig & 1))
1889 {
1890roundoff:
1891 while (*--s == '9')
1892 if (s == s0)
1893 {
1894 k++;
1895 *s++= '1';
1896 goto ret;
1897 }
1898 ++*s++;
1899 }
1900 else
1901 {
1902#ifdef Honor_FLT_ROUNDS
1903trimzeros:
1904#endif
1905 while (*--s == '0');
1906 s++;
1907 }
1908ret:
1909 Bfree(S, &alloc);
1910 if (mhi)
1911 {
1912 if (mlo && mlo != mhi)
1913 Bfree(mlo, &alloc);
1914 Bfree(mhi, &alloc);
1915 }
1916ret1:
1917 Bfree(b, &alloc);
1918 *s= 0;
1919 *decpt= k + 1;
1920 if (rve)
1921 *rve= s;
1922 return s0;
1923}
1924