1 /* $NetBSD: strtodg.c,v 1.12 2013/04/19 10:41:53 joerg Exp $ */
3 /****************************************************************
5 The author of this software is David M. Gay.
7 Copyright (C) 1998-2001 by Lucent Technologies
10 Permission to use, copy, modify, and distribute this software and
11 its documentation for any purpose and without fee is hereby
12 granted, provided that the above copyright notice appear in all
13 copies and that both that the copyright notice and this
14 permission notice and warranty disclaimer appear in supporting
15 documentation, and that the name of Lucent or any of its entities
16 not be used in advertising or publicity pertaining to
17 distribution of the software without specific, written prior
20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
29 ****************************************************************/
31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
32 * with " at " changed at "@" and " dot " changed to "."). */
42 fivesbits
[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21,
43 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
50 increment(b
) Bigint
*b
;
65 if (*x
< (ULong
)0xffffffffL
) {
82 if (b
->wds
>= b
->maxwds
) {
97 decrement(b
) Bigint
*b
;
121 borrow
= (y
& 0x10000) >> 16;
123 } while(borrow
&& x
< xe
);
129 all_on(b
, n
) CONST Bigint
*b
; int n
;
131 all_on(CONST Bigint
*b
, int n
)
137 xe
= x
+ ((unsigned int)n
>> kshift
);
139 if ((*x
++ & ALL_ON
) != ALL_ON
)
142 return ((*x
| (ALL_ON
<< n
)) & ALL_ON
) == ALL_ON
;
148 set_ones(b
, n
) Bigint
*b
; int n
;
150 set_ones(Bigint
*b
, int n
)
156 k
= (unsigned int)(n
+ ((1 << kshift
) - 1)) >> kshift
;
163 k
= (unsigned int)n
>> kshift
;
172 x
[-1] >>= ULbits
- n
;
179 (d
, fpi
, expt
, bits
, exact
, rd
, irv
)
180 U
*d
; CONST FPI
*fpi
; Long
*expt
; ULong
*bits
; int exact
, rd
, *irv
;
182 (U
*d
, CONST FPI
*fpi
, Long
*expt
, ULong
*bits
, int exact
, int rd
, int *irv
)
186 ULong carry
, inex
, lostbits
;
187 int bdif
, e
, j
, k
, k1
, nb
, rv
;
190 b
= d2b(dval(d
), &e
, &bdif
);
191 bdif
-= nb
= fpi
->nbits
;
200 #ifndef IMPRECISE_INEXACT
213 case 1: /* round down (toward -Infinity) */
215 case 2: /* round up (toward +Infinity) */
217 default: /* round near */
226 if (b
->x
[(unsigned int)k
>>kshift
] & ((ULong
)1 << (k
& kmask
)))
230 /* "break" cases: round up 1 bit, then truncate; bdif > 0 */
235 if ( (lostbits
= any_on(b
, bdif
)) !=0)
236 inex
= STRTOG_Inexlo
;
239 inex
= STRTOG_Inexhi
;
241 if ( (j
= nb
& kmask
) !=0)
243 if (hi0bits(b
->x
[b
->wds
- 1]) != j
) {
245 lostbits
= b
->x
[0] & 1;
252 b
= lshift(b
, -bdif
);
256 if (k
> nb
|| fpi
->sudden_underflow
) {
258 *irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
262 if (k1
> 0 && !lostbits
)
263 lostbits
= any_on(b
, k1
);
264 if (!lostbits
&& !exact
)
267 carry
= b
->x
[(unsigned int)k1
>>kshift
] &
270 *irv
= STRTOG_Denormal
;
273 inex
= STRTOG_Inexhi
| STRTOG_Underflow
;
276 inex
= STRTOG_Inexlo
| STRTOG_Underflow
;
279 else if (e
> fpi
->emax
) {
281 *irv
= STRTOG_Infinite
| STRTOG_Overflow
| STRTOG_Inexhi
;
288 copybits(bits
, nb
, b
);
306 L
= word1(d
) << 16 | word1(d
) >> 16;
309 if ( (L
= word1(d
)) !=0)
311 return P
- lo0bits(&L
);
313 L
= word0(d
) << 16 | word0(d
) >> 16 | Exp_msk11
;
315 L
= word0(d
) | Exp_msk1
;
317 return P
- 32 - lo0bits(&L
);
322 strtodg(CONST
char *s00
, char **se
, CONST FPI
*fpi
, Long
*expt
, ULong
*bits
,
325 int abe
, abits
, asub
;
329 int bb0
, bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, denorm
;
330 int dsign
, e
, e1
, e2
, emin
, esign
, finished
, i
, inex
, irv
;
331 int j
, k
, nbits
, nd
, nd0
, nf
, nz
, nz0
, rd
, rvbits
, rve
, rve1
, sign
;
332 int sudden_underflow
= 0; /* pacify gcc */
333 CONST
char *s
, *s0
, *s1
;
338 Bigint
*ab
, *bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
, *rvb
, *rvb0
;
339 #ifdef USE_LOCALE /*{{*/
340 char *decimalpoint
= localeconv_l(loc
)->decimal_point
;
341 size_t dplen
= strlen(decimalpoint
);
342 #endif /*USE_LOCALE}}*/
344 e2
= 0; /* XXX gcc */
347 denorm
= sign
= nz0
= nz
= 0;
351 for(s
= s00
;;s
++) switch(*s
) {
361 irv
= STRTOG_NoNumber
;
380 irv
= gethex(&s
, fpi
, expt
, &rvb
, sign
, loc
);
381 if (irv
== STRTOG_NoNumber
) {
393 sudden_underflow
= fpi
->sudden_underflow
;
399 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
406 if (c
== *decimalpoint
) {
407 for(i
= 1; decimalpoint
[i
]; ++i
)
408 if (s
[i
] != decimalpoint
[i
])
420 for(; c
== '0'; c
= *++s
)
422 if (c
> '0' && c
<= '9') {
430 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
435 for(i
= 1; i
< nz
; i
++)
438 else if (nd
<= DBL_DIG
+ 1)
442 else if (nd
<= DBL_DIG
+ 1)
450 if (c
== 'e' || c
== 'E') {
451 if (!nd
&& !nz
&& !nz0
) {
452 irv
= STRTOG_NoNumber
;
465 if (c
>= '0' && c
<= '9') {
468 if (c
> '0' && c
<= '9') {
471 while((c
= *++s
) >= '0' && c
<= '9')
473 if (s
- s1
> 8 || L
> 19999)
474 /* Avoid confusion from exponents
475 * so large that e might overflow.
477 e
= 19999; /* safe for 16 bit ints */
492 /* Check for Nan and Infinity */
497 if (match(&s
,"nf")) {
499 if (!match(&s
,"inity"))
501 irv
= STRTOG_Infinite
;
507 if (match(&s
, "an")) {
509 *expt
= fpi
->emax
+ 1;
512 irv
= hexnan(&s
, fpi
, bits
);
517 #endif /* INFNAN_CHECK */
518 irv
= STRTOG_NoNumber
;
527 switch(fpi
->rounding
& 3) {
538 /* Now we have nd0 digits, starting at s0, followed by a
539 * decimal point, followed by nd-nd0 digits. The number we're
540 * after is the integer represented by those digits times
545 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
548 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
550 if (nbits
<= P
&& nd
<= DBL_DIG
) {
552 if (rvOK(&rv
, fpi
, expt
, bits
, 1, rd
, &irv
))
560 i
= fivesbits
[e
] + mantbits(&rv
) <= P
;
561 /* rv = */ rounded_product(dval(&rv
), tens
[e
]);
562 if (rvOK(&rv
, fpi
, expt
, bits
, i
, rd
, &irv
))
569 if (e
<= Ten_pmax
+ i
) {
570 /* A fancier test would sometimes let us do
571 * this for larger i values.
575 dval(&rv
) *= tens
[i
];
577 /* VAX exponent range is so narrow we must
578 * worry about overflow here...
581 dval(&adj
) = dval(&rv
);
582 word0(&adj
) -= P
*Exp_msk1
;
583 /* adj = */ rounded_product(dval(&adj
), tens
[e2
]);
584 if ((word0(&adj
) & Exp_mask
)
585 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
))
587 word0(&adj
) += P
*Exp_msk1
;
588 dval(&rv
) = dval(&adj
);
590 /* rv = */ rounded_product(dval(&rv
), tens
[e2
]);
592 if (rvOK(&rv
, fpi
, expt
, bits
, 0, rd
, &irv
))
597 #ifndef Inaccurate_Divide
598 else if (e
>= -Ten_pmax
) {
599 /* rv = */ rounded_quotient(dval(&rv
), tens
[-e
]);
600 if (rvOK(&rv
, fpi
, expt
, bits
, 0, rd
, &irv
))
609 /* Get starting approximation = rv * 10**e1 */
613 if ( (i
= e1
& 15) !=0)
614 dval(&rv
) *= tens
[i
];
616 e1
= (unsigned int)e1
>> 4;
617 while(e1
>= (1 << (n_bigtens
-1))) {
618 e2
+= ((word0(&rv
) & Exp_mask
)
619 >> Exp_shift1
) - Bias
;
620 word0(&rv
) &= ~Exp_mask
;
621 word0(&rv
) |= Bias
<< Exp_shift1
;
622 dval(&rv
) *= bigtens
[n_bigtens
-1];
623 e1
-= 1 << (n_bigtens
-1);
625 e2
+= ((word0(&rv
) & Exp_mask
) >> Exp_shift1
) - Bias
;
626 word0(&rv
) &= ~Exp_mask
;
627 word0(&rv
) |= Bias
<< Exp_shift1
;
628 for(j
= 0; e1
> 0; j
++, e1
= (unsigned int)e1
>> 1)
630 dval(&rv
) *= bigtens
[j
];
635 if ( (i
= e1
& 15) !=0)
636 dval(&rv
) /= tens
[i
];
638 e1
= (unsigned int)e1
>> 4;
639 while(e1
>= (1 << (n_bigtens
-1))) {
640 e2
+= ((word0(&rv
) & Exp_mask
)
641 >> Exp_shift1
) - Bias
;
642 word0(&rv
) &= ~Exp_mask
;
643 word0(&rv
) |= Bias
<< Exp_shift1
;
644 dval(&rv
) *= tinytens
[n_bigtens
-1];
645 e1
-= 1 << (n_bigtens
-1);
647 e2
+= ((word0(&rv
) & Exp_mask
) >> Exp_shift1
) - Bias
;
648 word0(&rv
) &= ~Exp_mask
;
649 word0(&rv
) |= Bias
<< Exp_shift1
;
650 for(j
= 0; e1
> 0; j
++, e1
= (unsigned int)e1
>> 1)
652 dval(&rv
) *= tinytens
[j
];
656 /* e2 is a correction to the (base 2) exponent of the return
657 * value, reflecting adjustments above to avoid overflow in the
658 * native arithmetic. For native IBM (base 16) arithmetic, we
659 * must multiply e2 by 4 to change from base 16 to 2.
663 rvb
= d2b(dval(&rv
), &rve
, &rvbits
); /* rv = rvb * 2^rve */
665 return STRTOG_NoMemory
;
667 if ((j
= rvbits
- nbits
) > 0) {
672 bb0
= 0; /* trailing zero bits in rvb */
673 e2
= rve
+ rvbits
- nbits
;
674 if (e2
> fpi
->emax
+ 1)
676 rve1
= rve
+ rvbits
- nbits
;
677 if (e2
< (emin
= fpi
->emin
)) {
681 rvb
= lshift(rvb
, j
);
692 irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
695 rvb
->x
[0] = rvb
->wds
= rvbits
= 1;
701 if (sudden_underflow
&& e2
+ 1 < emin
)
705 /* Now the hard part -- adjusting rv to the correct value.*/
707 /* Put digits into bd: true value = bd * 10^e */
709 bd0
= s2b(s0
, nd0
, nd
, y
, dplen
);
714 return STRTOG_NoMemory
;
718 return STRTOG_NoMemory
;
720 bbbits
= rvbits
- bb0
;
724 return STRTOG_NoMemory
;
739 j
= nbits
+ 1 - bbbits
;
740 i
= bbe
+ bbbits
- nbits
;
741 if (i
< emin
) /* denormal */
745 i
= bb2
< bd2
? bb2
: bd2
;
754 bs
= pow5mult(bs
, bb5
);
756 return STRTOG_NoMemory
;
759 return STRTOG_NoMemory
;
765 bb
= lshift(bb
, bb2
);
767 return STRTOG_NoMemory
;
772 bd
= pow5mult(bd
, bd5
);
774 return STRTOG_NoMemory
;
777 bd
= lshift(bd
, bd2
);
779 return STRTOG_NoMemory
;
782 bs
= lshift(bs
, bs2
);
784 return STRTOG_NoMemory
;
787 inex
= STRTOG_Inexhi
;
788 delta
= diff(bb
, bd
);
790 return STRTOG_NoMemory
;
791 if (delta
->wds
<= 1 && !delta
->x
[0])
794 delta
->sign
= finished
= 0;
799 if ( (finished
= dsign
^ (rd
&1)) !=0) {
801 irv
|= STRTOG_Inexhi
;
804 irv
|= STRTOG_Inexlo
;
807 for(i
= 0, j
= nbits
; j
>= ULbits
;
809 if (rvb
->x
[i
] & ALL_ON
)
812 if (j
> 1 && lo0bits(rvb
->x
+ i
) < j
- 1)
815 rvb
= set_ones(rvb
, rvbits
= nbits
);
817 return STRTOG_NoMemory
;
820 irv
|= dsign
? STRTOG_Inexlo
: STRTOG_Inexhi
;
824 /* Error is less than half an ulp -- check for
825 * special case of mantissa a power of two.
828 ? STRTOG_Normal
| STRTOG_Inexlo
829 : STRTOG_Normal
| STRTOG_Inexhi
;
830 if (dsign
|| bbbits
> 1 || denorm
|| rve1
== emin
)
832 delta
= lshift(delta
,1);
834 return STRTOG_NoMemory
;
835 if (cmp(delta
, bs
) > 0) {
836 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
842 /* exactly half-way between */
844 if (denorm
&& all_on(rvb
, rvbits
)) {
845 /*boundary case -- increment exponent*/
848 rve
= emin
+ nbits
- (rvbits
= 1);
849 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
853 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
855 else if (bbbits
== 1) {
858 /* boundary case -- decrement exponent */
860 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
861 if (rvb
->wds
== 1 && rvb
->x
[0] == 1)
862 sudden_underflow
= 1;
866 rvb
= set_ones(rvb
, rvbits
= nbits
);
868 return STRTOG_NoMemory
;
872 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
873 if ((bbbits
< nbits
&& !denorm
) || !(rvb
->x
[0] & 1))
876 rvb
= increment(rvb
);
878 return STRTOG_NoMemory
;
879 j
= kmask
& (ULbits
- (rvbits
& kmask
));
880 if (hi0bits(rvb
->x
[rvb
->wds
- 1]) != j
)
882 irv
= STRTOG_Normal
| STRTOG_Inexhi
;
888 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
892 if ((dval(&adj
) = ratio(delta
, bs
)) <= 2.) {
894 inex
= STRTOG_Inexlo
;
897 inex
= STRTOG_Inexhi
;
899 else if (denorm
&& bbbits
<= 1) {
903 irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
906 adj0
= dval(&adj
) = 1.;
909 adj0
= dval(&adj
) *= 0.5;
912 inex
= STRTOG_Inexlo
;
914 if (dval(&adj
) < 2147483647.) {
923 if (asub
&& adj0
> 0.)
927 if (!asub
&& adj0
> 0.) {
930 inex
= STRTOG_Inexact
- inex
;
938 /* adj *= ulp(dval(&rv)); */
939 /* if (asub) rv -= adj; else rv += adj; */
941 if (!denorm
&& rvbits
< nbits
) {
942 rvb
= lshift(rvb
, j
= nbits
- rvbits
);
944 return STRTOG_NoMemory
;
948 ab
= d2b(dval(&adj
), &abe
, &abits
);
950 return STRTOG_NoMemory
;
954 ab
= lshift(ab
, abe
);
958 j
= hi0bits(rvb
->x
[rvb
->wds
-1]);
961 return STRTOG_NoMemory
;
965 else if (rvb
->wds
<= k
966 || hi0bits( rvb
->x
[k
]) >
967 hi0bits(rvb0
->x
[k
])) {
968 /* unlikely; can only have lost 1 high bit */
974 rvb
= lshift(rvb
, 1);
976 return STRTOG_NoMemory
;
986 return STRTOG_NoMemory
;
989 || hi0bits(rvb
->x
[k
]) < hi0bits(rvb0
->x
[k
])) {
991 if (++rvbits
== nbits
)
1009 /* Can we stop now? */
1010 tol
= dval(&adj
) * 5e-16; /* > max rel error */
1011 dval(&adj
) = adj0
- .5;
1012 if (dval(&adj
) < -tol
) {
1018 else if (dval(&adj
) > tol
&& adj0
< 1. - tol
) {
1023 bb0
= denorm
? 0 : trailz(rvb
);
1029 if (!denorm
&& (j
= nbits
- rvbits
)) {
1031 rvb
= lshift(rvb
, j
);
1042 if (rve
> fpi
->emax
) {
1043 switch(fpi
->rounding
& 3) {
1044 case FPI_Round_near
:
1050 case FPI_Round_down
:
1054 /* Round to largest representable magnitude */
1057 irv
= STRTOG_Normal
| STRTOG_Inexlo
;
1060 be
= b
+ ((unsigned int)(fpi
->nbits
+ 31) >> 5);
1062 *b
++ = (unsigned int)-1;
1063 if ((j
= fpi
->nbits
& 0x1f) != 0)
1068 irv
= STRTOG_Infinite
| STRTOG_Overflow
| STRTOG_Inexhi
;
1075 *expt
= fpi
->emax
+ 1;
1079 if (sudden_underflow
) {
1081 irv
= STRTOG_Underflow
| STRTOG_Inexlo
;
1087 irv
= (irv
& ~STRTOG_Retmask
) |
1088 (rvb
->wds
> 0 ? STRTOG_Denormal
: STRTOG_Zero
);
1089 if (irv
& STRTOG_Inexact
) {
1090 irv
|= STRTOG_Underflow
;
1102 copybits(bits
, nbits
, rvb
);