10 #include "floatscan.h"
12 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
15 #define LD_B1B_MAX 9007199, 254740991
18 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
21 #define LD_B1B_MAX 18, 446744073, 709551615
24 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
27 #define LD_B1B_MAX 10384593, 717069655, 257060992, 658440191
31 #error Unsupported long double representation
36 #define CONCAT2(x,y) x ## y
37 #define CONCAT(x,y) CONCAT2(x,y)
39 static long long scanexp(FILE *f
, int pok
)
47 if (c
=='+' || c
=='-') {
50 if (c
-'0'>=10U && pok
) shunget(f
);
56 for (x
=0; c
-'0'<10U && x
<INT_MAX
/10; c
= shgetc(f
))
58 for (y
=x
; c
-'0'<10U && y
<LLONG_MAX
/100; c
= shgetc(f
))
60 for (; c
-'0'<10U; c
= shgetc(f
));
66 static long double decfloat(FILE *f
, int c
, int bits
, int emin
, int sign
, int pok
)
69 static const uint32_t th
[] = { LD_B1B_MAX
};
71 long long lrp
=0, dc
=0;
74 int gotdig
= 0, gotrad
= 0;
77 int emax
= -emin
-bits
+3;
82 static const int p10s
[] = { 10, 100, 1000, 10000,
83 100000, 1000000, 10000000, 100000000 };
88 /* Don't let leading zeros consume buffer space */
89 for (; c
=='0'; c
= shgetc(f
)) gotdig
=1;
92 for (c
= shgetc(f
); c
=='0'; c
= shgetc(f
)) gotdig
=1, lrp
--;
96 for (; c
-'0'<10U || c
=='.'; c
= shgetc(f
)) {
101 } else if (k
< KMAX
-3) {
103 if (c
!='0') lnz
= dc
;
104 if (j
) x
[k
] = x
[k
]*10 + c
-'0';
121 if (gotdig
&& (c
|32)=='e') {
122 e10
= scanexp(f
, pok
);
123 if (e10
== LLONG_MIN
) {
142 /* Handle zero specially to avoid nasty special cases later */
143 if (!x
[0]) return sign
* 0.0;
145 /* Optimize small integers (w/no exponent) and over/under-flow */
146 if (lrp
==dc
&& dc
<10 && (bits
>30 || x
[0]>>bits
==0))
147 return sign
* (long double)x
[0];
150 return sign
* LDBL_MAX
* LDBL_MAX
;
152 if (lrp
< emin
-2*LDBL_MANT_DIG
) {
154 return sign
* LDBL_MIN
* LDBL_MIN
;
157 /* Align incomplete final B1B digit */
159 for (; j
<9; j
++) x
[k
]*=10;
169 /* Optimize small to mid-size integers (even in exp. notation) */
170 if (lnz
<9 && lnz
<=rp
&& rp
< 18) {
171 if (rp
== 9) return sign
* (long double)x
[0];
172 if (rp
< 9) return sign
* (long double)x
[0] / p10s
[8-rp
];
173 int bitlim
= bits
-3*(int)(rp
-9);
174 if (bitlim
>30 || x
[0]>>bitlim
==0)
175 return sign
* (long double)x
[0] * p10s
[rp
-10];
178 /* Drop trailing zeros */
179 for (; !x
[z
-1]; z
--);
181 /* Align radix point to B1B digit boundary */
183 int rpm9
= rp
>=0 ? rp
%9 : rp
%9+9;
184 int p10
= p10s
[8-rpm9
];
186 for (k
=a
; k
!=z
; k
++) {
187 uint32_t tmp
= x
[k
] % p10
;
188 x
[k
] = x
[k
]/p10
+ carry
;
189 carry
= 1000000000/p10
* tmp
;
195 if (carry
) x
[z
++] = carry
;
199 /* Upscale until desired number of bits are left of radix point */
200 while (rp
< 9*LD_B1B_DIG
|| (rp
== 9*LD_B1B_DIG
&& x
[a
]<th
[0])) {
203 for (k
=(z
-1 & MASK
); ; k
=(k
-1 & MASK
)) {
204 uint64_t tmp
= ((uint64_t)x
[k
] << 29) + carry
;
205 if (tmp
> 1000000000) {
206 carry
= tmp
/ 1000000000;
207 x
[k
] = tmp
% 1000000000;
212 if (k
==(z
-1 & MASK
) && k
!=a
&& !x
[k
]) z
= k
;
220 x
[z
-1 & MASK
] |= x
[z
];
226 /* Downscale until exactly number of bits are left of radix point */
230 for (i
=0; i
<LD_B1B_DIG
; i
++) {
232 if (k
== z
|| x
[k
] < th
[i
]) {
236 if (x
[a
+i
& MASK
] > th
[i
]) break;
238 if (i
==LD_B1B_DIG
&& rp
==9*LD_B1B_DIG
) break;
239 /* FIXME: find a way to compute optimal sh */
240 if (rp
> 9+9*LD_B1B_DIG
) sh
= 9;
242 for (k
=a
; k
!=z
; k
=(k
+1 & MASK
)) {
243 uint32_t tmp
= x
[k
] & (1<<sh
)-1;
244 x
[k
] = (x
[k
]>>sh
) + carry
;
245 carry
= (1000000000>>sh
) * tmp
;
253 if ((z
+1 & MASK
) != a
) {
256 } else x
[z
-1 & MASK
] |= 1;
260 /* Assemble desired bits into floating point variable */
261 for (y
=i
=0; i
<LD_B1B_DIG
; i
++) {
262 if ((a
+i
& MASK
)==z
) x
[(z
=(z
+1 & MASK
))-1] = 0;
263 y
= 1000000000.0L * y
+ x
[a
+i
& MASK
];
268 /* Limit precision for denormal results */
269 if (bits
> LDBL_MANT_DIG
+e2
-emin
) {
270 bits
= LDBL_MANT_DIG
+e2
-emin
;
275 /* Calculate bias term to force rounding, move out lower bits */
276 if (bits
< LDBL_MANT_DIG
) {
277 bias
= copysignl(scalbn(1, 2*LDBL_MANT_DIG
-bits
-1), y
);
278 frac
= fmodl(y
, scalbn(1, LDBL_MANT_DIG
-bits
));
283 /* Process tail of decimal input so it can affect rounding */
284 if ((a
+i
& MASK
) != z
) {
285 uint32_t t
= x
[a
+i
& MASK
];
286 if (t
< 500000000 && (t
|| (a
+i
+1 & MASK
) != z
))
288 else if (t
> 500000000)
290 else if (t
== 500000000) {
291 if ((a
+i
+1 & MASK
) == z
)
296 if (LDBL_MANT_DIG
-bits
>= 2 && !fmodl(frac
, 1))
303 if ((e2
+LDBL_MANT_DIG
& INT_MAX
) > emax
-5) {
304 if (fabs(y
) >= CONCAT(0x1p
, LDBL_MANT_DIG
)) {
305 if (denormal
&& bits
==LDBL_MANT_DIG
+e2
-emin
)
310 if (e2
+LDBL_MANT_DIG
>emax
|| (denormal
&& frac
))
314 return scalbnl(y
, e2
);
317 static long double hexfloat(FILE *f
, int bits
, int emin
, int sign
, int pok
)
321 long double scale
= 1;
322 long double bias
= 0;
323 int gottail
= 0, gotrad
= 0, gotdig
= 0;
332 /* Skip leading zeros */
333 for (; c
=='0'; c
= shgetc(f
)) gotdig
= 1;
338 /* Count zeros after the radix point before significand */
339 for (rp
=0; c
=='0'; c
= shgetc(f
), rp
--) gotdig
= 1;
342 for (; c
-'0'<10U || (c
|32)-'a'<6U || c
=='.'; c
= shgetc(f
)) {
349 if (c
> '9') d
= (c
|32)+10-'a';
353 } else if (dc
< LDBL_MANT_DIG
/4+1) {
355 } else if (d
&& !gottail
) {
366 if (gotrad
) shunget(f
);
372 if (!gotrad
) rp
= dc
;
373 while (dc
<8) x
*= 16, dc
++;
375 e2
= scanexp(f
, pok
);
376 if (e2
== LLONG_MIN
) {
390 if (!x
) return sign
* 0.0;
393 return sign
* LDBL_MAX
* LDBL_MAX
;
395 if (e2
< emin
-2*LDBL_MANT_DIG
) {
397 return sign
* LDBL_MIN
* LDBL_MIN
;
400 while (x
< 0x80000000) {
411 if (bits
> 32+e2
-emin
) {
416 if (bits
< LDBL_MANT_DIG
)
417 bias
= copysignl(scalbn(1, 32+LDBL_MANT_DIG
-bits
-1), sign
);
419 if (bits
<32 && y
&& !(x
&1)) x
++, y
=0;
421 y
= bias
+ sign
*(long double)x
+ sign
*y
;
424 if (!y
) errno
= ERANGE
;
426 return scalbnl(y
, e2
);
429 long double __floatscan(FILE *f
, int prec
, int pok
)
440 emin
= FLT_MIN_EXP
-bits
;
444 emin
= DBL_MIN_EXP
-bits
;
447 bits
= LDBL_MANT_DIG
;
448 emin
= LDBL_MIN_EXP
-bits
;
454 while (isspace((c
=shgetc(f
))));
456 if (c
=='+' || c
=='-') {
461 for (i
=0; i
<8 && (c
|32)=="infinity"[i
]; i
++)
462 if (i
<7) c
= shgetc(f
);
463 if (i
==3 || i
==8 || (i
>3 && pok
)) {
466 if (pok
) for (; i
>3; i
--) shunget(f
);
468 return sign
* INFINITY
;
470 if (!i
) for (i
=0; i
<3 && (c
|32)=="nan"[i
]; i
++)
471 if (i
<2) c
= shgetc(f
);
473 if (shgetc(f
) != '(') {
479 if (c
-'0'<10U || c
-'A'<26U || c
-'a'<26U || c
=='_')
481 if (c
==')') return NAN
;
488 while (i
--) shunget(f
);
504 return hexfloat(f
, bits
, emin
, sign
, pok
);
509 return decfloat(f
, c
, bits
, emin
, sign
, pok
);