SphinxBase 5prealpha
dtoa.c
1/****************************************************************
2 *
3 * The author of this software is David M. Gay.
4 *
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
6 *
7 * Permission to use, copy, modify, and distribute this software for any
8 * purpose without fee is hereby granted, provided that this entire notice
9 * is included in all copies of any software which is or includes a copy
10 * or modification of this software and in all copies of the supporting
11 * documentation for such software.
12 *
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17 *
18 ***************************************************************/
19
20/****************************************************************
21 * This is dtoa.c by David M. Gay, downloaded from
22 * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for
23 * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith.
24 * It was taken from Python distribution then and imported into sphinxbase.
25 * Python version is preferred due to cleanups, though original
26 * version at netlib is still maintained.
27 *
28 * Please remember to check http://www.netlib.org/fp regularly for bugfixes and updates.
29 *
30 * The major modifications from Gay's original code are as follows:
31 *
32 * 0. The original code has been specialized to Sphinxbase's needs by removing
33 * many of the #ifdef'd sections. In particular, code to support VAX and
34 * IBM floating-point formats, hex NaNs, hex floats, locale-aware
35 * treatment of the decimal point, and setting of the inexact flag have
36 * been removed.
37 *
38 * 1. We use cdk_calloc and ckd_free in place of malloc and free.
39 *
40 * 2. The public functions strtod, dtoa and freedtoa all now have
41 * a sb_ prefix.
42 *
43 * 3. Instead of assuming that malloc always succeeds, we thread
44 * malloc failures through the code. The functions
45 *
46 * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b
47 *
48 * of return type *Bigint all return NULL to indicate a malloc failure.
49 * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on
50 * failure. bigcomp now has return type int (it used to be void) and
51 * returns -1 on failure and 0 otherwise. sb_dtoa returns NULL
52 * on failure. sb_strtod indicates failure due to malloc failure
53 * by returning -1.0, setting errno=ENOMEM and *se to s00.
54 *
55 * 4. The static variable dtoa_result has been removed. Callers of
56 * sb_dtoa are expected to call sb_freedtoa to free the memory allocated
57 * by sb_dtoa.
58 *
59 * 5. The code has been reformatted to better fit with C style.
60 *
61 * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory
62 * that hasn't been MALLOC'ed, private_mem should only be used when k <=
63 * Kmax.
64 *
65 * 7. sb_strtod has been modified so that it doesn't accept strings with
66 * leading whitespace.
67 *
68 * 8. Global static variables are not used due to memory access issues. Fixes
69 * usage from multiple threads.
70 *
71 ***************************************************************/
72
73/* Please send bug reports for the original dtoa.c code to David M. Gay (dmg
74 * at acm dot org, with " at " changed at "@" and " dot " changed to ".").
75 */
76
77/* On a machine with IEEE extended-precision registers, it is
78 * necessary to specify double-precision (53-bit) rounding precision
79 * before invoking strtod or dtoa. If the machine uses (the equivalent
80 * of) Intel 80x87 arithmetic, the call
81 * _control87(PC_53, MCW_PC);
82 * does this with many compilers. Whether this or another call is
83 * appropriate depends on the compiler; for this to work, it may be
84 * necessary to #include "float.h" or another system-dependent header
85 * file.
86 */
87
88/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
89 *
90 * This strtod returns a nearest machine number to the input decimal
91 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
92 * broken by the IEEE round-even rule. Otherwise ties are broken by
93 * biased rounding (add half and chop).
94 *
95 * Inspired loosely by William D. Clinger's paper "How to Read Floating
96 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
97 *
98 * Modifications:
99 *
100 * 1. We only require IEEE, IBM, or VAX double-precision
101 * arithmetic (not IEEE double-extended).
102 * 2. We get by with floating-point arithmetic in a case that
103 * Clinger missed -- when we're computing d * 10^n
104 * for a small integer d and the integer n is not too
105 * much larger than 22 (the maximum integer k for which
106 * we can represent 10^k exactly), we may be able to
107 * compute (d*10^k) * 10^(e-k) with just one roundoff.
108 * 3. Rather than a bit-at-a-time adjustment of the binary
109 * result in the hard case, we use floating-point
110 * arithmetic to determine the adjustment to within
111 * one bit; only in really hard cases do we need to
112 * compute a second residual.
113 * 4. Because of 3., we don't need a large table of powers of 10
114 * for ten-to-e (just some small tables, e.g. of 10^k
115 * for 0 <= k <= 22).
116 */
117
118/* Linking of sphinxbase's #defines to Gay's #defines starts here. */
119
120#ifdef HAVE_CONFIG_H
121#include "config.h"
122#endif
123
124#include <errno.h>
125#include <string.h>
126#include <assert.h>
127#include <stdio.h>
128
129#include <sphinxbase/ckd_alloc.h>
130#include <sphinxbase/prim_type.h>
131
132#ifdef WORDS_BIGENDIAN
133#define IEEE_MC68k
134#else
135#define IEEE_8087
136#endif
137
138#define Long int32 /* ZOMG */
139#define ULong uint32 /* WTF */
140#ifdef HAVE_LONG_LONG
141#define ULLong uint64
142#endif
143
144#define MALLOC ckd_malloc
145#define FREE ckd_free
146
147#define DBL_DIG 15
148#define DBL_MAX_10_EXP 308
149#define DBL_MAX_EXP 1024
150#define FLT_RADIX 2
151
152/* maximum permitted exponent value for strtod; exponents larger than
153 MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP
154 should fit into an int. */
155#ifndef MAX_ABS_EXP
156#define MAX_ABS_EXP 1100000000U
157#endif
158/* Bound on length of pieces of input strings in sb_strtod; specifically,
159 this is used to bound the total number of digits ignoring leading zeros and
160 the number of digits that follow the decimal point. Ideally, MAX_DIGITS
161 should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the
162 exponent clipping in sb_strtod can't affect the value of the output. */
163#ifndef MAX_DIGITS
164#define MAX_DIGITS 1000000000U
165#endif
166
167/* End sphinxbase #define linking */
168
169#ifdef DEBUG
170#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
171#endif
172
173
174#ifdef __cplusplus
175extern "C" {
176#endif
177
178typedef union { double d; ULong L[2]; } U;
179
180#ifdef IEEE_8087
181#define word0(x) (x)->L[1]
182#define word1(x) (x)->L[0]
183#else
184#define word0(x) (x)->L[0]
185#define word1(x) (x)->L[1]
186#endif
187#define dval(x) (x)->d
188
189#ifndef STRTOD_DIGLIM
190#define STRTOD_DIGLIM 40
191#endif
192
193/* maximum permitted exponent value for strtod; exponents larger than
194 MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP
195 should fit into an int. */
196#ifndef MAX_ABS_EXP
197#define MAX_ABS_EXP 1100000000U
198#endif
199/* Bound on length of pieces of input strings in sb_strtod; specifically,
200 this is used to bound the total number of digits ignoring leading zeros and
201 the number of digits that follow the decimal point. Ideally, MAX_DIGITS
202 should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the
203 exponent clipping in sb_strtod can't affect the value of the output. */
204#ifndef MAX_DIGITS
205#define MAX_DIGITS 1000000000U
206#endif
207
208/* Guard against trying to use the above values on unusual platforms with ints
209 * of width less than 32 bits. */
210#if MAX_ABS_EXP > 0x7fffffff
211#error "MAX_ABS_EXP should fit in an int"
212#endif
213#if MAX_DIGITS > 0x7fffffff
214#error "MAX_DIGITS should fit in an int"
215#endif
216
217/* The following definition of Storeinc is appropriate for MIPS processors.
218 * An alternative that might be better on some machines is
219 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
220 */
221#if defined(IEEE_8087)
222#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
223 ((unsigned short *)a)[0] = (unsigned short)c, a++)
224#else
225#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
226 ((unsigned short *)a)[1] = (unsigned short)c, a++)
227#endif
228
229/* #define P DBL_MANT_DIG */
230/* Ten_pmax = floor(P*log(2)/log(5)) */
231/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
232/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
233/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
234
235#define Exp_shift 20
236#define Exp_shift1 20
237#define Exp_msk1 0x100000
238#define Exp_msk11 0x100000
239#define Exp_mask 0x7ff00000
240#define P 53
241#define Nbits 53
242#define Bias 1023
243#define Emax 1023
244#define Emin (-1022)
245#define Etiny (-1074) /* smallest denormal is 2**Etiny */
246#define Exp_1 0x3ff00000
247#define Exp_11 0x3ff00000
248#define Ebits 11
249#define Frac_mask 0xfffff
250#define Frac_mask1 0xfffff
251#define Ten_pmax 22
252#define Bletch 0x10
253#define Bndry_mask 0xfffff
254#define Bndry_mask1 0xfffff
255#define Sign_bit 0x80000000
256#define Log2P 1
257#define Tiny0 0
258#define Tiny1 1
259#define Quick_max 14
260#define Int_max 14
261
262#ifndef Flt_Rounds
263#ifdef FLT_ROUNDS
264#define Flt_Rounds FLT_ROUNDS
265#else
266#define Flt_Rounds 1
267#endif
268#endif /*Flt_Rounds*/
269
270#define Rounding Flt_Rounds
271
272#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
273#define Big1 0xffffffff
274
275/* Standard NaN used by sb_stdnan. */
276
277#define NAN_WORD0 0x7ff80000
278#define NAN_WORD1 0
279
280/* Bits of the representation of positive infinity. */
281
282#define POSINF_WORD0 0x7ff00000
283#define POSINF_WORD1 0
284
285/* struct BCinfo is used to pass information from sb_strtod to bigcomp */
286
287typedef struct BCinfo BCinfo;
288struct
289BCinfo {
290 int e0, nd, nd0, scale;
291};
292
293#define FFFFFFFF 0xffffffffUL
294
295#define Kmax 7
296
297/* struct Bigint is used to represent arbitrary-precision integers. These
298 integers are stored in sign-magnitude format, with the magnitude stored as
299 an array of base 2**32 digits. Bigints are always normalized: if x is a
300 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
301
302 The Bigint fields are as follows:
303
304 - next is a header used by Balloc and Bfree to keep track of lists
305 of freed Bigints; it's also used for the linked list of
306 powers of 5 of the form 5**2**i used by pow5mult.
307 - k indicates which pool this Bigint was allocated from
308 - maxwds is the maximum number of words space was allocated for
309 (usually maxwds == 2**k)
310 - sign is 1 for negative Bigints, 0 for positive. The sign is unused
311 (ignored on inputs, set to 0 on outputs) in almost all operations
312 involving Bigints: a notable exception is the diff function, which
313 ignores signs on inputs but sets the sign of the output correctly.
314 - wds is the actual number of significant words
315 - x contains the vector of words (digits) for this Bigint, from least
316 significant (x[0]) to most significant (x[wds-1]).
317*/
318
319struct
320Bigint {
321 struct Bigint *next;
322 int k, maxwds, sign, wds;
323 ULong x[1];
324};
325
326typedef struct Bigint Bigint;
327
328#define SPHINXBASE_USING_MEMORY_DEBUGGER 1
329
330#ifndef SPHINXBASE_USING_MEMORY_DEBUGGER
331
332#ifndef PRIVATE_MEM
333#define PRIVATE_MEM 2304
334#endif
335#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
336static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
337
338/* Memory management: memory is allocated from, and returned to, Kmax+1 pools
339 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
340 1 << k. These pools are maintained as linked lists, with freelist[k]
341 pointing to the head of the list for pool k.
342
343 On allocation, if there's no free slot in the appropriate pool, MALLOC is
344 called to get more memory. This memory is not returned to the system until
345 Python quits. There's also a private memory pool that's allocated from
346 in preference to using MALLOC.
347
348 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
349 decimal digits), memory is directly allocated using MALLOC, and freed using
350 FREE.
351
352 XXX: it would be easy to bypass this memory-management system and
353 translate each call to Balloc into a call to PyMem_Malloc, and each
354 Bfree to PyMem_Free. Investigate whether this has any significant
355 performance on impact. */
356
357static Bigint *freelist[Kmax+1];
358
359/* Allocate space for a Bigint with up to 1<<k digits */
360
361static Bigint *
362Balloc(int k)
363{
364 int x;
365 Bigint *rv;
366 unsigned int len;
367
368 if (k <= Kmax && (rv = freelist[k]))
369 freelist[k] = rv->next;
370 else {
371 x = 1 << k;
372 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
373 /sizeof(double);
374 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
375 rv = (Bigint*)pmem_next;
376 pmem_next += len;
377 }
378 else {
379 rv = (Bigint*)MALLOC(len*sizeof(double));
380 if (rv == NULL)
381 return NULL;
382 }
383 rv->k = k;
384 rv->maxwds = x;
385 }
386 rv->sign = rv->wds = 0;
387 return rv;
388}
389
390/* Free a Bigint allocated with Balloc */
391
392static void
393Bfree(Bigint *v)
394{
395 if (v) {
396 if (v->k > Kmax)
397 FREE((void*)v);
398 else {
399 v->next = freelist[v->k];
400 freelist[v->k] = v;
401 }
402 }
403}
404
405#else
406
407/* Alternative versions of Balloc and Bfree that use PyMem_Malloc and
408 PyMem_Free directly in place of the custom memory allocation scheme above.
409 These are provided for the benefit of memory debugging tools like
410 Valgrind. */
411
412/* Allocate space for a Bigint with up to 1<<k digits */
413
414static Bigint *
415Balloc(int k)
416{
417 int x;
418 Bigint *rv;
419 unsigned int len;
420
421 x = 1 << k;
422 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
423 /sizeof(double);
424
425 rv = (Bigint*)MALLOC(len*sizeof(double));
426 if (rv == NULL)
427 return NULL;
428
429 rv->k = k;
430 rv->maxwds = x;
431 rv->sign = rv->wds = 0;
432 return rv;
433}
434
435/* Free a Bigint allocated with Balloc */
436
437static void
438Bfree(Bigint *v)
439{
440 if (v) {
441 FREE((void*)v);
442 }
443}
444
445#endif /* SPHINXBASE_USING_MEMORY_DEBUGGER */
446
447#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
448 y->wds*sizeof(Long) + 2*sizeof(int))
449
450/* Multiply a Bigint b by m and add a. Either modifies b in place and returns
451 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
452 On failure, return NULL. In this case, b will have been already freed. */
453
454static Bigint *
455multadd(Bigint *b, int m, int a) /* multiply by m and add a */
456{
457 int i, wds;
458#ifdef ULLong
459 ULong *x;
460 ULLong carry, y;
461#else
462 ULong carry, *x, y;
463 ULong xi, z;
464#endif
465 Bigint *b1;
466
467 wds = b->wds;
468 x = b->x;
469 i = 0;
470 carry = a;
471 do {
472#ifdef ULLong
473 y = *x * (ULLong)m + carry;
474 carry = y >> 32;
475 *x++ = (ULong)(y & FFFFFFFF);
476#else
477 xi = *x;
478 y = (xi & 0xffff) * m + carry;
479 z = (xi >> 16) * m + (y >> 16);
480 carry = z >> 16;
481 *x++ = (z << 16) + (y & 0xffff);
482#endif
483 }
484 while(++i < wds);
485 if (carry) {
486 if (wds >= b->maxwds) {
487 b1 = Balloc(b->k+1);
488 if (b1 == NULL){
489 Bfree(b);
490 return NULL;
491 }
492 Bcopy(b1, b);
493 Bfree(b);
494 b = b1;
495 }
496 b->x[wds++] = (ULong)carry;
497 b->wds = wds;
498 }
499 return b;
500}
501
502/* convert a string s containing nd decimal digits (possibly containing a
503 decimal separator at position nd0, which is ignored) to a Bigint. This
504 function carries on where the parsing code in sb_strtod leaves off: on
505 entry, y9 contains the result of converting the first 9 digits. Returns
506 NULL on failure. */
507
508static Bigint *
509s2b(const char *s, int nd0, int nd, ULong y9)
510{
511 Bigint *b;
512 int i, k;
513 Long x, y;
514
515 x = (nd + 8) / 9;
516 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
517 b = Balloc(k);
518 if (b == NULL)
519 return NULL;
520 b->x[0] = y9;
521 b->wds = 1;
522
523 if (nd <= 9)
524 return b;
525
526 s += 9;
527 for (i = 9; i < nd0; i++) {
528 b = multadd(b, 10, *s++ - '0');
529 if (b == NULL)
530 return NULL;
531 }
532 s++;
533 for(; i < nd; i++) {
534 b = multadd(b, 10, *s++ - '0');
535 if (b == NULL)
536 return NULL;
537 }
538 return b;
539}
540
541/* count leading 0 bits in the 32-bit integer x. */
542
543static int
544hi0bits(ULong x)
545{
546 int k = 0;
547
548 if (!(x & 0xffff0000)) {
549 k = 16;
550 x <<= 16;
551 }
552 if (!(x & 0xff000000)) {
553 k += 8;
554 x <<= 8;
555 }
556 if (!(x & 0xf0000000)) {
557 k += 4;
558 x <<= 4;
559 }
560 if (!(x & 0xc0000000)) {
561 k += 2;
562 x <<= 2;
563 }
564 if (!(x & 0x80000000)) {
565 k++;
566 if (!(x & 0x40000000))
567 return 32;
568 }
569 return k;
570}
571
572/* count trailing 0 bits in the 32-bit integer y, and shift y right by that
573 number of bits. */
574
575static int
576lo0bits(ULong *y)
577{
578 int k;
579 ULong x = *y;
580
581 if (x & 7) {
582 if (x & 1)
583 return 0;
584 if (x & 2) {
585 *y = x >> 1;
586 return 1;
587 }
588 *y = x >> 2;
589 return 2;
590 }
591 k = 0;
592 if (!(x & 0xffff)) {
593 k = 16;
594 x >>= 16;
595 }
596 if (!(x & 0xff)) {
597 k += 8;
598 x >>= 8;
599 }
600 if (!(x & 0xf)) {
601 k += 4;
602 x >>= 4;
603 }
604 if (!(x & 0x3)) {
605 k += 2;
606 x >>= 2;
607 }
608 if (!(x & 1)) {
609 k++;
610 x >>= 1;
611 if (!x)
612 return 32;
613 }
614 *y = x;
615 return k;
616}
617
618/* convert a small nonnegative integer to a Bigint */
619
620static Bigint *
621i2b(int i)
622{
623 Bigint *b;
624
625 b = Balloc(1);
626 if (b == NULL)
627 return NULL;
628 b->x[0] = i;
629 b->wds = 1;
630 return b;
631}
632
633/* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
634 the signs of a and b. */
635
636static Bigint *
637mult(Bigint *a, Bigint *b)
638{
639 Bigint *c;
640 int k, wa, wb, wc;
641 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
642 ULong y;
643#ifdef ULLong
644 ULLong carry, z;
645#else
646 ULong carry, z;
647 ULong z2;
648#endif
649
650 if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) {
651 c = Balloc(0);
652 if (c == NULL)
653 return NULL;
654 c->wds = 1;
655 c->x[0] = 0;
656 return c;
657 }
658
659 if (a->wds < b->wds) {
660 c = a;
661 a = b;
662 b = c;
663 }
664 k = a->k;
665 wa = a->wds;
666 wb = b->wds;
667 wc = wa + wb;
668 if (wc > a->maxwds)
669 k++;
670 c = Balloc(k);
671 if (c == NULL)
672 return NULL;
673 for(x = c->x, xa = x + wc; x < xa; x++)
674 *x = 0;
675 xa = a->x;
676 xae = xa + wa;
677 xb = b->x;
678 xbe = xb + wb;
679 xc0 = c->x;
680#ifdef ULLong
681 for(; xb < xbe; xc0++) {
682 if ((y = *xb++)) {
683 x = xa;
684 xc = xc0;
685 carry = 0;
686 do {
687 z = *x++ * (ULLong)y + *xc + carry;
688 carry = z >> 32;
689 *xc++ = (ULong)(z & FFFFFFFF);
690 }
691 while(x < xae);
692 *xc = (ULong)carry;
693 }
694 }
695#else
696 for(; xb < xbe; xb++, xc0++) {
697 if (y = *xb & 0xffff) {
698 x = xa;
699 xc = xc0;
700 carry = 0;
701 do {
702 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
703 carry = z >> 16;
704 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
705 carry = z2 >> 16;
706 Storeinc(xc, z2, z);
707 }
708 while(x < xae);
709 *xc = carry;
710 }
711 if (y = *xb >> 16) {
712 x = xa;
713 xc = xc0;
714 carry = 0;
715 z2 = *xc;
716 do {
717 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
718 carry = z >> 16;
719 Storeinc(xc, z, z2);
720 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
721 carry = z2 >> 16;
722 }
723 while(x < xae);
724 *xc = z2;
725 }
726 }
727#endif
728 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
729 c->wds = wc;
730 return c;
731}
732
733#ifndef SPHINXBASE_USING_MEMORY_DEBUGGER
734
735/* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
736
737static Bigint *p5s;
738
739/* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
740 failure; if the returned pointer is distinct from b then the original
741 Bigint b will have been Bfree'd. Ignores the sign of b. */
742
743static Bigint *
744pow5mult(Bigint *b, int k)
745{
746 Bigint *b1, *p5, *p51;
747 int i;
748 static int p05[3] = { 5, 25, 125 };
749
750 if ((i = k & 3)) {
751 b = multadd(b, p05[i-1], 0);
752 if (b == NULL)
753 return NULL;
754 }
755
756 if (!(k >>= 2))
757 return b;
758 p5 = p5s;
759 if (!p5) {
760 /* first time */
761 p5 = i2b(625);
762 if (p5 == NULL) {
763 Bfree(b);
764 return NULL;
765 }
766 p5s = p5;
767 p5->next = 0;
768 }
769 for(;;) {
770 if (k & 1) {
771 b1 = mult(b, p5);
772 Bfree(b);
773 b = b1;
774 if (b == NULL)
775 return NULL;
776 }
777 if (!(k >>= 1))
778 break;
779 p51 = p5->next;
780 if (!p51) {
781 p51 = mult(p5,p5);
782 if (p51 == NULL) {
783 Bfree(b);
784 return NULL;
785 }
786 p51->next = 0;
787 p5->next = p51;
788 }
789 p5 = p51;
790 }
791 return b;
792}
793
794#else
795
796/* Version of pow5mult that doesn't cache powers of 5. Provided for
797 the benefit of memory debugging tools like Valgrind. */
798
799static Bigint *
800pow5mult(Bigint *b, int k)
801{
802 Bigint *b1, *p5, *p51;
803 int i;
804 static int p05[3] = { 5, 25, 125 };
805
806 if ((i = k & 3)) {
807 b = multadd(b, p05[i-1], 0);
808 if (b == NULL)
809 return NULL;
810 }
811
812 if (!(k >>= 2))
813 return b;
814 p5 = i2b(625);
815 if (p5 == NULL) {
816 Bfree(b);
817 return NULL;
818 }
819
820 for(;;) {
821 if (k & 1) {
822 b1 = mult(b, p5);
823 Bfree(b);
824 b = b1;
825 if (b == NULL) {
826 Bfree(p5);
827 return NULL;
828 }
829 }
830 if (!(k >>= 1))
831 break;
832 p51 = mult(p5, p5);
833 Bfree(p5);
834 p5 = p51;
835 if (p5 == NULL) {
836 Bfree(b);
837 return NULL;
838 }
839 }
840 Bfree(p5);
841 return b;
842}
843
844#endif /* SPHINXBASE_USING_MEMORY_DEBUGGER */
845
846/* shift a Bigint b left by k bits. Return a pointer to the shifted result,
847 or NULL on failure. If the returned pointer is distinct from b then the
848 original b will have been Bfree'd. Ignores the sign of b. */
849
850static Bigint *
851lshift(Bigint *b, int k)
852{
853 int i, k1, n, n1;
854 Bigint *b1;
855 ULong *x, *x1, *xe, z;
856
857 if (!k || (!b->x[0] && b->wds == 1))
858 return b;
859
860 n = k >> 5;
861 k1 = b->k;
862 n1 = n + b->wds + 1;
863 for(i = b->maxwds; n1 > i; i <<= 1)
864 k1++;
865 b1 = Balloc(k1);
866 if (b1 == NULL) {
867 Bfree(b);
868 return NULL;
869 }
870 x1 = b1->x;
871 for(i = 0; i < n; i++)
872 *x1++ = 0;
873 x = b->x;
874 xe = x + b->wds;
875 if (k &= 0x1f) {
876 k1 = 32 - k;
877 z = 0;
878 do {
879 *x1++ = *x << k | z;
880 z = *x++ >> k1;
881 }
882 while(x < xe);
883 if ((*x1 = z))
884 ++n1;
885 }
886 else do
887 *x1++ = *x++;
888 while(x < xe);
889 b1->wds = n1 - 1;
890 Bfree(b);
891 return b1;
892}
893
894/* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
895 1 if a > b. Ignores signs of a and b. */
896
897static int
898cmp(Bigint *a, Bigint *b)
899{
900 ULong *xa, *xa0, *xb, *xb0;
901 int i, j;
902
903 i = a->wds;
904 j = b->wds;
905#ifdef DEBUG
906 if (i > 1 && !a->x[i-1])
907 Bug("cmp called with a->x[a->wds-1] == 0");
908 if (j > 1 && !b->x[j-1])
909 Bug("cmp called with b->x[b->wds-1] == 0");
910#endif
911 if (i -= j)
912 return i;
913 xa0 = a->x;
914 xa = xa0 + j;
915 xb0 = b->x;
916 xb = xb0 + j;
917 for(;;) {
918 if (*--xa != *--xb)
919 return *xa < *xb ? -1 : 1;
920 if (xa <= xa0)
921 break;
922 }
923 return 0;
924}
925
926/* Take the difference of Bigints a and b, returning a new Bigint. Returns
927 NULL on failure. The signs of a and b are ignored, but the sign of the
928 result is set appropriately. */
929
930static Bigint *
931diff(Bigint *a, Bigint *b)
932{
933 Bigint *c;
934 int i, wa, wb;
935 ULong *xa, *xae, *xb, *xbe, *xc;
936#ifdef ULLong
937 ULLong borrow, y;
938#else
939 ULong borrow, y;
940 ULong z;
941#endif
942
943 i = cmp(a,b);
944 if (!i) {
945 c = Balloc(0);
946 if (c == NULL)
947 return NULL;
948 c->wds = 1;
949 c->x[0] = 0;
950 return c;
951 }
952 if (i < 0) {
953 c = a;
954 a = b;
955 b = c;
956 i = 1;
957 }
958 else
959 i = 0;
960 c = Balloc(a->k);
961 if (c == NULL)
962 return NULL;
963 c->sign = i;
964 wa = a->wds;
965 xa = a->x;
966 xae = xa + wa;
967 wb = b->wds;
968 xb = b->x;
969 xbe = xb + wb;
970 xc = c->x;
971 borrow = 0;
972#ifdef ULLong
973 do {
974 y = (ULLong)*xa++ - *xb++ - borrow;
975 borrow = y >> 32 & (ULong)1;
976 *xc++ = (ULong)(y & FFFFFFFF);
977 }
978 while(xb < xbe);
979 while(xa < xae) {
980 y = *xa++ - borrow;
981 borrow = y >> 32 & (ULong)1;
982 *xc++ = (ULong)(y & FFFFFFFF);
983 }
984#else
985 do {
986 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
987 borrow = (y & 0x10000) >> 16;
988 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
989 borrow = (z & 0x10000) >> 16;
990 Storeinc(xc, z, y);
991 }
992 while(xb < xbe);
993 while(xa < xae) {
994 y = (*xa & 0xffff) - borrow;
995 borrow = (y & 0x10000) >> 16;
996 z = (*xa++ >> 16) - borrow;
997 borrow = (z & 0x10000) >> 16;
998 Storeinc(xc, z, y);
999 }
1000#endif
1001 while(!*--xc)
1002 wa--;
1003 c->wds = wa;
1004 return c;
1005}
1006
1007/* Given a positive normal double x, return the difference between x and the
1008 next double up. Doesn't give correct results for subnormals. */
1009
1010static double
1011ulp(U *x)
1012{
1013 Long L;
1014 U u;
1015
1016 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1017 word0(&u) = L;
1018 word1(&u) = 0;
1019 return dval(&u);
1020}
1021
1022/* Convert a Bigint to a double plus an exponent */
1023
1024static double
1025b2d(Bigint *a, int *e)
1026{
1027 ULong *xa, *xa0, w, y, z;
1028 int k;
1029 U d;
1030
1031 xa0 = a->x;
1032 xa = xa0 + a->wds;
1033 y = *--xa;
1034#ifdef DEBUG
1035 if (!y) Bug("zero y in b2d");
1036#endif
1037 k = hi0bits(y);
1038 *e = 32 - k;
1039 if (k < Ebits) {
1040 word0(&d) = Exp_1 | y >> (Ebits - k);
1041 w = xa > xa0 ? *--xa : 0;
1042 word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
1043 goto ret_d;
1044 }
1045 z = xa > xa0 ? *--xa : 0;
1046 if (k -= Ebits) {
1047 word0(&d) = Exp_1 | y << k | z >> (32 - k);
1048 y = xa > xa0 ? *--xa : 0;
1049 word1(&d) = z << k | y >> (32 - k);
1050 }
1051 else {
1052 word0(&d) = Exp_1 | y;
1053 word1(&d) = z;
1054 }
1055 ret_d:
1056 return dval(&d);
1057}
1058
1059/* Convert a scaled double to a Bigint plus an exponent. Similar to d2b,
1060 except that it accepts the scale parameter used in sb_strtod (which
1061 should be either 0 or 2*P), and the normalization for the return value is
1062 different (see below). On input, d should be finite and nonnegative, and d
1063 / 2**scale should be exactly representable as an IEEE 754 double.
1064
1065 Returns a Bigint b and an integer e such that
1066
1067 dval(d) / 2**scale = b * 2**e.
1068
1069 Unlike d2b, b is not necessarily odd: b and e are normalized so
1070 that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P
1071 and e == Etiny. This applies equally to an input of 0.0: in that
1072 case the return values are b = 0 and e = Etiny.
1073
1074 The above normalization ensures that for all possible inputs d,
1075 2**e gives ulp(d/2**scale).
1076
1077 Returns NULL on failure.
1078*/
1079
1080static Bigint *
1081sd2b(U *d, int scale, int *e)
1082{
1083 Bigint *b;
1084
1085 b = Balloc(1);
1086 if (b == NULL)
1087 return NULL;
1088
1089 /* First construct b and e assuming that scale == 0. */
1090 b->wds = 2;
1091 b->x[0] = word1(d);
1092 b->x[1] = word0(d) & Frac_mask;
1093 *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift);
1094 if (*e < Etiny)
1095 *e = Etiny;
1096 else
1097 b->x[1] |= Exp_msk1;
1098
1099 /* Now adjust for scale, provided that b != 0. */
1100 if (scale && (b->x[0] || b->x[1])) {
1101 *e -= scale;
1102 if (*e < Etiny) {
1103 scale = Etiny - *e;
1104 *e = Etiny;
1105 /* We can't shift more than P-1 bits without shifting out a 1. */
1106 assert(0 < scale && scale <= P - 1);
1107 if (scale >= 32) {
1108 /* The bits shifted out should all be zero. */
1109 assert(b->x[0] == 0);
1110 b->x[0] = b->x[1];
1111 b->x[1] = 0;
1112 scale -= 32;
1113 }
1114 if (scale) {
1115 /* The bits shifted out should all be zero. */
1116 assert(b->x[0] << (32 - scale) == 0);
1117 b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale));
1118 b->x[1] >>= scale;
1119 }
1120 }
1121 }
1122 /* Ensure b is normalized. */
1123 if (!b->x[1])
1124 b->wds = 1;
1125
1126 return b;
1127}
1128
1129/* Convert a double to a Bigint plus an exponent. Return NULL on failure.
1130
1131 Given a finite nonzero double d, return an odd Bigint b and exponent *e
1132 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
1133 significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
1134
1135 If d is zero, then b == 0, *e == -1010, *bbits = 0.
1136 */
1137
1138static Bigint *
1139d2b(U *d, int *e, int *bits)
1140{
1141 Bigint *b;
1142 int de, k;
1143 ULong *x, y, z;
1144 int i;
1145
1146 b = Balloc(1);
1147 if (b == NULL)
1148 return NULL;
1149 x = b->x;
1150
1151 z = word0(d) & Frac_mask;
1152 word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */
1153 if ((de = (int)(word0(d) >> Exp_shift)))
1154 z |= Exp_msk1;
1155 if ((y = word1(d))) {
1156 if ((k = lo0bits(&y))) {
1157 x[0] = y | z << (32 - k);
1158 z >>= k;
1159 }
1160 else
1161 x[0] = y;
1162 i =
1163 b->wds = (x[1] = z) ? 2 : 1;
1164 }
1165 else {
1166 k = lo0bits(&z);
1167 x[0] = z;
1168 i =
1169 b->wds = 1;
1170 k += 32;
1171 }
1172 if (de) {
1173 *e = de - Bias - (P-1) + k;
1174 *bits = P - k;
1175 }
1176 else {
1177 *e = de - Bias - (P-1) + 1 + k;
1178 *bits = 32*i - hi0bits(x[i-1]);
1179 }
1180 return b;
1181}
1182
1183/* Compute the ratio of two Bigints, as a double. The result may have an
1184 error of up to 2.5 ulps. */
1185
1186static double
1187ratio(Bigint *a, Bigint *b)
1188{
1189 U da, db;
1190 int k, ka, kb;
1191
1192 dval(&da) = b2d(a, &ka);
1193 dval(&db) = b2d(b, &kb);
1194 k = ka - kb + 32*(a->wds - b->wds);
1195 if (k > 0)
1196 word0(&da) += k*Exp_msk1;
1197 else {
1198 k = -k;
1199 word0(&db) += k*Exp_msk1;
1200 }
1201 return dval(&da) / dval(&db);
1202}
1203
1204static const double
1205tens[] = {
1206 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1207 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1208 1e20, 1e21, 1e22
1209};
1210
1211static const double
1212bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1213static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1214 9007199254740992.*9007199254740992.e-256
1215 /* = 2^106 * 1e-256 */
1216};
1217/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1218/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1219#define Scale_Bit 0x10
1220#define n_bigtens 5
1221
1222#define ULbits 32
1223#define kshift 5
1224#define kmask 31
1225
1226
1227static int
1228dshift(Bigint *b, int p2)
1229{
1230 int rv = hi0bits(b->x[b->wds-1]) - 4;
1231 if (p2 > 0)
1232 rv -= p2;
1233 return rv & kmask;
1234}
1235
1236/* special case of Bigint division. The quotient is always in the range 0 <=
1237 quotient < 10, and on entry the divisor S is normalized so that its top 4
1238 bits (28--31) are zero and bit 27 is set. */
1239
1240static int
1241quorem(Bigint *b, Bigint *S)
1242{
1243 int n;
1244 ULong *bx, *bxe, q, *sx, *sxe;
1245#ifdef ULLong
1246 ULLong borrow, carry, y, ys;
1247#else
1248 ULong borrow, carry, y, ys;
1249 ULong si, z, zs;
1250#endif
1251
1252 n = S->wds;
1253#ifdef DEBUG
1254 /*debug*/ if (b->wds > n)
1255 /*debug*/ Bug("oversize b in quorem");
1256#endif
1257 if (b->wds < n)
1258 return 0;
1259 sx = S->x;
1260 sxe = sx + --n;
1261 bx = b->x;
1262 bxe = bx + n;
1263 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1264#ifdef DEBUG
1265 /*debug*/ if (q > 9)
1266 /*debug*/ Bug("oversized quotient in quorem");
1267#endif
1268 if (q) {
1269 borrow = 0;
1270 carry = 0;
1271 do {
1272#ifdef ULLong
1273 ys = *sx++ * (ULLong)q + carry;
1274 carry = ys >> 32;
1275 y = *bx - (ys & FFFFFFFF) - borrow;
1276 borrow = y >> 32 & (ULong)1;
1277 *bx++ = (ULong)(y & FFFFFFFF);
1278#else
1279 si = *sx++;
1280 ys = (si & 0xffff) * q + carry;
1281 zs = (si >> 16) * q + (ys >> 16);
1282 carry = zs >> 16;
1283 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1284 borrow = (y & 0x10000) >> 16;
1285 z = (*bx >> 16) - (zs & 0xffff) - borrow;
1286 borrow = (z & 0x10000) >> 16;
1287 Storeinc(bx, z, y);
1288#endif
1289 }
1290 while(sx <= sxe);
1291 if (!*bxe) {
1292 bx = b->x;
1293 while(--bxe > bx && !*bxe)
1294 --n;
1295 b->wds = n;
1296 }
1297 }
1298 if (cmp(b, S) >= 0) {
1299 q++;
1300 borrow = 0;
1301 carry = 0;
1302 bx = b->x;
1303 sx = S->x;
1304 do {
1305#ifdef ULLong
1306 ys = *sx++ + carry;
1307 carry = ys >> 32;
1308 y = *bx - (ys & FFFFFFFF) - borrow;
1309 borrow = y >> 32 & (ULong)1;
1310 *bx++ = (ULong)(y & FFFFFFFF);
1311#else
1312 si = *sx++;
1313 ys = (si & 0xffff) + carry;
1314 zs = (si >> 16) + (ys >> 16);
1315 carry = zs >> 16;
1316 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1317 borrow = (y & 0x10000) >> 16;
1318 z = (*bx >> 16) - (zs & 0xffff) - borrow;
1319 borrow = (z & 0x10000) >> 16;
1320 Storeinc(bx, z, y);
1321#endif
1322 }
1323 while(sx <= sxe);
1324 bx = b->x;
1325 bxe = bx + n;
1326 if (!*bxe) {
1327 while(--bxe > bx && !*bxe)
1328 --n;
1329 b->wds = n;
1330 }
1331 }
1332 return q;
1333}
1334
1335/* sulp(x) is a version of ulp(x) that takes bc.scale into account.
1336
1337 Assuming that x is finite and nonnegative (positive zero is fine
1338 here) and x / 2^bc.scale is exactly representable as a double,
1339 sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
1340
1341static double
1342sulp(U *x, BCinfo *bc)
1343{
1344 U u;
1345
1346 if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) {
1347 /* rv/2^bc->scale is subnormal */
1348 word0(&u) = (P+2)*Exp_msk1;
1349 word1(&u) = 0;
1350 return u.d;
1351 }
1352 else {
1353 assert(word0(x) || word1(x)); /* x != 0.0 */
1354 return ulp(x);
1355 }
1356}
1357
1358/* The bigcomp function handles some hard cases for strtod, for inputs
1359 with more than STRTOD_DIGLIM digits. It's called once an initial
1360 estimate for the double corresponding to the input string has
1361 already been obtained by the code in sb_strtod.
1362
1363 The bigcomp function is only called after sb_strtod has found a
1364 double value rv such that either rv or rv + 1ulp represents the
1365 correctly rounded value corresponding to the original string. It
1366 determines which of these two values is the correct one by
1367 computing the decimal digits of rv + 0.5ulp and comparing them with
1368 the corresponding digits of s0.
1369
1370 In the following, write dv for the absolute value of the number represented
1371 by the input string.
1372
1373 Inputs:
1374
1375 s0 points to the first significant digit of the input string.
1376
1377 rv is a (possibly scaled) estimate for the closest double value to the
1378 value represented by the original input to sb_strtod. If
1379 bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to
1380 the input value.
1381
1382 bc is a struct containing information gathered during the parsing and
1383 estimation steps of sb_strtod. Description of fields follows:
1384
1385 bc->e0 gives the exponent of the input value, such that dv = (integer
1386 given by the bd->nd digits of s0) * 10**e0
1387
1388 bc->nd gives the total number of significant digits of s0. It will
1389 be at least 1.
1390
1391 bc->nd0 gives the number of significant digits of s0 before the
1392 decimal separator. If there's no decimal separator, bc->nd0 ==
1393 bc->nd.
1394
1395 bc->scale is the value used to scale rv to avoid doing arithmetic with
1396 subnormal values. It's either 0 or 2*P (=106).
1397
1398 Outputs:
1399
1400 On successful exit, rv/2^(bc->scale) is the closest double to dv.
1401
1402 Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */
1403
1404static int
1405bigcomp(U *rv, const char *s0, BCinfo *bc)
1406{
1407 Bigint *b, *d;
1408 int b2, d2, dd, i, nd, nd0, odd, p2, p5;
1409
1410 nd = bc->nd;
1411 nd0 = bc->nd0;
1412 p5 = nd + bc->e0;
1413 b = sd2b(rv, bc->scale, &p2);
1414 if (b == NULL)
1415 return -1;
1416
1417 /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
1418 case, this is used for round to even. */
1419 odd = b->x[0] & 1;
1420
1421 /* left shift b by 1 bit and or a 1 into the least significant bit;
1422 this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */
1423 b = lshift(b, 1);
1424 if (b == NULL)
1425 return -1;
1426 b->x[0] |= 1;
1427 p2--;
1428
1429 p2 -= p5;
1430 d = i2b(1);
1431 if (d == NULL) {
1432 Bfree(b);
1433 return -1;
1434 }
1435 /* Arrange for convenient computation of quotients:
1436 * shift left if necessary so divisor has 4 leading 0 bits.
1437 */
1438 if (p5 > 0) {
1439 d = pow5mult(d, p5);
1440 if (d == NULL) {
1441 Bfree(b);
1442 return -1;
1443 }
1444 }
1445 else if (p5 < 0) {
1446 b = pow5mult(b, -p5);
1447 if (b == NULL) {
1448 Bfree(d);
1449 return -1;
1450 }
1451 }
1452 if (p2 > 0) {
1453 b2 = p2;
1454 d2 = 0;
1455 }
1456 else {
1457 b2 = 0;
1458 d2 = -p2;
1459 }
1460 i = dshift(d, d2);
1461 if ((b2 += i) > 0) {
1462 b = lshift(b, b2);
1463 if (b == NULL) {
1464 Bfree(d);
1465 return -1;
1466 }
1467 }
1468 if ((d2 += i) > 0) {
1469 d = lshift(d, d2);
1470 if (d == NULL) {
1471 Bfree(b);
1472 return -1;
1473 }
1474 }
1475
1476 /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==
1477 * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing
1478 * a number in the range [0.1, 1). */
1479 if (cmp(b, d) >= 0)
1480 /* b/d >= 1 */
1481 dd = -1;
1482 else {
1483 i = 0;
1484 for(;;) {
1485 b = multadd(b, 10, 0);
1486 if (b == NULL) {
1487 Bfree(d);
1488 return -1;
1489 }
1490 dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d);
1491 i++;
1492
1493 if (dd)
1494 break;
1495 if (!b->x[0] && b->wds == 1) {
1496 /* b/d == 0 */
1497 dd = i < nd;
1498 break;
1499 }
1500 if (!(i < nd)) {
1501 /* b/d != 0, but digits of s0 exhausted */
1502 dd = -1;
1503 break;
1504 }
1505 }
1506 }
1507 Bfree(b);
1508 Bfree(d);
1509 if (dd > 0 || (dd == 0 && odd))
1510 dval(rv) += sulp(rv, bc);
1511 return 0;
1512}
1513
1514/* Return a 'standard' NaN value.
1515
1516 There are exactly two quiet NaNs that don't arise by 'quieting' signaling
1517 NaNs (see IEEE 754-2008, section 6.2.1). If sign == 0, return the one whose
1518 sign bit is cleared. Otherwise, return the one whose sign bit is set.
1519*/
1520
1521double
1522sb_stdnan(int sign)
1523{
1524 U rv;
1525 word0(&rv) = NAN_WORD0;
1526 word1(&rv) = NAN_WORD1;
1527 if (sign)
1528 word0(&rv) |= Sign_bit;
1529 return dval(&rv);
1530}
1531
1532/* Return positive or negative infinity, according to the given sign (0 for
1533 * positive infinity, 1 for negative infinity). */
1534
1535double
1536sb_infinity(int sign)
1537{
1538 U rv;
1539 word0(&rv) = POSINF_WORD0;
1540 word1(&rv) = POSINF_WORD1;
1541 return sign ? -dval(&rv) : dval(&rv);
1542}
1543
1544double
1545sb_strtod(const char *s00, char **se)
1546{
1547 int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error;
1548 int esign, i, j, k, lz, nd, nd0, odd, sign;
1549 const char *s, *s0, *s1;
1550 double aadj, aadj1;
1551 U aadj2, adj, rv, rv0;
1552 ULong y, z, abs_exp;
1553 Long L;
1554 BCinfo bc;
1555 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1556 size_t ndigits, fraclen;
1557
1558 dval(&rv) = 0.;
1559
1560 /* Start parsing. */
1561 c = *(s = s00);
1562
1563 /* Parse optional sign, if present. */
1564 sign = 0;
1565 switch (c) {
1566 case '-':
1567 sign = 1;
1568 /* no break */
1569 case '+':
1570 c = *++s;
1571 }
1572
1573 /* Skip leading zeros: lz is true iff there were leading zeros. */
1574 s1 = s;
1575 while (c == '0')
1576 c = *++s;
1577 lz = s != s1;
1578
1579 /* Point s0 at the first nonzero digit (if any). fraclen will be the
1580 number of digits between the decimal point and the end of the
1581 digit string. ndigits will be the total number of digits ignoring
1582 leading zeros. */
1583 s0 = s1 = s;
1584 while ('0' <= c && c <= '9')
1585 c = *++s;
1586 ndigits = s - s1;
1587 fraclen = 0;
1588
1589 /* Parse decimal point and following digits. */
1590 if (c == '.') {
1591 c = *++s;
1592 if (!ndigits) {
1593 s1 = s;
1594 while (c == '0')
1595 c = *++s;
1596 lz = lz || s != s1;
1597 fraclen += (s - s1);
1598 s0 = s;
1599 }
1600 s1 = s;
1601 while ('0' <= c && c <= '9')
1602 c = *++s;
1603 ndigits += s - s1;
1604 fraclen += s - s1;
1605 }
1606
1607 /* Now lz is true if and only if there were leading zero digits, and
1608 ndigits gives the total number of digits ignoring leading zeros. A
1609 valid input must have at least one digit. */
1610 if (!ndigits && !lz) {
1611 if (se)
1612 *se = (char *)s00;
1613 goto parse_error;
1614 }
1615
1616 /* Range check ndigits and fraclen to make sure that they, and values
1617 computed with them, can safely fit in an int. */
1618 if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) {
1619 if (se)
1620 *se = (char *)s00;
1621 goto parse_error;
1622 }
1623 nd = (int)ndigits;
1624 nd0 = (int)ndigits - (int)fraclen;
1625
1626 /* Parse exponent. */
1627 e = 0;
1628 if (c == 'e' || c == 'E') {
1629 s00 = s;
1630 c = *++s;
1631
1632 /* Exponent sign. */
1633 esign = 0;
1634 switch (c) {
1635 case '-':
1636 esign = 1;
1637 /* no break */
1638 case '+':
1639 c = *++s;
1640 }
1641
1642 /* Skip zeros. lz is true iff there are leading zeros. */
1643 s1 = s;
1644 while (c == '0')
1645 c = *++s;
1646 lz = s != s1;
1647
1648 /* Get absolute value of the exponent. */
1649 s1 = s;
1650 abs_exp = 0;
1651 while ('0' <= c && c <= '9') {
1652 abs_exp = 10*abs_exp + (c - '0');
1653 c = *++s;
1654 }
1655
1656 /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
1657 there are at most 9 significant exponent digits then overflow is
1658 impossible. */
1659 if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)
1660 e = (int)MAX_ABS_EXP;
1661 else
1662 e = (int)abs_exp;
1663 if (esign)
1664 e = -e;
1665
1666 /* A valid exponent must have at least one digit. */
1667 if (s == s1 && !lz)
1668 s = s00;
1669 }
1670
1671 /* Adjust exponent to take into account position of the point. */
1672 e -= nd - nd0;
1673 if (nd0 <= 0)
1674 nd0 = nd;
1675
1676 /* Finished parsing. Set se to indicate how far we parsed */
1677 if (se)
1678 *se = (char *)s;
1679
1680 /* If all digits were zero, exit with return value +-0.0. Otherwise,
1681 strip trailing zeros: scan back until we hit a nonzero digit. */
1682 if (!nd)
1683 goto ret;
1684 for (i = nd; i > 0; ) {
1685 --i;
1686 if (s0[i < nd0 ? i : i+1] != '0') {
1687 ++i;
1688 break;
1689 }
1690 }
1691 e += nd - i;
1692 nd = i;
1693 if (nd0 > nd)
1694 nd0 = nd;
1695
1696 /* Summary of parsing results. After parsing, and dealing with zero
1697 * inputs, we have values s0, nd0, nd, e, sign, where:
1698 *
1699 * - s0 points to the first significant digit of the input string
1700 *
1701 * - nd is the total number of significant digits (here, and
1702 * below, 'significant digits' means the set of digits of the
1703 * significand of the input that remain after ignoring leading
1704 * and trailing zeros).
1705 *
1706 * - nd0 indicates the position of the decimal point, if present; it
1707 * satisfies 1 <= nd0 <= nd. The nd significant digits are in
1708 * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
1709 * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
1710 * nd0 == nd, then s0[nd0] could be any non-digit character.)
1711 *
1712 * - e is the adjusted exponent: the absolute value of the number
1713 * represented by the original input string is n * 10**e, where
1714 * n is the integer represented by the concatenation of
1715 * s0[0:nd0] and s0[nd0+1:nd+1]
1716 *
1717 * - sign gives the sign of the input: 1 for negative, 0 for positive
1718 *
1719 * - the first and last significant digits are nonzero
1720 */
1721
1722 /* put first DBL_DIG+1 digits into integer y and z.
1723 *
1724 * - y contains the value represented by the first min(9, nd)
1725 * significant digits
1726 *
1727 * - if nd > 9, z contains the value represented by significant digits
1728 * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
1729 * gives the value represented by the first min(16, nd) sig. digits.
1730 */
1731
1732 bc.e0 = e1 = e;
1733 y = z = 0;
1734 for (i = 0; i < nd; i++) {
1735 if (i < 9)
1736 y = 10*y + s0[i < nd0 ? i : i+1] - '0';
1737 else if (i < DBL_DIG+1)
1738 z = 10*z + s0[i < nd0 ? i : i+1] - '0';
1739 else
1740 break;
1741 }
1742
1743 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1744 dval(&rv) = y;
1745 if (k > 9) {
1746 dval(&rv) = tens[k - 9] * dval(&rv) + z;
1747 }
1748 bd0 = 0;
1749 if (nd <= DBL_DIG
1750 && Flt_Rounds == 1
1751 ) {
1752 if (!e)
1753 goto ret;
1754 if (e > 0) {
1755 if (e <= Ten_pmax) {
1756 dval(&rv) *= tens[e];
1757 goto ret;
1758 }
1759 i = DBL_DIG - nd;
1760 if (e <= Ten_pmax + i) {
1761 /* A fancier test would sometimes let us do
1762 * this for larger i values.
1763 */
1764 e -= i;
1765 dval(&rv) *= tens[i];
1766 dval(&rv) *= tens[e];
1767 goto ret;
1768 }
1769 }
1770 else if (e >= -Ten_pmax) {
1771 dval(&rv) /= tens[-e];
1772 goto ret;
1773 }
1774 }
1775 e1 += nd - k;
1776
1777 bc.scale = 0;
1778
1779 /* Get starting approximation = rv * 10**e1 */
1780
1781 if (e1 > 0) {
1782 if ((i = e1 & 15))
1783 dval(&rv) *= tens[i];
1784 if (e1 &= ~15) {
1785 if (e1 > DBL_MAX_10_EXP)
1786 goto ovfl;
1787 e1 >>= 4;
1788 for(j = 0; e1 > 1; j++, e1 >>= 1)
1789 if (e1 & 1)
1790 dval(&rv) *= bigtens[j];
1791 /* The last multiplication could overflow. */
1792 word0(&rv) -= P*Exp_msk1;
1793 dval(&rv) *= bigtens[j];
1794 if ((z = word0(&rv) & Exp_mask)
1795 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1796 goto ovfl;
1797 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1798 /* set to largest number */
1799 /* (Can't trust DBL_MAX) */
1800 word0(&rv) = Big0;
1801 word1(&rv) = Big1;
1802 }
1803 else
1804 word0(&rv) += P*Exp_msk1;
1805 }
1806 }
1807 else if (e1 < 0) {
1808 /* The input decimal value lies in [10**e1, 10**(e1+16)).
1809
1810 If e1 <= -512, underflow immediately.
1811 If e1 <= -256, set bc.scale to 2*P.
1812
1813 So for input value < 1e-256, bc.scale is always set;
1814 for input value >= 1e-240, bc.scale is never set.
1815 For input values in [1e-256, 1e-240), bc.scale may or may
1816 not be set. */
1817
1818 e1 = -e1;
1819 if ((i = e1 & 15))
1820 dval(&rv) /= tens[i];
1821 if (e1 >>= 4) {
1822 if (e1 >= 1 << n_bigtens)
1823 goto undfl;
1824 if (e1 & Scale_Bit)
1825 bc.scale = 2*P;
1826 for(j = 0; e1 > 0; j++, e1 >>= 1)
1827 if (e1 & 1)
1828 dval(&rv) *= tinytens[j];
1829 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
1830 >> Exp_shift)) > 0) {
1831 /* scaled rv is denormal; clear j low bits */
1832 if (j >= 32) {
1833 word1(&rv) = 0;
1834 if (j >= 53)
1835 word0(&rv) = (P+2)*Exp_msk1;
1836 else
1837 word0(&rv) &= 0xffffffff << (j-32);
1838 }
1839 else
1840 word1(&rv) &= 0xffffffff << j;
1841 }
1842 if (!dval(&rv))
1843 goto undfl;
1844 }
1845 }
1846
1847 /* Now the hard part -- adjusting rv to the correct value.*/
1848
1849 /* Put digits into bd: true value = bd * 10^e */
1850
1851 bc.nd = nd;
1852 bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */
1853 /* to silence an erroneous warning about bc.nd0 */
1854 /* possibly not being initialized. */
1855 if (nd > STRTOD_DIGLIM) {
1856 /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
1857 /* minimum number of decimal digits to distinguish double values */
1858 /* in IEEE arithmetic. */
1859
1860 /* Truncate input to 18 significant digits, then discard any trailing
1861 zeros on the result by updating nd, nd0, e and y suitably. (There's
1862 no need to update z; it's not reused beyond this point.) */
1863 for (i = 18; i > 0; ) {
1864 /* scan back until we hit a nonzero digit. significant digit 'i'
1865 is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
1866 --i;
1867 if (s0[i < nd0 ? i : i+1] != '0') {
1868 ++i;
1869 break;
1870 }
1871 }
1872 e += nd - i;
1873 nd = i;
1874 if (nd0 > nd)
1875 nd0 = nd;
1876 if (nd < 9) { /* must recompute y */
1877 y = 0;
1878 for(i = 0; i < nd0; ++i)
1879 y = 10*y + s0[i] - '0';
1880 for(; i < nd; ++i)
1881 y = 10*y + s0[i+1] - '0';
1882 }
1883 }
1884 bd0 = s2b(s0, nd0, nd, y);
1885 if (bd0 == NULL)
1886 goto failed_malloc;
1887
1888 /* Notation for the comments below. Write:
1889
1890 - dv for the absolute value of the number represented by the original
1891 decimal input string.
1892
1893 - if we've truncated dv, write tdv for the truncated value.
1894 Otherwise, set tdv == dv.
1895
1896 - srv for the quantity rv/2^bc.scale; so srv is the current binary
1897 approximation to tdv (and dv). It should be exactly representable
1898 in an IEEE 754 double.
1899 */
1900
1901 for(;;) {
1902
1903 /* This is the main correction loop for sb_strtod.
1904
1905 We've got a decimal value tdv, and a floating-point approximation
1906 srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is
1907 close enough (i.e., within 0.5 ulps) to tdv, and to compute a new
1908 approximation if not.
1909
1910 To determine whether srv is close enough to tdv, compute integers
1911 bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv)
1912 respectively, and then use integer arithmetic to determine whether
1913 |tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv).
1914 */
1915
1916 bd = Balloc(bd0->k);
1917 if (bd == NULL) {
1918 Bfree(bd0);
1919 goto failed_malloc;
1920 }
1921 Bcopy(bd, bd0);
1922 bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */
1923 if (bb == NULL) {
1924 Bfree(bd);
1925 Bfree(bd0);
1926 goto failed_malloc;
1927 }
1928 /* Record whether lsb of bb is odd, in case we need this
1929 for the round-to-even step later. */
1930 odd = bb->x[0] & 1;
1931
1932 /* tdv = bd * 10**e; srv = bb * 2**bbe */
1933 bs = i2b(1);
1934 if (bs == NULL) {
1935 Bfree(bb);
1936 Bfree(bd);
1937 Bfree(bd0);
1938 goto failed_malloc;
1939 }
1940
1941 if (e >= 0) {
1942 bb2 = bb5 = 0;
1943 bd2 = bd5 = e;
1944 }
1945 else {
1946 bb2 = bb5 = -e;
1947 bd2 = bd5 = 0;
1948 }
1949 if (bbe >= 0)
1950 bb2 += bbe;
1951 else
1952 bd2 -= bbe;
1953 bs2 = bb2;
1954 bb2++;
1955 bd2++;
1956
1957 /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1,
1958 and bs == 1, so:
1959
1960 tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5)
1961 srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2)
1962 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2)
1963
1964 It follows that:
1965
1966 M * tdv = bd * 2**bd2 * 5**bd5
1967 M * srv = bb * 2**bb2 * 5**bb5
1968 M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5
1969
1970 for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but
1971 this fact is not needed below.)
1972 */
1973
1974 /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */
1975 i = bb2 < bd2 ? bb2 : bd2;
1976 if (i > bs2)
1977 i = bs2;
1978 if (i > 0) {
1979 bb2 -= i;
1980 bd2 -= i;
1981 bs2 -= i;
1982 }
1983
1984 /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
1985 if (bb5 > 0) {
1986 bs = pow5mult(bs, bb5);
1987 if (bs == NULL) {
1988 Bfree(bb);
1989 Bfree(bd);
1990 Bfree(bd0);
1991 goto failed_malloc;
1992 }
1993 bb1 = mult(bs, bb);
1994 Bfree(bb);
1995 bb = bb1;
1996 if (bb == NULL) {
1997 Bfree(bs);
1998 Bfree(bd);
1999 Bfree(bd0);
2000 goto failed_malloc;
2001 }
2002 }
2003 if (bb2 > 0) {
2004 bb = lshift(bb, bb2);
2005 if (bb == NULL) {
2006 Bfree(bs);
2007 Bfree(bd);
2008 Bfree(bd0);
2009 goto failed_malloc;
2010 }
2011 }
2012 if (bd5 > 0) {
2013 bd = pow5mult(bd, bd5);
2014 if (bd == NULL) {
2015 Bfree(bb);
2016 Bfree(bs);
2017 Bfree(bd0);
2018 goto failed_malloc;
2019 }
2020 }
2021 if (bd2 > 0) {
2022 bd = lshift(bd, bd2);
2023 if (bd == NULL) {
2024 Bfree(bb);
2025 Bfree(bs);
2026 Bfree(bd0);
2027 goto failed_malloc;
2028 }
2029 }
2030 if (bs2 > 0) {
2031 bs = lshift(bs, bs2);
2032 if (bs == NULL) {
2033 Bfree(bb);
2034 Bfree(bd);
2035 Bfree(bd0);
2036 goto failed_malloc;
2037 }
2038 }
2039
2040 /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),
2041 respectively. Compute the difference |tdv - srv|, and compare
2042 with 0.5 ulp(srv). */
2043
2044 delta = diff(bb, bd);
2045 if (delta == NULL) {
2046 Bfree(bb);
2047 Bfree(bs);
2048 Bfree(bd);
2049 Bfree(bd0);
2050 goto failed_malloc;
2051 }
2052 dsign = delta->sign;
2053 delta->sign = 0;
2054 i = cmp(delta, bs);
2055 if (bc.nd > nd && i <= 0) {
2056 if (dsign)
2057 break; /* Must use bigcomp(). */
2058
2059 /* Here rv overestimates the truncated decimal value by at most
2060 0.5 ulp(rv). Hence rv either overestimates the true decimal
2061 value by <= 0.5 ulp(rv), or underestimates it by some small
2062 amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of
2063 the true decimal value, so it's possible to exit.
2064
2065 Exception: if scaled rv is a normal exact power of 2, but not
2066 DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the
2067 next double, so the correctly rounded result is either rv - 0.5
2068 ulp(rv) or rv; in this case, use bigcomp to distinguish. */
2069
2070 if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) {
2071 /* rv can't be 0, since it's an overestimate for some
2072 nonzero value. So rv is a normal power of 2. */
2073 j = (int)(word0(&rv) & Exp_mask) >> Exp_shift;
2074 /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if
2075 rv / 2^bc.scale >= 2^-1021. */
2076 if (j - bc.scale >= 2) {
2077 dval(&rv) -= 0.5 * sulp(&rv, &bc);
2078 break; /* Use bigcomp. */
2079 }
2080 }
2081
2082 {
2083 bc.nd = nd;
2084 i = -1; /* Discarded digits make delta smaller. */
2085 }
2086 }
2087
2088 if (i < 0) {
2089 /* Error is less than half an ulp -- check for
2090 * special case of mantissa a power of two.
2091 */
2092 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
2093 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2094 ) {
2095 break;
2096 }
2097 if (!delta->x[0] && delta->wds <= 1) {
2098 /* exact result */
2099 break;
2100 }
2101 delta = lshift(delta,Log2P);
2102 if (delta == NULL) {
2103 Bfree(bb);
2104 Bfree(bs);
2105 Bfree(bd);
2106 Bfree(bd0);
2107 goto failed_malloc;
2108 }
2109 if (cmp(delta, bs) > 0)
2110 goto drop_down;
2111 break;
2112 }
2113 if (i == 0) {
2114 /* exactly half-way between */
2115 if (dsign) {
2116 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
2117 && word1(&rv) == (
2118 (bc.scale &&
2119 (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ?
2120 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2121 0xffffffff)) {
2122 /*boundary case -- increment exponent*/
2123 word0(&rv) = (word0(&rv) & Exp_mask)
2124 + Exp_msk1
2125 ;
2126 word1(&rv) = 0;
2127 /* dsign = 0; */
2128 break;
2129 }
2130 }
2131 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
2132 drop_down:
2133 /* boundary case -- decrement exponent */
2134 if (bc.scale) {
2135 L = word0(&rv) & Exp_mask;
2136 if (L <= (2*P+1)*Exp_msk1) {
2137 if (L > (P+2)*Exp_msk1)
2138 /* round even ==> */
2139 /* accept rv */
2140 break;
2141 /* rv = smallest denormal */
2142 if (bc.nd > nd)
2143 break;
2144 goto undfl;
2145 }
2146 }
2147 L = (word0(&rv) & Exp_mask) - Exp_msk1;
2148 word0(&rv) = L | Bndry_mask1;
2149 word1(&rv) = 0xffffffff;
2150 break;
2151 }
2152 if (!odd)
2153 break;
2154 if (dsign)
2155 dval(&rv) += sulp(&rv, &bc);
2156 else {
2157 dval(&rv) -= sulp(&rv, &bc);
2158 if (!dval(&rv)) {
2159 if (bc.nd >nd)
2160 break;
2161 goto undfl;
2162 }
2163 }
2164 /* dsign = 1 - dsign; */
2165 break;
2166 }
2167 if ((aadj = ratio(delta, bs)) <= 2.) {
2168 if (dsign)
2169 aadj = aadj1 = 1.;
2170 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
2171 if (word1(&rv) == Tiny1 && !word0(&rv)) {
2172 if (bc.nd >nd)
2173 break;
2174 goto undfl;
2175 }
2176 aadj = 1.;
2177 aadj1 = -1.;
2178 }
2179 else {
2180 /* special case -- power of FLT_RADIX to be */
2181 /* rounded down... */
2182
2183 if (aadj < 2./FLT_RADIX)
2184 aadj = 1./FLT_RADIX;
2185 else
2186 aadj *= 0.5;
2187 aadj1 = -aadj;
2188 }
2189 }
2190 else {
2191 aadj *= 0.5;
2192 aadj1 = dsign ? aadj : -aadj;
2193 if (Flt_Rounds == 0)
2194 aadj1 += 0.5;
2195 }
2196 y = word0(&rv) & Exp_mask;
2197
2198 /* Check for overflow */
2199
2200 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2201 dval(&rv0) = dval(&rv);
2202 word0(&rv) -= P*Exp_msk1;
2203 adj.d = aadj1 * ulp(&rv);
2204 dval(&rv) += adj.d;
2205 if ((word0(&rv) & Exp_mask) >=
2206 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2207 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) {
2208 Bfree(bb);
2209 Bfree(bd);
2210 Bfree(bs);
2211 Bfree(bd0);
2212 Bfree(delta);
2213 goto ovfl;
2214 }
2215 word0(&rv) = Big0;
2216 word1(&rv) = Big1;
2217 goto cont;
2218 }
2219 else
2220 word0(&rv) += P*Exp_msk1;
2221 }
2222 else {
2223 if (bc.scale && y <= 2*P*Exp_msk1) {
2224 if (aadj <= 0x7fffffff) {
2225 if ((z = (ULong)aadj) <= 0)
2226 z = 1;
2227 aadj = z;
2228 aadj1 = dsign ? aadj : -aadj;
2229 }
2230 dval(&aadj2) = aadj1;
2231 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
2232 aadj1 = dval(&aadj2);
2233 }
2234 adj.d = aadj1 * ulp(&rv);
2235 dval(&rv) += adj.d;
2236 }
2237 z = word0(&rv) & Exp_mask;
2238 if (bc.nd == nd) {
2239 if (!bc.scale)
2240 if (y == z) {
2241 /* Can we stop now? */
2242 L = (Long)aadj;
2243 aadj -= L;
2244 /* The tolerances below are conservative. */
2245 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
2246 if (aadj < .4999999 || aadj > .5000001)
2247 break;
2248 }
2249 else if (aadj < .4999999/FLT_RADIX)
2250 break;
2251 }
2252 }
2253 cont:
2254 Bfree(bb);
2255 Bfree(bd);
2256 Bfree(bs);
2257 Bfree(delta);
2258 }
2259 Bfree(bb);
2260 Bfree(bd);
2261 Bfree(bs);
2262 Bfree(bd0);
2263 Bfree(delta);
2264 if (bc.nd > nd) {
2265 error = bigcomp(&rv, s0, &bc);
2266 if (error)
2267 goto failed_malloc;
2268 }
2269
2270 if (bc.scale) {
2271 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
2272 word1(&rv0) = 0;
2273 dval(&rv) *= dval(&rv0);
2274 }
2275
2276 ret:
2277 return sign ? -dval(&rv) : dval(&rv);
2278
2279 parse_error:
2280 return 0.0;
2281
2282 failed_malloc:
2283 errno = ENOMEM;
2284 return -1.0;
2285
2286 undfl:
2287 return sign ? -0.0 : 0.0;
2288
2289 ovfl:
2290 errno = ERANGE;
2291 /* Can't trust HUGE_VAL */
2292 word0(&rv) = Exp_mask;
2293 word1(&rv) = 0;
2294 return sign ? -dval(&rv) : dval(&rv);
2295
2296}
2297
2298static char *
2299rv_alloc(int i)
2300{
2301 int j, k, *r;
2302
2303 j = sizeof(ULong);
2304 for(k = 0;
2305 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
2306 j <<= 1)
2307 k++;
2308 r = (int*)Balloc(k);
2309 if (r == NULL)
2310 return NULL;
2311 *r = k;
2312 return (char *)(r+1);
2313}
2314
2315static char *
2316nrv_alloc(char *s, char **rve, int n)
2317{
2318 char *rv, *t;
2319
2320 rv = rv_alloc(n);
2321 if (rv == NULL)
2322 return NULL;
2323 t = rv;
2324 while((*t = *s++)) t++;
2325 if (rve)
2326 *rve = t;
2327 return rv;
2328}
2329
2330/* freedtoa(s) must be used to free values s returned by dtoa
2331 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2332 * but for consistency with earlier versions of dtoa, it is optional
2333 * when MULTIPLE_THREADS is not defined.
2334 */
2335
2336void
2337sb_freedtoa(char *s)
2338{
2339 Bigint *b = (Bigint *)((int *)s - 1);
2340 b->maxwds = 1 << (b->k = *(int*)b);
2341 Bfree(b);
2342}
2343
2344/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2345 *
2346 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2347 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2348 *
2349 * Modifications:
2350 * 1. Rather than iterating, we use a simple numeric overestimate
2351 * to determine k = floor(log10(d)). We scale relevant
2352 * quantities using O(log2(k)) rather than O(k) multiplications.
2353 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2354 * try to generate digits strictly left to right. Instead, we
2355 * compute with fewer bits and propagate the carry if necessary
2356 * when rounding the final digit up. This is often faster.
2357 * 3. Under the assumption that input will be rounded nearest,
2358 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2359 * That is, we allow equality in stopping tests when the
2360 * round-nearest rule will give the same floating-point value
2361 * as would satisfaction of the stopping test with strict
2362 * inequality.
2363 * 4. We remove common factors of powers of 2 from relevant
2364 * quantities.
2365 * 5. When converting floating-point integers less than 1e16,
2366 * we use floating-point arithmetic rather than resorting
2367 * to multiple-precision integers.
2368 * 6. When asked to produce fewer than 15 digits, we first try
2369 * to get by with floating-point arithmetic; we resort to
2370 * multiple-precision integer arithmetic only if we cannot
2371 * guarantee that the floating-point calculation has given
2372 * the correctly rounded result. For k requested digits and
2373 * "uniformly" distributed input, the probability is
2374 * something like 10^(k-15) that we must resort to the Long
2375 * calculation.
2376 */
2377
2378/* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2379 leakage, a successful call to sb_dtoa should always be matched by a
2380 call to sb_freedtoa. */
2381
2382char *
2383sb_dtoa(double dd, int mode, int ndigits,
2384 int *decpt, int *sign, char **rve)
2385{
2386 /* Arguments ndigits, decpt, sign are similar to those
2387 of ecvt and fcvt; trailing zeros are suppressed from
2388 the returned string. If not null, *rve is set to point
2389 to the end of the return value. If d is +-Infinity or NaN,
2390 then *decpt is set to 9999.
2391
2392 mode:
2393 0 ==> shortest string that yields d when read in
2394 and rounded to nearest.
2395 1 ==> like 0, but with Steele & White stopping rule;
2396 e.g. with IEEE P754 arithmetic , mode 0 gives
2397 1e23 whereas mode 1 gives 9.999999999999999e22.
2398 2 ==> max(1,ndigits) significant digits. This gives a
2399 return value similar to that of ecvt, except
2400 that trailing zeros are suppressed.
2401 3 ==> through ndigits past the decimal point. This
2402 gives a return value similar to that from fcvt,
2403 except that trailing zeros are suppressed, and
2404 ndigits can be negative.
2405 4,5 ==> similar to 2 and 3, respectively, but (in
2406 round-nearest mode) with the tests of mode 0 to
2407 possibly return a shorter string that rounds to d.
2408 With IEEE arithmetic and compilation with
2409 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2410 as modes 2 and 3 when FLT_ROUNDS != 1.
2411 6-9 ==> Debugging modes similar to mode - 4: don't try
2412 fast floating-point estimate (if applicable).
2413
2414 Values of mode other than 0-9 are treated as mode 0.
2415
2416 Sufficient space is allocated to the return value
2417 to hold the suppressed trailing zeros.
2418 */
2419
2420 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2421 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2422 spec_case, try_quick;
2423 Long L;
2424 int denorm;
2425 ULong x;
2426 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2427 U d2, eps, u;
2428 double ds;
2429 char *s, *s0;
2430
2431 /* set pointers to NULL, to silence gcc compiler warnings and make
2432 cleanup easier on error */
2433 mlo = mhi = S = 0;
2434 s0 = 0;
2435
2436 u.d = dd;
2437 if (word0(&u) & Sign_bit) {
2438 /* set sign for everything, including 0's and NaNs */
2439 *sign = 1;
2440 word0(&u) &= ~Sign_bit; /* clear sign bit */
2441 }
2442 else
2443 *sign = 0;
2444
2445 /* quick return for Infinities, NaNs and zeros */
2446 if ((word0(&u) & Exp_mask) == Exp_mask)
2447 {
2448 /* Infinity or NaN */
2449 *decpt = 9999;
2450 if (!word1(&u) && !(word0(&u) & 0xfffff))
2451 return nrv_alloc("Infinity", rve, 8);
2452 return nrv_alloc("NaN", rve, 3);
2453 }
2454 if (!dval(&u)) {
2455 *decpt = 1;
2456 return nrv_alloc("0", rve, 1);
2457 }
2458
2459 /* compute k = floor(log10(d)). The computation may leave k
2460 one too large, but should never leave k too small. */
2461 b = d2b(&u, &be, &bbits);
2462 if (b == NULL)
2463 goto failed_malloc;
2464 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
2465 dval(&d2) = dval(&u);
2466 word0(&d2) &= Frac_mask1;
2467 word0(&d2) |= Exp_11;
2468
2469 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2470 * log10(x) = log(x) / log(10)
2471 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2472 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2473 *
2474 * This suggests computing an approximation k to log10(d) by
2475 *
2476 * k = (i - Bias)*0.301029995663981
2477 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2478 *
2479 * We want k to be too large rather than too small.
2480 * The error in the first-order Taylor series approximation
2481 * is in our favor, so we just round up the constant enough
2482 * to compensate for any error in the multiplication of
2483 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2484 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2485 * adding 1e-13 to the constant term more than suffices.
2486 * Hence we adjust the constant term to 0.1760912590558.
2487 * (We could get a more accurate k by invoking log10,
2488 * but this is probably not worthwhile.)
2489 */
2490
2491 i -= Bias;
2492 denorm = 0;
2493 }
2494 else {
2495 /* d is denormalized */
2496
2497 i = bbits + be + (Bias + (P-1) - 1);
2498 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
2499 : word1(&u) << (32 - i);
2500 dval(&d2) = x;
2501 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
2502 i -= (Bias + (P-1) - 1) + 1;
2503 denorm = 1;
2504 }
2505 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2506 i*0.301029995663981;
2507 k = (int)ds;
2508 if (ds < 0. && ds != k)
2509 k--; /* want k = floor(ds) */
2510 k_check = 1;
2511 if (k >= 0 && k <= Ten_pmax) {
2512 if (dval(&u) < tens[k])
2513 k--;
2514 k_check = 0;
2515 }
2516 j = bbits - i - 1;
2517 if (j >= 0) {
2518 b2 = 0;
2519 s2 = j;
2520 }
2521 else {
2522 b2 = -j;
2523 s2 = 0;
2524 }
2525 if (k >= 0) {
2526 b5 = 0;
2527 s5 = k;
2528 s2 += k;
2529 }
2530 else {
2531 b2 -= k;
2532 b5 = -k;
2533 s5 = 0;
2534 }
2535 if (mode < 0 || mode > 9)
2536 mode = 0;
2537
2538 try_quick = 1;
2539
2540 if (mode > 5) {
2541 mode -= 4;
2542 try_quick = 0;
2543 }
2544 leftright = 1;
2545 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
2546 /* silence erroneous "gcc -Wall" warning. */
2547 switch(mode) {
2548 case 0:
2549 case 1:
2550 i = 18;
2551 ndigits = 0;
2552 break;
2553 case 2:
2554 leftright = 0;
2555 /* no break */
2556 case 4:
2557 if (ndigits <= 0)
2558 ndigits = 1;
2559 ilim = ilim1 = i = ndigits;
2560 break;
2561 case 3:
2562 leftright = 0;
2563 /* no break */
2564 case 5:
2565 i = ndigits + k + 1;
2566 ilim = i;
2567 ilim1 = i - 1;
2568 if (i <= 0)
2569 i = 1;
2570 }
2571 s0 = rv_alloc(i);
2572 if (s0 == NULL)
2573 goto failed_malloc;
2574 s = s0;
2575
2576
2577 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2578
2579 /* Try to get by with floating-point arithmetic. */
2580
2581 i = 0;
2582 dval(&d2) = dval(&u);
2583 k0 = k;
2584 ilim0 = ilim;
2585 ieps = 2; /* conservative */
2586 if (k > 0) {
2587 ds = tens[k&0xf];
2588 j = k >> 4;
2589 if (j & Bletch) {
2590 /* prevent overflows */
2591 j &= Bletch - 1;
2592 dval(&u) /= bigtens[n_bigtens-1];
2593 ieps++;
2594 }
2595 for(; j; j >>= 1, i++)
2596 if (j & 1) {
2597 ieps++;
2598 ds *= bigtens[i];
2599 }
2600 dval(&u) /= ds;
2601 }
2602 else if ((j1 = -k)) {
2603 dval(&u) *= tens[j1 & 0xf];
2604 for(j = j1 >> 4; j; j >>= 1, i++)
2605 if (j & 1) {
2606 ieps++;
2607 dval(&u) *= bigtens[i];
2608 }
2609 }
2610 if (k_check && dval(&u) < 1. && ilim > 0) {
2611 if (ilim1 <= 0)
2612 goto fast_failed;
2613 ilim = ilim1;
2614 k--;
2615 dval(&u) *= 10.;
2616 ieps++;
2617 }
2618 dval(&eps) = ieps*dval(&u) + 7.;
2619 word0(&eps) -= (P-1)*Exp_msk1;
2620 if (ilim == 0) {
2621 S = mhi = 0;
2622 dval(&u) -= 5.;
2623 if (dval(&u) > dval(&eps))
2624 goto one_digit;
2625 if (dval(&u) < -dval(&eps))
2626 goto no_digits;
2627 goto fast_failed;
2628 }
2629 if (leftright) {
2630 /* Use Steele & White method of only
2631 * generating digits needed.
2632 */
2633 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
2634 for(i = 0;;) {
2635 L = (Long)dval(&u);
2636 dval(&u) -= L;
2637 *s++ = '0' + (int)L;
2638 if (dval(&u) < dval(&eps))
2639 goto ret1;
2640 if (1. - dval(&u) < dval(&eps))
2641 goto bump_up;
2642 if (++i >= ilim)
2643 break;
2644 dval(&eps) *= 10.;
2645 dval(&u) *= 10.;
2646 }
2647 }
2648 else {
2649 /* Generate ilim digits, then fix them up. */
2650 dval(&eps) *= tens[ilim-1];
2651 for(i = 1;; i++, dval(&u) *= 10.) {
2652 L = (Long)(dval(&u));
2653 if (!(dval(&u) -= L))
2654 ilim = i;
2655 *s++ = '0' + (int)L;
2656 if (i == ilim) {
2657 if (dval(&u) > 0.5 + dval(&eps))
2658 goto bump_up;
2659 else if (dval(&u) < 0.5 - dval(&eps)) {
2660 while(*--s == '0');
2661 s++;
2662 goto ret1;
2663 }
2664 break;
2665 }
2666 }
2667 }
2668 fast_failed:
2669 s = s0;
2670 dval(&u) = dval(&d2);
2671 k = k0;
2672 ilim = ilim0;
2673 }
2674
2675 /* Do we have a "small" integer? */
2676
2677 if (be >= 0 && k <= Int_max) {
2678 /* Yes. */
2679 ds = tens[k];
2680 if (ndigits < 0 && ilim <= 0) {
2681 S = mhi = 0;
2682 if (ilim < 0 || dval(&u) <= 5*ds)
2683 goto no_digits;
2684 goto one_digit;
2685 }
2686 for(i = 1;; i++, dval(&u) *= 10.) {
2687 L = (Long)(dval(&u) / ds);
2688 dval(&u) -= L*ds;
2689 *s++ = '0' + (int)L;
2690 if (!dval(&u)) {
2691 break;
2692 }
2693 if (i == ilim) {
2694 dval(&u) += dval(&u);
2695 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
2696 bump_up:
2697 while(*--s == '9')
2698 if (s == s0) {
2699 k++;
2700 *s = '0';
2701 break;
2702 }
2703 ++*s++;
2704 }
2705 break;
2706 }
2707 }
2708 goto ret1;
2709 }
2710
2711 m2 = b2;
2712 m5 = b5;
2713 if (leftright) {
2714 i =
2715 denorm ? be + (Bias + (P-1) - 1 + 1) :
2716 1 + P - bbits;
2717 b2 += i;
2718 s2 += i;
2719 mhi = i2b(1);
2720 if (mhi == NULL)
2721 goto failed_malloc;
2722 }
2723 if (m2 > 0 && s2 > 0) {
2724 i = m2 < s2 ? m2 : s2;
2725 b2 -= i;
2726 m2 -= i;
2727 s2 -= i;
2728 }
2729 if (b5 > 0) {
2730 if (leftright) {
2731 if (m5 > 0) {
2732 mhi = pow5mult(mhi, m5);
2733 if (mhi == NULL)
2734 goto failed_malloc;
2735 b1 = mult(mhi, b);
2736 Bfree(b);
2737 b = b1;
2738 if (b == NULL)
2739 goto failed_malloc;
2740 }
2741 if ((j = b5 - m5)) {
2742 b = pow5mult(b, j);
2743 if (b == NULL)
2744 goto failed_malloc;
2745 }
2746 }
2747 else {
2748 b = pow5mult(b, b5);
2749 if (b == NULL)
2750 goto failed_malloc;
2751 }
2752 }
2753 S = i2b(1);
2754 if (S == NULL)
2755 goto failed_malloc;
2756 if (s5 > 0) {
2757 S = pow5mult(S, s5);
2758 if (S == NULL)
2759 goto failed_malloc;
2760 }
2761
2762 /* Check for special case that d is a normalized power of 2. */
2763
2764 spec_case = 0;
2765 if ((mode < 2 || leftright)
2766 ) {
2767 if (!word1(&u) && !(word0(&u) & Bndry_mask)
2768 && word0(&u) & (Exp_mask & ~Exp_msk1)
2769 ) {
2770 /* The special case */
2771 b2 += Log2P;
2772 s2 += Log2P;
2773 spec_case = 1;
2774 }
2775 }
2776
2777 /* Arrange for convenient computation of quotients:
2778 * shift left if necessary so divisor has 4 leading 0 bits.
2779 *
2780 * Perhaps we should just compute leading 28 bits of S once
2781 * and for all and pass them and a shift to quorem, so it
2782 * can do shifts and ors to compute the numerator for q.
2783 */
2784#define iInc 28
2785 i = dshift(S, s2);
2786 b2 += i;
2787 m2 += i;
2788 s2 += i;
2789 if (b2 > 0) {
2790 b = lshift(b, b2);
2791 if (b == NULL)
2792 goto failed_malloc;
2793 }
2794 if (s2 > 0) {
2795 S = lshift(S, s2);
2796 if (S == NULL)
2797 goto failed_malloc;
2798 }
2799 if (k_check) {
2800 if (cmp(b,S) < 0) {
2801 k--;
2802 b = multadd(b, 10, 0); /* we botched the k estimate */
2803 if (b == NULL)
2804 goto failed_malloc;
2805 if (leftright) {
2806 mhi = multadd(mhi, 10, 0);
2807 if (mhi == NULL)
2808 goto failed_malloc;
2809 }
2810 ilim = ilim1;
2811 }
2812 }
2813 if (ilim <= 0 && (mode == 3 || mode == 5)) {
2814 if (ilim < 0) {
2815 /* no digits, fcvt style */
2816 no_digits:
2817 k = -1 - ndigits;
2818 goto ret;
2819 }
2820 else {
2821 S = multadd(S, 5, 0);
2822 if (S == NULL)
2823 goto failed_malloc;
2824 if (cmp(b, S) <= 0)
2825 goto no_digits;
2826 }
2827 one_digit:
2828 *s++ = '1';
2829 k++;
2830 goto ret;
2831 }
2832 if (leftright) {
2833 if (m2 > 0) {
2834 mhi = lshift(mhi, m2);
2835 if (mhi == NULL)
2836 goto failed_malloc;
2837 }
2838
2839 /* Compute mlo -- check for special case
2840 * that d is a normalized power of 2.
2841 */
2842
2843 mlo = mhi;
2844 if (spec_case) {
2845 mhi = Balloc(mhi->k);
2846 if (mhi == NULL)
2847 goto failed_malloc;
2848 Bcopy(mhi, mlo);
2849 mhi = lshift(mhi, Log2P);
2850 if (mhi == NULL)
2851 goto failed_malloc;
2852 }
2853
2854 for(i = 1;;i++) {
2855 dig = quorem(b,S) + '0';
2856 /* Do we yet have the shortest decimal string
2857 * that will round to d?
2858 */
2859 j = cmp(b, mlo);
2860 delta = diff(S, mhi);
2861 if (delta == NULL)
2862 goto failed_malloc;
2863 j1 = delta->sign ? 1 : cmp(b, delta);
2864 Bfree(delta);
2865 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
2866 ) {
2867 if (dig == '9')
2868 goto round_9_up;
2869 if (j > 0)
2870 dig++;
2871 *s++ = dig;
2872 goto ret;
2873 }
2874 if (j < 0 || (j == 0 && mode != 1
2875 && !(word1(&u) & 1)
2876 )) {
2877 if (!b->x[0] && b->wds <= 1) {
2878 goto accept_dig;
2879 }
2880 if (j1 > 0) {
2881 b = lshift(b, 1);
2882 if (b == NULL)
2883 goto failed_malloc;
2884 j1 = cmp(b, S);
2885 if ((j1 > 0 || (j1 == 0 && dig & 1))
2886 && dig++ == '9')
2887 goto round_9_up;
2888 }
2889 accept_dig:
2890 *s++ = dig;
2891 goto ret;
2892 }
2893 if (j1 > 0) {
2894 if (dig == '9') { /* possible if i == 1 */
2895 round_9_up:
2896 *s++ = '9';
2897 goto roundoff;
2898 }
2899 *s++ = dig + 1;
2900 goto ret;
2901 }
2902 *s++ = dig;
2903 if (i == ilim)
2904 break;
2905 b = multadd(b, 10, 0);
2906 if (b == NULL)
2907 goto failed_malloc;
2908 if (mlo == mhi) {
2909 mlo = mhi = multadd(mhi, 10, 0);
2910 if (mlo == NULL)
2911 goto failed_malloc;
2912 }
2913 else {
2914 mlo = multadd(mlo, 10, 0);
2915 if (mlo == NULL)
2916 goto failed_malloc;
2917 mhi = multadd(mhi, 10, 0);
2918 if (mhi == NULL)
2919 goto failed_malloc;
2920 }
2921 }
2922 }
2923 else
2924 for(i = 1;; i++) {
2925 *s++ = dig = quorem(b,S) + '0';
2926 if (!b->x[0] && b->wds <= 1) {
2927 goto ret;
2928 }
2929 if (i >= ilim)
2930 break;
2931 b = multadd(b, 10, 0);
2932 if (b == NULL)
2933 goto failed_malloc;
2934 }
2935
2936 /* Round off last digit */
2937
2938 b = lshift(b, 1);
2939 if (b == NULL)
2940 goto failed_malloc;
2941 j = cmp(b, S);
2942 if (j > 0 || (j == 0 && dig & 1)) {
2943 roundoff:
2944 while(*--s == '9')
2945 if (s == s0) {
2946 k++;
2947 *s++ = '1';
2948 goto ret;
2949 }
2950 ++*s++;
2951 }
2952 else {
2953 while(*--s == '0');
2954 s++;
2955 }
2956 ret:
2957 Bfree(S);
2958 if (mhi) {
2959 if (mlo && mlo != mhi)
2960 Bfree(mlo);
2961 Bfree(mhi);
2962 }
2963 ret1:
2964 Bfree(b);
2965 *s = 0;
2966 *decpt = k + 1;
2967 if (rve)
2968 *rve = s;
2969 return s0;
2970 failed_malloc:
2971 if (S)
2972 Bfree(S);
2973 if (mlo && mlo != mhi)
2974 Bfree(mlo);
2975 if (mhi)
2976 Bfree(mhi);
2977 if (b)
2978 Bfree(b);
2979 if (s0)
2980 sb_freedtoa(s0);
2981 return NULL;
2982}
2983#ifdef __cplusplus
2984}
2985#endif
Sphinx's memory allocation/deallocation routines.
Basic type definitions used in Sphinx.
Definition dtoa.c:289
Definition dtoa.c:320
Definition dtoa.c:178