1/* Copyright (c) 2007, 2012, Oracle and/or its affiliates. All rights reserved.
2 Copyright (c) 2017, MariaDB Corporation.
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
43/**
44 Appears to suffice to not call malloc() in most cases.
45 @todo
46 see if it is possible to get rid of malloc().
47 this constant is sufficient to avoid malloc() on all inputs I have tried.
48*/
49#define DTOA_BUFF_SIZE (460 * sizeof(void *))
50
51/* Magic value returned by dtoa() to indicate overflow */
52#define DTOA_OVERFLOW 9999
53
54static double my_strtod_int(const char *, char **, int *, char *, size_t);
55static char *dtoa(double, int, int, int *, int *, char **, char *, size_t);
56static void dtoa_free(char *, char *, size_t);
57
58/**
59 @brief
60 Converts a given floating point number to a zero-terminated string
61 representation using the 'f' format.
62
63 @details
64 This function is a wrapper around dtoa() to do the same as
65 sprintf(to, "%-.*f", precision, x), though the conversion is usually more
66 precise. The only difference is in handling [-,+]infinity and nan values,
67 in which case we print '0\0' to the output string and indicate an overflow.
68
69 @param x the input floating point number.
70 @param precision the number of digits after the decimal point.
71 All properties of sprintf() apply:
72 - if the number of significant digits after the decimal
73 point is less than precision, the resulting string is
74 right-padded with zeros
75 - if the precision is 0, no decimal point appears
76 - if a decimal point appears, at least one digit appears
77 before it
78 @param to pointer to the output buffer. The longest string which
79 my_fcvt() can return is FLOATING_POINT_BUFFER bytes
80 (including the terminating '\0').
81 @param error if not NULL, points to a location where the status of
82 conversion is stored upon return.
83 FALSE successful conversion
84 TRUE the input number is [-,+]infinity or nan.
85 The output string in this case is always '0'.
86 @return number of written characters (excluding terminating '\0')
87*/
88
89size_t my_fcvt(double x, int precision, char *to, my_bool *error)
90{
91 int decpt, sign, len, i;
92 char *res, *src, *end, *dst= to;
93 char buf[DTOA_BUFF_SIZE];
94 DBUG_ASSERT(precision >= 0 && precision < DECIMAL_NOT_SPECIFIED && to != NULL);
95
96 res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
97
98 if (decpt == DTOA_OVERFLOW)
99 {
100 dtoa_free(res, buf, sizeof(buf));
101 *to++= '0';
102 *to= '\0';
103 if (error != NULL)
104 *error= TRUE;
105 return 1;
106 }
107
108 src= res;
109 len= (int)(end - src);
110
111 if (sign)
112 *dst++= '-';
113
114 if (decpt <= 0)
115 {
116 *dst++= '0';
117 *dst++= '.';
118 for (i= decpt; i < 0; i++)
119 *dst++= '0';
120 }
121
122 for (i= 1; i <= len; i++)
123 {
124 *dst++= *src++;
125 if (i == decpt && i < len)
126 *dst++= '.';
127 }
128 while (i++ <= decpt)
129 *dst++= '0';
130
131 if (precision > 0)
132 {
133 if (len <= decpt)
134 *dst++= '.';
135
136 for (i= precision - MY_MAX(0, (len - decpt)); i > 0; i--)
137 *dst++= '0';
138 }
139
140 *dst= '\0';
141 if (error != NULL)
142 *error= FALSE;
143
144 dtoa_free(res, buf, sizeof(buf));
145
146 return dst - to;
147}
148
149/**
150 @brief
151 Converts a given floating point number to a zero-terminated string
152 representation with a given field width using the 'e' format
153 (aka scientific notation) or the 'f' one.
154
155 @details
156 The format is chosen automatically to provide the most number of significant
157 digits (and thus, precision) with a given field width. In many cases, the
158 result is similar to that of sprintf(to, "%g", x) with a few notable
159 differences:
160 - the conversion is usually more precise than C library functions.
161 - there is no 'precision' argument. instead, we specify the number of
162 characters available for conversion (i.e. a field width).
163 - the result never exceeds the specified field width. If the field is too
164 short to contain even a rounded decimal representation, my_gcvt()
165 indicates overflow and truncates the output string to the specified width.
166 - float-type arguments are handled differently than double ones. For a
167 float input number (i.e. when the 'type' argument is MY_GCVT_ARG_FLOAT)
168 we deliberately limit the precision of conversion by FLT_DIG digits to
169 avoid garbage past the significant digits.
170 - unlike sprintf(), in cases where the 'e' format is preferred, we don't
171 zero-pad the exponent to save space for significant digits. The '+' sign
172 for a positive exponent does not appear for the same reason.
173
174 @param x the input floating point number.
175 @param type is either MY_GCVT_ARG_FLOAT or MY_GCVT_ARG_DOUBLE.
176 Specifies the type of the input number (see notes above).
177 @param width field width in characters. The minimal field width to
178 hold any number representation (albeit rounded) is 7
179 characters ("-Ne-NNN").
180 @param to pointer to the output buffer. The result is always
181 zero-terminated, and the longest returned string is thus
182 'width + 1' bytes.
183 @param error if not NULL, points to a location where the status of
184 conversion is stored upon return.
185 FALSE successful conversion
186 TRUE the input number is [-,+]infinity or nan.
187 The output string in this case is always '0'.
188 @return number of written characters (excluding terminating '\0')
189
190 @todo
191 Check if it is possible and makes sense to do our own rounding on top of
192 dtoa() instead of calling dtoa() twice in (rare) cases when the resulting
193 string representation does not fit in the specified field width and we want
194 to re-round the input number with fewer significant digits. Examples:
195
196 my_gcvt(-9e-3, ..., 4, ...);
197 my_gcvt(-9e-3, ..., 2, ...);
198 my_gcvt(1.87e-3, ..., 4, ...);
199 my_gcvt(55, ..., 1, ...);
200
201 We do our best to minimize such cases by:
202
203 - passing to dtoa() the field width as the number of significant digits
204
205 - removing the sign of the number early (and decreasing the width before
206 passing it to dtoa())
207
208 - choosing the proper format to preserve the most number of significant
209 digits.
210*/
211
212size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
213 my_bool *error)
214{
215 int decpt, sign, len, exp_len;
216 char *res, *src, *end, *dst= to, *dend= dst + width;
217 char buf[DTOA_BUFF_SIZE];
218 my_bool have_space, force_e_format;
219 DBUG_ASSERT(width > 0 && to != NULL);
220
221 /* We want to remove '-' from equations early */
222 if (x < 0.)
223 width--;
224
225 res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : MY_MIN(width, FLT_DIG),
226 &decpt, &sign, &end, buf, sizeof(buf));
227 if (decpt == DTOA_OVERFLOW)
228 {
229 dtoa_free(res, buf, sizeof(buf));
230 *to++= '0';
231 *to= '\0';
232 if (error != NULL)
233 *error= TRUE;
234 return 1;
235 }
236
237 if (error != NULL)
238 *error= FALSE;
239
240 src= res;
241 len= (int)(end - res);
242
243 /*
244 Number of digits in the exponent from the 'e' conversion.
245 The sign of the exponent is taken into account separetely, we don't need
246 to count it here.
247 */
248 exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
249
250 /*
251 Do we have enough space for all digits in the 'f' format?
252 Let 'len' be the number of significant digits returned by dtoa,
253 and F be the length of the resulting decimal representation.
254 Consider the following cases:
255 1. decpt <= 0, i.e. we have "0.NNN" => F = len - decpt + 2
256 2. 0 < decpt < len, i.e. we have "NNN.NNN" => F = len + 1
257 3. len <= decpt, i.e. we have "NNN00" => F = decpt
258 */
259 have_space= (decpt <= 0 ? len - decpt + 2 :
260 decpt > 0 && decpt < len ? len + 1 :
261 decpt) <= width;
262 /*
263 The following is true when no significant digits can be placed with the
264 specified field width using the 'f' format, and the 'e' format
265 will not be truncated.
266 */
267 force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len);
268 /*
269 Assume that we don't have enough space to place all significant digits in
270 the 'f' format. We have to choose between the 'e' format and the 'f' one
271 to keep as many significant digits as possible.
272 Let E and F be the lengths of decimal representation in the 'e' and 'f'
273 formats, respectively. We want to use the 'f' format if, and only if F <= E.
274 Consider the following cases:
275 1. decpt <= 0.
276 F = len - decpt + 2 (see above)
277 E = len + (len > 1) + 1 + 1 (decpt <= -99) + (decpt <= -9) + 1
278 ("N.NNe-MMM")
279 (F <= E) <=> (len == 1 && decpt >= -1) || (len > 1 && decpt >= -2)
280 We also need to ensure that if the 'f' format is chosen,
281 the field width allows us to place at least one significant digit
282 (i.e. width > 2 - decpt). If not, we prefer the 'e' format.
283 2. 0 < decpt < len
284 F = len + 1 (see above)
285 E = len + 1 + 1 + ... ("N.NNeMMM")
286 F is always less than E.
287 3. len <= decpt <= width
288 In this case we have enough space to represent the number in the 'f'
289 format, so we prefer it with some exceptions.
290 4. width < decpt
291 The number cannot be represented in the 'f' format at all, always use
292 the 'e' 'one.
293 */
294 if ((have_space ||
295 /*
296 Not enough space, let's see if the 'f' format provides the most number
297 of significant digits.
298 */
299 ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
300 (len > 1 || !force_e_format)))) &&
301 !force_e_format)) &&
302
303 /*
304 Use the 'e' format in some cases even if we have enough space for the
305 'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
306 */
307 (!have_space || (decpt >= -MAX_DECPT_FOR_F_FORMAT + 1 &&
308 (decpt <= MAX_DECPT_FOR_F_FORMAT || len > decpt))))
309 {
310 /* 'f' format */
311 int i;
312
313 width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0);
314
315 /* Do we have to truncate any digits? */
316 if (width < len)
317 {
318 if (width < decpt)
319 {
320 if (error != NULL)
321 *error= TRUE;
322 width= decpt;
323 }
324
325 /*
326 We want to truncate (len - width) least significant digits after the
327 decimal point. For this we are calling dtoa with mode=5, passing the
328 number of significant digits = (len-decpt) - (len-width) = width-decpt
329 */
330 dtoa_free(res, buf, sizeof(buf));
331 res= dtoa(x, 5, width - decpt, &decpt, &sign, &end, buf, sizeof(buf));
332 src= res;
333 len= (int)(end - res);
334 }
335
336 if (len == 0)
337 {
338 /* Underflow. Just print '0' and exit */
339 *dst++= '0';
340 goto end;
341 }
342
343 /*
344 At this point we are sure we have enough space to put all digits
345 returned by dtoa
346 */
347 if (sign && dst < dend)
348 *dst++= '-';
349 if (decpt <= 0)
350 {
351 if (dst < dend)
352 *dst++= '0';
353 if (len > 0 && dst < dend)
354 *dst++= '.';
355 for (; decpt < 0 && dst < dend; decpt++)
356 *dst++= '0';
357 }
358
359 for (i= 1; i <= len && dst < dend; i++)
360 {
361 *dst++= *src++;
362 if (i == decpt && i < len && dst < dend)
363 *dst++= '.';
364 }
365 while (i++ <= decpt && dst < dend)
366 *dst++= '0';
367 }
368 else
369 {
370 /* 'e' format */
371 int decpt_sign= 0;
372
373 if (--decpt < 0)
374 {
375 decpt= -decpt;
376 width--;
377 decpt_sign= 1;
378 }
379 width-= 1 + exp_len; /* eNNN */
380
381 if (len > 1)
382 width--;
383
384 if (width <= 0)
385 {
386 /* Overflow */
387 if (error != NULL)
388 *error= TRUE;
389 width= 0;
390 }
391
392 /* Do we have to truncate any digits? */
393 if (width < len)
394 {
395 /* Yes, re-convert with a smaller width */
396 dtoa_free(res, buf, sizeof(buf));
397 res= dtoa(x, 4, width, &decpt, &sign, &end, buf, sizeof(buf));
398 src= res;
399 len= (int)(end - res);
400 if (--decpt < 0)
401 decpt= -decpt;
402 }
403 /*
404 At this point we are sure we have enough space to put all digits
405 returned by dtoa
406 */
407 if (sign && dst < dend)
408 *dst++= '-';
409 if (dst < dend)
410 *dst++= *src++;
411 if (len > 1 && dst < dend)
412 {
413 *dst++= '.';
414 while (src < end && dst < dend)
415 *dst++= *src++;
416 }
417 if (dst < dend)
418 *dst++= 'e';
419 if (decpt_sign && dst < dend)
420 *dst++= '-';
421
422 if (decpt >= 100 && dst < dend)
423 {
424 *dst++= decpt / 100 + '0';
425 decpt%= 100;
426 if (dst < dend)
427 *dst++= decpt / 10 + '0';
428 }
429 else if (decpt >= 10 && dst < dend)
430 *dst++= decpt / 10 + '0';
431 if (dst < dend)
432 *dst++= decpt % 10 + '0';
433
434 }
435
436end:
437 dtoa_free(res, buf, sizeof(buf));
438 *dst= '\0';
439
440 return dst - to;
441}
442
443/**
444 @brief
445 Converts string to double (string does not have to be zero-terminated)
446
447 @details
448 This is a wrapper around dtoa's version of strtod().
449
450 @param str input string
451 @param end address of a pointer to the first character after the input
452 string. Upon return the pointer is set to point to the first
453 rejected character.
454 @param error Upon return is set to EOVERFLOW in case of underflow or
455 overflow.
456
457 @return The resulting double value. In case of underflow, 0.0 is
458 returned. In case overflow, signed DBL_MAX is returned.
459*/
460
461double my_strtod(const char *str, char **end, int *error)
462{
463 char buf[DTOA_BUFF_SIZE];
464 double res;
465 DBUG_ASSERT(end != NULL && ((str != NULL && *end != NULL) ||
466 (str == NULL && *end == NULL)) &&
467 error != NULL);
468
469 res= my_strtod_int(str, end, error, buf, sizeof(buf));
470 return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
471}
472
473
474double my_atof(const char *nptr)
475{
476 int error;
477 const char *end= nptr+65535; /* Should be enough */
478 return (my_strtod(nptr, (char**) &end, &error));
479}
480
481
482/****************************************************************
483 *
484 * The author of this software is David M. Gay.
485 *
486 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
487 *
488 * Permission to use, copy, modify, and distribute this software for any
489 * purpose without fee is hereby granted, provided that this entire notice
490 * is included in all copies of any software which is or includes a copy
491 * or modification of this software and in all copies of the supporting
492 * documentation for such software.
493 *
494 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
495 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
496 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
497 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
498 *
499 ***************************************************************/
500/* Please send bug reports to David M. Gay (dmg at acm dot org,
501 * with " at " changed at "@" and " dot " changed to "."). */
502
503/*
504 Original copy of the software is located at http://www.netlib.org/fp/dtoa.c
505 It was adjusted to serve MySQL server needs:
506 * strtod() was modified to not expect a zero-terminated string.
507 It now honors 'se' (end of string) argument as the input parameter,
508 not just as the output one.
509 * in dtoa(), in case of overflow/underflow/NaN result string now contains "0";
510 decpt is set to DTOA_OVERFLOW to indicate overflow.
511 * support for VAX, IBM mainframe and 16-bit hardware removed
512 * we always assume that 64-bit integer type is available
513 * support for Kernigan-Ritchie style headers (pre-ANSI compilers)
514 removed
515 * all gcc warnings ironed out
516 * we always assume multithreaded environment, so we had to change
517 memory allocation procedures to use stack in most cases;
518 malloc is used as the last resort.
519 * pow5mult rewritten to use pre-calculated pow5 list instead of
520 the one generated on the fly.
521*/
522
523
524/*
525 On a machine with IEEE extended-precision registers, it is
526 necessary to specify double-precision (53-bit) rounding precision
527 before invoking strtod or dtoa. If the machine uses (the equivalent
528 of) Intel 80x87 arithmetic, the call
529 _control87(PC_53, MCW_PC);
530 does this with many compilers. Whether this or another call is
531 appropriate depends on the compiler; for this to work, it may be
532 necessary to #include "float.h" or another system-dependent header
533 file.
534*/
535
536/*
537 #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
538 and dtoa should round accordingly.
539 #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
540 and Honor_FLT_ROUNDS is not #defined.
541
542 TODO: check if we can get rid of the above two
543*/
544
545typedef int32 Long;
546typedef uint32 ULong;
547typedef int64 LLong;
548typedef uint64 ULLong;
549
550typedef union { double d; ULong L[2]; } U;
551
552#if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) && \
553 (__FLOAT_WORD_ORDER == __BIG_ENDIAN))
554#define word0(x) (x)->L[0]
555#define word1(x) (x)->L[1]
556#else
557#define word0(x) (x)->L[1]
558#define word1(x) (x)->L[0]
559#endif
560
561#define dval(x) (x)->d
562
563/* #define P DBL_MANT_DIG */
564/* Ten_pmax= floor(P*log(2)/log(5)) */
565/* Bletch= (highest power of 2 < DBL_MAX_10_EXP) / 16 */
566/* Quick_max= floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
567/* Int_max= floor(P*log(FLT_RADIX)/log(10) - 1) */
568
569#define Exp_shift 20
570#define Exp_shift1 20
571#define Exp_msk1 0x100000
572#define Exp_mask 0x7ff00000
573#define P 53
574#define Bias 1023
575#define Emin (-1022)
576#define Exp_1 0x3ff00000
577#define Exp_11 0x3ff00000
578#define Ebits 11
579#define Frac_mask 0xfffff
580#define Frac_mask1 0xfffff
581#define Ten_pmax 22
582#define Bletch 0x10
583#define Bndry_mask 0xfffff
584#define Bndry_mask1 0xfffff
585#define LSB 1
586#define Sign_bit 0x80000000
587#define Log2P 1
588#define Tiny1 1
589#define Quick_max 14
590#define Int_max 14
591
592#ifndef Flt_Rounds
593#ifdef FLT_ROUNDS
594#define Flt_Rounds FLT_ROUNDS
595#else
596#define Flt_Rounds 1
597#endif
598#endif /*Flt_Rounds*/
599
600#ifdef Honor_FLT_ROUNDS
601#define Rounding rounding
602#undef Check_FLT_ROUNDS
603#define Check_FLT_ROUNDS
604#else
605#define Rounding Flt_Rounds
606#endif
607
608#define rounded_product(a,b) a*= b
609#define rounded_quotient(a,b) a/= b
610
611#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
612#define Big1 0xffffffff
613#define FFFFFFFF 0xffffffffUL
614
615/* This is tested to be enough for dtoa */
616
617#define Kmax 15
618
619#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
620 2*sizeof(int) + y->wds*sizeof(ULong))
621
622/* Arbitrary-length integer */
623
624typedef struct Bigint
625{
626 union {
627 ULong *x; /* points right after this Bigint object */
628 struct Bigint *next; /* to maintain free lists */
629 } p;
630 int k; /* 2^k = maxwds */
631 int maxwds; /* maximum length in 32-bit words */
632 int sign; /* not zero if number is negative */
633 int wds; /* current length in 32-bit words */
634} Bigint;
635
636
637/* A simple stack-memory based allocator for Bigints */
638
639typedef struct Stack_alloc
640{
641 char *begin;
642 char *free;
643 char *end;
644 /*
645 Having list of free blocks lets us reduce maximum required amount
646 of memory from ~4000 bytes to < 1680 (tested on x86).
647 */
648 Bigint *freelist[Kmax+1];
649} Stack_alloc;
650
651
652/*
653 Try to allocate object on stack, and resort to malloc if all
654 stack memory is used. Ensure allocated objects to be aligned by the pointer
655 size in order to not break the alignment rules when storing a pointer to a
656 Bigint.
657*/
658
659static Bigint *Balloc(int k, Stack_alloc *alloc)
660{
661 Bigint *rv;
662 DBUG_ASSERT(k <= Kmax);
663 if (k <= Kmax && alloc->freelist[k])
664 {
665 rv= alloc->freelist[k];
666 alloc->freelist[k]= rv->p.next;
667 }
668 else
669 {
670 int x, len;
671
672 x= 1 << k;
673 len= MY_ALIGN(sizeof(Bigint) + x * sizeof(ULong), SIZEOF_CHARP);
674
675 if (alloc->free + len <= alloc->end)
676 {
677 rv= (Bigint*) alloc->free;
678 alloc->free+= len;
679 }
680 else
681 rv= (Bigint*) malloc(len);
682
683 rv->k= k;
684 rv->maxwds= x;
685 }
686 rv->sign= rv->wds= 0;
687 rv->p.x= (ULong*) (rv + 1);
688 return rv;
689}
690
691
692/*
693 If object was allocated on stack, try putting it to the free
694 list. Otherwise call free().
695*/
696
697static void Bfree(Bigint *v, Stack_alloc *alloc)
698{
699 char *gptr= (char*) v; /* generic pointer */
700 if (gptr < alloc->begin || gptr >= alloc->end)
701 free(gptr);
702 else if (v->k <= Kmax)
703 {
704 /*
705 Maintain free lists only for stack objects: this way we don't
706 have to bother with freeing lists in the end of dtoa;
707 heap should not be used normally anyway.
708 */
709 v->p.next= alloc->freelist[v->k];
710 alloc->freelist[v->k]= v;
711 }
712}
713
714
715/*
716 This is to place return value of dtoa in: tries to use stack
717 as well, but passes by free lists management and just aligns len by
718 the pointer size in order to not break the alignment rules when storing a
719 pointer to a Bigint.
720*/
721
722static char *dtoa_alloc(int i, Stack_alloc *alloc)
723{
724 char *rv;
725 int aligned_size= MY_ALIGN(i, SIZEOF_CHARP);
726 if (alloc->free + aligned_size <= alloc->end)
727 {
728 rv= alloc->free;
729 alloc->free+= aligned_size;
730 }
731 else
732 rv= malloc(i);
733 return rv;
734}
735
736
737/*
738 dtoa_free() must be used to free values s returned by dtoa()
739 This is the counterpart of dtoa_alloc()
740*/
741
742static void dtoa_free(char *gptr, char *buf, size_t buf_size)
743{
744 if (gptr < buf || gptr >= buf + buf_size)
745 free(gptr);
746}
747
748
749/* Bigint arithmetic functions */
750
751/* Multiply by m and add a */
752
753static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc)
754{
755 int i, wds;
756 ULong *x;
757 ULLong carry, y;
758 Bigint *b1;
759
760 wds= b->wds;
761 x= b->p.x;
762 i= 0;
763 carry= a;
764 do
765 {
766 y= *x * (ULLong)m + carry;
767 carry= y >> 32;
768 *x++= (ULong)(y & FFFFFFFF);
769 }
770 while (++i < wds);
771 if (carry)
772 {
773 if (wds >= b->maxwds)
774 {
775 b1= Balloc(b->k+1, alloc);
776 Bcopy(b1, b);
777 Bfree(b, alloc);
778 b= b1;
779 }
780 b->p.x[wds++]= (ULong) carry;
781 b->wds= wds;
782 }
783 return b;
784}
785
786/**
787 Converts a string to Bigint.
788
789 Now we have nd0 digits, starting at s, followed by a
790 decimal point, followed by nd-nd0 digits.
791 Unless nd0 == nd, in which case we have a number of the form:
792 ".xxxxxx" or "xxxxxx."
793
794 @param s Input string, already partially parsed by my_strtod_int().
795 @param nd0 Number of digits before decimal point.
796 @param nd Total number of digits.
797 @param y9 Pre-computed value of the first nine digits.
798 @param alloc Stack allocator for Bigints.
799 */
800static Bigint *s2b(const char *s, int nd0, int nd, ULong y9, Stack_alloc *alloc)
801{
802 Bigint *b;
803 int i, k;
804 Long x, y;
805
806 x= (nd + 8) / 9;
807 for (k= 0, y= 1; x > y; y <<= 1, k++) ;
808 b= Balloc(k, alloc);
809 b->p.x[0]= y9;
810 b->wds= 1;
811
812 i= 9;
813 if (9 < nd0)
814 {
815 s+= 9;
816 do
817 b= multadd(b, 10, *s++ - '0', alloc);
818 while (++i < nd0);
819 s++; /* skip '.' */
820 }
821 else
822 s+= 10;
823 /* now do the fractional part */
824 for(; i < nd; i++)
825 b= multadd(b, 10, *s++ - '0', alloc);
826 return b;
827}
828
829
830static int hi0bits(register ULong x)
831{
832 register int k= 0;
833
834 if (!(x & 0xffff0000))
835 {
836 k= 16;
837 x<<= 16;
838 }
839 if (!(x & 0xff000000))
840 {
841 k+= 8;
842 x<<= 8;
843 }
844 if (!(x & 0xf0000000))
845 {
846 k+= 4;
847 x<<= 4;
848 }
849 if (!(x & 0xc0000000))
850 {
851 k+= 2;
852 x<<= 2;
853 }
854 if (!(x & 0x80000000))
855 {
856 k++;
857 if (!(x & 0x40000000))
858 return 32;
859 }
860 return k;
861}
862
863
864static int lo0bits(ULong *y)
865{
866 register int k;
867 register ULong x= *y;
868
869 if (x & 7)
870 {
871 if (x & 1)
872 return 0;
873 if (x & 2)
874 {
875 *y= x >> 1;
876 return 1;
877 }
878 *y= x >> 2;
879 return 2;
880 }
881 k= 0;
882 if (!(x & 0xffff))
883 {
884 k= 16;
885 x>>= 16;
886 }
887 if (!(x & 0xff))
888 {
889 k+= 8;
890 x>>= 8;
891 }
892 if (!(x & 0xf))
893 {
894 k+= 4;
895 x>>= 4;
896 }
897 if (!(x & 0x3))
898 {
899 k+= 2;
900 x>>= 2;
901 }
902 if (!(x & 1))
903 {
904 k++;
905 x>>= 1;
906 if (!x)
907 return 32;
908 }
909 *y= x;
910 return k;
911}
912
913
914/* Convert integer to Bigint number */
915
916static Bigint *i2b(int i, Stack_alloc *alloc)
917{
918 Bigint *b;
919
920 b= Balloc(1, alloc);
921 b->p.x[0]= i;
922 b->wds= 1;
923 return b;
924}
925
926
927/* Multiply two Bigint numbers */
928
929static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc)
930{
931 Bigint *c;
932 int k, wa, wb, wc;
933 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
934 ULong y;
935 ULLong carry, z;
936
937 if (a->wds < b->wds)
938 {
939 c= a;
940 a= b;
941 b= c;
942 }
943 k= a->k;
944 wa= a->wds;
945 wb= b->wds;
946 wc= wa + wb;
947 if (wc > a->maxwds)
948 k++;
949 c= Balloc(k, alloc);
950 for (x= c->p.x, xa= x + wc; x < xa; x++)
951 *x= 0;
952 xa= a->p.x;
953 xae= xa + wa;
954 xb= b->p.x;
955 xbe= xb + wb;
956 xc0= c->p.x;
957 for (; xb < xbe; xc0++)
958 {
959 if ((y= *xb++))
960 {
961 x= xa;
962 xc= xc0;
963 carry= 0;
964 do
965 {
966 z= *x++ * (ULLong)y + *xc + carry;
967 carry= z >> 32;
968 *xc++= (ULong) (z & FFFFFFFF);
969 }
970 while (x < xae);
971 *xc= (ULong) carry;
972 }
973 }
974 for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ;
975 c->wds= wc;
976 return c;
977}
978
979
980/*
981 Precalculated array of powers of 5: tested to be enough for
982 vasting majority of dtoa_r cases.
983*/
984
985static ULong powers5[]=
986{
987 625UL,
988
989 390625UL,
990
991 2264035265UL, 35UL,
992
993 2242703233UL, 762134875UL, 1262UL,
994
995 3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL,
996
997 781532673UL, 64985353UL, 253049085UL, 594863151UL, 3553621484UL,
998 3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL,
999
1000 2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL,
1001 3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL,
1002 1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL,
1003 2161952759UL, 4100910556UL, 1608314830UL, 349175UL
1004};
1005
1006
1007static Bigint p5_a[]=
1008{
1009 /* { x } - k - maxwds - sign - wds */
1010 { { powers5 }, 1, 1, 0, 1 },
1011 { { powers5 + 1 }, 1, 1, 0, 1 },
1012 { { powers5 + 2 }, 1, 2, 0, 2 },
1013 { { powers5 + 4 }, 2, 3, 0, 3 },
1014 { { powers5 + 7 }, 3, 5, 0, 5 },
1015 { { powers5 + 12 }, 4, 10, 0, 10 },
1016 { { powers5 + 22 }, 5, 19, 0, 19 }
1017};
1018
1019#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
1020
1021static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc)
1022{
1023 Bigint *b1, *p5, *p51=NULL;
1024 int i;
1025 static int p05[3]= { 5, 25, 125 };
1026 my_bool overflow= FALSE;
1027
1028 if ((i= k & 3))
1029 b= multadd(b, p05[i-1], 0, alloc);
1030
1031 if (!(k>>= 2))
1032 return b;
1033 p5= p5_a;
1034 for (;;)
1035 {
1036 if (k & 1)
1037 {
1038 b1= mult(b, p5, alloc);
1039 Bfree(b, alloc);
1040 b= b1;
1041 }
1042 if (!(k>>= 1))
1043 break;
1044 /* Calculate next power of 5 */
1045 if (overflow)
1046 {
1047 p51= mult(p5, p5, alloc);
1048 Bfree(p5, alloc);
1049 p5= p51;
1050 }
1051 else if (p5 < p5_a + P5A_MAX)
1052 ++p5;
1053 else if (p5 == p5_a + P5A_MAX)
1054 {
1055 p5= mult(p5, p5, alloc);
1056 overflow= TRUE;
1057 }
1058 }
1059 if (p51)
1060 Bfree(p51, alloc);
1061 return b;
1062}
1063
1064
1065static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc)
1066{
1067 int i, k1, n, n1;
1068 Bigint *b1;
1069 ULong *x, *x1, *xe, z;
1070
1071 n= k >> 5;
1072 k1= b->k;
1073 n1= n + b->wds + 1;
1074 for (i= b->maxwds; n1 > i; i<<= 1)
1075 k1++;
1076 b1= Balloc(k1, alloc);
1077 x1= b1->p.x;
1078 for (i= 0; i < n; i++)
1079 *x1++= 0;
1080 x= b->p.x;
1081 xe= x + b->wds;
1082 if (k&= 0x1f)
1083 {
1084 k1= 32 - k;
1085 z= 0;
1086 do
1087 {
1088 *x1++= *x << k | z;
1089 z= *x++ >> k1;
1090 }
1091 while (x < xe);
1092 if ((*x1= z))
1093 ++n1;
1094 }
1095 else
1096 do
1097 *x1++= *x++;
1098 while (x < xe);
1099 b1->wds= n1 - 1;
1100 Bfree(b, alloc);
1101 return b1;
1102}
1103
1104
1105static int cmp(Bigint *a, Bigint *b)
1106{
1107 ULong *xa, *xa0, *xb, *xb0;
1108 int i, j;
1109
1110 i= a->wds;
1111 j= b->wds;
1112 if (i-= j)
1113 return i;
1114 xa0= a->p.x;
1115 xa= xa0 + j;
1116 xb0= b->p.x;
1117 xb= xb0 + j;
1118 for (;;)
1119 {
1120 if (*--xa != *--xb)
1121 return *xa < *xb ? -1 : 1;
1122 if (xa <= xa0)
1123 break;
1124 }
1125 return 0;
1126}
1127
1128
1129static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc)
1130{
1131 Bigint *c;
1132 int i, wa, wb;
1133 ULong *xa, *xae, *xb, *xbe, *xc;
1134 ULLong borrow, y;
1135
1136 i= cmp(a,b);
1137 if (!i)
1138 {
1139 c= Balloc(0, alloc);
1140 c->wds= 1;
1141 c->p.x[0]= 0;
1142 return c;
1143 }
1144 if (i < 0)
1145 {
1146 c= a;
1147 a= b;
1148 b= c;
1149 i= 1;
1150 }
1151 else
1152 i= 0;
1153 c= Balloc(a->k, alloc);
1154 c->sign= i;
1155 wa= a->wds;
1156 xa= a->p.x;
1157 xae= xa + wa;
1158 wb= b->wds;
1159 xb= b->p.x;
1160 xbe= xb + wb;
1161 xc= c->p.x;
1162 borrow= 0;
1163 do
1164 {
1165 y= (ULLong)*xa++ - *xb++ - borrow;
1166 borrow= y >> 32 & (ULong)1;
1167 *xc++= (ULong) (y & FFFFFFFF);
1168 }
1169 while (xb < xbe);
1170 while (xa < xae)
1171 {
1172 y= *xa++ - borrow;
1173 borrow= y >> 32 & (ULong)1;
1174 *xc++= (ULong) (y & FFFFFFFF);
1175 }
1176 while (!*--xc)
1177 wa--;
1178 c->wds= wa;
1179 return c;
1180}
1181
1182
1183static double ulp(U *x)
1184{
1185 register Long L;
1186 U u;
1187
1188 L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1;
1189 word0(&u) = L;
1190 word1(&u) = 0;
1191 return dval(&u);
1192}
1193
1194
1195static double b2d(Bigint *a, int *e)
1196{
1197 ULong *xa, *xa0, w, y, z;
1198 int k;
1199 U d;
1200#define d0 word0(&d)
1201#define d1 word1(&d)
1202
1203 xa0= a->p.x;
1204 xa= xa0 + a->wds;
1205 y= *--xa;
1206 k= hi0bits(y);
1207 *e= 32 - k;
1208 if (k < Ebits)
1209 {
1210 d0= Exp_1 | y >> (Ebits - k);
1211 w= xa > xa0 ? *--xa : 0;
1212 d1= y << ((32-Ebits) + k) | w >> (Ebits - k);
1213 goto ret_d;
1214 }
1215 z= xa > xa0 ? *--xa : 0;
1216 if (k-= Ebits)
1217 {
1218 d0= Exp_1 | y << k | z >> (32 - k);
1219 y= xa > xa0 ? *--xa : 0;
1220 d1= z << k | y >> (32 - k);
1221 }
1222 else
1223 {
1224 d0= Exp_1 | y;
1225 d1= z;
1226 }
1227 ret_d:
1228#undef d0
1229#undef d1
1230 return dval(&d);
1231}
1232
1233
1234static Bigint *d2b(U *d, int *e, int *bits, Stack_alloc *alloc)
1235{
1236 Bigint *b;
1237 int de, k;
1238 ULong *x, y, z;
1239 int i;
1240#define d0 word0(d)
1241#define d1 word1(d)
1242
1243 b= Balloc(1, alloc);
1244 x= b->p.x;
1245
1246 z= d0 & Frac_mask;
1247 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1248 if ((de= (int)(d0 >> Exp_shift)))
1249 z|= Exp_msk1;
1250 if ((y= d1))
1251 {
1252 if ((k= lo0bits(&y)))
1253 {
1254 x[0]= y | z << (32 - k);
1255 z>>= k;
1256 }
1257 else
1258 x[0]= y;
1259 i= b->wds= (x[1]= z) ? 2 : 1;
1260 }
1261 else
1262 {
1263 k= lo0bits(&z);
1264 x[0]= z;
1265 i= b->wds= 1;
1266 k+= 32;
1267 }
1268 if (de)
1269 {
1270 *e= de - Bias - (P-1) + k;
1271 *bits= P - k;
1272 }
1273 else
1274 {
1275 *e= de - Bias - (P-1) + 1 + k;
1276 *bits= 32*i - hi0bits(x[i-1]);
1277 }
1278 return b;
1279#undef d0
1280#undef d1
1281}
1282
1283
1284static double ratio(Bigint *a, Bigint *b)
1285{
1286 U da, db;
1287 int k, ka, kb;
1288
1289 dval(&da)= b2d(a, &ka);
1290 dval(&db)= b2d(b, &kb);
1291 k= ka - kb + 32*(a->wds - b->wds);
1292 if (k > 0)
1293 word0(&da)+= (ULong)(k*Exp_msk1 * 1.0);
1294 else
1295 {
1296 k= -k;
1297 word0(&db)+= k*Exp_msk1;
1298 }
1299 return dval(&da) / dval(&db);
1300}
1301
1302static const double tens[] =
1303{
1304 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1305 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1306 1e20, 1e21, 1e22
1307};
1308
1309static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 };
1310static const double tinytens[]=
1311{ 1e-16, 1e-32, 1e-64, 1e-128,
1312 9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1313};
1314/*
1315 The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
1316 flag unnecessarily. It leads to a song and dance at the end of strtod.
1317*/
1318#define Scale_Bit 0x10
1319#define n_bigtens 5
1320
1321/*
1322 strtod for IEEE--arithmetic machines.
1323
1324 This strtod returns a nearest machine number to the input decimal
1325 string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1326 rule.
1327
1328 Inspired loosely by William D. Clinger's paper "How to Read Floating
1329 Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1330
1331 Modifications:
1332
1333 1. We only require IEEE (not IEEE double-extended).
1334 2. We get by with floating-point arithmetic in a case that
1335 Clinger missed -- when we're computing d * 10^n
1336 for a small integer d and the integer n is not too
1337 much larger than 22 (the maximum integer k for which
1338 we can represent 10^k exactly), we may be able to
1339 compute (d*10^k) * 10^(e-k) with just one roundoff.
1340 3. Rather than a bit-at-a-time adjustment of the binary
1341 result in the hard case, we use floating-point
1342 arithmetic to determine the adjustment to within
1343 one bit; only in really hard cases do we need to
1344 compute a second residual.
1345 4. Because of 3., we don't need a large table of powers of 10
1346 for ten-to-e (just some small tables, e.g. of 10^k
1347 for 0 <= k <= 22).
1348*/
1349
1350static double my_strtod_int(const char *s00, char **se, int *error, char *buf, size_t buf_size)
1351{
1352 int scale;
1353 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, UNINIT_VAR(c), dsign,
1354 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1355 const char *s, *s0, *s1, *end = *se;
1356 double aadj, aadj1;
1357 U aadj2, adj, rv, rv0;
1358 Long L;
1359 ULong y, z;
1360 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1361#ifdef SET_INEXACT
1362 int inexact, oldinexact;
1363#endif
1364#ifdef Honor_FLT_ROUNDS
1365 int rounding;
1366#endif
1367 Stack_alloc alloc;
1368
1369 *error= 0;
1370
1371 alloc.begin= alloc.free= buf;
1372 alloc.end= buf + buf_size;
1373 memset(alloc.freelist, 0, sizeof(alloc.freelist));
1374
1375 sign= nz0= nz= 0;
1376 dval(&rv)= 0.;
1377 for (s= s00; s < end; s++)
1378 switch (*s) {
1379 case '-':
1380 sign= 1;
1381 /* fall through */
1382 case '+':
1383 s++;
1384 goto break2;
1385 case '\t':
1386 case '\n':
1387 case '\v':
1388 case '\f':
1389 case '\r':
1390 case ' ':
1391 continue;
1392 default:
1393 goto break2;
1394 }
1395 break2:
1396 if (s >= end)
1397 goto ret0;
1398
1399 if (*s == '0')
1400 {
1401 nz0= 1;
1402 while (++s < end && *s == '0') ;
1403 if (s >= end)
1404 goto ret;
1405 }
1406 s0= s;
1407 y= z= 0;
1408 for (nd= nf= 0; s < end && (c= *s) >= '0' && c <= '9'; nd++, s++)
1409 if (nd < 9)
1410 y= 10*y + c - '0';
1411 else if (nd < 16)
1412 z= 10*z + c - '0';
1413 nd0= nd;
1414 if (s < end && c == '.')
1415 {
1416 ++s;
1417 if (!nd)
1418 {
1419 for (; s < end && (c= *s) == '0'; ++s)
1420 nz++;
1421 if (s < end && (c= *s) > '0' && c <= '9')
1422 {
1423 s0= s;
1424 nf+= nz;
1425 nz= 0;
1426 goto have_dig;
1427 }
1428 goto dig_done;
1429 }
1430 for (; s < end && (c= *s) >= '0' && c <= '9'; ++s)
1431 {
1432 have_dig:
1433 /*
1434 Here we are parsing the fractional part.
1435 We can stop counting digits after a while: the extra digits
1436 will not contribute to the actual result produced by s2b().
1437 We have to continue scanning, in case there is an exponent part.
1438 */
1439 if (nd < 2 * DBL_DIG)
1440 {
1441 nz++;
1442 if (c-= '0')
1443 {
1444 nf+= nz;
1445 for (i= 1; i < nz; i++)
1446 if (nd++ < 9)
1447 y*= 10;
1448 else if (nd <= DBL_DIG + 1)
1449 z*= 10;
1450 if (nd++ < 9)
1451 y= 10*y + c;
1452 else if (nd <= DBL_DIG + 1)
1453 z= 10*z + c;
1454 nz= 0;
1455 }
1456 }
1457 }
1458 }
1459 dig_done:
1460 e= 0;
1461 if (s < end && (c == 'e' || c == 'E'))
1462 {
1463 if (!nd && !nz && !nz0)
1464 goto ret0;
1465 s00= s;
1466 esign= 0;
1467 if (++s < end)
1468 switch (c= *s) {
1469 case '-': esign= 1;
1470 /* fall through */
1471 case '+': c= *++s;
1472 }
1473 if (s < end && c >= '0' && c <= '9')
1474 {
1475 while (s < end && c == '0')
1476 c= *++s;
1477 if (s < end && c > '0' && c <= '9') {
1478 L= c - '0';
1479 s1= s;
1480 while (++s < end && (c= *s) >= '0' && c <= '9')
1481 L= 10*L + c - '0';
1482 if (s - s1 > 8 || L > 19999)
1483 /* Avoid confusion from exponents
1484 * so large that e might overflow.
1485 */
1486 e= 19999; /* safe for 16 bit ints */
1487 else
1488 e= (int)L;
1489 if (esign)
1490 e= -e;
1491 }
1492 else
1493 e= 0;
1494 }
1495 else
1496 s= s00;
1497 }
1498 if (!nd)
1499 {
1500 if (!nz && !nz0)
1501 {
1502 ret0:
1503 s= s00;
1504 sign= 0;
1505 }
1506 goto ret;
1507 }
1508 e1= e -= nf;
1509
1510 /*
1511 Now we have nd0 digits, starting at s0, followed by a
1512 decimal point, followed by nd-nd0 digits. The number we're
1513 after is the integer represented by those digits times
1514 10**e
1515 */
1516
1517 if (!nd0)
1518 nd0= nd;
1519 k= nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1520 dval(&rv)= y;
1521 if (k > 9)
1522 {
1523#ifdef SET_INEXACT
1524 if (k > DBL_DIG)
1525 oldinexact = get_inexact();
1526#endif
1527 dval(&rv)= tens[k - 9] * dval(&rv) + z;
1528 }
1529 bd0= 0;
1530 if (nd <= DBL_DIG
1531#ifndef Honor_FLT_ROUNDS
1532 && Flt_Rounds == 1
1533#endif
1534 )
1535 {
1536 if (!e)
1537 goto ret;
1538 if (e > 0)
1539 {
1540 if (e <= Ten_pmax)
1541 {
1542#ifdef Honor_FLT_ROUNDS
1543 /* round correctly FLT_ROUNDS = 2 or 3 */
1544 if (sign)
1545 {
1546 rv.d= -rv.d;
1547 sign= 0;
1548 }
1549#endif
1550 /* rv = */ rounded_product(dval(&rv), tens[e]);
1551 goto ret;
1552 }
1553 i= DBL_DIG - nd;
1554 if (e <= Ten_pmax + i)
1555 {
1556 /*
1557 A fancier test would sometimes let us do
1558 this for larger i values.
1559 */
1560#ifdef Honor_FLT_ROUNDS
1561 /* round correctly FLT_ROUNDS = 2 or 3 */
1562 if (sign)
1563 {
1564 rv.d= -rv.d;
1565 sign= 0;
1566 }
1567#endif
1568 e-= i;
1569 dval(&rv)*= tens[i];
1570 /* rv = */ rounded_product(dval(&rv), tens[e]);
1571 goto ret;
1572 }
1573 }
1574#ifndef Inaccurate_Divide
1575 else if (e >= -Ten_pmax)
1576 {
1577#ifdef Honor_FLT_ROUNDS
1578 /* round correctly FLT_ROUNDS = 2 or 3 */
1579 if (sign)
1580 {
1581 rv.d= -rv.d;
1582 sign= 0;
1583 }
1584#endif
1585 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
1586 goto ret;
1587 }
1588#endif
1589 }
1590 e1+= nd - k;
1591
1592#ifdef SET_INEXACT
1593 inexact= 1;
1594 if (k <= DBL_DIG)
1595 oldinexact= get_inexact();
1596#endif
1597 scale= 0;
1598#ifdef Honor_FLT_ROUNDS
1599 if ((rounding= Flt_Rounds) >= 2)
1600 {
1601 if (sign)
1602 rounding= rounding == 2 ? 0 : 2;
1603 else
1604 if (rounding != 2)
1605 rounding= 0;
1606 }
1607#endif
1608
1609 /* Get starting approximation = rv * 10**e1 */
1610
1611 if (e1 > 0)
1612 {
1613 if ((i= e1 & 15))
1614 dval(&rv)*= tens[i];
1615 if (e1&= ~15)
1616 {
1617 if (e1 > DBL_MAX_10_EXP)
1618 {
1619 ovfl:
1620 *error= EOVERFLOW;
1621 /* Can't trust HUGE_VAL */
1622#ifdef Honor_FLT_ROUNDS
1623 switch (rounding)
1624 {
1625 case 0: /* toward 0 */
1626 case 3: /* toward -infinity */
1627 word0(&rv)= Big0;
1628 word1(&rv)= Big1;
1629 break;
1630 default:
1631 word0(&rv)= Exp_mask;
1632 word1(&rv)= 0;
1633 }
1634#else /*Honor_FLT_ROUNDS*/
1635 word0(&rv)= Exp_mask;
1636 word1(&rv)= 0;
1637#endif /*Honor_FLT_ROUNDS*/
1638#ifdef SET_INEXACT
1639 /* set overflow bit */
1640 dval(&rv0)= 1e300;
1641 dval(&rv0)*= dval(&rv0);
1642#endif
1643 if (bd0)
1644 goto retfree;
1645 goto ret;
1646 }
1647 e1>>= 4;
1648 for(j= 0; e1 > 1; j++, e1>>= 1)
1649 if (e1 & 1)
1650 dval(&rv)*= bigtens[j];
1651 /* The last multiplication could overflow. */
1652 word0(&rv)-= P*Exp_msk1;
1653 dval(&rv)*= bigtens[j];
1654 if ((z= word0(&rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
1655 goto ovfl;
1656 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
1657 {
1658 /* set to largest number (Can't trust DBL_MAX) */
1659 word0(&rv)= Big0;
1660 word1(&rv)= Big1;
1661 }
1662 else
1663 word0(&rv)+= P*Exp_msk1;
1664 }
1665 }
1666 else if (e1 < 0)
1667 {
1668 e1= -e1;
1669 if ((i= e1 & 15))
1670 dval(&rv)/= tens[i];
1671 if ((e1>>= 4))
1672 {
1673 if (e1 >= 1 << n_bigtens)
1674 goto undfl;
1675 if (e1 & Scale_Bit)
1676 scale= 2 * P;
1677 for(j= 0; e1 > 0; j++, e1>>= 1)
1678 if (e1 & 1)
1679 dval(&rv)*= tinytens[j];
1680 if (scale && (j = 2 * P + 1 - ((word0(&rv) & Exp_mask) >> Exp_shift)) > 0)
1681 {
1682 /* scaled rv is denormal; zap j low bits */
1683 if (j >= 32)
1684 {
1685 word1(&rv)= 0;
1686 if (j >= 53)
1687 word0(&rv)= (P + 2) * Exp_msk1;
1688 else
1689 word0(&rv)&= 0xffffffff << (j - 32);
1690 }
1691 else
1692 word1(&rv)&= 0xffffffff << j;
1693 }
1694 if (!dval(&rv))
1695 {
1696 undfl:
1697 dval(&rv)= 0.;
1698 if (bd0)
1699 goto retfree;
1700 goto ret;
1701 }
1702 }
1703 }
1704
1705 /* Now the hard part -- adjusting rv to the correct value.*/
1706
1707 /* Put digits into bd: true value = bd * 10^e */
1708
1709 bd0= s2b(s0, nd0, nd, y, &alloc);
1710
1711 for(;;)
1712 {
1713 bd= Balloc(bd0->k, &alloc);
1714 Bcopy(bd, bd0);
1715 bb= d2b(&rv, &bbe, &bbbits, &alloc); /* rv = bb * 2^bbe */
1716 bs= i2b(1, &alloc);
1717
1718 if (e >= 0)
1719 {
1720 bb2= bb5= 0;
1721 bd2= bd5= e;
1722 }
1723 else
1724 {
1725 bb2= bb5= -e;
1726 bd2= bd5= 0;
1727 }
1728 if (bbe >= 0)
1729 bb2+= bbe;
1730 else
1731 bd2-= bbe;
1732 bs2= bb2;
1733#ifdef Honor_FLT_ROUNDS
1734 if (rounding != 1)
1735 bs2++;
1736#endif
1737 j= bbe - scale;
1738 i= j + bbbits - 1; /* logb(rv) */
1739 if (i < Emin) /* denormal */
1740 j+= P - Emin;
1741 else
1742 j= P + 1 - bbbits;
1743 bb2+= j;
1744 bd2+= j;
1745 bd2+= scale;
1746 i= bb2 < bd2 ? bb2 : bd2;
1747 if (i > bs2)
1748 i= bs2;
1749 if (i > 0)
1750 {
1751 bb2-= i;
1752 bd2-= i;
1753 bs2-= i;
1754 }
1755 if (bb5 > 0)
1756 {
1757 bs= pow5mult(bs, bb5, &alloc);
1758 bb1= mult(bs, bb, &alloc);
1759 Bfree(bb, &alloc);
1760 bb= bb1;
1761 }
1762 if (bb2 > 0)
1763 bb= lshift(bb, bb2, &alloc);
1764 if (bd5 > 0)
1765 bd= pow5mult(bd, bd5, &alloc);
1766 if (bd2 > 0)
1767 bd= lshift(bd, bd2, &alloc);
1768 if (bs2 > 0)
1769 bs= lshift(bs, bs2, &alloc);
1770 delta= diff(bb, bd, &alloc);
1771 dsign= delta->sign;
1772 delta->sign= 0;
1773 i= cmp(delta, bs);
1774#ifdef Honor_FLT_ROUNDS
1775 if (rounding != 1)
1776 {
1777 if (i < 0)
1778 {
1779 /* Error is less than an ulp */
1780 if (!delta->p.x[0] && delta->wds <= 1)
1781 {
1782 /* exact */
1783#ifdef SET_INEXACT
1784 inexact= 0;
1785#endif
1786 break;
1787 }
1788 if (rounding)
1789 {
1790 if (dsign)
1791 {
1792 adj.d= 1.;
1793 goto apply_adj;
1794 }
1795 }
1796 else if (!dsign)
1797 {
1798 adj.d= -1.;
1799 if (!word1(&rv) && !(word0(&rv) & Frac_mask))
1800 {
1801 y= word0(&rv) & Exp_mask;
1802 if (!scale || y > 2*P*Exp_msk1)
1803 {
1804 delta= lshift(delta, Log2P, &alloc);
1805 if (cmp(delta, bs) <= 0)
1806 adj.d= -0.5;
1807 }
1808 }
1809 apply_adj:
1810 if (scale && (y= word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1)
1811 word0(&adj)+= (2 * P + 1) * Exp_msk1 - y;
1812 dval(&rv)+= adj.d * ulp(&rv);
1813 }
1814 break;
1815 }
1816 adj.d= ratio(delta, bs);
1817 if (adj.d < 1.)
1818 adj.d= 1.;
1819 if (adj.d <= 0x7ffffffe)
1820 {
1821 /* adj = rounding ? ceil(adj) : floor(adj); */
1822 y= adj.d;
1823 if (y != adj.d)
1824 {
1825 if (!((rounding >> 1) ^ dsign))
1826 y++;
1827 adj.d= y;
1828 }
1829 }
1830 if (scale && (y= word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1)
1831 word0(&adj)+= (2 * P + 1) * Exp_msk1 - y;
1832 adj.d*= ulp(&rv);
1833 if (dsign)
1834 dval(&rv)+= adj.d;
1835 else
1836 dval(&rv)-= adj.d;
1837 goto cont;
1838 }
1839#endif /*Honor_FLT_ROUNDS*/
1840
1841 if (i < 0)
1842 {
1843 /*
1844 Error is less than half an ulp -- check for special case of mantissa
1845 a power of two.
1846 */
1847 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask ||
1848 (word0(&rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1)
1849 {
1850#ifdef SET_INEXACT
1851 if (!delta->x[0] && delta->wds <= 1)
1852 inexact= 0;
1853#endif
1854 break;
1855 }
1856 if (!delta->p.x[0] && delta->wds <= 1)
1857 {
1858 /* exact result */
1859#ifdef SET_INEXACT
1860 inexact= 0;
1861#endif
1862 break;
1863 }
1864 delta= lshift(delta, Log2P, &alloc);
1865 if (cmp(delta, bs) > 0)
1866 goto drop_down;
1867 break;
1868 }
1869 if (i == 0)
1870 {
1871 /* exactly half-way between */
1872 if (dsign)
1873 {
1874 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 &&
1875 word1(&rv) ==
1876 ((scale && (y = word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1) ?
1877 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
1878 0xffffffff))
1879 {
1880 /*boundary case -- increment exponent*/
1881 word0(&rv)= (word0(&rv) & Exp_mask) + Exp_msk1;
1882 word1(&rv) = 0;
1883 dsign = 0;
1884 break;
1885 }
1886 }
1887 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv))
1888 {
1889 drop_down:
1890 /* boundary case -- decrement exponent */
1891 if (scale)
1892 {
1893 L= word0(&rv) & Exp_mask;
1894 if (L <= (2 *P + 1) * Exp_msk1)
1895 {
1896 if (L > (P + 2) * Exp_msk1)
1897 /* round even ==> accept rv */
1898 break;
1899 /* rv = smallest denormal */
1900 goto undfl;
1901 }
1902 }
1903 L= (word0(&rv) & Exp_mask) - Exp_msk1;
1904 word0(&rv)= L | Bndry_mask1;
1905 word1(&rv)= 0xffffffff;
1906 break;
1907 }
1908 if (!(word1(&rv) & LSB))
1909 break;
1910 if (dsign)
1911 dval(&rv)+= ulp(&rv);
1912 else
1913 {
1914 dval(&rv)-= ulp(&rv);
1915 if (!dval(&rv))
1916 goto undfl;
1917 }
1918 dsign= 1 - dsign;
1919 break;
1920 }
1921 if ((aadj= ratio(delta, bs)) <= 2.)
1922 {
1923 if (dsign)
1924 aadj= aadj1= 1.;
1925 else if (word1(&rv) || word0(&rv) & Bndry_mask)
1926 {
1927 if (word1(&rv) == Tiny1 && !word0(&rv))
1928 goto undfl;
1929 aadj= 1.;
1930 aadj1= -1.;
1931 }
1932 else
1933 {
1934 /* special case -- power of FLT_RADIX to be rounded down... */
1935 if (aadj < 2. / FLT_RADIX)
1936 aadj= 1. / FLT_RADIX;
1937 else
1938 aadj*= 0.5;
1939 aadj1= -aadj;
1940 }
1941 }
1942 else
1943 {
1944 aadj*= 0.5;
1945 aadj1= dsign ? aadj : -aadj;
1946#ifdef Check_FLT_ROUNDS
1947 switch (Rounding)
1948 {
1949 case 2: /* towards +infinity */
1950 aadj1-= 0.5;
1951 break;
1952 case 0: /* towards 0 */
1953 case 3: /* towards -infinity */
1954 aadj1+= 0.5;
1955 }
1956#else
1957 if (Flt_Rounds == 0)
1958 aadj1+= 0.5;
1959#endif /*Check_FLT_ROUNDS*/
1960 }
1961 y= word0(&rv) & Exp_mask;
1962
1963 /* Check for overflow */
1964
1965 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
1966 {
1967 dval(&rv0)= dval(&rv);
1968 word0(&rv)-= P * Exp_msk1;
1969 adj.d= aadj1 * ulp(&rv);
1970 dval(&rv)+= adj.d;
1971 if ((word0(&rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P))
1972 {
1973 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
1974 goto ovfl;
1975 word0(&rv)= Big0;
1976 word1(&rv)= Big1;
1977 goto cont;
1978 }
1979 else
1980 word0(&rv)+= P * Exp_msk1;
1981 }
1982 else
1983 {
1984 if (scale && y <= 2 * P * Exp_msk1)
1985 {
1986 if (aadj <= 0x7fffffff)
1987 {
1988 if ((z= (ULong) aadj) <= 0)
1989 z= 1;
1990 aadj= z;
1991 aadj1= dsign ? aadj : -aadj;
1992 }
1993 dval(&aadj2) = aadj1;
1994 word0(&aadj2)+= (2 * P + 1) * Exp_msk1 - y;
1995 aadj1= dval(&aadj2);
1996 adj.d= aadj1 * ulp(&rv);
1997 dval(&rv)+= adj.d;
1998 if (rv.d == 0.)
1999 goto undfl;
2000 }
2001 else
2002 {
2003 adj.d= aadj1 * ulp(&rv);
2004 dval(&rv)+= adj.d;
2005 }
2006 }
2007 z= word0(&rv) & Exp_mask;
2008#ifndef SET_INEXACT
2009 if (!scale)
2010 if (y == z)
2011 {
2012 /* Can we stop now? */
2013 L= (Long)aadj;
2014 aadj-= L;
2015 /* The tolerances below are conservative. */
2016 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask)
2017 {
2018 if (aadj < .4999999 || aadj > .5000001)
2019 break;
2020 }
2021 else if (aadj < .4999999 / FLT_RADIX)
2022 break;
2023 }
2024#endif
2025 cont:
2026 Bfree(bb, &alloc);
2027 Bfree(bd, &alloc);
2028 Bfree(bs, &alloc);
2029 Bfree(delta, &alloc);
2030 }
2031#ifdef SET_INEXACT
2032 if (inexact)
2033 {
2034 if (!oldinexact)
2035 {
2036 word0(&rv0)= Exp_1 + (70 << Exp_shift);
2037 word1(&rv0)= 0;
2038 dval(&rv0)+= 1.;
2039 }
2040 }
2041 else if (!oldinexact)
2042 clear_inexact();
2043#endif
2044 if (scale)
2045 {
2046 word0(&rv0)= Exp_1 - 2 * P * Exp_msk1;
2047 word1(&rv0)= 0;
2048 dval(&rv)*= dval(&rv0);
2049 }
2050#ifdef SET_INEXACT
2051 if (inexact && !(word0(&rv) & Exp_mask))
2052 {
2053 /* set underflow bit */
2054 dval(&rv0)= 1e-300;
2055 dval(&rv0)*= dval(&rv0);
2056 }
2057#endif
2058 retfree:
2059 Bfree(bb, &alloc);
2060 Bfree(bd, &alloc);
2061 Bfree(bs, &alloc);
2062 Bfree(bd0, &alloc);
2063 Bfree(delta, &alloc);
2064 ret:
2065 *se= (char *)s;
2066 return sign ? -dval(&rv) : dval(&rv);
2067}
2068
2069
2070static int quorem(Bigint *b, Bigint *S)
2071{
2072 int n;
2073 ULong *bx, *bxe, q, *sx, *sxe;
2074 ULLong borrow, carry, y, ys;
2075
2076 n= S->wds;
2077 if (b->wds < n)
2078 return 0;
2079 sx= S->p.x;
2080 sxe= sx + --n;
2081 bx= b->p.x;
2082 bxe= bx + n;
2083 q= *bxe / (*sxe + 1); /* ensure q <= true quotient */
2084 if (q)
2085 {
2086 borrow= 0;
2087 carry= 0;
2088 do
2089 {
2090 ys= *sx++ * (ULLong)q + carry;
2091 carry= ys >> 32;
2092 y= *bx - (ys & FFFFFFFF) - borrow;
2093 borrow= y >> 32 & (ULong)1;
2094 *bx++= (ULong) (y & FFFFFFFF);
2095 }
2096 while (sx <= sxe);
2097 if (!*bxe)
2098 {
2099 bx= b->p.x;
2100 while (--bxe > bx && !*bxe)
2101 --n;
2102 b->wds= n;
2103 }
2104 }
2105 if (cmp(b, S) >= 0)
2106 {
2107 q++;
2108 borrow= 0;
2109 carry= 0;
2110 bx= b->p.x;
2111 sx= S->p.x;
2112 do
2113 {
2114 ys= *sx++ + carry;
2115 carry= ys >> 32;
2116 y= *bx - (ys & FFFFFFFF) - borrow;
2117 borrow= y >> 32 & (ULong)1;
2118 *bx++= (ULong) (y & FFFFFFFF);
2119 }
2120 while (sx <= sxe);
2121 bx= b->p.x;
2122 bxe= bx + n;
2123 if (!*bxe)
2124 {
2125 while (--bxe > bx && !*bxe)
2126 --n;
2127 b->wds= n;
2128 }
2129 }
2130 return q;
2131}
2132
2133
2134/*
2135 dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2136
2137 Inspired by "How to Print Floating-Point Numbers Accurately" by
2138 Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2139
2140 Modifications:
2141 1. Rather than iterating, we use a simple numeric overestimate
2142 to determine k= floor(log10(d)). We scale relevant
2143 quantities using O(log2(k)) rather than O(k) multiplications.
2144 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2145 try to generate digits strictly left to right. Instead, we
2146 compute with fewer bits and propagate the carry if necessary
2147 when rounding the final digit up. This is often faster.
2148 3. Under the assumption that input will be rounded nearest,
2149 mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2150 That is, we allow equality in stopping tests when the
2151 round-nearest rule will give the same floating-point value
2152 as would satisfaction of the stopping test with strict
2153 inequality.
2154 4. We remove common factors of powers of 2 from relevant
2155 quantities.
2156 5. When converting floating-point integers less than 1e16,
2157 we use floating-point arithmetic rather than resorting
2158 to multiple-precision integers.
2159 6. When asked to produce fewer than 15 digits, we first try
2160 to get by with floating-point arithmetic; we resort to
2161 multiple-precision integer arithmetic only if we cannot
2162 guarantee that the floating-point calculation has given
2163 the correctly rounded result. For k requested digits and
2164 "uniformly" distributed input, the probability is
2165 something like 10^(k-15) that we must resort to the Long
2166 calculation.
2167 */
2168
2169static char *dtoa(double dd, int mode, int ndigits, int *decpt, int *sign,
2170 char **rve, char *buf, size_t buf_size)
2171{
2172 /*
2173 Arguments ndigits, decpt, sign are similar to those
2174 of ecvt and fcvt; trailing zeros are suppressed from
2175 the returned string. If not null, *rve is set to point
2176 to the end of the return value. If d is +-Infinity or NaN,
2177 then *decpt is set to DTOA_OVERFLOW.
2178
2179 mode:
2180 0 ==> shortest string that yields d when read in
2181 and rounded to nearest.
2182 1 ==> like 0, but with Steele & White stopping rule;
2183 e.g. with IEEE P754 arithmetic , mode 0 gives
2184 1e23 whereas mode 1 gives 9.999999999999999e22.
2185 2 ==> MY_MAX(1,ndigits) significant digits. This gives a
2186 return value similar to that of ecvt, except
2187 that trailing zeros are suppressed.
2188 3 ==> through ndigits past the decimal point. This
2189 gives a return value similar to that from fcvt,
2190 except that trailing zeros are suppressed, and
2191 ndigits can be negative.
2192 4,5 ==> similar to 2 and 3, respectively, but (in
2193 round-nearest mode) with the tests of mode 0 to
2194 possibly return a shorter string that rounds to d.
2195 With IEEE arithmetic and compilation with
2196 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2197 as modes 2 and 3 when FLT_ROUNDS != 1.
2198 6-9 ==> Debugging modes similar to mode - 4: don't try
2199 fast floating-point estimate (if applicable).
2200
2201 Values of mode other than 0-9 are treated as mode 0.
2202
2203 Sufficient space is allocated to the return value
2204 to hold the suppressed trailing zeros.
2205 */
2206
2207 int bbits, b2, b5, be, dig, i, ieps, UNINIT_VAR(ilim), ilim0,
2208 UNINIT_VAR(ilim1), j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2209 spec_case, try_quick;
2210 Long L;
2211 int denorm;
2212 ULong x;
2213 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2214 U d2, eps, u;
2215 double ds;
2216 char *s, *s0;
2217#ifdef Honor_FLT_ROUNDS
2218 int rounding;
2219#endif
2220 Stack_alloc alloc;
2221
2222 alloc.begin= alloc.free= buf;
2223 alloc.end= buf + buf_size;
2224 memset(alloc.freelist, 0, sizeof(alloc.freelist));
2225
2226 u.d= dd;
2227 if (word0(&u) & Sign_bit)
2228 {
2229 /* set sign for everything, including 0's and NaNs */
2230 *sign= 1;
2231 word0(&u) &= ~Sign_bit; /* clear sign bit */
2232 }
2233 else
2234 *sign= 0;
2235
2236 /* If infinity, set decpt to DTOA_OVERFLOW, if 0 set it to 1 */
2237 if (((word0(&u) & Exp_mask) == Exp_mask && (*decpt= DTOA_OVERFLOW)) ||
2238 (!dval(&u) && (*decpt= 1)))
2239 {
2240 /* Infinity, NaN, 0 */
2241 char *res= (char*) dtoa_alloc(2, &alloc);
2242 res[0]= '0';
2243 res[1]= '\0';
2244 if (rve)
2245 *rve= res + 1;
2246 return res;
2247 }
2248
2249#ifdef Honor_FLT_ROUNDS
2250 if ((rounding= Flt_Rounds) >= 2)
2251 {
2252 if (*sign)
2253 rounding= rounding == 2 ? 0 : 2;
2254 else
2255 if (rounding != 2)
2256 rounding= 0;
2257 }
2258#endif
2259
2260 b= d2b(&u, &be, &bbits, &alloc);
2261 if ((i= (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
2262 {
2263 dval(&d2)= dval(&u);
2264 word0(&d2) &= Frac_mask1;
2265 word0(&d2) |= Exp_11;
2266
2267 /*
2268 log(x) ~=~ log(1.5) + (x-1.5)/1.5
2269 log10(x) = log(x) / log(10)
2270 ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2271 log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
2272
2273 This suggests computing an approximation k to log10(d) by
2274
2275 k= (i - Bias)*0.301029995663981
2276 + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2277
2278 We want k to be too large rather than too small.
2279 The error in the first-order Taylor series approximation
2280 is in our favor, so we just round up the constant enough
2281 to compensate for any error in the multiplication of
2282 (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2283 and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2284 adding 1e-13 to the constant term more than suffices.
2285 Hence we adjust the constant term to 0.1760912590558.
2286 (We could get a more accurate k by invoking log10,
2287 but this is probably not worthwhile.)
2288 */
2289
2290 i-= Bias;
2291 denorm= 0;
2292 }
2293 else
2294 {
2295 /* d is denormalized */
2296
2297 i= bbits + be + (Bias + (P-1) - 1);
2298 x= i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
2299 : word1(&u) << (32 - i);
2300 dval(&d2)= x;
2301 word0(&d2)-= 31*Exp_msk1; /* adjust exponent */
2302 i-= (Bias + (P-1) - 1) + 1;
2303 denorm= 1;
2304 }
2305 ds= (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2306 k= (int)ds;
2307 if (ds < 0. && ds != k)
2308 k--; /* want k= floor(ds) */
2309 k_check= 1;
2310 if (k >= 0 && k <= Ten_pmax)
2311 {
2312 if (dval(&u) < tens[k])
2313 k--;
2314 k_check= 0;
2315 }
2316 j= bbits - i - 1;
2317 if (j >= 0)
2318 {
2319 b2= 0;
2320 s2= j;
2321 }
2322 else
2323 {
2324 b2= -j;
2325 s2= 0;
2326 }
2327 if (k >= 0)
2328 {
2329 b5= 0;
2330 s5= k;
2331 s2+= k;
2332 }
2333 else
2334 {
2335 b2-= k;
2336 b5= -k;
2337 s5= 0;
2338 }
2339 if (mode < 0 || mode > 9)
2340 mode= 0;
2341
2342#ifdef Check_FLT_ROUNDS
2343 try_quick= Rounding == 1;
2344#else
2345 try_quick= 1;
2346#endif
2347
2348 if (mode > 5)
2349 {
2350 mode-= 4;
2351 try_quick= 0;
2352 }
2353 leftright= 1;
2354 switch (mode) {
2355 case 0:
2356 case 1:
2357 ilim= ilim1= -1;
2358 i= 18;
2359 ndigits= 0;
2360 break;
2361 case 2:
2362 leftright= 0;
2363 /* fall through */
2364 case 4:
2365 if (ndigits <= 0)
2366 ndigits= 1;
2367 ilim= ilim1= i= ndigits;
2368 break;
2369 case 3:
2370 leftright= 0;
2371 /* fall through */
2372 case 5:
2373 i= ndigits + k + 1;
2374 ilim= i;
2375 ilim1= i - 1;
2376 if (i <= 0)
2377 i= 1;
2378 }
2379 s= s0= dtoa_alloc(i, &alloc);
2380
2381#ifdef Honor_FLT_ROUNDS
2382 if (mode > 1 && rounding != 1)
2383 leftright= 0;
2384#endif
2385
2386 if (ilim >= 0 && ilim <= Quick_max && try_quick)
2387 {
2388 /* Try to get by with floating-point arithmetic. */
2389 i= 0;
2390 dval(&d2)= dval(&u);
2391 k0= k;
2392 ilim0= ilim;
2393 ieps= 2; /* conservative */
2394 if (k > 0)
2395 {
2396 ds= tens[k&0xf];
2397 j= k >> 4;
2398 if (j & Bletch)
2399 {
2400 /* prevent overflows */
2401 j&= Bletch - 1;
2402 dval(&u)/= bigtens[n_bigtens-1];
2403 ieps++;
2404 }
2405 for (; j; j>>= 1, i++)
2406 {
2407 if (j & 1)
2408 {
2409 ieps++;
2410 ds*= bigtens[i];
2411 }
2412 }
2413 dval(&u)/= ds;
2414 }
2415 else if ((j1= -k))
2416 {
2417 dval(&u)*= tens[j1 & 0xf];
2418 for (j= j1 >> 4; j; j>>= 1, i++)
2419 {
2420 if (j & 1)
2421 {
2422 ieps++;
2423 dval(&u)*= bigtens[i];
2424 }
2425 }
2426 }
2427 if (k_check && dval(&u) < 1. && ilim > 0)
2428 {
2429 if (ilim1 <= 0)
2430 goto fast_failed;
2431 ilim= ilim1;
2432 k--;
2433 dval(&u)*= 10.;
2434 ieps++;
2435 }
2436 dval(&eps)= ieps*dval(&u) + 7.;
2437 word0(&eps)-= (P-1)*Exp_msk1;
2438 if (ilim == 0)
2439 {
2440 S= mhi= 0;
2441 dval(&u)-= 5.;
2442 if (dval(&u) > dval(&eps))
2443 goto one_digit;
2444 if (dval(&u) < -dval(&eps))
2445 goto no_digits;
2446 goto fast_failed;
2447 }
2448 if (leftright)
2449 {
2450 /* Use Steele & White method of only generating digits needed. */
2451 dval(&eps)= 0.5/tens[ilim-1] - dval(&eps);
2452 for (i= 0;;)
2453 {
2454 L= (Long) dval(&u);
2455 dval(&u)-= L;
2456 *s++= '0' + (int)L;
2457 if (dval(&u) < dval(&eps))
2458 goto ret1;
2459 if (1. - dval(&u) < dval(&eps))
2460 goto bump_up;
2461 if (++i >= ilim)
2462 break;
2463 dval(&eps)*= 10.;
2464 dval(&u)*= 10.;
2465 }
2466 }
2467 else
2468 {
2469 /* Generate ilim digits, then fix them up. */
2470 dval(&eps)*= tens[ilim-1];
2471 for (i= 1;; i++, dval(&u)*= 10.)
2472 {
2473 L= (Long)(dval(&u));
2474 if (!(dval(&u)-= L))
2475 ilim= i;
2476 *s++= '0' + (int)L;
2477 if (i == ilim)
2478 {
2479 if (dval(&u) > 0.5 + dval(&eps))
2480 goto bump_up;
2481 else if (dval(&u) < 0.5 - dval(&eps))
2482 {
2483 while (*--s == '0');
2484 s++;
2485 goto ret1;
2486 }
2487 break;
2488 }
2489 }
2490 }
2491 fast_failed:
2492 s= s0;
2493 dval(&u)= dval(&d2);
2494 k= k0;
2495 ilim= ilim0;
2496 }
2497
2498 /* Do we have a "small" integer? */
2499
2500 if (be >= 0 && k <= Int_max)
2501 {
2502 /* Yes. */
2503 ds= tens[k];
2504 if (ndigits < 0 && ilim <= 0)
2505 {
2506 S= mhi= 0;
2507 if (ilim < 0 || dval(&u) <= 5*ds)
2508 goto no_digits;
2509 goto one_digit;
2510 }
2511 for (i= 1;; i++, dval(&u)*= 10.)
2512 {
2513 L= (Long)(dval(&u) / ds);
2514 dval(&u)-= L*ds;
2515#ifdef Check_FLT_ROUNDS
2516 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2517 if (dval(&u) < 0)
2518 {
2519 L--;
2520 dval(&u)+= ds;
2521 }
2522#endif
2523 *s++= '0' + (int)L;
2524 if (!dval(&u))
2525 {
2526 break;
2527 }
2528 if (i == ilim)
2529 {
2530#ifdef Honor_FLT_ROUNDS
2531 if (mode > 1)
2532 {
2533 switch (rounding) {
2534 case 0: goto ret1;
2535 case 2: goto bump_up;
2536 }
2537 }
2538#endif
2539 dval(&u)+= dval(&u);
2540 if (dval(&u) > ds || (dval(&u) == ds && L & 1))
2541 {
2542bump_up:
2543 while (*--s == '9')
2544 if (s == s0)
2545 {
2546 k++;
2547 *s= '0';
2548 break;
2549 }
2550 ++*s++;
2551 }
2552 break;
2553 }
2554 }
2555 goto ret1;
2556 }
2557
2558 m2= b2;
2559 m5= b5;
2560 mhi= mlo= 0;
2561 if (leftright)
2562 {
2563 i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2564 b2+= i;
2565 s2+= i;
2566 mhi= i2b(1, &alloc);
2567 }
2568 if (m2 > 0 && s2 > 0)
2569 {
2570 i= m2 < s2 ? m2 : s2;
2571 b2-= i;
2572 m2-= i;
2573 s2-= i;
2574 }
2575 if (b5 > 0)
2576 {
2577 if (leftright)
2578 {
2579 if (m5 > 0)
2580 {
2581 mhi= pow5mult(mhi, m5, &alloc);
2582 b1= mult(mhi, b, &alloc);
2583 Bfree(b, &alloc);
2584 b= b1;
2585 }
2586 if ((j= b5 - m5))
2587 b= pow5mult(b, j, &alloc);
2588 }
2589 else
2590 b= pow5mult(b, b5, &alloc);
2591 }
2592 S= i2b(1, &alloc);
2593 if (s5 > 0)
2594 S= pow5mult(S, s5, &alloc);
2595
2596 /* Check for special case that d is a normalized power of 2. */
2597
2598 spec_case= 0;
2599 if ((mode < 2 || leftright)
2600#ifdef Honor_FLT_ROUNDS
2601 && rounding == 1
2602#endif
2603 )
2604 {
2605 if (!word1(&u) && !(word0(&u) & Bndry_mask) &&
2606 word0(&u) & (Exp_mask & ~Exp_msk1)
2607 )
2608 {
2609 /* The special case */
2610 b2+= Log2P;
2611 s2+= Log2P;
2612 spec_case= 1;
2613 }
2614 }
2615
2616 /*
2617 Arrange for convenient computation of quotients:
2618 shift left if necessary so divisor has 4 leading 0 bits.
2619
2620 Perhaps we should just compute leading 28 bits of S once
2621 a nd for all and pass them and a shift to quorem, so it
2622 can do shifts and ors to compute the numerator for q.
2623 */
2624 if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
2625 i= 32 - i;
2626 if (i > 4)
2627 {
2628 i-= 4;
2629 b2+= i;
2630 m2+= i;
2631 s2+= i;
2632 }
2633 else if (i < 4)
2634 {
2635 i+= 28;
2636 b2+= i;
2637 m2+= i;
2638 s2+= i;
2639 }
2640 if (b2 > 0)
2641 b= lshift(b, b2, &alloc);
2642 if (s2 > 0)
2643 S= lshift(S, s2, &alloc);
2644 if (k_check)
2645 {
2646 if (cmp(b,S) < 0)
2647 {
2648 k--;
2649 /* we botched the k estimate */
2650 b= multadd(b, 10, 0, &alloc);
2651 if (leftright)
2652 mhi= multadd(mhi, 10, 0, &alloc);
2653 ilim= ilim1;
2654 }
2655 }
2656 if (ilim <= 0 && (mode == 3 || mode == 5))
2657 {
2658 if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
2659 {
2660 /* no digits, fcvt style */
2661no_digits:
2662 k= -1 - ndigits;
2663 goto ret;
2664 }
2665one_digit:
2666 *s++= '1';
2667 k++;
2668 goto ret;
2669 }
2670 if (leftright)
2671 {
2672 if (m2 > 0)
2673 mhi= lshift(mhi, m2, &alloc);
2674
2675 /*
2676 Compute mlo -- check for special case that d is a normalized power of 2.
2677 */
2678
2679 mlo= mhi;
2680 if (spec_case)
2681 {
2682 mhi= Balloc(mhi->k, &alloc);
2683 Bcopy(mhi, mlo);
2684 mhi= lshift(mhi, Log2P, &alloc);
2685 }
2686
2687 for (i= 1;;i++)
2688 {
2689 dig= quorem(b,S) + '0';
2690 /* Do we yet have the shortest decimal string that will round to d? */
2691 j= cmp(b, mlo);
2692 delta= diff(S, mhi, &alloc);
2693 j1= delta->sign ? 1 : cmp(b, delta);
2694 Bfree(delta, &alloc);
2695 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
2696#ifdef Honor_FLT_ROUNDS
2697 && rounding >= 1
2698#endif
2699 )
2700 {
2701 if (dig == '9')
2702 goto round_9_up;
2703 if (j > 0)
2704 dig++;
2705 *s++= dig;
2706 goto ret;
2707 }
2708 if (j < 0 || (j == 0 && mode != 1 && !(word1(&u) & 1)))
2709 {
2710 if (!b->p.x[0] && b->wds <= 1)
2711 {
2712 goto accept_dig;
2713 }
2714#ifdef Honor_FLT_ROUNDS
2715 if (mode > 1)
2716 switch (rounding) {
2717 case 0: goto accept_dig;
2718 case 2: goto keep_dig;
2719 }
2720#endif /*Honor_FLT_ROUNDS*/
2721 if (j1 > 0)
2722 {
2723 b= lshift(b, 1, &alloc);
2724 j1= cmp(b, S);
2725 if ((j1 > 0 || (j1 == 0 && dig & 1))
2726 && dig++ == '9')
2727 goto round_9_up;
2728 }
2729accept_dig:
2730 *s++= dig;
2731 goto ret;
2732 }
2733 if (j1 > 0)
2734 {
2735#ifdef Honor_FLT_ROUNDS
2736 if (!rounding)
2737 goto accept_dig;
2738#endif
2739 if (dig == '9')
2740 { /* possible if i == 1 */
2741round_9_up:
2742 *s++= '9';
2743 goto roundoff;
2744 }
2745 *s++= dig + 1;
2746 goto ret;
2747 }
2748#ifdef Honor_FLT_ROUNDS
2749keep_dig:
2750#endif
2751 *s++= dig;
2752 if (i == ilim)
2753 break;
2754 b= multadd(b, 10, 0, &alloc);
2755 if (mlo == mhi)
2756 mlo= mhi= multadd(mhi, 10, 0, &alloc);
2757 else
2758 {
2759 mlo= multadd(mlo, 10, 0, &alloc);
2760 mhi= multadd(mhi, 10, 0, &alloc);
2761 }
2762 }
2763 }
2764 else
2765 for (i= 1;; i++)
2766 {
2767 *s++= dig= quorem(b,S) + '0';
2768 if (!b->p.x[0] && b->wds <= 1)
2769 {
2770 goto ret;
2771 }
2772 if (i >= ilim)
2773 break;
2774 b= multadd(b, 10, 0, &alloc);
2775 }
2776
2777 /* Round off last digit */
2778
2779#ifdef Honor_FLT_ROUNDS
2780 switch (rounding) {
2781 case 0: goto trimzeros;
2782 case 2: goto roundoff;
2783 }
2784#endif
2785 b= lshift(b, 1, &alloc);
2786 j= cmp(b, S);
2787 if (j > 0 || (j == 0 && dig & 1))
2788 {
2789roundoff:
2790 while (*--s == '9')
2791 if (s == s0)
2792 {
2793 k++;
2794 *s++= '1';
2795 goto ret;
2796 }
2797 ++*s++;
2798 }
2799 else
2800 {
2801#ifdef Honor_FLT_ROUNDS
2802trimzeros:
2803#endif
2804 while (*--s == '0');
2805 s++;
2806 }
2807ret:
2808 Bfree(S, &alloc);
2809 if (mhi)
2810 {
2811 if (mlo && mlo != mhi)
2812 Bfree(mlo, &alloc);
2813 Bfree(mhi, &alloc);
2814 }
2815ret1:
2816 Bfree(b, &alloc);
2817 *s= 0;
2818 *decpt= k + 1;
2819 if (rve)
2820 *rve= s;
2821 return s0;
2822}
2823