SphinxBase  0.6
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 /* Please send bug reports to David M. Gay (dmg at acm dot org,
21  * with " at " changed at "@" and " dot " changed to "."). */
22 
23 /* On a machine with IEEE extended-precision registers, it is
24  * necessary to specify double-precision (53-bit) rounding precision
25  * before invoking strtod or dtoa. If the machine uses (the equivalent
26  * of) Intel 80x87 arithmetic, the call
27  * _control87(PC_53, MCW_PC);
28  * does this with many compilers. Whether this or another call is
29  * appropriate depends on the compiler; for this to work, it may be
30  * necessary to #include "float.h" or another system-dependent header
31  * file.
32  */
33 
34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
35  *
36  * This strtod returns a nearest machine number to the input decimal
37  * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
38  * broken by the IEEE round-even rule. Otherwise ties are broken by
39  * biased rounding (add half and chop).
40  *
41  * Inspired loosely by William D. Clinger's paper "How to Read Floating
42  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
43  *
44  * Modifications:
45  *
46  * 1. We only require IEEE, IBM, or VAX double-precision
47  * arithmetic (not IEEE double-extended).
48  * 2. We get by with floating-point arithmetic in a case that
49  * Clinger missed -- when we're computing d * 10^n
50  * for a small integer d and the integer n is not too
51  * much larger than 22 (the maximum integer k for which
52  * we can represent 10^k exactly), we may be able to
53  * compute (d*10^k) * 10^(e-k) with just one roundoff.
54  * 3. Rather than a bit-at-a-time adjustment of the binary
55  * result in the hard case, we use floating-point
56  * arithmetic to determine the adjustment to within
57  * one bit; only in really hard cases do we need to
58  * compute a second residual.
59  * 4. Because of 3., we don't need a large table of powers of 10
60  * for ten-to-e (just some small tables, e.g. of 10^k
61  * for 0 <= k <= 22).
62  */
63 
64 /*
65  * This file has been modified to remove dtoa() and all
66  * non-reentrancy. This makes it slower, but it also makes life a lot
67  * easier on Windows and other platforms without static lock
68  * initializers (grumble).
69  */
70 
71 /* Added by dhuggins@cs.cmu.edu to use autoconf results. */
72 /* We do not care about the VAX. */
73 #include "config.h"
74 #ifdef WORDS_BIGENDIAN
75 #define IEEE_MC68k
76 #else
77 #define IEEE_8087
78 #endif
79 #ifndef HAVE_LONG_LONG
80 #define NO_LONG_LONG
81 #endif
82 #define Omit_Private_Memory
83 #include "sphinxbase/ckd_alloc.h"
84 #undef USE_LOCALE
85 
86 /* Correct totally bogus typedefs in this code. */
87 #include "sphinxbase/prim_type.h"
88 #define Long int32 /* ZOMG */
89 #define ULong uint32 /* WTF */
90 
91 /*
92  * #define IEEE_8087 for IEEE-arithmetic machines where the least
93  * significant byte has the lowest address.
94  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
95  * significant byte has the lowest address.
96  * #define Long int on machines with 32-bit ints and 64-bit longs.
97  * #define IBM for IBM mainframe-style floating-point arithmetic.
98  * #define VAX for VAX-style floating-point arithmetic (D_floating).
99  * #define No_leftright to omit left-right logic in fast floating-point
100  * computation of dtoa.
101  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
102  * and strtod and dtoa should round accordingly.
103  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
104  * and Honor_FLT_ROUNDS is not #defined.
105  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
106  * that use extended-precision instructions to compute rounded
107  * products and quotients) with IBM.
108  * #define ROUND_BIASED for IEEE-format with biased rounding.
109  * #define Inaccurate_Divide for IEEE-format with correctly rounded
110  * products but inaccurate quotients, e.g., for Intel i860.
111  * #define NO_LONG_LONG on machines that do not have a "long long"
112  * integer type (of >= 64 bits). On such machines, you can
113  * #define Just_16 to store 16 bits per 32-bit Long when doing
114  * high-precision integer arithmetic. Whether this speeds things
115  * up or slows things down depends on the machine and the number
116  * being converted. If long long is available and the name is
117  * something other than "long long", #define Llong to be the name,
118  * and if "unsigned Llong" does not work as an unsigned version of
119  * Llong, #define #ULLong to be the corresponding unsigned type.
120  * #define KR_headers for old-style C function headers.
121  * #define Bad_float_h if your system lacks a float.h or if it does not
122  * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
123  * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
124  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
125  * if memory is available and otherwise does something you deem
126  * appropriate. If MALLOC is undefined, malloc will be invoked
127  * directly -- and assumed always to succeed.
128  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
129  * memory allocations from a private pool of memory when possible.
130  * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
131  * unless #defined to be a different length. This default length
132  * suffices to get rid of MALLOC calls except for unusual cases,
133  * such as decimal-to-binary conversion of a very long string of
134  * digits. The longest string dtoa can return is about 751 bytes
135  * long. For conversions by strtod of strings of 800 digits and
136  * all dtoa conversions in single-threaded executions with 8-byte
137  * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
138  * pointers, PRIVATE_MEM >= 7112 appears adequate.
139  * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
140  * #defined automatically on IEEE systems. On such systems,
141  * when INFNAN_CHECK is #defined, strtod checks
142  * for Infinity and NaN (case insensitively). On some systems
143  * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
144  * appropriately -- to the most significant word of a quiet NaN.
145  * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
146  * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
147  * strtod also accepts (case insensitively) strings of the form
148  * NaN(x), where x is a string of hexadecimal digits and spaces;
149  * if there is only one string of hexadecimal digits, it is taken
150  * for the 52 fraction bits of the resulting NaN; if there are two
151  * or more strings of hex digits, the first is for the high 20 bits,
152  * the second and subsequent for the low 32 bits, with intervening
153  * white space ignored; but if this results in none of the 52
154  * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
155  * and NAN_WORD1 are used instead.
156  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
157  * multiple threads. In this case, you must provide (or suitably
158  * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
159  * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
160  * in pow5mult, ensures lazy evaluation of only one copy of high
161  * powers of 5; omitting this lock would introduce a small
162  * probability of wasting memory, but would otherwise be harmless.)
163  * You must also invoke freedtoa(s) to free the value s returned by
164  * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
165  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
166  * avoids underflows on inputs whose result does not underflow.
167  * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
168  * floating-point numbers and flushes underflows to zero rather
169  * than implementing gradual underflow, then you must also #define
170  * Sudden_Underflow.
171  * #define YES_ALIAS to permit aliasing certain double values with
172  * arrays of ULongs. This leads to slightly better code with
173  * some compilers and was always used prior to 19990916, but it
174  * is not strictly legal and can cause trouble with aggressively
175  * optimizing compilers (e.g., gcc 2.95.1 under -O2).
176  * #define USE_LOCALE to use the current locale's decimal_point value.
177  * #define SET_INEXACT if IEEE arithmetic is being used and extra
178  * computation should be done to set the inexact flag when the
179  * result is inexact and avoid setting inexact when the result
180  * is exact. In this case, dtoa.c must be compiled in
181  * an environment, perhaps provided by #include "dtoa.c" in a
182  * suitable wrapper, that defines two functions,
183  * int get_inexact(void);
184  * void clear_inexact(void);
185  * such that get_inexact() returns a nonzero value if the
186  * inexact bit is already set, and clear_inexact() sets the
187  * inexact bit to 0. When SET_INEXACT is #defined, strtod
188  * also does extra computations to set the underflow and overflow
189  * flags when appropriate (i.e., when the result is tiny and
190  * inexact or when it is a numeric value rounded to +-infinity).
191  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
192  * the result overflows to +-Infinity or underflows to 0.
193  */
194 
195 #ifndef Long
196 #define Long long
197 #endif
198 #ifndef ULong
199 typedef unsigned Long ULong;
200 #endif
201 
202 #ifdef DEBUG
203 #include "stdio.h"
204 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
205 #endif
206 
207 #include "stdlib.h"
208 #include "string.h"
209 
210 #ifdef USE_LOCALE
211 #include "locale.h"
212 #endif
213 
214 /* Private memory and other non-reentrant stuff removed. */
215 
216 #undef IEEE_Arith
217 #undef Avoid_Underflow
218 #ifdef IEEE_MC68k
219 #define IEEE_Arith
220 #endif
221 #ifdef IEEE_8087
222 #define IEEE_Arith
223 #endif
224 
225 #ifdef IEEE_Arith
226 #ifndef NO_INFNAN_CHECK
227 #undef INFNAN_CHECK
228 #define INFNAN_CHECK
229 #endif
230 #else
231 #undef INFNAN_CHECK
232 #endif
233 
234 #include "errno.h"
235 
236 #ifdef Bad_float_h
237 
238 #ifdef IEEE_Arith
239 #define DBL_DIG 15
240 #define DBL_MAX_10_EXP 308
241 #define DBL_MAX_EXP 1024
242 #define FLT_RADIX 2
243 #endif /*IEEE_Arith*/
244 
245 #ifdef IBM
246 #define DBL_DIG 16
247 #define DBL_MAX_10_EXP 75
248 #define DBL_MAX_EXP 63
249 #define FLT_RADIX 16
250 #define DBL_MAX 7.2370055773322621e+75
251 #endif
252 
253 #ifdef VAX
254 #define DBL_DIG 16
255 #define DBL_MAX_10_EXP 38
256 #define DBL_MAX_EXP 127
257 #define FLT_RADIX 2
258 #define DBL_MAX 1.7014118346046923e+38
259 #endif
260 
261 #ifndef LONG_MAX
262 #define LONG_MAX 2147483647
263 #endif
264 
265 #else /* ifndef Bad_float_h */
266 #include "float.h"
267 #endif /* Bad_float_h */
268 
269 #ifndef __MATH_H__
270 #include "math.h"
271 #endif
272 
273 #ifdef __cplusplus
274 extern "C" {
275 #endif
276 
277 #ifndef CONST
278 #ifdef KR_headers
279 #define CONST /* blank */
280 #else
281 #define CONST const
282 #endif
283 #endif
284 
285 
286 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
287 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
288 #endif
289 
291 typedef union { double d; ULong L[2]; } U;
292 
293 #ifdef YES_ALIAS
294 #define dval(x) x
295 #ifdef IEEE_8087
296 #define word0(x) ((ULong *)&x)[1]
297 #define word1(x) ((ULong *)&x)[0]
298 #else
299 #define word0(x) ((ULong *)&x)[0]
300 #define word1(x) ((ULong *)&x)[1]
301 #endif
302 #else
303 #ifdef IEEE_8087
304 #define word0(x) ((U*)&x)->L[1]
305 #define word1(x) ((U*)&x)->L[0]
306 #else
307 #define word0(x) ((U*)&x)->L[0]
308 #define word1(x) ((U*)&x)->L[1]
309 #endif
310 #define dval(x) ((U*)&x)->d
311 #endif
312 
313 /* The following definition of Storeinc is appropriate for MIPS processors.
314  * An alternative that might be better on some machines is
315  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
316  */
317 #if defined(IEEE_8087) + defined(VAX)
318 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
319 ((unsigned short *)a)[0] = (unsigned short)c, a++)
320 #else
321 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
322 ((unsigned short *)a)[1] = (unsigned short)c, a++)
323 #endif
324 
325 /* #define P DBL_MANT_DIG */
326 /* Ten_pmax = floor(P*log(2)/log(5)) */
327 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
328 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
329 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
330 
331 #ifdef IEEE_Arith
332 #define Exp_shift 20
333 #define Exp_shift1 20
334 #define Exp_msk1 0x100000
335 #define Exp_msk11 0x100000
336 #define Exp_mask 0x7ff00000
337 #define P 53
338 #define Bias 1023
339 #define Emin (-1022)
340 #define Exp_1 0x3ff00000
341 #define Exp_11 0x3ff00000
342 #define Ebits 11
343 #define Frac_mask 0xfffff
344 #define Frac_mask1 0xfffff
345 #define Ten_pmax 22
346 #define Bletch 0x10
347 #define Bndry_mask 0xfffff
348 #define Bndry_mask1 0xfffff
349 #define LSB 1
350 #define Sign_bit 0x80000000
351 #define Log2P 1
352 #define Tiny0 0
353 #define Tiny1 1
354 #define Quick_max 14
355 #define Int_max 14
356 #ifndef NO_IEEE_Scale
357 #define Avoid_Underflow
358 #ifdef Flush_Denorm /* debugging option */
359 #undef Sudden_Underflow
360 #endif
361 #endif
362 
363 #ifndef Flt_Rounds
364 #ifdef FLT_ROUNDS
365 #define Flt_Rounds FLT_ROUNDS
366 #else
367 #define Flt_Rounds 1
368 #endif
369 #endif /*Flt_Rounds*/
370 
371 #ifdef Honor_FLT_ROUNDS
372 #define Rounding rounding
373 #undef Check_FLT_ROUNDS
374 #define Check_FLT_ROUNDS
375 #else
376 #define Rounding Flt_Rounds
377 #endif
378 
379 #else /* ifndef IEEE_Arith */
380 #undef Check_FLT_ROUNDS
381 #undef Honor_FLT_ROUNDS
382 #undef SET_INEXACT
383 #undef Sudden_Underflow
384 #define Sudden_Underflow
385 #ifdef IBM
386 #undef Flt_Rounds
387 #define Flt_Rounds 0
388 #define Exp_shift 24
389 #define Exp_shift1 24
390 #define Exp_msk1 0x1000000
391 #define Exp_msk11 0x1000000
392 #define Exp_mask 0x7f000000
393 #define P 14
394 #define Bias 65
395 #define Exp_1 0x41000000
396 #define Exp_11 0x41000000
397 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
398 #define Frac_mask 0xffffff
399 #define Frac_mask1 0xffffff
400 #define Bletch 4
401 #define Ten_pmax 22
402 #define Bndry_mask 0xefffff
403 #define Bndry_mask1 0xffffff
404 #define LSB 1
405 #define Sign_bit 0x80000000
406 #define Log2P 4
407 #define Tiny0 0x100000
408 #define Tiny1 0
409 #define Quick_max 14
410 #define Int_max 15
411 #else /* VAX */
412 #undef Flt_Rounds
413 #define Flt_Rounds 1
414 #define Exp_shift 23
415 #define Exp_shift1 7
416 #define Exp_msk1 0x80
417 #define Exp_msk11 0x800000
418 #define Exp_mask 0x7f80
419 #define P 56
420 #define Bias 129
421 #define Exp_1 0x40800000
422 #define Exp_11 0x4080
423 #define Ebits 8
424 #define Frac_mask 0x7fffff
425 #define Frac_mask1 0xffff007f
426 #define Ten_pmax 24
427 #define Bletch 2
428 #define Bndry_mask 0xffff007f
429 #define Bndry_mask1 0xffff007f
430 #define LSB 0x10000
431 #define Sign_bit 0x8000
432 #define Log2P 1
433 #define Tiny0 0x80
434 #define Tiny1 0
435 #define Quick_max 15
436 #define Int_max 15
437 #endif /* IBM, VAX */
438 #endif /* IEEE_Arith */
439 
440 #ifndef IEEE_Arith
441 #define ROUND_BIASED
442 #endif
443 
444 #ifdef RND_PRODQUOT
445 #define rounded_product(a,b) a = rnd_prod(a, b)
446 #define rounded_quotient(a,b) a = rnd_quot(a, b)
447 #ifdef KR_headers
448 extern double rnd_prod(), rnd_quot();
449 #else
450 extern double rnd_prod(double, double), rnd_quot(double, double);
451 #endif
452 #else
453 #define rounded_product(a,b) a *= b
454 #define rounded_quotient(a,b) a /= b
455 #endif
456 
457 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
458 #define Big1 0xffffffff
459 
460 #ifndef Pack_32
461 #define Pack_32
462 #endif
463 
464 #ifdef KR_headers
465 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
466 #else
467 #define FFFFFFFF 0xffffffffUL
468 #endif
469 
470 #ifdef NO_LONG_LONG
471 #undef ULLong
472 #ifdef Just_16
473 #undef Pack_32
474 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
475  * This makes some inner loops simpler and sometimes saves work
476  * during multiplications, but it often seems to make things slightly
477  * slower. Hence the default is now to store 32 bits per Long.
478  */
479 #endif
480 #else /* long long available */
481 #ifndef Llong
482 #define Llong long long
483 #endif
484 #ifndef ULLong
485 #define ULLong unsigned Llong
486 #endif
487 #endif /* NO_LONG_LONG */
488 
489 #ifndef MULTIPLE_THREADS
490 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
491 #define FREE_DTOA_LOCK(n) /*nothing*/
492 #endif
493 
494 #define Kmax 15
495 
496 #ifdef __cplusplus
497 extern "C" double sb_strtod(const char *s00, char **se);
498 #endif
499 
500  struct
501 Bigint {
502  struct Bigint *next;
503  int k, maxwds, sign, wds;
504  ULong x[1];
505  };
506 
507  typedef struct Bigint Bigint;
508 
509  static Bigint *
510 Balloc
511 #ifdef KR_headers
512  (k) int k;
513 #else
514  (int k)
515 #endif
516 {
517  int x;
518  size_t len;
519  Bigint *rv;
520 
521  x = 1 << k;
522  len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
523  /sizeof(double);
524  rv = ckd_malloc(len*sizeof(double));
525  rv->k = k;
526  rv->maxwds = x;
527  rv->sign = rv->wds = 0;
528  return rv;
529 }
530 
531  static void
532 Bfree
533 #ifdef KR_headers
534  (v) Bigint *v;
535 #else
536  (Bigint *v)
537 #endif
538 {
539  ckd_free(v);
540 }
541 
542 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
543 y->wds*sizeof(Long) + 2*sizeof(int))
544 
545  static Bigint *
546 multadd
547 #ifdef KR_headers
548  (b, m, a) Bigint *b; int m, a;
549 #else
550  (Bigint *b, int m, int a) /* multiply by m and add a */
551 #endif
552 {
553  int i, wds;
554 #ifdef ULLong
555  ULong *x;
556  ULLong carry, y;
557 #else
558  ULong carry, *x, y;
559 #ifdef Pack_32
560  ULong xi, z;
561 #endif
562 #endif
563  Bigint *b1;
564 
565  wds = b->wds;
566  x = b->x;
567  i = 0;
568  carry = a;
569  do {
570 #ifdef ULLong
571  y = *x * (ULLong)m + carry;
572  carry = y >> 32;
573  *x++ = y & FFFFFFFF;
574 #else
575 #ifdef Pack_32
576  xi = *x;
577  y = (xi & 0xffff) * m + carry;
578  z = (xi >> 16) * m + (y >> 16);
579  carry = z >> 16;
580  *x++ = (z << 16) + (y & 0xffff);
581 #else
582  y = *x * m + carry;
583  carry = y >> 16;
584  *x++ = y & 0xffff;
585 #endif
586 #endif
587  }
588  while(++i < wds);
589  if (carry) {
590  if (wds >= b->maxwds) {
591  b1 = Balloc(b->k+1);
592  Bcopy(b1, b);
593  Bfree(b);
594  b = b1;
595  }
596  b->x[wds++] = carry;
597  b->wds = wds;
598  }
599  return b;
600  }
601 
602  static Bigint *
603 s2b
604 #ifdef KR_headers
605  (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
606 #else
607  (CONST char *s, int nd0, int nd, ULong y9)
608 #endif
609 {
610  Bigint *b;
611  int i, k;
612  Long x, y;
613 
614  x = (nd + 8) / 9;
615  for(k = 0, y = 1; x > y; y <<= 1, k++) ;
616 #ifdef Pack_32
617  b = Balloc(k);
618  b->x[0] = y9;
619  b->wds = 1;
620 #else
621  b = Balloc(k+1);
622  b->x[0] = y9 & 0xffff;
623  b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
624 #endif
625 
626  i = 9;
627  if (9 < nd0) {
628  s += 9;
629  do b = multadd(b, 10, *s++ - '0');
630  while(++i < nd0);
631  s++;
632  }
633  else
634  s += 10;
635  for(; i < nd; i++)
636  b = multadd(b, 10, *s++ - '0');
637  return b;
638  }
639 
640  static int
641 hi0bits
642 #ifdef KR_headers
643  (x) register ULong x;
644 #else
645  (register ULong x)
646 #endif
647 {
648  register int k = 0;
649 
650  if (!(x & 0xffff0000)) {
651  k = 16;
652  x <<= 16;
653  }
654  if (!(x & 0xff000000)) {
655  k += 8;
656  x <<= 8;
657  }
658  if (!(x & 0xf0000000)) {
659  k += 4;
660  x <<= 4;
661  }
662  if (!(x & 0xc0000000)) {
663  k += 2;
664  x <<= 2;
665  }
666  if (!(x & 0x80000000)) {
667  k++;
668  if (!(x & 0x40000000))
669  return 32;
670  }
671  return k;
672  }
673 
674  static int
675 lo0bits
676 #ifdef KR_headers
677  (y) ULong *y;
678 #else
679  (ULong *y)
680 #endif
681 {
682  register int k;
683  register ULong x = *y;
684 
685  if (x & 7) {
686  if (x & 1)
687  return 0;
688  if (x & 2) {
689  *y = x >> 1;
690  return 1;
691  }
692  *y = x >> 2;
693  return 2;
694  }
695  k = 0;
696  if (!(x & 0xffff)) {
697  k = 16;
698  x >>= 16;
699  }
700  if (!(x & 0xff)) {
701  k += 8;
702  x >>= 8;
703  }
704  if (!(x & 0xf)) {
705  k += 4;
706  x >>= 4;
707  }
708  if (!(x & 0x3)) {
709  k += 2;
710  x >>= 2;
711  }
712  if (!(x & 1)) {
713  k++;
714  x >>= 1;
715  if (!x)
716  return 32;
717  }
718  *y = x;
719  return k;
720  }
721 
722  static Bigint *
723 i2b
724 #ifdef KR_headers
725  (i) int i;
726 #else
727  (int i)
728 #endif
729 {
730  Bigint *b;
731 
732  b = Balloc(1);
733  b->x[0] = i;
734  b->wds = 1;
735  return b;
736  }
737 
738  static Bigint *
739 mult
740 #ifdef KR_headers
741  (a, b) Bigint *a, *b;
742 #else
743  (Bigint *a, Bigint *b)
744 #endif
745 {
746  Bigint *c;
747  int k, wa, wb, wc;
748  ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
749  ULong y;
750 #ifdef ULLong
751  ULLong carry, z;
752 #else
753  ULong carry, z;
754 #ifdef Pack_32
755  ULong z2;
756 #endif
757 #endif
758 
759  if (a->wds < b->wds) {
760  c = a;
761  a = b;
762  b = c;
763  }
764  k = a->k;
765  wa = a->wds;
766  wb = b->wds;
767  wc = wa + wb;
768  if (wc > a->maxwds)
769  k++;
770  c = Balloc(k);
771  for(x = c->x, xa = x + wc; x < xa; x++)
772  *x = 0;
773  xa = a->x;
774  xae = xa + wa;
775  xb = b->x;
776  xbe = xb + wb;
777  xc0 = c->x;
778 #ifdef ULLong
779  for(; xb < xbe; xc0++) {
780  if ((y = *xb++)) {
781  x = xa;
782  xc = xc0;
783  carry = 0;
784  do {
785  z = *x++ * (ULLong)y + *xc + carry;
786  carry = z >> 32;
787  *xc++ = z & FFFFFFFF;
788  }
789  while(x < xae);
790  *xc = carry;
791  }
792  }
793 #else
794 #ifdef Pack_32
795  for(; xb < xbe; xb++, xc0++) {
796  if (y = *xb & 0xffff) {
797  x = xa;
798  xc = xc0;
799  carry = 0;
800  do {
801  z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
802  carry = z >> 16;
803  z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
804  carry = z2 >> 16;
805  Storeinc(xc, z2, z);
806  }
807  while(x < xae);
808  *xc = carry;
809  }
810  if (y = *xb >> 16) {
811  x = xa;
812  xc = xc0;
813  carry = 0;
814  z2 = *xc;
815  do {
816  z = (*x & 0xffff) * y + (*xc >> 16) + carry;
817  carry = z >> 16;
818  Storeinc(xc, z, z2);
819  z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
820  carry = z2 >> 16;
821  }
822  while(x < xae);
823  *xc = z2;
824  }
825  }
826 #else
827  for(; xb < xbe; xc0++) {
828  if (y = *xb++) {
829  x = xa;
830  xc = xc0;
831  carry = 0;
832  do {
833  z = *x++ * y + *xc + carry;
834  carry = z >> 16;
835  *xc++ = z & 0xffff;
836  }
837  while(x < xae);
838  *xc = carry;
839  }
840  }
841 #endif
842 #endif
843  for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
844  c->wds = wc;
845  return c;
846  }
847 
848  static Bigint *
849 pow5mult
850 #ifdef KR_headers
851  (b, k) Bigint *b; int k;
852 #else
853  (Bigint *b, int k)
854 #endif
855 {
856  Bigint *b1, *p5, *p51;
857  int i;
858  static int CONST p05[3] = { 5, 25, 125 };
859 
860  if ((i = k & 3))
861  b = multadd(b, p05[i-1], 0);
862 
863  if (!(k >>= 2))
864  return b;
865 
866  p5 = i2b(625);
867  for(;;) {
868  if (k & 1) {
869  b1 = mult(b, p5);
870  Bfree(b);
871  b = b1;
872  }
873  if (!(k >>= 1))
874  break;
875  p51 = mult(p5,p5);
876  Bfree(p5);
877  p5 = p51;
878  }
879  Bfree(p5);
880  return b;
881 }
882 
883  static Bigint *
884 lshift
885 #ifdef KR_headers
886  (b, k) Bigint *b; int k;
887 #else
888  (Bigint *b, int k)
889 #endif
890 {
891  int i, k1, n, n1;
892  Bigint *b1;
893  ULong *x, *x1, *xe, z;
894 
895 #ifdef Pack_32
896  n = k >> 5;
897 #else
898  n = k >> 4;
899 #endif
900  k1 = b->k;
901  n1 = n + b->wds + 1;
902  for(i = b->maxwds; n1 > i; i <<= 1)
903  k1++;
904  b1 = Balloc(k1);
905  x1 = b1->x;
906  for(i = 0; i < n; i++)
907  *x1++ = 0;
908  x = b->x;
909  xe = x + b->wds;
910 #ifdef Pack_32
911  if (k &= 0x1f) {
912  k1 = 32 - k;
913  z = 0;
914  do {
915  *x1++ = *x << k | z;
916  z = *x++ >> k1;
917  }
918  while(x < xe);
919  if ((*x1 = z))
920  ++n1;
921  }
922 #else
923  if (k &= 0xf) {
924  k1 = 16 - k;
925  z = 0;
926  do {
927  *x1++ = *x << k & 0xffff | z;
928  z = *x++ >> k1;
929  }
930  while(x < xe);
931  if (*x1 = z)
932  ++n1;
933  }
934 #endif
935  else do
936  *x1++ = *x++;
937  while(x < xe);
938  b1->wds = n1 - 1;
939  Bfree(b);
940  return b1;
941  }
942 
943  static int
944 cmp
945 #ifdef KR_headers
946  (a, b) Bigint *a, *b;
947 #else
948  (Bigint *a, Bigint *b)
949 #endif
950 {
951  ULong *xa, *xa0, *xb, *xb0;
952  int i, j;
953 
954  i = a->wds;
955  j = b->wds;
956 #ifdef DEBUG
957  if (i > 1 && !a->x[i-1])
958  Bug("cmp called with a->x[a->wds-1] == 0");
959  if (j > 1 && !b->x[j-1])
960  Bug("cmp called with b->x[b->wds-1] == 0");
961 #endif
962  if (i -= j)
963  return i;
964  xa0 = a->x;
965  xa = xa0 + j;
966  xb0 = b->x;
967  xb = xb0 + j;
968  for(;;) {
969  if (*--xa != *--xb)
970  return *xa < *xb ? -1 : 1;
971  if (xa <= xa0)
972  break;
973  }
974  return 0;
975  }
976 
977  static Bigint *
978 diff
979 #ifdef KR_headers
980  (a, b) Bigint *a, *b;
981 #else
982  (Bigint *a, Bigint *b)
983 #endif
984 {
985  Bigint *c;
986  int i, wa, wb;
987  ULong *xa, *xae, *xb, *xbe, *xc;
988 #ifdef ULLong
989  ULLong borrow, y;
990 #else
991  ULong borrow, y;
992 #ifdef Pack_32
993  ULong z;
994 #endif
995 #endif
996 
997  i = cmp(a,b);
998  if (!i) {
999  c = Balloc(0);
1000  c->wds = 1;
1001  c->x[0] = 0;
1002  return c;
1003  }
1004  if (i < 0) {
1005  c = a;
1006  a = b;
1007  b = c;
1008  i = 1;
1009  }
1010  else
1011  i = 0;
1012  c = Balloc(a->k);
1013  c->sign = i;
1014  wa = a->wds;
1015  xa = a->x;
1016  xae = xa + wa;
1017  wb = b->wds;
1018  xb = b->x;
1019  xbe = xb + wb;
1020  xc = c->x;
1021  borrow = 0;
1022 #ifdef ULLong
1023  do {
1024  y = (ULLong)*xa++ - *xb++ - borrow;
1025  borrow = y >> 32 & (ULong)1;
1026  *xc++ = y & FFFFFFFF;
1027  }
1028  while(xb < xbe);
1029  while(xa < xae) {
1030  y = *xa++ - borrow;
1031  borrow = y >> 32 & (ULong)1;
1032  *xc++ = y & FFFFFFFF;
1033  }
1034 #else
1035 #ifdef Pack_32
1036  do {
1037  y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1038  borrow = (y & 0x10000) >> 16;
1039  z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1040  borrow = (z & 0x10000) >> 16;
1041  Storeinc(xc, z, y);
1042  }
1043  while(xb < xbe);
1044  while(xa < xae) {
1045  y = (*xa & 0xffff) - borrow;
1046  borrow = (y & 0x10000) >> 16;
1047  z = (*xa++ >> 16) - borrow;
1048  borrow = (z & 0x10000) >> 16;
1049  Storeinc(xc, z, y);
1050  }
1051 #else
1052  do {
1053  y = *xa++ - *xb++ - borrow;
1054  borrow = (y & 0x10000) >> 16;
1055  *xc++ = y & 0xffff;
1056  }
1057  while(xb < xbe);
1058  while(xa < xae) {
1059  y = *xa++ - borrow;
1060  borrow = (y & 0x10000) >> 16;
1061  *xc++ = y & 0xffff;
1062  }
1063 #endif
1064 #endif
1065  while(!*--xc)
1066  wa--;
1067  c->wds = wa;
1068  return c;
1069  }
1070 
1071  static double
1072 ulp
1073 #ifdef KR_headers
1074  (x) double x;
1075 #else
1076  (double x)
1077 #endif
1078 {
1079  register Long L;
1080  U a;
1081 
1082  L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1083 #ifndef Avoid_Underflow
1084 #ifndef Sudden_Underflow
1085  if (L > 0) {
1086 #endif
1087 #endif
1088 #ifdef IBM
1089  L |= Exp_msk1 >> 4;
1090 #endif
1091  word0(a) = L;
1092  word1(a) = 0;
1093 #ifndef Avoid_Underflow
1094 #ifndef Sudden_Underflow
1095  }
1096  else {
1097  L = -L >> Exp_shift;
1098  if (L < Exp_shift) {
1099  word0(a) = 0x80000 >> L;
1100  word1(a) = 0;
1101  }
1102  else {
1103  word0(a) = 0;
1104  L -= Exp_shift;
1105  word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1106  }
1107  }
1108 #endif
1109 #endif
1110  return dval(a);
1111  }
1112 
1113  static double
1114 b2d
1115 #ifdef KR_headers
1116  (a, e) Bigint *a; int *e;
1117 #else
1118  (Bigint *a, int *e)
1119 #endif
1120 {
1121  ULong *xa, *xa0, w, y, z;
1122  int k;
1123  U d;
1124 #ifdef VAX
1125  ULong d0, d1;
1126 #else
1127 #define d0 word0(d)
1128 #define d1 word1(d)
1129 #endif
1130 
1131  xa0 = a->x;
1132  xa = xa0 + a->wds;
1133  y = *--xa;
1134 #ifdef DEBUG
1135  if (!y) Bug("zero y in b2d");
1136 #endif
1137  k = hi0bits(y);
1138  *e = 32 - k;
1139 #ifdef Pack_32
1140  if (k < Ebits) {
1141  d0 = Exp_1 | y >> (Ebits - k);
1142  w = xa > xa0 ? *--xa : 0;
1143  d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1144  goto ret_d;
1145  }
1146  z = xa > xa0 ? *--xa : 0;
1147  if (k -= Ebits) {
1148  d0 = Exp_1 | y << k | z >> (32 - k);
1149  y = xa > xa0 ? *--xa : 0;
1150  d1 = z << k | y >> (32 - k);
1151  }
1152  else {
1153  d0 = Exp_1 | y;
1154  d1 = z;
1155  }
1156 #else
1157  if (k < Ebits + 16) {
1158  z = xa > xa0 ? *--xa : 0;
1159  d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1160  w = xa > xa0 ? *--xa : 0;
1161  y = xa > xa0 ? *--xa : 0;
1162  d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1163  goto ret_d;
1164  }
1165  z = xa > xa0 ? *--xa : 0;
1166  w = xa > xa0 ? *--xa : 0;
1167  k -= Ebits + 16;
1168  d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1169  y = xa > xa0 ? *--xa : 0;
1170  d1 = w << k + 16 | y << k;
1171 #endif
1172  ret_d:
1173 #ifdef VAX
1174  word0(d) = d0 >> 16 | d0 << 16;
1175  word1(d) = d1 >> 16 | d1 << 16;
1176 #else
1177 #undef d0
1178 #undef d1
1179 #endif
1180  return dval(d);
1181  }
1182 
1183  static Bigint *
1184 d2b
1185 #ifdef KR_headers
1186  (d, e, bits) double d; int *e, *bits;
1187 #else
1188  (double _d, int *e, int *bits)
1189 #endif
1190 {
1191  Bigint *b;
1192  int de, k;
1193  ULong *x, y, z;
1194  U d;
1195 #ifndef Sudden_Underflow
1196  int i;
1197 #endif
1198 #ifdef VAX
1199  ULong d0, d1;
1200  d0 = word0(d) >> 16 | word0(d) << 16;
1201  d1 = word1(d) >> 16 | word1(d) << 16;
1202 #else
1203 #define d0 word0(d)
1204 #define d1 word1(d)
1205 #endif
1206  dval(d) = _d;
1207 
1208 #ifdef Pack_32
1209  b = Balloc(1);
1210 #else
1211  b = Balloc(2);
1212 #endif
1213  x = b->x;
1214 
1215  z = d0 & Frac_mask;
1216  d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1217 #ifdef Sudden_Underflow
1218  de = (int)(d0 >> Exp_shift);
1219 #ifndef IBM
1220  z |= Exp_msk11;
1221 #endif
1222 #else
1223  if ((de = (int)(d0 >> Exp_shift)))
1224  z |= Exp_msk1;
1225 #endif
1226 #ifdef Pack_32
1227  if ((y = d1)) {
1228  if ((k = lo0bits(&y))) {
1229  x[0] = y | z << (32 - k);
1230  z >>= k;
1231  }
1232  else
1233  x[0] = y;
1234 #ifndef Sudden_Underflow
1235  i =
1236 #endif
1237  b->wds = (x[1] = z) ? 2 : 1;
1238  }
1239  else {
1240 #ifdef DEBUG
1241  if (!z)
1242  Bug("Zero passed to d2b");
1243 #endif
1244  k = lo0bits(&z);
1245  x[0] = z;
1246 #ifndef Sudden_Underflow
1247  i =
1248 #endif
1249  b->wds = 1;
1250  k += 32;
1251  }
1252 #else
1253  if (y = d1) {
1254  if (k = lo0bits(&y))
1255  if (k >= 16) {
1256  x[0] = y | z << 32 - k & 0xffff;
1257  x[1] = z >> k - 16 & 0xffff;
1258  x[2] = z >> k;
1259  i = 2;
1260  }
1261  else {
1262  x[0] = y & 0xffff;
1263  x[1] = y >> 16 | z << 16 - k & 0xffff;
1264  x[2] = z >> k & 0xffff;
1265  x[3] = z >> k+16;
1266  i = 3;
1267  }
1268  else {
1269  x[0] = y & 0xffff;
1270  x[1] = y >> 16;
1271  x[2] = z & 0xffff;
1272  x[3] = z >> 16;
1273  i = 3;
1274  }
1275  }
1276  else {
1277 #ifdef DEBUG
1278  if (!z)
1279  Bug("Zero passed to d2b");
1280 #endif
1281  k = lo0bits(&z);
1282  if (k >= 16) {
1283  x[0] = z;
1284  i = 0;
1285  }
1286  else {
1287  x[0] = z & 0xffff;
1288  x[1] = z >> 16;
1289  i = 1;
1290  }
1291  k += 32;
1292  }
1293  while(!x[i])
1294  --i;
1295  b->wds = i + 1;
1296 #endif
1297 #ifndef Sudden_Underflow
1298  if (de) {
1299 #endif
1300 #ifdef IBM
1301  *e = (de - Bias - (P-1) << 2) + k;
1302  *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1303 #else
1304  *e = de - Bias - (P-1) + k;
1305  *bits = P - k;
1306 #endif
1307 #ifndef Sudden_Underflow
1308  }
1309  else {
1310  *e = de - Bias - (P-1) + 1 + k;
1311 #ifdef Pack_32
1312  *bits = 32*i - hi0bits(x[i-1]);
1313 #else
1314  *bits = (i+2)*16 - hi0bits(x[i]);
1315 #endif
1316  }
1317 #endif
1318  return b;
1319  }
1320 #undef d0
1321 #undef d1
1322 
1323  static double
1324 ratio
1325 #ifdef KR_headers
1326  (a, b) Bigint *a, *b;
1327 #else
1328  (Bigint *a, Bigint *b)
1329 #endif
1330 {
1331  U da, db;
1332  int k, ka, kb;
1333 
1334  dval(da) = b2d(a, &ka);
1335  dval(db) = b2d(b, &kb);
1336 #ifdef Pack_32
1337  k = ka - kb + 32*(a->wds - b->wds);
1338 #else
1339  k = ka - kb + 16*(a->wds - b->wds);
1340 #endif
1341 #ifdef IBM
1342  if (k > 0) {
1343  word0(da) += (k >> 2)*Exp_msk1;
1344  if (k &= 3)
1345  dval(da) *= 1 << k;
1346  }
1347  else {
1348  k = -k;
1349  word0(db) += (k >> 2)*Exp_msk1;
1350  if (k &= 3)
1351  dval(db) *= 1 << k;
1352  }
1353 #else
1354  if (k > 0)
1355  word0(da) += k*Exp_msk1;
1356  else {
1357  k = -k;
1358  word0(db) += k*Exp_msk1;
1359  }
1360 #endif
1361  return dval(da) / dval(db);
1362  }
1363 
1364  static CONST double
1365 tens[] = {
1366  1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1367  1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1368  1e20, 1e21, 1e22
1369 #ifdef VAX
1370  , 1e23, 1e24
1371 #endif
1372  };
1373 
1374  static CONST double
1375 #ifdef IEEE_Arith
1376 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1377 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1378 #ifdef Avoid_Underflow
1379  9007199254740992.*9007199254740992.e-256
1380  /* = 2^106 * 1e-53 */
1381 #else
1382  1e-256
1383 #endif
1384  };
1385 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1386 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1387 #define Scale_Bit 0x10
1388 #define n_bigtens 5
1389 #else
1390 #ifdef IBM
1391 bigtens[] = { 1e16, 1e32, 1e64 };
1392 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1393 #define n_bigtens 3
1394 #else
1395 bigtens[] = { 1e16, 1e32 };
1396 static CONST double tinytens[] = { 1e-16, 1e-32 };
1397 #define n_bigtens 2
1398 #endif
1399 #endif
1400 
1401 #ifdef INFNAN_CHECK
1402 
1403 #ifndef NAN_WORD0
1404 #define NAN_WORD0 0x7ff80000
1405 #endif
1406 
1407 #ifndef NAN_WORD1
1408 #define NAN_WORD1 0
1409 #endif
1410 
1411  static int
1412 match
1413 #ifdef KR_headers
1414  (sp, t) char **sp, *t;
1415 #else
1416  (CONST char **sp, char *t)
1417 #endif
1418 {
1419  int c, d;
1420  CONST char *s = *sp;
1421 
1422  while((d = *t++)) {
1423  if ((c = *++s) >= 'A' && c <= 'Z')
1424  c += 'a' - 'A';
1425  if (c != d)
1426  return 0;
1427  }
1428  *sp = s + 1;
1429  return 1;
1430  }
1431 
1432 #ifndef No_Hex_NaN
1433  static void
1434 hexnan
1435 #ifdef KR_headers
1436  (rvp, sp) double *rvp; CONST char **sp;
1437 #else
1438  (U *rvp, CONST char **sp)
1439 #endif
1440 {
1441  ULong c, x[2];
1442  CONST char *s;
1443  int havedig, udx0, xshift;
1444 
1445  x[0] = x[1] = 0;
1446  havedig = xshift = 0;
1447  udx0 = 1;
1448  s = *sp;
1449  /* allow optional initial 0x or 0X */
1450  while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
1451  ++s;
1452  if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
1453  s += 2;
1454  while((c = *(CONST unsigned char*)++s)) {
1455  if (c >= '0' && c <= '9')
1456  c -= '0';
1457  else if (c >= 'a' && c <= 'f')
1458  c += 10 - 'a';
1459  else if (c >= 'A' && c <= 'F')
1460  c += 10 - 'A';
1461  else if (c <= ' ') {
1462  if (udx0 && havedig) {
1463  udx0 = 0;
1464  xshift = 1;
1465  }
1466  continue;
1467  }
1468 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1469  else if (/*(*/ c == ')' && havedig) {
1470  *sp = s + 1;
1471  break;
1472  }
1473  else
1474  return; /* invalid form: don't change *sp */
1475 #else
1476  else {
1477  do {
1478  if (/*(*/ c == ')') {
1479  *sp = s + 1;
1480  break;
1481  }
1482  } while((c = *++s));
1483  break;
1484  }
1485 #endif
1486  havedig = 1;
1487  if (xshift) {
1488  xshift = 0;
1489  x[0] = x[1];
1490  x[1] = 0;
1491  }
1492  if (udx0)
1493  x[0] = (x[0] << 4) | (x[1] >> 28);
1494  x[1] = (x[1] << 4) | c;
1495  }
1496  if ((x[0] &= 0xfffff) || x[1]) {
1497  word0(*rvp) = Exp_mask | x[0];
1498  word1(*rvp) = x[1];
1499  }
1500  }
1501 #endif /*No_Hex_NaN*/
1502 #endif /* INFNAN_CHECK */
1503 
1504  double
1505 sb_strtod
1506 #ifdef KR_headers
1507  (s00, se) CONST char *s00; char **se;
1508 #else
1509  (CONST char *s00, char **se)
1510 #endif
1511 {
1512 #ifdef Avoid_Underflow
1513  int scale;
1514 #endif
1515  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1516  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1517  CONST char *s, *s0, *s1;
1518  double aadj, adj;
1519  U rv, rv0, aadj1;
1520  Long L;
1521  ULong y, z;
1522  Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1523 #ifdef SET_INEXACT
1524  int inexact, oldinexact;
1525 #endif
1526 #ifdef Honor_FLT_ROUNDS
1527  int rounding;
1528 #endif
1529 #ifdef USE_LOCALE
1530  CONST char *s2;
1531 #endif
1532 
1533  sign = nz0 = nz = 0;
1534  dval(rv) = 0.;
1535  for(s = s00;;s++) switch(*s) {
1536  case '-':
1537  sign = 1;
1538  /* no break */
1539  case '+':
1540  if (*++s)
1541  goto break2;
1542  /* no break */
1543  case 0:
1544  goto ret0;
1545  case '\t':
1546  case '\n':
1547  case '\v':
1548  case '\f':
1549  case '\r':
1550  case ' ':
1551  continue;
1552  default:
1553  goto break2;
1554  }
1555  break2:
1556  if (*s == '0') {
1557  nz0 = 1;
1558  while(*++s == '0') ;
1559  if (!*s)
1560  goto ret;
1561  }
1562  s0 = s;
1563  y = z = 0;
1564  for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1565  if (nd < 9)
1566  y = 10*y + c - '0';
1567  else if (nd < 16)
1568  z = 10*z + c - '0';
1569  nd0 = nd;
1570 #ifdef USE_LOCALE
1571  s1 = localeconv()->decimal_point;
1572  if (c == *s1) {
1573  c = '.';
1574  if (*++s1) {
1575  s2 = s;
1576  for(;;) {
1577  if (*++s2 != *s1) {
1578  c = 0;
1579  break;
1580  }
1581  if (!*++s1) {
1582  s = s2;
1583  break;
1584  }
1585  }
1586  }
1587  }
1588 #endif
1589  if (c == '.') {
1590  c = *++s;
1591  if (!nd) {
1592  for(; c == '0'; c = *++s)
1593  nz++;
1594  if (c > '0' && c <= '9') {
1595  s0 = s;
1596  nf += nz;
1597  nz = 0;
1598  goto have_dig;
1599  }
1600  goto dig_done;
1601  }
1602  for(; c >= '0' && c <= '9'; c = *++s) {
1603  have_dig:
1604  nz++;
1605  if (c -= '0') {
1606  nf += nz;
1607  for(i = 1; i < nz; i++)
1608  if (nd++ < 9)
1609  y *= 10;
1610  else if (nd <= DBL_DIG + 1)
1611  z *= 10;
1612  if (nd++ < 9)
1613  y = 10*y + c;
1614  else if (nd <= DBL_DIG + 1)
1615  z = 10*z + c;
1616  nz = 0;
1617  }
1618  }
1619  }
1620  dig_done:
1621  e = 0;
1622  if (c == 'e' || c == 'E') {
1623  if (!nd && !nz && !nz0) {
1624  goto ret0;
1625  }
1626  s00 = s;
1627  esign = 0;
1628  switch(c = *++s) {
1629  case '-':
1630  esign = 1;
1631  case '+':
1632  c = *++s;
1633  }
1634  if (c >= '0' && c <= '9') {
1635  while(c == '0')
1636  c = *++s;
1637  if (c > '0' && c <= '9') {
1638  L = c - '0';
1639  s1 = s;
1640  while((c = *++s) >= '0' && c <= '9')
1641  L = 10*L + c - '0';
1642  if (s - s1 > 8 || L > 19999)
1643  /* Avoid confusion from exponents
1644  * so large that e might overflow.
1645  */
1646  e = 19999; /* safe for 16 bit ints */
1647  else
1648  e = (int)L;
1649  if (esign)
1650  e = -e;
1651  }
1652  else
1653  e = 0;
1654  }
1655  else
1656  s = s00;
1657  }
1658  if (!nd) {
1659  if (!nz && !nz0) {
1660 #ifdef INFNAN_CHECK
1661  /* Check for Nan and Infinity */
1662  switch(c) {
1663  case 'i':
1664  case 'I':
1665  if (match(&s,"nf")) {
1666  --s;
1667  if (!match(&s,"inity"))
1668  ++s;
1669  word0(rv) = 0x7ff00000;
1670  word1(rv) = 0;
1671  goto ret;
1672  }
1673  break;
1674  case 'n':
1675  case 'N':
1676  if (match(&s, "an")) {
1677  word0(rv) = NAN_WORD0;
1678  word1(rv) = NAN_WORD1;
1679 #ifndef No_Hex_NaN
1680  if (*s == '(') /*)*/
1681  hexnan(&rv, &s);
1682 #endif
1683  goto ret;
1684  }
1685  }
1686 #endif /* INFNAN_CHECK */
1687  ret0:
1688  s = s00;
1689  sign = 0;
1690  }
1691  goto ret;
1692  }
1693  e1 = e -= nf;
1694 
1695  /* Now we have nd0 digits, starting at s0, followed by a
1696  * decimal point, followed by nd-nd0 digits. The number we're
1697  * after is the integer represented by those digits times
1698  * 10**e */
1699 
1700  if (!nd0)
1701  nd0 = nd;
1702  k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1703  dval(rv) = y;
1704  if (k > 9) {
1705 #ifdef SET_INEXACT
1706  if (k > DBL_DIG)
1707  oldinexact = get_inexact();
1708 #endif
1709  dval(rv) = tens[k - 9] * dval(rv) + z;
1710  }
1711  bd0 = 0;
1712  if (nd <= DBL_DIG
1713 #ifndef RND_PRODQUOT
1714 #ifndef Honor_FLT_ROUNDS
1715  && Flt_Rounds == 1
1716 #endif
1717 #endif
1718  ) {
1719  if (!e)
1720  goto ret;
1721  if (e > 0) {
1722  if (e <= Ten_pmax) {
1723 #ifdef VAX
1724  goto vax_ovfl_check;
1725 #else
1726 #ifdef Honor_FLT_ROUNDS
1727  /* round correctly FLT_ROUNDS = 2 or 3 */
1728  if (sign) {
1729  rv = -rv;
1730  sign = 0;
1731  }
1732 #endif
1733  /* rv = */ rounded_product(dval(rv), tens[e]);
1734  goto ret;
1735 #endif
1736  }
1737  i = DBL_DIG - nd;
1738  if (e <= Ten_pmax + i) {
1739  /* A fancier test would sometimes let us do
1740  * this for larger i values.
1741  */
1742 #ifdef Honor_FLT_ROUNDS
1743  /* round correctly FLT_ROUNDS = 2 or 3 */
1744  if (sign) {
1745  rv = -rv;
1746  sign = 0;
1747  }
1748 #endif
1749  e -= i;
1750  dval(rv) *= tens[i];
1751 #ifdef VAX
1752  /* VAX exponent range is so narrow we must
1753  * worry about overflow here...
1754  */
1755  vax_ovfl_check:
1756  word0(rv) -= P*Exp_msk1;
1757  /* rv = */ rounded_product(dval(rv), tens[e]);
1758  if ((word0(rv) & Exp_mask)
1759  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1760  goto ovfl;
1761  word0(rv) += P*Exp_msk1;
1762 #else
1763  /* rv = */ rounded_product(dval(rv), tens[e]);
1764 #endif
1765  goto ret;
1766  }
1767  }
1768 #ifndef Inaccurate_Divide
1769  else if (e >= -Ten_pmax) {
1770 #ifdef Honor_FLT_ROUNDS
1771  /* round correctly FLT_ROUNDS = 2 or 3 */
1772  if (sign) {
1773  rv = -rv;
1774  sign = 0;
1775  }
1776 #endif
1777  /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1778  goto ret;
1779  }
1780 #endif
1781  }
1782  e1 += nd - k;
1783 
1784 #ifdef IEEE_Arith
1785 #ifdef SET_INEXACT
1786  inexact = 1;
1787  if (k <= DBL_DIG)
1788  oldinexact = get_inexact();
1789 #endif
1790 #ifdef Avoid_Underflow
1791  scale = 0;
1792 #endif
1793 #ifdef Honor_FLT_ROUNDS
1794  if ((rounding = Flt_Rounds) >= 2) {
1795  if (sign)
1796  rounding = rounding == 2 ? 0 : 2;
1797  else
1798  if (rounding != 2)
1799  rounding = 0;
1800  }
1801 #endif
1802 #endif /*IEEE_Arith*/
1803 
1804  /* Get starting approximation = rv * 10**e1 */
1805 
1806  if (e1 > 0) {
1807  if ((i = e1 & 15))
1808  dval(rv) *= tens[i];
1809  if (e1 &= ~15) {
1810  if (e1 > DBL_MAX_10_EXP) {
1811  ovfl:
1812 #ifndef NO_ERRNO
1813  errno = ERANGE;
1814 #endif
1815  /* Can't trust HUGE_VAL */
1816 #ifdef IEEE_Arith
1817 #ifdef Honor_FLT_ROUNDS
1818  switch(rounding) {
1819  case 0: /* toward 0 */
1820  case 3: /* toward -infinity */
1821  word0(rv) = Big0;
1822  word1(rv) = Big1;
1823  break;
1824  default:
1825  word0(rv) = Exp_mask;
1826  word1(rv) = 0;
1827  }
1828 #else /*Honor_FLT_ROUNDS*/
1829  word0(rv) = Exp_mask;
1830  word1(rv) = 0;
1831 #endif /*Honor_FLT_ROUNDS*/
1832 #ifdef SET_INEXACT
1833  /* set overflow bit */
1834  dval(rv0) = 1e300;
1835  dval(rv0) *= dval(rv0);
1836 #endif
1837 #else /*IEEE_Arith*/
1838  word0(rv) = Big0;
1839  word1(rv) = Big1;
1840 #endif /*IEEE_Arith*/
1841  if (bd0)
1842  goto retfree;
1843  goto ret;
1844  }
1845  e1 >>= 4;
1846  for(j = 0; e1 > 1; j++, e1 >>= 1)
1847  if (e1 & 1)
1848  dval(rv) *= bigtens[j];
1849  /* The last multiplication could overflow. */
1850  word0(rv) -= P*Exp_msk1;
1851  dval(rv) *= bigtens[j];
1852  if ((z = word0(rv) & Exp_mask)
1853  > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1854  goto ovfl;
1855  if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1856  /* set to largest number */
1857  /* (Can't trust DBL_MAX) */
1858  word0(rv) = Big0;
1859  word1(rv) = Big1;
1860  }
1861  else
1862  word0(rv) += P*Exp_msk1;
1863  }
1864  }
1865  else if (e1 < 0) {
1866  e1 = -e1;
1867  if ((i = e1 & 15))
1868  dval(rv) /= tens[i];
1869  if (e1 >>= 4) {
1870  if (e1 >= 1 << n_bigtens)
1871  goto undfl;
1872 #ifdef Avoid_Underflow
1873  if (e1 & Scale_Bit)
1874  scale = 2*P;
1875  for(j = 0; e1 > 0; j++, e1 >>= 1)
1876  if (e1 & 1)
1877  dval(rv) *= tinytens[j];
1878  if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1879  >> Exp_shift)) > 0) {
1880  /* scaled rv is denormal; zap j low bits */
1881  if (j >= 32) {
1882  word1(rv) = 0;
1883  if (j >= 53)
1884  word0(rv) = (P+2)*Exp_msk1;
1885  else
1886  word0(rv) &= 0xffffffff << (j-32);
1887  }
1888  else
1889  word1(rv) &= 0xffffffff << j;
1890  }
1891 #else
1892  for(j = 0; e1 > 1; j++, e1 >>= 1)
1893  if (e1 & 1)
1894  dval(rv) *= tinytens[j];
1895  /* The last multiplication could underflow. */
1896  dval(rv0) = dval(rv);
1897  dval(rv) *= tinytens[j];
1898  if (!dval(rv)) {
1899  dval(rv) = 2.*dval(rv0);
1900  dval(rv) *= tinytens[j];
1901 #endif
1902  if (!dval(rv)) {
1903  undfl:
1904  dval(rv) = 0.;
1905 #ifndef NO_ERRNO
1906  errno = ERANGE;
1907 #endif
1908  if (bd0)
1909  goto retfree;
1910  goto ret;
1911  }
1912 #ifndef Avoid_Underflow
1913  word0(rv) = Tiny0;
1914  word1(rv) = Tiny1;
1915  /* The refinement below will clean
1916  * this approximation up.
1917  */
1918  }
1919 #endif
1920  }
1921  }
1922 
1923  /* Now the hard part -- adjusting rv to the correct value.*/
1924 
1925  /* Put digits into bd: true value = bd * 10^e */
1926 
1927  bd0 = s2b(s0, nd0, nd, y);
1928 
1929  for(;;) {
1930  bd = Balloc(bd0->k);
1931  Bcopy(bd, bd0);
1932  bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
1933  bs = i2b(1);
1934 
1935  if (e >= 0) {
1936  bb2 = bb5 = 0;
1937  bd2 = bd5 = e;
1938  }
1939  else {
1940  bb2 = bb5 = -e;
1941  bd2 = bd5 = 0;
1942  }
1943  if (bbe >= 0)
1944  bb2 += bbe;
1945  else
1946  bd2 -= bbe;
1947  bs2 = bb2;
1948 #ifdef Honor_FLT_ROUNDS
1949  if (rounding != 1)
1950  bs2++;
1951 #endif
1952 #ifdef Avoid_Underflow
1953  j = bbe - scale;
1954  i = j + bbbits - 1; /* logb(rv) */
1955  if (i < Emin) /* denormal */
1956  j += P - Emin;
1957  else
1958  j = P + 1 - bbbits;
1959 #else /*Avoid_Underflow*/
1960 #ifdef Sudden_Underflow
1961 #ifdef IBM
1962  j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1963 #else
1964  j = P + 1 - bbbits;
1965 #endif
1966 #else /*Sudden_Underflow*/
1967  j = bbe;
1968  i = j + bbbits - 1; /* logb(rv) */
1969  if (i < Emin) /* denormal */
1970  j += P - Emin;
1971  else
1972  j = P + 1 - bbbits;
1973 #endif /*Sudden_Underflow*/
1974 #endif /*Avoid_Underflow*/
1975  bb2 += j;
1976  bd2 += j;
1977 #ifdef Avoid_Underflow
1978  bd2 += scale;
1979 #endif
1980  i = bb2 < bd2 ? bb2 : bd2;
1981  if (i > bs2)
1982  i = bs2;
1983  if (i > 0) {
1984  bb2 -= i;
1985  bd2 -= i;
1986  bs2 -= i;
1987  }
1988  if (bb5 > 0) {
1989  bs = pow5mult(bs, bb5);
1990  bb1 = mult(bs, bb);
1991  Bfree(bb);
1992  bb = bb1;
1993  }
1994  if (bb2 > 0)
1995  bb = lshift(bb, bb2);
1996  if (bd5 > 0)
1997  bd = pow5mult(bd, bd5);
1998  if (bd2 > 0)
1999  bd = lshift(bd, bd2);
2000  if (bs2 > 0)
2001  bs = lshift(bs, bs2);
2002  delta = diff(bb, bd);
2003  dsign = delta->sign;
2004  delta->sign = 0;
2005  i = cmp(delta, bs);
2006 #ifdef Honor_FLT_ROUNDS
2007  if (rounding != 1) {
2008  if (i < 0) {
2009  /* Error is less than an ulp */
2010  if (!delta->x[0] && delta->wds <= 1) {
2011  /* exact */
2012 #ifdef SET_INEXACT
2013  inexact = 0;
2014 #endif
2015  break;
2016  }
2017  if (rounding) {
2018  if (dsign) {
2019  adj = 1.;
2020  goto apply_adj;
2021  }
2022  }
2023  else if (!dsign) {
2024  adj = -1.;
2025  if (!word1(rv)
2026  && !(word0(rv) & Frac_mask)) {
2027  y = word0(rv) & Exp_mask;
2028 #ifdef Avoid_Underflow
2029  if (!scale || y > 2*P*Exp_msk1)
2030 #else
2031  if (y)
2032 #endif
2033  {
2034  delta = lshift(delta,Log2P);
2035  if (cmp(delta, bs) <= 0)
2036  adj = -0.5;
2037  }
2038  }
2039  apply_adj:
2040 #ifdef Avoid_Underflow
2041  if (scale && (y = word0(rv) & Exp_mask)
2042  <= 2*P*Exp_msk1)
2043  word0(adj) += (2*P+1)*Exp_msk1 - y;
2044 #else
2045 #ifdef Sudden_Underflow
2046  if ((word0(rv) & Exp_mask) <=
2047  P*Exp_msk1) {
2048  word0(rv) += P*Exp_msk1;
2049  dval(rv) += adj*ulp(dval(rv));
2050  word0(rv) -= P*Exp_msk1;
2051  }
2052  else
2053 #endif /*Sudden_Underflow*/
2054 #endif /*Avoid_Underflow*/
2055  dval(rv) += adj*ulp(dval(rv));
2056  }
2057  break;
2058  }
2059  adj = ratio(delta, bs);
2060  if (adj < 1.)
2061  adj = 1.;
2062  if (adj <= 0x7ffffffe) {
2063  /* adj = rounding ? ceil(adj) : floor(adj); */
2064  y = adj;
2065  if (y != adj) {
2066  if (!((rounding>>1) ^ dsign))
2067  y++;
2068  adj = y;
2069  }
2070  }
2071 #ifdef Avoid_Underflow
2072  if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2073  word0(adj) += (2*P+1)*Exp_msk1 - y;
2074 #else
2075 #ifdef Sudden_Underflow
2076  if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2077  word0(rv) += P*Exp_msk1;
2078  adj *= ulp(dval(rv));
2079  if (dsign)
2080  dval(rv) += adj;
2081  else
2082  dval(rv) -= adj;
2083  word0(rv) -= P*Exp_msk1;
2084  goto cont;
2085  }
2086 #endif /*Sudden_Underflow*/
2087 #endif /*Avoid_Underflow*/
2088  adj *= ulp(dval(rv));
2089  if (dsign)
2090  dval(rv) += adj;
2091  else
2092  dval(rv) -= adj;
2093  goto cont;
2094  }
2095 #endif /*Honor_FLT_ROUNDS*/
2096 
2097  if (i < 0) {
2098  /* Error is less than half an ulp -- check for
2099  * special case of mantissa a power of two.
2100  */
2101  if (dsign || word1(rv) || word0(rv) & Bndry_mask
2102 #ifdef IEEE_Arith
2103 #ifdef Avoid_Underflow
2104  || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2105 #else
2106  || (word0(rv) & Exp_mask) <= Exp_msk1
2107 #endif
2108 #endif
2109  ) {
2110 #ifdef SET_INEXACT
2111  if (!delta->x[0] && delta->wds <= 1)
2112  inexact = 0;
2113 #endif
2114  break;
2115  }
2116  if (!delta->x[0] && delta->wds <= 1) {
2117  /* exact result */
2118 #ifdef SET_INEXACT
2119  inexact = 0;
2120 #endif
2121  break;
2122  }
2123  delta = lshift(delta,Log2P);
2124  if (cmp(delta, bs) > 0)
2125  goto drop_down;
2126  break;
2127  }
2128  if (i == 0) {
2129  /* exactly half-way between */
2130  if (dsign) {
2131  if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2132  && word1(rv) == (
2133 #ifdef Avoid_Underflow
2134  (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2135  ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2136 #endif
2137  0xffffffff)) {
2138  /*boundary case -- increment exponent*/
2139  word0(rv) = (word0(rv) & Exp_mask)
2140  + Exp_msk1
2141 #ifdef IBM
2142  | Exp_msk1 >> 4
2143 #endif
2144  ;
2145  word1(rv) = 0;
2146 #ifdef Avoid_Underflow
2147  dsign = 0;
2148 #endif
2149  break;
2150  }
2151  }
2152  else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2153  drop_down:
2154  /* boundary case -- decrement exponent */
2155 #ifdef Sudden_Underflow /*{{*/
2156  L = word0(rv) & Exp_mask;
2157 #ifdef IBM
2158  if (L < Exp_msk1)
2159 #else
2160 #ifdef Avoid_Underflow
2161  if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2162 #else
2163  if (L <= Exp_msk1)
2164 #endif /*Avoid_Underflow*/
2165 #endif /*IBM*/
2166  goto undfl;
2167  L -= Exp_msk1;
2168 #else /*Sudden_Underflow}{*/
2169 #ifdef Avoid_Underflow
2170  if (scale) {
2171  L = word0(rv) & Exp_mask;
2172  if (L <= (2*P+1)*Exp_msk1) {
2173  if (L > (P+2)*Exp_msk1)
2174  /* round even ==> */
2175  /* accept rv */
2176  break;
2177  /* rv = smallest denormal */
2178  goto undfl;
2179  }
2180  }
2181 #endif /*Avoid_Underflow*/
2182  L = (word0(rv) & Exp_mask) - Exp_msk1;
2183 #endif /*Sudden_Underflow}}*/
2184  word0(rv) = L | Bndry_mask1;
2185  word1(rv) = 0xffffffff;
2186 #ifdef IBM
2187  goto cont;
2188 #else
2189  break;
2190 #endif
2191  }
2192 #ifndef ROUND_BIASED
2193  if (!(word1(rv) & LSB))
2194  break;
2195 #endif
2196  if (dsign)
2197  dval(rv) += ulp(dval(rv));
2198 #ifndef ROUND_BIASED
2199  else {
2200  dval(rv) -= ulp(dval(rv));
2201 #ifndef Sudden_Underflow
2202  if (!dval(rv))
2203  goto undfl;
2204 #endif
2205  }
2206 #ifdef Avoid_Underflow
2207  dsign = 1 - dsign;
2208 #endif
2209 #endif
2210  break;
2211  }
2212  if ((aadj = ratio(delta, bs)) <= 2.) {
2213  if (dsign)
2214  aadj = dval(aadj1) = 1.;
2215  else if (word1(rv) || word0(rv) & Bndry_mask) {
2216 #ifndef Sudden_Underflow
2217  if (word1(rv) == Tiny1 && !word0(rv))
2218  goto undfl;
2219 #endif
2220  aadj = 1.;
2221  dval(aadj1) = -1.;
2222  }
2223  else {
2224  /* special case -- power of FLT_RADIX to be */
2225  /* rounded down... */
2226 
2227  if (aadj < 2./FLT_RADIX)
2228  aadj = 1./FLT_RADIX;
2229  else
2230  aadj *= 0.5;
2231  dval(aadj1) = -aadj;
2232  }
2233  }
2234  else {
2235  aadj *= 0.5;
2236  dval(aadj1) = dsign ? aadj : -aadj;
2237 #ifdef Check_FLT_ROUNDS
2238  switch(Rounding) {
2239  case 2: /* towards +infinity */
2240  dval(aadj1) -= 0.5;
2241  break;
2242  case 0: /* towards 0 */
2243  case 3: /* towards -infinity */
2244  dval(aadj1) += 0.5;
2245  }
2246 #else
2247  if (Flt_Rounds == 0)
2248  dval(aadj1) += 0.5;
2249 #endif /*Check_FLT_ROUNDS*/
2250  }
2251  y = word0(rv) & Exp_mask;
2252 
2253  /* Check for overflow */
2254 
2255  if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2256  dval(rv0) = dval(rv);
2257  word0(rv) -= P*Exp_msk1;
2258  adj = dval(aadj1) * ulp(dval(rv));
2259  dval(rv) += adj;
2260  if ((word0(rv) & Exp_mask) >=
2261  Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2262  if (word0(rv0) == Big0 && word1(rv0) == Big1)
2263  goto ovfl;
2264  word0(rv) = Big0;
2265  word1(rv) = Big1;
2266  goto cont;
2267  }
2268  else
2269  word0(rv) += P*Exp_msk1;
2270  }
2271  else {
2272 #ifdef Avoid_Underflow
2273  if (scale && y <= 2*P*Exp_msk1) {
2274  if (aadj <= 0x7fffffff) {
2275  if ((z = (uint32)aadj) <= 0)
2276  z = 1;
2277  aadj = z;
2278  dval(aadj1) = dsign ? aadj : -aadj;
2279  }
2280  word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2281  }
2282  adj = dval(aadj1) * ulp(dval(rv));
2283  dval(rv) += adj;
2284 #else
2285 #ifdef Sudden_Underflow
2286  if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2287  dval(rv0) = dval(rv);
2288  word0(rv) += P*Exp_msk1;
2289  adj = aadj1 * ulp(dval(rv));
2290  dval(rv) += adj;
2291 #ifdef IBM
2292  if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2293 #else
2294  if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2295 #endif
2296  {
2297  if (word0(rv0) == Tiny0
2298  && word1(rv0) == Tiny1)
2299  goto undfl;
2300  word0(rv) = Tiny0;
2301  word1(rv) = Tiny1;
2302  goto cont;
2303  }
2304  else
2305  word0(rv) -= P*Exp_msk1;
2306  }
2307  else {
2308  adj = aadj1 * ulp(dval(rv));
2309  dval(rv) += adj;
2310  }
2311 #else /*Sudden_Underflow*/
2312  /* Compute adj so that the IEEE rounding rules will
2313  * correctly round rv + adj in some half-way cases.
2314  * If rv * ulp(rv) is denormalized (i.e.,
2315  * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2316  * trouble from bits lost to denormalization;
2317  * example: 1.2e-307 .
2318  */
2319  if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2320  aadj1 = (double)(int)(aadj + 0.5);
2321  if (!dsign)
2322  aadj1 = -aadj1;
2323  }
2324  adj = aadj1 * ulp(dval(rv));
2325  dval(rv) += adj;
2326 #endif /*Sudden_Underflow*/
2327 #endif /*Avoid_Underflow*/
2328  }
2329  z = word0(rv) & Exp_mask;
2330 #ifndef SET_INEXACT
2331 #ifdef Avoid_Underflow
2332  if (!scale)
2333 #endif
2334  if (y == z) {
2335  /* Can we stop now? */
2336  L = (Long)aadj;
2337  aadj -= L;
2338  /* The tolerances below are conservative. */
2339  if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2340  if (aadj < .4999999 || aadj > .5000001)
2341  break;
2342  }
2343  else if (aadj < .4999999/FLT_RADIX)
2344  break;
2345  }
2346 #endif
2347  cont:
2348  Bfree(bb);
2349  Bfree(bd);
2350  Bfree(bs);
2351  Bfree(delta);
2352  }
2353 #ifdef SET_INEXACT
2354  if (inexact) {
2355  if (!oldinexact) {
2356  word0(rv0) = Exp_1 + (70 << Exp_shift);
2357  word1(rv0) = 0;
2358  dval(rv0) += 1.;
2359  }
2360  }
2361  else if (!oldinexact)
2362  clear_inexact();
2363 #endif
2364 #ifdef Avoid_Underflow
2365  if (scale) {
2366  word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2367  word1(rv0) = 0;
2368  dval(rv) *= dval(rv0);
2369 #ifndef NO_ERRNO
2370  /* try to avoid the bug of testing an 8087 register value */
2371  if (word0(rv) == 0 && word1(rv) == 0)
2372  errno = ERANGE;
2373 #endif
2374  }
2375 #endif /* Avoid_Underflow */
2376 #ifdef SET_INEXACT
2377  if (inexact && !(word0(rv) & Exp_mask)) {
2378  /* set underflow bit */
2379  dval(rv0) = 1e-300;
2380  dval(rv0) *= dval(rv0);
2381  }
2382 #endif
2383  retfree:
2384  Bfree(bb);
2385  Bfree(bd);
2386  Bfree(bs);
2387  Bfree(bd0);
2388  Bfree(delta);
2389  ret:
2390  if (se)
2391  *se = (char *)s;
2392  return sign ? -dval(rv) : dval(rv);
2393  }
2394 
2395 #ifdef __cplusplus
2396 }
2397 #endif
Sphinx's memory allocation/deallocation routines.
Basic type definitions used in Sphinx.
SPHINXBASE_EXPORT void ckd_free(void *ptr)
Test and free a 1-D array.
Definition: ckd_alloc.c:241
Definition: dtoa.c:500
#define ckd_malloc(sz)
Macro for ckd_malloc
Definition: ckd_alloc.h:253
Union to extract the bytes of a double.
Definition: dtoa.c:291