1 /****************************************************************
3 The author of this software is David M. Gay.
5 Copyright (C) 1998-2001 by Lucent Technologies
8 Permission to use, copy, modify, and distribute this software and
9 its documentation for any purpose and without fee is hereby
10 granted, provided that the above copyright notice appear in all
11 copies and that both that the copyright notice and this
12 permission notice and warranty disclaimer appear in supporting
13 documentation, and that the name of Lucent or any of its entities
14 not be used in advertising or publicity pertaining to
15 distribution of the software without specific, written prior
18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
27 ****************************************************************/
29 /* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
39 fivesbits
[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21,
40 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
49 increment(b
) Bigint
*b
;
64 if (*x
< (ULong
)0xffffffffL
) {
81 if (b
->wds
>= b
->maxwds
) {
94 decrement(b
) Bigint
*b
;
118 borrow
= (y
& 0x10000) >> 16;
120 } while(borrow
&& x
< xe
);
126 all_on(b
, n
) Bigint
*b
; int n
;
128 all_on(Bigint
*b
, int n
)
134 xe
= x
+ (n
>> kshift
);
136 if ((*x
++ & ALL_ON
) != ALL_ON
)
139 return ((*x
| (ALL_ON
<< n
)) & ALL_ON
) == ALL_ON
;
145 set_ones(b
, n
) Bigint
*b
; int n
;
147 set_ones(Bigint
*b
, int n
)
153 k
= (n
+ ((1 << kshift
) - 1)) >> kshift
;
167 x
[-1] >>= ULbits
- n
;
174 (d
, fpi
, exp
, bits
, exact
, rd
, irv
)
175 U
*d
; FPI
*fpi
; Long
*exp
; ULong
*bits
; int exact
, rd
, *irv
;
177 (U
*d
, FPI
*fpi
, Long
*exp
, ULong
*bits
, int exact
, int rd
, int *irv
)
181 ULong carry
, inex
, lostbits
;
182 int bdif
, e
, j
, k
, k1
, nb
, rv
;
185 b
= d2b(dval(d
), &e
, &bdif
);
186 bdif
-= nb
= fpi
->nbits
;
195 #ifndef IMPRECISE_INEXACT
208 case 1: /* round down (toward -Infinity) */
210 case 2: /* round up (toward +Infinity) */
212 default: /* round near */
223 if (b
->x
[k
>>kshift
] & ((ULong
)1 << (k
& kmask
)))
227 /* "break" cases: round up 1 bit, then truncate; bdif > 0 */
232 if ( (lostbits
= any_on(b
, bdif
)) !=0)
233 inex
= STRTOG_Inexlo
;
236 inex
= STRTOG_Inexhi
;
238 if ( (j
= nb
& kmask
) !=0)
240 if (hi0bits(b
->x
[b
->wds
- 1]) != j
) {
242 lostbits
= b
->x
[0] & 1;
249 b
= lshift(b
, -bdif
);
253 if (k
> nb
|| fpi
->sudden_underflow
) {
255 *irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
259 if (k1
> 0 && !lostbits
)
260 lostbits
= any_on(b
, k1
);
261 if (!lostbits
&& !exact
)
264 carry
= b
->x
[k1
>>kshift
] & (1 << (k1
& kmask
));
266 *irv
= STRTOG_Denormal
;
269 inex
= STRTOG_Inexhi
| STRTOG_Underflow
;
272 inex
= STRTOG_Inexlo
| STRTOG_Underflow
;
275 else if (e
> fpi
->emax
) {
277 *irv
= STRTOG_Infinite
| STRTOG_Overflow
| STRTOG_Inexhi
;
284 copybits(bits
, nb
, b
);
301 L
= word1(d
) << 16 | word1(d
) >> 16;
304 if ( (L
= word1(d
)) !=0)
306 return P
- lo0bits(&L
);
308 L
= word0(d
) << 16 | word0(d
) >> 16 | Exp_msk11
;
310 L
= word0(d
) | Exp_msk1
;
312 return P
- 32 - lo0bits(&L
);
318 (s00
, se
, fpi
, exp
, bits
)
319 CONST
char *s00
; char **se
; FPI
*fpi
; Long
*exp
; ULong
*bits
;
321 (CONST
char *s00
, char **se
, FPI
*fpi
, Long
*exp
, ULong
*bits
)
324 int abe
, abits
, asub
;
325 int bb0
, bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, decpt
, denorm
;
326 int dsign
, e
, e1
, e2
, emin
, esign
, finished
, i
, inex
, irv
;
327 int j
, k
, nbits
, nd
, nd0
, nf
, nz
, nz0
, rd
, rvbits
, rve
, rve1
, sign
;
328 int sudden_underflow
;
329 CONST
char *s
, *s0
, *s1
;
334 Bigint
*ab
, *bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
, *rvb
, *rvb0
;
335 #ifdef USE_LOCALE /*{{*/
336 #ifdef NO_LOCALE_CACHE
337 char *decimalpoint
= localeconv()->decimal_point
;
338 int dplen
= strlen(decimalpoint
);
341 static char *decimalpoint_cache
;
343 if (!(s0
= decimalpoint_cache
)) {
344 s0
= localeconv()->decimal_point
;
345 if ((decimalpoint_cache
= (char*)MALLOC(strlen(s0
) + 1))) {
346 strcpy(decimalpoint_cache
, s0
);
347 s0
= decimalpoint_cache
;
351 decimalpoint
= (char*)s0
;
352 #endif /*NO_LOCALE_CACHE*/
353 #else /*USE_LOCALE}{*/
355 #endif /*USE_LOCALE}}*/
358 denorm
= sign
= nz0
= nz
= 0;
362 for(s
= s00
;;s
++) switch(*s
) {
372 irv
= STRTOG_NoNumber
;
391 irv
= gethex(&s
, fpi
, exp
, &rvb
, sign
);
392 if (irv
== STRTOG_NoNumber
) {
404 sudden_underflow
= fpi
->sudden_underflow
;
407 for(decpt
= nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
414 if (c
== *decimalpoint
) {
415 for(i
= 1; decimalpoint
[i
]; ++i
)
416 if (s
[i
] != decimalpoint
[i
])
426 for(; c
== '0'; c
= *++s
)
428 if (c
> '0' && c
<= '9') {
436 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
441 for(i
= 1; i
< nz
; i
++)
444 else if (nd
<= DBL_DIG
+ 1)
448 else if (nd
<= DBL_DIG
+ 1)
456 if (c
== 'e' || c
== 'E') {
457 if (!nd
&& !nz
&& !nz0
) {
458 irv
= STRTOG_NoNumber
;
470 if (c
>= '0' && c
<= '9') {
473 if (c
> '0' && c
<= '9') {
476 while((c
= *++s
) >= '0' && c
<= '9')
478 if (s
- s1
> 8 || L
> 19999)
479 /* Avoid confusion from exponents
480 * so large that e might overflow.
482 e
= 19999; /* safe for 16 bit ints */
497 /* Check for Nan and Infinity */
502 if (match(&s
,"nf")) {
504 if (!match(&s
,"inity"))
506 irv
= STRTOG_Infinite
;
512 if (match(&s
, "an")) {
514 *exp
= fpi
->emax
+ 1;
517 irv
= hexnan(&s
, fpi
, bits
);
522 #endif /* INFNAN_CHECK */
523 irv
= STRTOG_NoNumber
;
532 switch(fpi
->rounding
& 3) {
543 /* Now we have nd0 digits, starting at s0, followed by a
544 * decimal point, followed by nd-nd0 digits. The number we're
545 * after is the integer represented by those digits times
550 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
553 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
555 if (nbits
<= P
&& nd
<= DBL_DIG
) {
557 if (rvOK(&rv
, fpi
, exp
, bits
, 1, rd
, &irv
))
565 i
= fivesbits
[e
] + mantbits(&rv
) <= P
;
566 /* rv = */ rounded_product(dval(&rv
), tens
[e
]);
567 if (rvOK(&rv
, fpi
, exp
, bits
, i
, rd
, &irv
))
574 if (e
<= Ten_pmax
+ i
) {
575 /* A fancier test would sometimes let us do
576 * this for larger i values.
580 dval(&rv
) *= tens
[i
];
582 /* VAX exponent range is so narrow we must
583 * worry about overflow here...
586 dval(&adj
) = dval(&rv
);
587 word0(&adj
) -= P
*Exp_msk1
;
588 /* adj = */ rounded_product(dval(&adj
), tens
[e2
]);
589 if ((word0(&adj
) & Exp_mask
)
590 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
))
592 word0(&adj
) += P
*Exp_msk1
;
593 dval(&rv
) = dval(&adj
);
595 /* rv = */ rounded_product(dval(&rv
), tens
[e2
]);
597 if (rvOK(&rv
, fpi
, exp
, bits
, 0, rd
, &irv
))
602 #ifndef Inaccurate_Divide
603 else if (e
>= -Ten_pmax
) {
604 /* rv = */ rounded_quotient(dval(&rv
), tens
[-e
]);
605 if (rvOK(&rv
, fpi
, exp
, bits
, 0, rd
, &irv
))
614 /* Get starting approximation = rv * 10**e1 */
618 if ( (i
= e1
& 15) !=0)
619 dval(&rv
) *= tens
[i
];
622 while(e1
>= (1 << (n_bigtens
-1))) {
623 e2
+= ((word0(&rv
) & Exp_mask
)
624 >> Exp_shift1
) - Bias
;
625 word0(&rv
) &= ~Exp_mask
;
626 word0(&rv
) |= Bias
<< Exp_shift1
;
627 dval(&rv
) *= bigtens
[n_bigtens
-1];
628 e1
-= 1 << (n_bigtens
-1);
630 e2
+= ((word0(&rv
) & Exp_mask
) >> Exp_shift1
) - Bias
;
631 word0(&rv
) &= ~Exp_mask
;
632 word0(&rv
) |= Bias
<< Exp_shift1
;
633 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
635 dval(&rv
) *= bigtens
[j
];
640 if ( (i
= e1
& 15) !=0)
641 dval(&rv
) /= tens
[i
];
644 while(e1
>= (1 << (n_bigtens
-1))) {
645 e2
+= ((word0(&rv
) & Exp_mask
)
646 >> Exp_shift1
) - Bias
;
647 word0(&rv
) &= ~Exp_mask
;
648 word0(&rv
) |= Bias
<< Exp_shift1
;
649 dval(&rv
) *= tinytens
[n_bigtens
-1];
650 e1
-= 1 << (n_bigtens
-1);
652 e2
+= ((word0(&rv
) & Exp_mask
) >> Exp_shift1
) - Bias
;
653 word0(&rv
) &= ~Exp_mask
;
654 word0(&rv
) |= Bias
<< Exp_shift1
;
655 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
657 dval(&rv
) *= tinytens
[j
];
661 /* e2 is a correction to the (base 2) exponent of the return
662 * value, reflecting adjustments above to avoid overflow in the
663 * native arithmetic. For native IBM (base 16) arithmetic, we
664 * must multiply e2 by 4 to change from base 16 to 2.
668 rvb
= d2b(dval(&rv
), &rve
, &rvbits
); /* rv = rvb * 2^rve */
670 if ((j
= rvbits
- nbits
) > 0) {
675 bb0
= 0; /* trailing zero bits in rvb */
676 e2
= rve
+ rvbits
- nbits
;
677 if (e2
> fpi
->emax
+ 1)
679 rve1
= rve
+ rvbits
- nbits
;
680 if (e2
< (emin
= fpi
->emin
)) {
684 rvb
= lshift(rvb
, j
);
695 irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
698 rvb
->x
[0] = rvb
->wds
= rvbits
= 1;
704 if (sudden_underflow
&& e2
+ 1 < emin
)
708 /* Now the hard part -- adjusting rv to the correct value.*/
710 /* Put digits into bd: true value = bd * 10^e */
712 bd0
= s2b(s0
, nd0
, nd
, y
, dplen
);
719 bbbits
= rvbits
- bb0
;
736 j
= nbits
+ 1 - bbbits
;
737 i
= bbe
+ bbbits
- nbits
;
738 if (i
< emin
) /* denormal */
742 i
= bb2
< bd2
? bb2
: bd2
;
751 bs
= pow5mult(bs
, bb5
);
758 bb
= lshift(bb
, bb2
);
762 bd
= pow5mult(bd
, bd5
);
764 bd
= lshift(bd
, bd2
);
766 bs
= lshift(bs
, bs2
);
768 inex
= STRTOG_Inexhi
;
769 delta
= diff(bb
, bd
);
770 if (delta
->wds
<= 1 && !delta
->x
[0])
773 delta
->sign
= finished
= 0;
778 if ( (finished
= dsign
^ (rd
&1)) !=0) {
780 irv
|= STRTOG_Inexhi
;
783 irv
|= STRTOG_Inexlo
;
786 for(i
= 0, j
= nbits
; j
>= ULbits
;
788 if (rvb
->x
[i
] & ALL_ON
)
791 if (j
> 1 && lo0bits(rvb
->x
+ i
) < j
- 1)
794 rvb
= set_ones(rvb
, rvbits
= nbits
);
797 irv
|= dsign
? STRTOG_Inexlo
: STRTOG_Inexhi
;
801 /* Error is less than half an ulp -- check for
802 * special case of mantissa a power of two.
805 ? STRTOG_Normal
| STRTOG_Inexlo
806 : STRTOG_Normal
| STRTOG_Inexhi
;
807 if (dsign
|| bbbits
> 1 || denorm
|| rve1
== emin
)
809 delta
= lshift(delta
,1);
810 if (cmp(delta
, bs
) > 0) {
811 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
817 /* exactly half-way between */
819 if (denorm
&& all_on(rvb
, rvbits
)) {
820 /*boundary case -- increment exponent*/
823 rve
= emin
+ nbits
- (rvbits
= 1);
824 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
828 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
830 else if (bbbits
== 1) {
833 /* boundary case -- decrement exponent */
835 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
836 if (rvb
->wds
== 1 && rvb
->x
[0] == 1)
837 sudden_underflow
= 1;
841 rvb
= set_ones(rvb
, rvbits
= nbits
);
845 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
846 if ((bbbits
< nbits
&& !denorm
) || !(rvb
->x
[0] & 1))
849 rvb
= increment(rvb
);
850 j
= kmask
& (ULbits
- (rvbits
& kmask
));
851 if (hi0bits(rvb
->x
[rvb
->wds
- 1]) != j
)
853 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
859 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
863 if ((dval(&adj
) = ratio(delta
, bs
)) <= 2.) {
865 inex
= STRTOG_Inexlo
;
868 inex
= STRTOG_Inexhi
;
870 else if (denorm
&& bbbits
<= 1) {
874 irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
877 adj0
= dval(&adj
) = 1.;
880 adj0
= dval(&adj
) *= 0.5;
883 inex
= STRTOG_Inexlo
;
885 if (dval(&adj
) < 2147483647.) {
894 if (asub
&& adj0
> 0.)
898 if (!asub
&& adj0
> 0.) {
901 inex
= STRTOG_Inexact
- inex
;
909 /* adj *= ulp(dval(&rv)); */
910 /* if (asub) rv -= adj; else rv += adj; */
912 if (!denorm
&& rvbits
< nbits
) {
913 rvb
= lshift(rvb
, j
= nbits
- rvbits
);
917 ab
= d2b(dval(&adj
), &abe
, &abits
);
921 ab
= lshift(ab
, abe
);
925 j
= hi0bits(rvb
->x
[rvb
->wds
-1]);
930 else if (rvb
->wds
<= k
931 || hi0bits( rvb
->x
[k
]) >
932 hi0bits(rvb0
->x
[k
])) {
933 /* unlikely; can only have lost 1 high bit */
939 rvb
= lshift(rvb
, 1);
950 || hi0bits(rvb
->x
[k
]) < hi0bits(rvb0
->x
[k
])) {
952 if (++rvbits
== nbits
)
970 /* Can we stop now? */
971 tol
= dval(&adj
) * 5e-16; /* > max rel error */
972 dval(&adj
) = adj0
- .5;
973 if (dval(&adj
) < -tol
) {
979 else if (dval(&adj
) > tol
&& adj0
< 1. - tol
) {
984 bb0
= denorm
? 0 : trailz(rvb
);
990 if (!denorm
&& (j
= nbits
- rvbits
)) {
992 rvb
= lshift(rvb
, j
);
1003 if (rve
> fpi
->emax
) {
1004 switch(fpi
->rounding
& 3) {
1005 case FPI_Round_near
:
1011 case FPI_Round_down
:
1015 /* Round to largest representable magnitude */
1018 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
1021 be
= b
+ ((fpi
->nbits
+ 31) >> 5);
1024 if ((j
= fpi
->nbits
& 0x1f))
1029 irv
= STRTOG_Infinite
| STRTOG_Overflow
| STRTOG_Inexhi
;
1034 *exp
= fpi
->emax
+ 1;
1038 if (sudden_underflow
) {
1040 irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
1046 irv
= (irv
& ~STRTOG_Retmask
) |
1047 (rvb
->wds
> 0 ? STRTOG_Denormal
: STRTOG_Zero
);
1048 if (irv
& STRTOG_Inexact
) {
1049 irv
|= STRTOG_Underflow
;
1061 copybits(bits
, nbits
, rvb
);