1 /*-
2 * Copyright (c) 2004 - 2011 CTPP Team
3 *
4 * Redistribution and use in source and binary forms, with or without
5 * modification, are permitted provided that the following conditions
6 * are met:
7 * 1. Redistributions of source code must retain the above copyright
8 * notice, this list of conditions and the following disclaimer.
9 * 2. Redistributions in binary form must reproduce the above copyright
10 * notice, this list of conditions and the following disclaimer in the
11 * documentation and/or other materials provided with the distribution.
12 * 4. Neither the name of the CTPP Team nor the names of its contributors
13 * may be used to endorse or promote products derived from this software
14 * without specific prior written permission.
15 *
16 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
20 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
21 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
22 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
23 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
24 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
25 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
26 * SUCH DAMAGE.
27 *
28 * CTPP2Dtoa.cpp
29 *
30 * $CTPP$
31 */
32
33 /****************************************************************
34 *
35 * The author of this software is David M. Gay.
36 *
37 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
38 *
39 * Permission to use, copy, modify, and distribute this software for any
40 * purpose without fee is hereby granted, provided that this entire notice
41 * is included in all copies of any software which is or includes a copy
42 * or modification of this software and in all copies of the supporting
43 * documentation for such software.
44 *
45 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
46 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
47 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
48 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
49 *
50 ***************************************************************/
51
52 /* Please send bug reports to
53 David M. Gay
54 Bell Laboratories, Room 2C-463
55 600 Mountain Avenue
56 Murray Hill, NJ 07974-0636
57 U.S.A.
58 dmg@bell-labs.com
59 */
60
61 /* On a machine with IEEE extended-precision registers, it is
62 * necessary to specify double-precision (53-bit) rounding precision
63 * before invoking strtod or dtoa. If the machine uses (the equivalent
64 * of) Intel 80x87 arithmetic, the call
65 * _control87(PC_53, MCW_PC);
66 * does this with many compilers. Whether this or another call is
67 * appropriate depends on the compiler; for this to work, it may be
68 * necessary to #include "float.h" or another system-dependent header
69 * file.
70 */
71
72 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
73 *
74 * This strtod returns a nearest machine number to the input decimal
75 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
76 * broken by the IEEE round-even rule. Otherwise ties are broken by
77 * biased rounding (add half and chop).
78 *
79 * Inspired loosely by William D. Clinger's paper "How to Read Floating
80 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
81 *
82 * Modifications:
83 *
84 * 1. We only require IEEE, IBM, or VAX double-precision
85 * arithmetic (not IEEE double-extended).
86 * 2. We get by with floating-point arithmetic in a case that
87 * Clinger missed -- when we're computing d * 10^n
88 * for a small integer d and the integer n is not too
89 * much larger than 22 (the maximum integer k for which
90 * we can represent 10^k exactly), we may be able to
91 * compute (d*10^k) * 10^(e-k) with just one roundoff.
92 * 3. Rather than a bit-at-a-time adjustment of the binary
93 * result in the hard case, we use floating-point
94 * arithmetic to determine the adjustment to within
95 * one bit; only in really hard cases do we need to
96 * compute a second residual.
97 * 4. Because of 3., we don't need a large table of powers of 10
98 * for ten-to-e (just some small tables, e.g. of 10^k
99 * for 0 <= k <= 22).
100 */
101 #define Honor_FLT_ROUNDS
102 #define FLT_ROUNDS 3
103 /*
104 * #define IEEE_8087 for IEEE-arithmetic machines where the least
105 * significant byte has the lowest address.
106 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
107 * significant byte has the lowest address.
108 * #define INT_32 int on machines with 32-bit ints and 64-bit longs.
109 * #define IBM for IBM mainframe-style floating-point arithmetic.
110 * #define VAX for VAX-style floating-point arithmetic (D_floating).
111 * #define No_leftright to omit left-right logic in fast floating-point
112 * computation of dtoa.
113 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
114 * and strtod and dtoa should round accordingly.
115 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
116 * and Honor_FLT_ROUNDS is not #defined.
117 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
118 * that use extended-precision instructions to compute rounded
119 * products and quotients) with IBM.
120 * #define ROUND_BIASED for IEEE-format with biased rounding.
121 * #define Inaccurate_Divide for IEEE-format with correctly rounded
122 * products but inaccurate quotients, e.g., for Intel i860.
123 * #define NO_LONG_LONG on machines that do not have a "INT_32 INT_32"
124 * integer type (of >= 64 bits). On such machines, you can
125 * #define Just_16 to store 16 bits per 32-bit INT_32 when doing
126 * high-precision integer arithmetic. Whether this speeds things
127 * up or slows things down depends on the machine and the number
128 * being converted. If INT_32 INT_32 is available and the name is
129 * something other than "INT_32 INT_32", #define INT_64 to be the name,
130 * and if "unsigned INT_64" does not work as an unsigned version of
131 * INT_64, #define #UINT_64 to be the corresponding unsigned type.
132 * #define Bad_float_h if your system lacks a float.h or if it does not
133 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
134 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
135 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
136 * if memory is available and otherwise does something you deem
137 * appropriate. If MALLOC is undefined, malloc will be invoked
138 * directly -- and assumed always to succeed.
139 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
140 * avoids underflows on inputs whose result does not underflow.
141 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
142 * floating-point numbers and flushes underflows to zero rather
143 * than implementing gradual underflow, then you must also #define
144 * #define YES_ALIAS to permit aliasing certain double values with
145 * arrays of ULongs. This leads to slightly better code with
146 * some compilers and was always used prior to 19990916, but it
147 * is not strictly legal and can cause trouble with aggressively
148 * optimizing compilers (e.g., gcc 2.95.1 under -O2).
149 * #define USE_LOCALE to use the current locale's decimal_point value.
150 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
151 * the result overflows to +-Infinity or underflows to 0.
152 */
153
154 #include "CTPP2DTOA.hpp"
155
156 #include <stdlib.h>
157 #include <string.h>
158 #include <stdio.h>
159 #include <errno.h>
160
161 #ifdef _MSC_VER
162 #include <WinSock2.h>
163 #ifndef BIG_ENDIAN
164 #define BIG_ENDIAN BIGENDIAN
165 #endif
166 #ifndef LITTLE_ENDIAN
167 #define LITTLE_ENDIAN LITTLEENDIAN
168 #endif
169 #ifndef BYTE_ORDER
170 #if (LITTLEENDIAN > BIGENDIAN)
171 #define BYTE_ORDER LITTLE_ENDIAN
172 #else
173 #define BYTE_ORDER BIG_ENDIAN
174 #endif
175 #endif
176 #endif
177
178 #ifndef BYTE_ORDER
179 #error "BYTE_ORDER is not defined!"
180 #endif
181
182 #define PLATFORM(x) (BYTE_ORDER == (x))
183
184 #if PLATFORM(BIG_ENDIAN)
185 #define IEEE_MC68k
186 #else
187 #define IEEE_8087
188 #endif
189
190 #ifdef IEEE_Arith
191 #undef IEEE_Arith
192 #endif
193
194 #ifdef Avoid_Underflow
195 #undef Avoid_Underflow
196 #endif
197
198 #ifdef IEEE_MC68k
199 #define IEEE_Arith
200 #endif
201
202 #ifdef IEEE_8087
203 #define IEEE_Arith
204 #endif
205
206 #ifdef Bad_float_h
207
208 #ifdef IEEE_Arith
209 #define DBL_DIG 15
210 #define DBL_MAX_10_EXP 308
211 #define DBL_MAX_EXP 1024
212 #define FLT_RADIX 2
213 #endif /*IEEE_Arith*/
214
215 #ifdef IBM
216 #define DBL_DIG 16
217 #define DBL_MAX_10_EXP 75
218 #define DBL_MAX_EXP 63
219 #define FLT_RADIX 16
220 #define DBL_MAX 7.2370055773322621e+75
221 #endif
222
223 #ifdef VAX
224 #define DBL_DIG 16
225 #define DBL_MAX_10_EXP 38
226 #define DBL_MAX_EXP 127
227 #define FLT_RADIX 2
228 #define DBL_MAX 1.7014118346046923e+38
229 #endif
230
231 #ifndef LONG_MAX
232 #define LONG_MAX 2147483647
233 #endif
234
235 #else /* ifndef Bad_float_h */
236 #include <float.h>
237 #endif /* Bad_float_h */
238
239 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
240 #error "Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined."
241 #endif
242
243 typedef union
244 {
245 double d;
246 UINT_32 L[2];
247 } U;
248
249 #define dval(x) (x).d
250 #ifdef IEEE_8087
251 #define word0(x) (x).L[1]
252 #define word1(x) (x).L[0]
253 #else
254 #define word0(x) (x).L[0]
255 #define word1(x) (x).L[1]
256 #endif
257
258 /* The following definition of Storeinc is appropriate for MIPS processors.
259 * An alternative that might be better on some machines is
260 */
261 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
262
263 /* #define P DBL_MANT_DIG */
264 /* Ten_pmax = floor(P*log(2)/log(5)) */
265 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
266 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
267 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
268
269 #ifdef IEEE_Arith
270 #define Exp_shift 20
271 #define Exp_shift1 20
272 #define Exp_msk1 0x100000
273 #define Exp_msk11 0x100000
274 #define Exp_mask 0x7ff00000
275 #define P 53
276 #define Bias 1023
277 #define Emin (-1022)
278 #define Exp_1 0x3ff00000
279 #define Exp_11 0x3ff00000
280 #define Ebits 11
281 #define Frac_mask 0xfffff
282 #define Frac_mask1 0xfffff
283 #define Ten_pmax 22
284 #define Bletch 0x10
285 #define Bndry_mask 0xfffff
286 #define Bndry_mask1 0xfffff
287 #define LSB 1
288 #define Sign_bit 0x80000000
289 #define Log2P 1
290 #define Tiny0 0
291 #define Tiny1 1
292 #define Quick_max 14
293 #define Int_max 14
294 #ifndef NO_IEEE_Scale
295 #define Avoid_Underflow
296 #endif
297
298 #ifndef Flt_Rounds
299 #ifdef FLT_ROUNDS
300 #define Flt_Rounds FLT_ROUNDS
301 #else
302 #define Flt_Rounds 1
303 #endif
304 #endif /*Flt_Rounds*/
305
306 #ifdef Honor_FLT_ROUNDS
307 #define Rounding rounding
308 #undef Check_FLT_ROUNDS
309 #define Check_FLT_ROUNDS
310 #else
311 #define Rounding Flt_Rounds
312 #endif
313
314 #else /* ifndef IEEE_Arith */
315 #undef Check_FLT_ROUNDS
316 #undef Honor_FLT_ROUNDS
317 #ifdef IBM
318 #undef Flt_Rounds
319 #define Flt_Rounds 0
320 #define Exp_shift 24
321 #define Exp_shift1 24
322 #define Exp_msk1 0x1000000
323 #define Exp_msk11 0x1000000
324 #define Exp_mask 0x7f000000
325 #define P 14
326 #define Bias 65
327 #define Exp_1 0x41000000
328 #define Exp_11 0x41000000
329 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
330 #define Frac_mask 0xffffff
331 #define Frac_mask1 0xffffff
332 #define Bletch 4
333 #define Ten_pmax 22
334 #define Bndry_mask 0xefffff
335 #define Bndry_mask1 0xffffff
336 #define LSB 1
337 #define Sign_bit 0x80000000
338 #define Log2P 4
339 #define Tiny0 0x100000
340 #define Tiny1 0
341 #define Quick_max 14
342 #define Int_max 15
343 #else /* VAX */
344 #undef Flt_Rounds
345 #define Flt_Rounds 1
346 #define Exp_shift 23
347 #define Exp_shift1 7
348 #define Exp_msk1 0x80
349 #define Exp_msk11 0x800000
350 #define Exp_mask 0x7f80
351 #define P 56
352 #define Bias 129
353 #define Exp_1 0x40800000
354 #define Exp_11 0x4080
355 #define Ebits 8
356 #define Frac_mask 0x7fffff
357 #define Frac_mask1 0xffff007f
358 #define Ten_pmax 24
359 #define Bletch 2
360 #define Bndry_mask 0xffff007f
361 #define Bndry_mask1 0xffff007f
362 #define LSB 0x10000
363 #define Sign_bit 0x8000
364 #define Log2P 1
365 #define Tiny0 0x80
366 #define Tiny1 0
367 #define Quick_max 15
368 #define Int_max 15
369 #endif /* IBM, VAX */
370 #endif /* IEEE_Arith */
371
372 #ifndef IEEE_Arith
373 #define ROUND_BIASED
374 #endif
375
376 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
377 #define Big1 0xffffffff
378 #define FFFFFFFF 0xffffffffUL
379
380 #define ACQUIRE_DTOA_LOCK(n)
381 #define FREE_DTOA_LOCK(n)
382
safe_malloc(AllocatedBlock ** aBlocks,const UINT_32 iSize)383 void * safe_malloc(AllocatedBlock ** aBlocks, const UINT_32 iSize)
384 {
385 // Allocate memory block placeholder
386 AllocatedBlock * pNewBlock = (AllocatedBlock *)malloc(sizeof(AllocatedBlock));
387 // Allocate memory block
388 pNewBlock -> data = malloc(iSize);
389 pNewBlock -> prev = *aBlocks;
390
391 // Ouch!
392 *aBlocks = pNewBlock;
393 //fprintf(stderr, "safe_malloc: %p(%p, %d)\n", (void*)pNewBlock, (void*)(pNewBlock -> data), iSize);
394
395 return pNewBlock -> data;
396 }
397
safe_free(AllocatedBlock ** aBlocks)398 void safe_free(AllocatedBlock ** aBlocks)
399 {
400 // Ouch!
401 if (aBlocks == NULL || *aBlocks == NULL) { return; }
402
403 for(;;)
404 {
405 AllocatedBlock * pPrevBlock = (*aBlocks) -> prev;
406
407 //fprintf(stderr, "safe_free: %p(%p)\n", (void*)(*aBlocks), (*aBlocks) -> data);
408
409 free((*aBlocks) -> data);
410 free(*aBlocks);
411
412 *aBlocks = pPrevBlock;
413 if (pPrevBlock == NULL) { break; }
414 }
415 }
416
Balloc(AllocatedBlock ** aBlocks,Bigint ** freelist,int k)417 static Bigint * Balloc(AllocatedBlock ** aBlocks, Bigint ** freelist, int k)
418 {
419 int x;
420 Bigint *rv;
421
422 ACQUIRE_DTOA_LOCK(0);
423 if ((rv = freelist[k]))
424 {
425 freelist[k] = rv->next;
426 }
427 else
428 {
429 x = 1 << k;
430 rv = (Bigint *)safe_malloc(aBlocks, sizeof(Bigint) + (x-1)*sizeof(UINT_32));
431 rv->k = k;
432 rv->maxwds = x;
433 }
434 FREE_DTOA_LOCK(0);
435 rv->sign = rv->wds = 0;
436
437 return rv;
438 }
439
440 #define Kmax (sizeof(size_t) << 3)
441
Bfree(Bigint ** freelist,Bigint * v)442 static void Bfree(Bigint ** freelist, Bigint *v)
443 {
444 if (v)
445 {
446 // if (v->k > Kmax)
447 // {
448 // free((void*)v);
449 // }
450 // else
451 {
452 ACQUIRE_DTOA_LOCK(0);
453 v->next = freelist[v->k];
454 freelist[v->k] = v;
455 FREE_DTOA_LOCK(0);
456 }
457 }
458 }
459
460 /*static void Bfree(Bigint ** freelist, Bigint *v)
461 {
462 if (v)
463 {
464 ACQUIRE_DTOA_LOCK(0);
465 v->next = freelist[v->k];
466 freelist[v->k] = v;
467 FREE_DTOA_LOCK(0);
468 }
469 }
470 */
471 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, y->wds*sizeof(INT_32) + 2*sizeof(int))
472
473 //
474 // multiply by m and add a
475 //
multadd(AllocatedBlock ** aBlocks,Bigint ** freelist,Bigint * b,int m,int a)476 static Bigint * multadd(AllocatedBlock ** aBlocks, Bigint ** freelist, Bigint *b, int m, int a)
477 {
478 int i, wds;
479
480 UINT_32 *x;
481 UINT_64 carry, y;
482 Bigint *b1;
483
484 wds = b->wds;
485 x = b->x;
486 i = 0;
487 carry = a;
488 do
489 {
490 y = *x * (UINT_64)m + carry;
491 carry = y >> 32;
492 *x++ = (UINT_32)y & FFFFFFFF;
493 }
494 while(++i < wds);
495
496 if (carry)
497 {
498 if (wds >= b->maxwds)
499 {
500 b1 = Balloc(aBlocks, freelist, b->k+1);
501 Bcopy(b1, b);
502 Bfree(freelist, b);
503 b = b1;
504 }
505 b->x[wds++] = (UINT_32)carry;
506 b->wds = wds;
507 }
508
509 return b;
510 }
511
512
hi0bits(register UINT_32 x)513 static int hi0bits(register UINT_32 x)
514 {
515 register int k = 0;
516
517 if (!(x & 0xffff0000))
518 {
519 k = 16;
520 x <<= 16;
521 }
522 if (!(x & 0xff000000))
523 {
524 k += 8;
525 x <<= 8;
526 }
527 if (!(x & 0xf0000000))
528 {
529 k += 4;
530 x <<= 4;
531 }
532 if (!(x & 0xc0000000))
533 {
534 k += 2;
535 x <<= 2;
536 }
537 if (!(x & 0x80000000))
538 {
539 k++;
540 if (!(x & 0x40000000))
541 {
542 return 32;
543 }
544 }
545
546 return k;
547 }
548
lo0bits(UINT_32 * y)549 static int lo0bits(UINT_32 *y)
550 {
551 register int k;
552 register UINT_32 x = *y;
553
554 if (x & 7)
555 {
556 if (x & 1)
557 {
558 return 0;
559 }
560 if (x & 2)
561 {
562 *y = x >> 1;
563 return 1;
564 }
565 *y = x >> 2;
566 return 2;
567 }
568
569 k = 0;
570 if (!(x & 0xffff))
571 {
572 k = 16;
573 x >>= 16;
574 }
575 if (!(x & 0xff))
576 {
577 k += 8;
578 x >>= 8;
579 }
580 if (!(x & 0xf))
581 {
582 k += 4;
583 x >>= 4;
584 }
585 if (!(x & 0x3))
586 {
587 k += 2;
588 x >>= 2;
589 }
590 if (!(x & 1))
591 {
592 k++;
593 x >>= 1;
594 if (!x & 1)
595 {
596 return 32;
597 }
598 }
599 *y = x;
600 return k;
601 }
602
i2b(AllocatedBlock ** aBlocks,Bigint ** freelist,int i)603 static Bigint * i2b(AllocatedBlock ** aBlocks, Bigint ** freelist, int i)
604 {
605 Bigint *b;
606
607 b = Balloc(aBlocks, freelist, 1);
608 b->x[0] = i;
609 b->wds = 1;
610
611 return b;
612 }
613
mult(AllocatedBlock ** aBlocks,Bigint ** freelist,Bigint * a,Bigint * b)614 static Bigint * mult(AllocatedBlock ** aBlocks, Bigint ** freelist, Bigint *a, Bigint *b)
615 {
616 Bigint *c;
617 int k, wa, wb, wc;
618 UINT_32 *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
619 UINT_32 y;
620 UINT_64 carry, z;
621
622 if (a->wds < b->wds)
623 {
624 c = a;
625 a = b;
626 b = c;
627 }
628
629 k = a->k;
630 wa = a->wds;
631 wb = b->wds;
632 wc = wa + wb;
633 if (wc > a->maxwds) { k++; }
634 c = Balloc(aBlocks, freelist, k);
635 for(x = c->x, xa = x + wc; x < xa; x++) { *x = 0; }
636
637 xa = a->x;
638 xae = xa + wa;
639
640 xb = b->x;
641 xbe = xb + wb;
642
643 xc0 = c->x;
644
645 for(; xb < xbe; xc0++)
646 {
647 if ((y = *xb++))
648 {
649 x = xa;
650 xc = xc0;
651 carry = 0;
652 do
653 {
654 z = *x++ * (UINT_64)y + *xc + carry;
655 carry = z >> 32;
656 *xc++ = (UINT_32)z & FFFFFFFF;
657 }
658 while(x < xae);
659 *xc = (UINT_32)carry;
660 }
661 }
662
663 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) { ;; }
664
665 c->wds = wc;
666
667 return c;
668 }
669
670 //static Bigint *p5s;
671
pow5mult(AllocatedBlock ** aBlocks,Bigint ** freelist,Bigint * b,int k)672 static Bigint * pow5mult(AllocatedBlock ** aBlocks, Bigint ** freelist, Bigint *b, int k)
673 {
674 Bigint *b1, *p5, *p51;
675 int i;
676 static int p05[3] = { 5, 25, 125 };
677
678 if ((i = k & 3)) { b = multadd(aBlocks, freelist, b, p05[i-1], 0); }
679
680 if (!(k >>= 2)) { return b; }
681 // if (!(p5 = p5s))
682 {
683 /* first time */
684 ACQUIRE_DTOA_LOCK(1);
685 // if (!(p5 = p5s))
686 {
687 p5 /* = p5s*/ = i2b(aBlocks, freelist, 625);
688 p5->next = 0;
689 }
690 FREE_DTOA_LOCK(1);
691 }
692
693 for(;;)
694 {
695 if (k & 1)
696 {
697 b1 = mult(aBlocks, freelist, b, p5);
698 Bfree(freelist, b);
699 b = b1;
700 }
701 if (!(k >>= 1)) { break; }
702
703 if (!(p51 = p5->next))
704 {
705 ACQUIRE_DTOA_LOCK(1);
706 if (!(p51 = p5->next))
707 {
708 p51 = p5->next = mult(aBlocks, freelist, p5,p5);
709 p51->next = 0;
710 }
711 FREE_DTOA_LOCK(1);
712 }
713 p5 = p51;
714 }
715 return b;
716 }
717
lshift(AllocatedBlock ** aBlocks,Bigint ** freelist,Bigint * b,int k)718 static Bigint * lshift(AllocatedBlock ** aBlocks, Bigint ** freelist, Bigint *b, int k)
719 {
720 int i, k1, n, n1;
721 Bigint *b1;
722 UINT_32 *x, *x1, *xe, z;
723
724 n = k >> 5;
725 k1 = b->k;
726 n1 = n + b->wds + 1;
727 for(i = b->maxwds; n1 > i; i <<= 1) { k1++; }
728
729 b1 = Balloc(aBlocks, freelist, k1);
730 x1 = b1->x;
731 for(i = 0; i < n; i++) { *x1++ = 0; }
732
733 x = b->x;
734 xe = x + b->wds;
735
736 if (k &= 0x1f)
737 {
738 k1 = 32 - k;
739 z = 0;
740 do
741 {
742 *x1++ = *x << k | z;
743 z = *x++ >> k1;
744 }
745 while(x < xe);
746 if ((*x1 = z)) { ++n1; }
747 }
748 else
749 {
750 do { *x1++ = *x++; } while(x < xe);
751 }
752 b1->wds = n1 - 1;
753 Bfree(freelist, b);
754
755 return b1;
756 }
757
cmp(Bigint * a,Bigint * b)758 static int cmp(Bigint *a, Bigint *b)
759 {
760 UINT_32 *xa, *xa0, *xb, *xb0;
761 int i, j;
762
763 i = a->wds;
764 j = b->wds;
765 #ifdef DEBUG
766 if (i > 1 && !a->x[i-1])
767 Bug("cmp called with a->x[a->wds-1] == 0");
768 if (j > 1 && !b->x[j-1])
769 Bug("cmp called with b->x[b->wds-1] == 0");
770 #endif
771 if (i -= j) { return i; }
772 xa0 = a->x;
773 xa = xa0 + j;
774 xb0 = b->x;
775 xb = xb0 + j;
776 for(;;)
777 {
778 if (*--xa != *--xb)
779 {
780 return *xa < *xb ? -1 : 1;
781 }
782 if (xa <= xa0)
783 {
784 break;
785 }
786 }
787
788 return 0;
789 }
790
diff(AllocatedBlock ** aBlocks,Bigint ** freelist,Bigint * a,Bigint * b)791 static Bigint * diff(AllocatedBlock ** aBlocks, Bigint ** freelist, Bigint *a, Bigint *b)
792 {
793 Bigint *c;
794 int i, wa, wb;
795 UINT_32 *xa, *xae, *xb, *xbe, *xc;
796
797 UINT_64 borrow, y;
798
799 i = cmp(a,b);
800 if (!i) {
801 c = Balloc(aBlocks, freelist, 0);
802 c->wds = 1;
803 c->x[0] = 0;
804 return c;
805 }
806 if (i < 0) {
807 c = a;
808 a = b;
809 b = c;
810 i = 1;
811 }
812 else
813 i = 0;
814 c = Balloc(aBlocks, freelist, a->k);
815 c->sign = i;
816 wa = a->wds;
817 xa = a->x;
818 xae = xa + wa;
819 wb = b->wds;
820 xb = b->x;
821 xbe = xb + wb;
822 xc = c->x;
823 borrow = 0;
824
825 do
826 {
827 y = (UINT_64)*xa++ - *xb++ - borrow;
828 borrow = y >> 32 & (UINT_32)1;
829 *xc++ = (UINT_32)y & FFFFFFFF;
830 }
831 while(xb < xbe);
832 while(xa < xae)
833 {
834 y = *xa++ - borrow;
835 borrow = y >> 32 & (UINT_32)1;
836 *xc++ = (UINT_32)y & FFFFFFFF;
837 }
838
839 while(!*--xc)
840 {
841 wa--;
842 }
843 c->wds = wa;
844
845 return c;
846 }
847
d2b(AllocatedBlock ** aBlocks,Bigint ** freelist,double dd,int * e,int * bits)848 static Bigint * d2b(AllocatedBlock ** aBlocks, Bigint ** freelist, double dd, int *e, int *bits)
849 {
850 U d;
851 Bigint *b;
852 int de, k;
853 UINT_32 *x, y, z;
854 int i;
855 #ifdef VAX
856 UINT_32 d0, d1;
857 #endif
858 dval(d) = dd;
859 #ifdef VAX
860 d0 = word0(d) >> 16 | word0(d) << 16;
861 d1 = word1(d) >> 16 | word1(d) << 16;
862 #else
863 #define d0 word0(d)
864 #define d1 word1(d)
865 #endif
866 b = Balloc(aBlocks, freelist, 1);
867 x = b->x;
868
869 z = d0 & Frac_mask;
870 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
871
872 if ((de = (int)(d0 >> Exp_shift))) { z |= Exp_msk1; }
873
874 if ((y = d1))
875 {
876 if ((k = lo0bits(&y)))
877 {
878 x[0] = y | (z << (32 - k));
879 z >>= k;
880 }
881 else
882 {
883 x[0] = y;
884 }
885 i = b->wds = (x[1] = z) ? 2 : 1;
886 }
887 else
888 {
889 #ifdef DEBUG
890 if (!z)
891 Bug("Zero passed to d2b");
892 #endif
893 k = lo0bits(&z);
894 x[0] = z;
895 i = b->wds = 1;
896 k += 32;
897 }
898
899 if (de)
900 {
901 #ifdef IBM
902 *e = (de - Bias - (P-1) << 2) + k;
903 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
904 #else
905 *e = de - Bias - (P-1) + k;
906 *bits = P - k;
907 #endif
908 }
909 else
910 {
911 *e = de - Bias - (P-1) + 1 + k;
912
913 *bits = 32*i - hi0bits(x[i-1]);
914 }
915
916 return b;
917 }
918
919 #undef d0
920 #undef d1
921
922
923 static const double tens[] =
924 {
925 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7,
926 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15,
927 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22
928 #ifdef VAX
929 , 1e23, 1e24
930 #endif
931 };
932
933 static const double
934 #ifdef IEEE_Arith
935 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
936 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
937 #ifdef Avoid_Underflow
938 9007199254740992.*9007199254740992.e-256
939 /* = 2^106 * 1e-53 */
940 #else
941 1e-256
942 #endif
943 };
944 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
945 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
946 #define Scale_Bit 0x10
947 #define n_bigtens 5
948 #else
949 #ifdef IBM
950 bigtens[] = { 1e16, 1e32, 1e64 };
951 static const double tinytens[] = { 1e-16, 1e-32, 1e-64 };
952 #define n_bigtens 3
953 #else
954 bigtens[] = { 1e16, 1e32 };
955 static const double tinytens[] = { 1e-16, 1e-32 };
956 #define n_bigtens 2
957 #endif
958 #endif
959
quorem(Bigint * b,Bigint * S)960 static int quorem(Bigint *b, Bigint *S)
961 {
962 int n;
963 UINT_32 *bx, *bxe, q, *sx, *sxe;
964
965 UINT_64 borrow, carry, y, ys;
966
967 n = S->wds;
968 #ifdef DEBUG
969 /*debug*/ if (b->wds > n)
970 /*debug*/ Bug("oversize b in quorem");
971 #endif
972 if (b->wds < n)
973 return 0;
974 sx = S->x;
975 sxe = sx + --n;
976 bx = b->x;
977 bxe = bx + n;
978 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
979 #ifdef DEBUG
980 /*debug*/ if (q > 9)
981 /*debug*/ Bug("oversized quotient in quorem");
982 #endif
983 if (q) {
984 borrow = 0;
985 carry = 0;
986 do
987 {
988 ys = *sx++ * (UINT_64)q + carry;
989 carry = ys >> 32;
990 y = *bx - (ys & FFFFFFFF) - borrow;
991 borrow = y >> 32 & (UINT_32)1;
992 *bx++ = (UINT_32)y & FFFFFFFF;
993 }
994 while(sx <= sxe);
995
996 if (!*bxe)
997 {
998 bx = b->x;
999 while(--bxe > bx && !*bxe)
1000 --n;
1001 b->wds = n;
1002 }
1003 }
1004 if (cmp(b, S) >= 0)
1005 {
1006 q++;
1007 borrow = 0;
1008 carry = 0;
1009 bx = b->x;
1010 sx = S->x;
1011 do
1012 {
1013 ys = *sx++ + carry;
1014 carry = ys >> 32;
1015 y = *bx - (ys & FFFFFFFF) - borrow;
1016 borrow = y >> 32 & (UINT_32)1;
1017 *bx++ = (UINT_32)y & FFFFFFFF;
1018 }
1019 while(sx <= sxe);
1020
1021 bx = b->x;
1022 bxe = bx + n;
1023 if (!*bxe)
1024 {
1025 while(--bxe > bx && !*bxe) { --n; }
1026 b->wds = n;
1027 }
1028 }
1029
1030 return q;
1031 }
1032
rv_alloc(AllocatedBlock ** aBlocks,Bigint ** freelist,int i)1033 static char * rv_alloc(AllocatedBlock ** aBlocks, Bigint ** freelist, int i)
1034 {
1035 int j, k, *r;
1036
1037 j = sizeof(UINT_32);
1038 for(k = 0; sizeof(Bigint) - sizeof(UINT_32) - sizeof(int) + j <= (unsigned)i; j <<= 1) { k++; }
1039 r = (int*)Balloc(aBlocks, freelist, k);
1040 *r = k;
1041 return (char *)(r+1);
1042 }
1043
nrv_alloc(AllocatedBlock ** aBlocks,Bigint ** freelist,const char * s,char ** rve,int n)1044 static char * nrv_alloc(AllocatedBlock ** aBlocks, Bigint ** freelist, const char *s, char **rve, int n)
1045 {
1046 char *rv, *t;
1047
1048 t = rv = rv_alloc(aBlocks, freelist, n);
1049 while((*t = *s++)) { t++; }
1050 if (rve) { *rve = t; }
1051 return rv;
1052 }
1053
freedtoa(AllocatedBlock ** aBlocks)1054 void freedtoa(AllocatedBlock ** aBlocks)
1055 {
1056 safe_free(aBlocks);
1057 }
1058
1059 /*void freedtoa(Bigint ** freelist, char *s)
1060 {
1061 Bigint *b = (Bigint *)((int *)s - 1);
1062 b->maxwds = 1 << (b->k = *(int*)b);
1063 Bfree(freelist, b);
1064 }
1065 */
1066 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1067 *
1068 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1069 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1070 *
1071 * Modifications:
1072 * 1. Rather than iterating, we use a simple numeric overestimate
1073 * to determine k = floor(log10(d)). We scale relevant
1074 * quantities using O(log2(k)) rather than O(k) multiplications.
1075 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1076 * try to generate digits strictly left to right. Instead, we
1077 * compute with fewer bits and propagate the carry if necessary
1078 * when rounding the final digit up. This is often faster.
1079 * 3. Under the assumption that input will be rounded nearest,
1080 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1081 * That is, we allow equality in stopping tests when the
1082 * round-nearest rule will give the same floating-point value
1083 * as would satisfaction of the stopping test with strict
1084 * inequality.
1085 * 4. We remove common factors of powers of 2 from relevant
1086 * quantities.
1087 * 5. When converting floating-point integers less than 1e16,
1088 * we use floating-point arithmetic rather than resorting
1089 * to multiple-precision integers.
1090 * 6. When asked to produce fewer than 15 digits, we first try
1091 * to get by with floating-point arithmetic; we resort to
1092 * multiple-precision integer arithmetic only if we cannot
1093 * guarantee that the floating-point calculation has given
1094 * the correctly rounded result. For k requested digits and
1095 * "uniformly" distributed input, the probability is
1096 * something like 10^(k-15) that we must resort to the INT_32
1097 * calculation.
1098 */
1099
ctpp_dtoa(AllocatedBlock ** aBlocks,Bigint ** freelist,double dd,int mode,int ndigits,int * decpt,int * sign,char ** rve)1100 char * ctpp_dtoa(AllocatedBlock ** aBlocks, Bigint ** freelist, double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
1101 {
1102 /* Arguments ndigits, decpt, sign are similar to those
1103 of ecvt and fcvt; trailing zeros are suppressed from
1104 the returned string. If not null, *rve is set to point
1105 to the end of the return value. If d is +-Infinity or NaN,
1106 then *decpt is set to 9999.
1107
1108 mode:
1109 0 ==> shortest string that yields d when read in
1110 and rounded to nearest.
1111 1 ==> like 0, but with Steele & White stopping rule;
1112 e.g. with IEEE P754 arithmetic , mode 0 gives
1113 1e23 whereas mode 1 gives 9.999999999999999e22.
1114 2 ==> max(1,ndigits) significant digits. This gives a
1115 return value similar to that of ecvt, except
1116 that trailing zeros are suppressed.
1117 3 ==> through ndigits past the decimal point. This
1118 gives a return value similar to that from fcvt,
1119 except that trailing zeros are suppressed, and
1120 ndigits can be negative.
1121 4,5 ==> similar to 2 and 3, respectively, but (in
1122 round-nearest mode) with the tests of mode 0 to
1123 possibly return a shorter string that rounds to d.
1124 With IEEE arithmetic and compilation with
1125 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
1126 as modes 2 and 3 when FLT_ROUNDS != 1.
1127 6-9 ==> Debugging modes similar to mode - 4: don't try
1128 fast floating-point estimate (if applicable).
1129
1130 Values of mode other than 0-9 are treated as mode 0.
1131
1132 Sufficient space is allocated to the return value
1133 to hold the suppressed trailing zeros.
1134 */
1135
1136 INT_32 bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, spec_case, try_quick;
1137 INT_32 L;
1138
1139 INT_32 denorm;
1140 UINT_32 x;
1141
1142 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
1143 U d, d2, eps;
1144 double ds;
1145 char *s, *s0;
1146 #ifdef Honor_FLT_ROUNDS
1147 int rounding;
1148 #endif
1149
1150 dval(d) = dd;
1151 if (word0(d) & Sign_bit)
1152 {
1153 /* set sign for everything, including 0's and NaNs */
1154 *sign = 1;
1155 word0(d) &= ~Sign_bit; /* clear sign bit */
1156 }
1157 else
1158 {
1159 *sign = 0;
1160 }
1161
1162 #if defined(IEEE_Arith) + defined(VAX)
1163 #ifdef IEEE_Arith
1164 if ((word0(d) & Exp_mask) == Exp_mask)
1165 #else
1166 if (word0(d) == 0x8000)
1167 #endif
1168 {
1169 /* Infinity or NaN */
1170 *decpt = 9999;
1171 #ifdef IEEE_Arith
1172 if (!word1(d) && !(word0(d) & 0xfffff)) { return nrv_alloc(aBlocks, freelist, "Infinity", rve, 8); }
1173 #endif
1174 return nrv_alloc(aBlocks, freelist, "NaN", rve, 3);
1175 }
1176 #endif
1177
1178 #ifdef IBM
1179 dval(d) += 0; /* normalize */
1180 #endif
1181 if (!dval(d))
1182 {
1183 *decpt = 1;
1184 return nrv_alloc(aBlocks, freelist, "0", rve, 1);
1185 }
1186
1187 #ifdef Honor_FLT_ROUNDS
1188 if ((rounding = Flt_Rounds) >= 2)
1189 {
1190 if (*sign)
1191 {
1192 rounding = (rounding == 2) ? 0 : 2;
1193 }
1194 else
1195 {
1196 if (rounding != 2)
1197 {
1198 rounding = 0;
1199 }
1200 }
1201 }
1202 #endif
1203
1204 b = d2b(aBlocks, freelist, dval(d), &be, &bbits);
1205
1206 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1207 {
1208 dval(d2) = dval(d);
1209 word0(d2) &= Frac_mask1;
1210 word0(d2) |= Exp_11;
1211 #ifdef IBM
1212 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) { dval(d2) /= 1 << j; }
1213 #endif
1214
1215 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1216 * log10(x) = log(x) / log(10)
1217 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1218 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1219 *
1220 * This suggests computing an approximation k to log10(d) by
1221 *
1222 * k = (i - Bias)*0.301029995663981
1223 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1224 *
1225 * We want k to be too large rather than too small.
1226 * The error in the first-order Taylor series approximation
1227 * is in our favor, so we just round up the constant enough
1228 * to compensate for any error in the multiplication of
1229 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1230 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1231 * adding 1e-13 to the constant term more than suffices.
1232 * Hence we adjust the constant term to 0.1760912590558.
1233 * (We could get a more accurate k by invoking log10,
1234 * but this is probably not worthwhile.)
1235 */
1236
1237 i -= Bias;
1238 #ifdef IBM
1239 i <<= 2;
1240 i += j;
1241 #endif
1242 denorm = 0;
1243 }
1244 else
1245 {
1246 /* d is denormalized */
1247
1248 i = bbits + be + (Bias + (P-1) - 1);
1249 x = i > 32 ? (word0(d) << (64 - i)) | (word1(d) >> (i - 32))
1250 : word1(d) << (32 - i);
1251 dval(d2) = x;
1252 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1253 i -= (Bias + (P-1) - 1) + 1;
1254 denorm = 1;
1255 }
1256
1257 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1258 k = (int)ds;
1259 if (ds < 0. && ds != k)
1260 k--; /* want k = floor(ds) */
1261 k_check = 1;
1262 if (k >= 0 && k <= Ten_pmax)
1263 {
1264 if (dval(d) < tens[k]) { k--; }
1265 k_check = 0;
1266 }
1267 j = bbits - i - 1;
1268
1269 if (j >= 0)
1270 {
1271 b2 = 0;
1272 s2 = j;
1273 }
1274 else
1275 {
1276 b2 = -j;
1277 s2 = 0;
1278 }
1279
1280 if (k >= 0)
1281 {
1282 b5 = 0;
1283 s5 = k;
1284 s2 += k;
1285 }
1286 else
1287 {
1288 b2 -= k;
1289 b5 = -k;
1290 s5 = 0;
1291 }
1292
1293 if (mode < 0 || mode > 9) { mode = 0; }
1294
1295 #ifdef Check_FLT_ROUNDS
1296 try_quick = Rounding == 1;
1297 #else
1298 try_quick = 1;
1299 #endif
1300
1301 if (mode > 5)
1302 {
1303 mode -= 4;
1304 try_quick = 0;
1305 }
1306 leftright = 1;
1307 switch(mode)
1308 {
1309 case 0:
1310 case 1:
1311 ilim = ilim1 = -1;
1312 i = 18;
1313 ndigits = 0;
1314 break;
1315 case 2:
1316 leftright = 0;
1317 /* no break */
1318 case 4:
1319 if (ndigits <= 0)
1320 ndigits = 1;
1321 ilim = ilim1 = i = ndigits;
1322 break;
1323 case 3:
1324 leftright = 0;
1325 /* no break */
1326 case 5:
1327 i = ndigits + k + 1;
1328 ilim = i;
1329 ilim1 = i - 1;
1330 if (i <= 0)
1331 i = 1;
1332 }
1333 s = s0 = rv_alloc(aBlocks, freelist, i);
1334
1335 #ifdef Honor_FLT_ROUNDS
1336 if (mode > 1 && rounding != 1) { leftright = 0; }
1337 #endif
1338
1339 if (ilim >= 0 && ilim <= Quick_max && try_quick)
1340 {
1341 /* Try to get by with floating-point arithmetic. */
1342
1343 i = 0;
1344 dval(d2) = dval(d);
1345 k0 = k;
1346 ilim0 = ilim;
1347 ieps = 2; /* conservative */
1348 if (k > 0)
1349 {
1350 ds = tens[k&0xf];
1351 j = k >> 4;
1352 if (j & Bletch)
1353 {
1354 /* prevent overflows */
1355 j &= Bletch - 1;
1356 dval(d) /= bigtens[n_bigtens-1];
1357 ieps++;
1358 }
1359
1360 for(; j; j >>= 1, i++)
1361 {
1362 if (j & 1)
1363 {
1364 ieps++;
1365 ds *= bigtens[i];
1366 }
1367 }
1368
1369 dval(d) /= ds;
1370 }
1371 else if ((j1 = -k))
1372 {
1373 dval(d) *= tens[j1 & 0xf];
1374 for(j = j1 >> 4; j; j >>= 1, i++)
1375 {
1376 if (j & 1)
1377 {
1378 ieps++;
1379 dval(d) *= bigtens[i];
1380 }
1381 }
1382 }
1383 if (k_check && dval(d) < 1. && ilim > 0)
1384 {
1385 if (ilim1 <= 0) { goto fast_failed; }
1386 ilim = ilim1;
1387 k--;
1388 dval(d) *= 10.;
1389 ieps++;
1390 }
1391 dval(eps) = ieps*dval(d) + 7.;
1392 word0(eps) -= (P-1)*Exp_msk1;
1393 if (ilim == 0)
1394 {
1395 S = mhi = 0;
1396 dval(d) -= 5.;
1397 if (dval(d) > dval(eps)) { goto one_digit; }
1398 if (dval(d) < -dval(eps)) { goto no_digits; }
1399 goto fast_failed;
1400 }
1401 #ifndef No_leftright
1402 if (leftright)
1403 {
1404 /* Use Steele & White method of only
1405 * generating digits needed.
1406 */
1407 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
1408 for(i = 0;;) {
1409 L = (INT_32)dval(d);
1410 dval(d) -= L;
1411 *s++ = '0' + (int)L;
1412 if (dval(d) < dval(eps))
1413 goto ret1;
1414 if (1. - dval(d) < dval(eps))
1415 goto bump_up;
1416 if (++i >= ilim)
1417 break;
1418 dval(eps) *= 10.;
1419 dval(d) *= 10.;
1420 }
1421 }
1422 else
1423 {
1424 #endif
1425 /* Generate ilim digits, then fix them up. */
1426 dval(eps) *= tens[ilim-1];
1427 for(i = 1;; i++, dval(d) *= 10.)
1428 {
1429 L = (INT_32)(dval(d));
1430 if (!(dval(d) -= L))
1431 {
1432 ilim = i;
1433 }
1434 *s++ = '0' + (int)L;
1435 if (i == ilim)
1436 {
1437 if (dval(d) > 0.5 + dval(eps))
1438 {
1439 goto bump_up;
1440 }
1441 else if (dval(d) < 0.5 - dval(eps))
1442 {
1443 while(*--s == '0')
1444 {
1445 ;
1446 }
1447 s++;
1448 goto ret1;
1449 }
1450 break;
1451 }
1452 }
1453 #ifndef No_leftright
1454 }
1455 #endif
1456 fast_failed:
1457 s = s0;
1458 dval(d) = dval(d2);
1459 k = k0;
1460 ilim = ilim0;
1461 }
1462
1463 /* Do we have a "small" integer? */
1464
1465 if (be >= 0 && k <= Int_max)
1466 {
1467 /* Yes. */
1468 ds = tens[k];
1469 if (ndigits < 0 && ilim <= 0)
1470 {
1471 S = mhi = 0;
1472 if (ilim < 0 || dval(d) <= 5*ds)
1473 {
1474 goto no_digits;
1475 }
1476 goto one_digit;
1477 }
1478 for(i = 1;; i++, dval(d) *= 10.)
1479 {
1480 L = (INT_32)(dval(d) / ds);
1481 dval(d) -= L*ds;
1482 #ifdef Check_FLT_ROUNDS
1483 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1484 if (dval(d) < 0) {
1485 L--;
1486 dval(d) += ds;
1487 }
1488 #endif
1489 *s++ = '0' + (int)L;
1490 if (!dval(d)) {
1491 break;
1492 }
1493 if (i == ilim) {
1494 #ifdef Honor_FLT_ROUNDS
1495 if (mode > 1)
1496 switch(rounding) {
1497 case 0: goto ret1;
1498 case 2: goto bump_up;
1499 }
1500 #endif
1501 dval(d) += dval(d);
1502 if (dval(d) > ds || ((dval(d) == ds) && (L & 1)))
1503 {
1504 bump_up:
1505 while(*--s == '9')
1506 {
1507 if (s == s0)
1508 {
1509 k++;
1510 *s = '0';
1511 break;
1512 }
1513 }
1514 ++*s++;
1515 }
1516 break;
1517 }
1518 }
1519 goto ret1;
1520 }
1521
1522 m2 = b2;
1523 m5 = b5;
1524 mhi = mlo = 0;
1525 if (leftright)
1526 {
1527 i =
1528 denorm ? be + (Bias + (P-1) - 1 + 1) :
1529 #ifdef IBM
1530 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
1531 #else
1532 1 + P - bbits;
1533 #endif
1534 b2 += i;
1535 s2 += i;
1536 mhi = i2b(aBlocks, freelist, 1);
1537 }
1538
1539 if (m2 > 0 && s2 > 0)
1540 {
1541 i = m2 < s2 ? m2 : s2;
1542 b2 -= i;
1543 m2 -= i;
1544 s2 -= i;
1545 }
1546 if (b5 > 0)
1547 {
1548 if (leftright)
1549 {
1550 if (m5 > 0)
1551 {
1552 mhi = pow5mult(aBlocks, freelist, mhi, m5);
1553 b1 = mult(aBlocks, freelist, mhi, b);
1554 Bfree(freelist, b);
1555 b = b1;
1556 }
1557 if ((j = b5 - m5))
1558 {
1559 b = pow5mult(aBlocks, freelist, b, j);
1560 }
1561 }
1562 else
1563 {
1564 b = pow5mult(aBlocks, freelist, b, b5);
1565 }
1566 }
1567 S = i2b(aBlocks, freelist, 1);
1568 if (s5 > 0)
1569 {
1570 S = pow5mult(aBlocks, freelist, S, s5);
1571 }
1572
1573 /* Check for special case that d is a normalized power of 2. */
1574
1575 spec_case = 0;
1576 if ((mode < 2 || leftright)
1577 #ifdef Honor_FLT_ROUNDS
1578 && rounding == 1
1579 #endif
1580 )
1581 {
1582 if (!word1(d) && !(word0(d) & Bndry_mask) && word0(d) & (Exp_mask & ~Exp_msk1)
1583 )
1584 {
1585 /* The special case */
1586 b2 += Log2P;
1587 s2 += Log2P;
1588 spec_case = 1;
1589 }
1590 }
1591
1592 /* Arrange for convenient computation of quotients:
1593 * shift left if necessary so divisor has 4 leading 0 bits.
1594 *
1595 * Perhaps we should just compute leading 28 bits of S once
1596 * and for all and pass them and a shift to quorem, so it
1597 * can do shifts and ors to compute the numerator for q.
1598 */
1599 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
1600 {
1601 i = 32 - i;
1602 }
1603
1604 if (i > 4)
1605 {
1606 i -= 4;
1607 b2 += i;
1608 m2 += i;
1609 s2 += i;
1610 }
1611 else if (i < 4)
1612 {
1613 i += 28;
1614 b2 += i;
1615 m2 += i;
1616 s2 += i;
1617 }
1618
1619 if (b2 > 0)
1620 {
1621 b = lshift(aBlocks, freelist, b, b2);
1622 }
1623
1624 if (s2 > 0)
1625 {
1626 S = lshift(aBlocks, freelist, S, s2);
1627 }
1628
1629 if (k_check)
1630 {
1631 if (cmp(b,S) < 0)
1632 {
1633 k--;
1634 b = multadd(aBlocks, freelist, b, 10, 0); /* we botched the k estimate */
1635 if (leftright) { mhi = multadd(aBlocks, freelist, mhi, 10, 0); }
1636 ilim = ilim1;
1637 }
1638 }
1639 if (ilim <= 0 && (mode == 3 || mode == 5))
1640 {
1641 if (ilim < 0 || cmp(b,S = multadd(aBlocks, freelist, S,5,0)) <= 0)
1642 {
1643 /* no digits, fcvt style */
1644 no_digits:
1645 k = -1 - ndigits;
1646 goto ret;
1647 }
1648 one_digit:
1649 *s++ = '1';
1650 k++;
1651 goto ret;
1652 }
1653 if (leftright)
1654 {
1655 if (m2 > 0)
1656 {
1657 mhi = lshift(aBlocks, freelist, mhi, m2);
1658 }
1659
1660 /* Compute mlo -- check for special case
1661 * that d is a normalized power of 2.
1662 */
1663
1664 mlo = mhi;
1665 if (spec_case)
1666 {
1667 mhi = Balloc(aBlocks, freelist, mhi->k);
1668 Bcopy(mhi, mlo);
1669 mhi = lshift(aBlocks, freelist, mhi, Log2P);
1670 }
1671
1672 for(i = 1;;i++)
1673 {
1674 dig = quorem(b,S) + '0';
1675 /* Do we yet have the shortest decimal string
1676 * that will round to d?
1677 */
1678 j = cmp(b, mlo);
1679 delta = diff(aBlocks, freelist, S, mhi);
1680 j1 = delta->sign ? 1 : cmp(b, delta);
1681 Bfree(freelist, delta);
1682 #ifndef ROUND_BIASED
1683 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
1684 #ifdef Honor_FLT_ROUNDS
1685 && rounding >= 1
1686 #endif
1687 )
1688 {
1689 if (dig == '9')
1690 goto round_9_up;
1691 if (j > 0)
1692 dig++;
1693 *s++ = dig;
1694 goto ret;
1695 }
1696 #endif
1697 if (j < 0 || ((j == 0) && (mode != 1)
1698 #ifndef ROUND_BIASED
1699 && !(word1(d) & 1))
1700 #else
1701 )
1702 #endif
1703 )
1704 {
1705 if (!b->x[0] && b->wds <= 1)
1706 {
1707 goto accept_dig;
1708 }
1709 #ifdef Honor_FLT_ROUNDS
1710 if (mode > 1)
1711 {
1712 switch(rounding)
1713 {
1714 case 0: goto accept_dig;
1715 case 2: goto keep_dig;
1716 }
1717 }
1718 #endif /*Honor_FLT_ROUNDS*/
1719 if (j1 > 0)
1720 {
1721 b = lshift(aBlocks, freelist, b, 1);
1722 j1 = cmp(b, S);
1723 if ((j1 > 0 || ((j1 == 0) && (dig & 1))) && dig++ == '9')
1724 {
1725 goto round_9_up;
1726 }
1727 }
1728 accept_dig:
1729 *s++ = dig;
1730 goto ret;
1731 }
1732 if (j1 > 0)
1733 {
1734 #ifdef Honor_FLT_ROUNDS
1735 if (!rounding) { goto accept_dig; }
1736 #endif
1737 if (dig == '9')
1738 {
1739 /* possible if i == 1 */
1740 round_9_up:
1741 *s++ = '9';
1742 goto roundoff;
1743 }
1744
1745 *s++ = dig + 1;
1746 goto ret;
1747 }
1748 #ifdef Honor_FLT_ROUNDS
1749 keep_dig:
1750 #endif
1751 *s++ = dig;
1752 if (i == ilim)
1753 {
1754 break;
1755 }
1756
1757 b = multadd(aBlocks, freelist, b, 10, 0);
1758
1759 if (mlo == mhi)
1760 {
1761 mlo = mhi = multadd(aBlocks, freelist, mhi, 10, 0);
1762 }
1763 else
1764 {
1765 mlo = multadd(aBlocks, freelist, mlo, 10, 0);
1766 mhi = multadd(aBlocks, freelist, mhi, 10, 0);
1767 }
1768 }
1769 }
1770 else
1771 {
1772 for(i = 1;; i++)
1773 {
1774 *s++ = dig = quorem(b,S) + '0';
1775 if (!b->x[0] && b->wds <= 1)
1776 {
1777 goto ret;
1778 }
1779 if (i >= ilim)
1780 {
1781 break;
1782 }
1783 b = multadd(aBlocks, freelist, b, 10, 0);
1784 }
1785 }
1786 /* Round off last digit */
1787 #ifdef Honor_FLT_ROUNDS
1788 switch(rounding)
1789 {
1790 case 0: goto trimzeros;
1791 case 2: goto roundoff;
1792 }
1793 #endif
1794 b = lshift(aBlocks, freelist, b, 1);
1795 j = cmp(b, S);
1796 if (j > 0 || ((j == 0) && (dig & 1)))
1797 {
1798 roundoff:
1799 while(*--s == '9')
1800 {
1801 if (s == s0)
1802 {
1803 k++;
1804 *s++ = '1';
1805 goto ret;
1806 }
1807 }
1808 ++*s++;
1809 }
1810 else
1811 {
1812 #ifdef Honor_FLT_ROUNDS
1813 trimzeros:
1814 #endif
1815 while(*--s == '0') { ;; }
1816 s++;
1817 }
1818 ret:
1819 Bfree(freelist, S);
1820 if (mhi)
1821 {
1822 if (mlo && mlo != mhi) { Bfree(freelist, mlo); }
1823 Bfree(freelist, mhi);
1824 }
1825 ret1:
1826 Bfree(freelist, b);
1827 *s = 0;
1828 *decpt = k + 1;
1829 if (rve) { *rve = s; }
1830 return s0;
1831 }
1832
1833 // End.
1834