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
175 extern "C" {
176 #endif
177 
178 typedef 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 
287 typedef struct BCinfo BCinfo;
288 struct
289 BCinfo {
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 
319 struct
320 Bigint {
321  struct Bigint *next;
322  int k, maxwds, sign, wds;
323  ULong x[1];
324 };
325 
326 typedef 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))
336 static 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 
357 static Bigint *freelist[Kmax+1];
358 
359 /* Allocate space for a Bigint with up to 1<<k digits */
360 
361 static Bigint *
362 Balloc(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 
392 static void
393 Bfree(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 
414 static Bigint *
415 Balloc(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 
437 static void
438 Bfree(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 
454 static Bigint *
455 multadd(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 
508 static Bigint *
509 s2b(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 
543 static int
544 hi0bits(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 
575 static int
576 lo0bits(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 
620 static Bigint *
621 i2b(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 
636 static Bigint *
637 mult(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 
737 static 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 
743 static Bigint *
744 pow5mult(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 
799 static Bigint *
800 pow5mult(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 
850 static Bigint *
851 lshift(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 
897 static int
898 cmp(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 
930 static Bigint *
931 diff(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 
1010 static double
1011 ulp(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 
1024 static double
1025 b2d(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 
1080 static Bigint *
1081 sd2b(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 
1138 static Bigint *
1139 d2b(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 
1186 static double
1187 ratio(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 
1204 static const double
1205 tens[] = {
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 
1211 static const double
1212 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1213 static 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 
1227 static int
1228 dshift(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 
1240 static int
1241 quorem(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 
1341 static double
1342 sulp(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 
1404 static int
1405 bigcomp(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 
1521 double
1522 sb_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 
1535 double
1536 sb_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 
1544 double
1545 sb_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 
2298 static char *
2299 rv_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 
2315 static char *
2316 nrv_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 
2336 void
2337 sb_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 
2382 char *
2383 sb_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&#39;s memory allocation/deallocation routines.
Basic type definitions used in Sphinx.
Definition: dtoa.c:319
Definition: dtoa.c:288
Definition: dtoa.c:178