1
1
using Base. MPFR, Base. GMP, Base. GMP. MPZ
2
2
3
3
_widen (x:: UInt64 ) = UInt128 (x)
4
+ _widen (x:: Int64 ) = Int128 (x)
4
5
5
6
const BIGINT = BigInt[]
6
7
@@ -17,10 +18,10 @@ function access_threaded(f, v::Vector)
17
18
end
18
19
@noinline _length_assert () = @assert false " 0 < tid <= v"
19
20
20
- function _widen (v:: UInt128 )
21
+ function _widen (v:: T ) where {T <: Union{Int128, UInt128} }
21
22
x = access_threaded (() -> (@static VERSION > v " 1.5" ? BigInt (; nbits= 256 ) : BigInt ()), BIGINT)
22
23
ccall ((:__gmpz_import , :libgmp ), Int32,
23
- (Ref{BigInt}, Csize_t, Cint, Csize_t, Cint, Csize_t, Ref{UInt128 }),
24
+ (Ref{BigInt}, Csize_t, Cint, Csize_t, Cint, Csize_t, Ref{T }),
24
25
x, 1 , 1 , 16 , 0 , 0 , v)
25
26
return x
26
27
end
29
30
maxdigits (:: Type{Float64} ) = 1079
30
31
maxdigits (:: Type{Float32} ) = 154
31
32
maxdigits (:: Type{Float16} ) = 29
32
- maxdigits (:: Type{BigFloat} ) = typemax (Int64)
33
+ maxdigits (T ) = typemax (Int64)
33
34
34
35
ten (:: Type{T} ) where {T} = T (10 )
35
36
const TEN = BigInt (10 )
@@ -43,12 +44,46 @@ function _muladd(ten, digits::BigInt, b)
43
44
return digits
44
45
end
45
46
47
+ @enum FloatType FLOAT16 FLOAT32 FLOAT64 BIGFLOAT
48
+ float_type (:: Type{T} , FT:: FloatType ) where {T <: AbstractFloat } = T
49
+ float_type (T, FT:: FloatType ) = FT === FLOAT16 ? Float16 :
50
+ FT === FLOAT32 ? Float32 :
51
+ FT === FLOAT64 ? Float64 :
52
+ BigFloat
53
+
46
54
# for non SupportedFloat Reals, parse as Float64, then convert
47
55
@inline function typeparser (:: Type{T} , source, pos, len, b, code, pl, options) where {T <: Real }
48
56
pos, code, pl, x = typeparser (Float64, source, pos, len, b, code, pl, options)
49
57
return pos, code, pl, T (x)
50
58
end
51
59
60
+ function typeparser (:: Type{BigFloat} , source, pos, len, b, code, pl, options)
61
+ base = 0
62
+ rounding = Base. MPFR. ROUNDING_MODE[]
63
+ z = BigFloat (precision= Base. MPFR. DEFAULT_PRECISION[])
64
+ if source isa AbstractVector{UInt8}
65
+ str = source
66
+ strpos = pos
67
+ else
68
+ _, _, _pl, _ = typeparser (String, source, pos, len, b, code, pl, options)
69
+ _pos = position (source)
70
+ vpos, vlen = _pl. pos, _pl. len
71
+ fastseek! (source, vpos - 1 )
72
+ str = Base. StringVector (vlen)
73
+ strpos = 1
74
+ readbytes! (source, str, vlen)
75
+ fastseek! (source, _pos) # reset IO to earlier position
76
+ end
77
+ GC. @preserve str begin
78
+ ptr = pointer (str, strpos)
79
+ endptr = Ref {Ptr{UInt8}} ()
80
+ err = ccall ((:mpfr_strtofr , :libmpfr ), Int32, (Ref{BigFloat}, Ptr{UInt8}, Ref{Ptr{UInt8}}, Int32, Base. MPFR. MPFRRoundingMode), z, ptr, endptr, base, rounding)
81
+ code |= endptr[] == ptr ? INVALID : OK
82
+ pos += Int (endptr[] - ptr)
83
+ return pos, code, PosLen (pl. pos, pos - pl. pos), z
84
+ end
85
+ end
86
+
52
87
@inline function typeparser (:: Type{T} , source, pos, len, b, code, pl, options) where {T <: SupportedFloats }
53
88
# keep track of starting pos in case of invalid, we can rewind to start of parsing
54
89
startpos = pos
89
124
if eof (source, pos, len)
90
125
code |= EOF
91
126
end
92
- code |= OK
127
+ code |= OK | SPECIAL_VALUE
93
128
@goto done
94
129
end
95
130
end
111
146
b = peekbyte (source, pos)
112
147
if b == UInt8 (' f' ) || b == UInt8 (' F' )
113
148
x = ifelse (neg, T (- Inf ), T (Inf ))
114
- code |= OK
149
+ code |= OK | SPECIAL_VALUE
115
150
pos += 1
116
151
incr! (source)
117
152
if eof (source, pos, len)
@@ -183,10 +218,10 @@ end
183
218
end
184
219
185
220
# if we need to _widen the type due to `digits` overflow, we want a non-inlined version so base case compilation doesn't get out of control
186
- @noinline _parsedigits (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos) where {T <: SupportedFloats , IntType} =
221
+ @noinline _parsedigits (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos) where {T, IntType} =
187
222
parsedigits (T, source, pos, len, b, code, options, digits, neg, startpos)
188
223
189
- @inline function parsedigits (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos) where {T <: SupportedFloats , IntType}
224
+ @inline function parsedigits (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos) where {T, IntType}
190
225
x = zero (T)
191
226
ndigits = 0
192
227
has_groupmark = options. groupmark != = nothing
@@ -267,12 +302,13 @@ end
267
302
# same as above; if digits overflows, we want a non-inlined version to call with a wider type
268
303
# note that we never expect `frac` to overflow, since it's just keep track of the # of digits
269
304
# we parse post-decimal point
270
- @noinline _parsefrac (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos, frac) where {T <: SupportedFloats , IntType} =
305
+ @noinline _parsefrac (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos, frac) where {T, IntType} =
271
306
parsefrac (T, source, pos, len, b, code, options, digits, neg, startpos, frac)
272
307
273
- @inline function parsefrac (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos, frac) where {T <: SupportedFloats , IntType}
308
+ @inline function parsefrac (:: Type{T} , source, pos, len, b, code, options, digits:: IntType , neg:: Bool , startpos, frac) where {T, IntType}
274
309
x = zero (T)
275
310
parsedanyfrac = false
311
+ FT = FLOAT64
276
312
# check if `b` is a digit
277
313
if b - UInt8 (' 0' ) < 0x0a
278
314
b -= UInt8 (' 0' )
285
321
frac += UInt64 (1 )
286
322
if eof (source, pos, len)
287
323
# input is simple non-scientific-notation floating number, like "1.1"
288
- x = scale (T, digits, - signed (frac), neg)
324
+ x = scale (T, FT, digits, - signed (frac), neg)
289
325
code |= OK | EOF
290
326
@goto done
291
327
end
299
335
end
300
336
# check for exponent notation
301
337
if b == UInt8 (' e' ) || b == UInt8 (' E' ) || b == UInt8 (' f' ) || b == UInt8 (' F' )
338
+ if b == UInt8 (' f' ) || b == UInt8 (' F' )
339
+ FT = FLOAT32
340
+ end
302
341
pos += 1
303
342
incr! (source)
304
343
if eof (source, pos, len)
@@ -327,11 +366,11 @@ end
327
366
328
367
# at this point, we've parsed X and Y in "X.YeZ", but not Z in a scientific notation exponent number
329
368
# we start our exponent number at UInt64(0)
330
- return parseexp (T, source, pos, len, b, code, options, digits, neg, startpos, frac, UInt64 (0 ), negexp)
369
+ return parseexp (T, source, pos, len, b, code, options, digits, neg, startpos, frac, UInt64 (0 ), negexp, FT )
331
370
else
332
371
# if no scientific notation, we're done, so scale digits + frac and return
333
372
if parsedanyfrac
334
- x = scale (T, digits, - signed (frac), neg)
373
+ x = scale (T, FT, digits, - signed (frac), neg)
335
374
else
336
375
x = ifelse (neg, - T (digits), T (digits))
337
376
end
@@ -344,10 +383,10 @@ end
344
383
345
384
# same no-inline story, but this time for exponent number; probably even more rare to overflow the exponent number
346
385
# compared to pre/post decimal digits, but we account for it all the same (a lot of float parsers don't account for this)
347
- @noinline _parseexp (:: Type{T} , source, pos, len, b, code, options, digits, neg:: Bool , startpos, frac, exp:: ExpType , negexp) where {T <: SupportedFloats , ExpType} =
348
- parseexp (T, source, pos, len, b, code, options, digits, neg, startpos, frac, exp, negexp)
386
+ @noinline _parseexp (:: Type{T} , source, pos, len, b, code, options, digits, neg:: Bool , startpos, frac, exp:: ExpType , negexp, FT ) where {T, ExpType} =
387
+ parseexp (T, source, pos, len, b, code, options, digits, neg, startpos, frac, exp, negexp, FT )
349
388
350
- @inline function parseexp (:: Type{T} , source, pos, len, b, code, options, digits, neg:: Bool , startpos, frac, exp:: ExpType , negexp) where {T <: SupportedFloats , ExpType}
389
+ @inline function parseexp (:: Type{T} , source, pos, len, b, code, options, digits, neg:: Bool , startpos, frac, exp:: ExpType , negexp, FT ) where {T, ExpType}
351
390
x = zero (T)
352
391
# note that `b` has already had `b - UInt8('0')` applied to it for parseexp
353
392
while true
@@ -356,19 +395,19 @@ end
356
395
incr! (source)
357
396
if eof (source, pos, len)
358
397
# we finished parsing input like "1.1e1"
359
- x = scale (T, digits, ifelse (negexp, - signed (exp), signed (exp)) - signed (frac), neg)
398
+ x = scale (T, FT, digits, ifelse (negexp, - signed (exp), signed (exp)) - signed (frac), neg)
360
399
code |= OK | EOF
361
400
@goto done
362
401
end
363
402
b = peekbyte (source, pos) - UInt8 (' 0' )
364
403
# if we encounter a non-digit, that must mean we're done
365
404
if b > 0x09
366
- x = scale (T, digits, ifelse (negexp, - signed (exp), signed (exp)) - signed (frac), neg)
405
+ x = scale (T, FT, digits, ifelse (negexp, - signed (exp), signed (exp)) - signed (frac), neg)
367
406
code |= OK
368
407
@goto done
369
408
end
370
409
if overflows (ExpType) && exp > overflowval (ExpType)
371
- return _parseexp (T, source, pos, len, b, code, options, digits, neg, startpos, frac, _widen (exp), negexp)
410
+ return _parseexp (T, source, pos, len, b, code, options, digits, neg, startpos, frac, _widen (exp), negexp, FT )
372
411
end
373
412
end
374
413
@label done
@@ -395,7 +434,14 @@ pow10(::Type{Float32}, e) = (@inbounds v = F32_SHORT_POWERS[e+1]; return v)
395
434
pow10 (:: Type{Float64} , e) = (@inbounds v = F64_SHORT_POWERS[e+ 1 ]; return v)
396
435
pow10 (:: Type{BigFloat} , e) = (@inbounds v = F64_SHORT_POWERS[e+ 1 ]; return v)
397
436
398
- function scale (:: Type{T} , v, exp, neg) where {T}
437
+ _unsigned (x:: BigInt ) = x
438
+ _unsigned (x) = unsigned (x)
439
+
440
+ function scale (:: Type{T} , FT:: FloatType , v, exp, neg) where {T}
441
+ return __scale (float_type (T, FT), _unsigned (v), exp, neg)
442
+ end
443
+
444
+ function __scale (:: Type{T} , v, exp, neg) where {T}
399
445
ms = maxsig (T)
400
446
cl = ceillog5 (T)
401
447
if v < ms
@@ -409,7 +455,7 @@ function scale(::Type{T}, v, exp, neg) where {T}
409
455
end
410
456
end
411
457
v == 0 && return zero (T)
412
- if exp > 308
458
+ if exp > 308 && T != BigFloat
413
459
return T (neg ? - Inf : Inf )
414
460
elseif exp < - 326
415
461
# https://github.com/JuliaData/Parsers.jl/issues/83
@@ -485,9 +531,9 @@ function _scale(::Type{T}, v::V, exp, neg) where {T, V <: UInt128}
485
531
if exp == 23
486
532
# special-case concluded from https://github.com/JuliaLang/julia/issues/38509
487
533
x = v * V (1e23 )
488
- elseif exp >= 0
534
+ elseif 0 <= exp < 290
489
535
x = v * exp10 (exp)
490
- elseif exp < - 308 || v > maxsig (T)
536
+ elseif exp < - 308 || exp > 308 || v > maxsig (T)
491
537
# if v is too large, we lose precision by just doing
492
538
# v / exp10(-exp) since that only promotes to Float64
493
539
# so detect and re-route to this branch where we widen v
@@ -500,15 +546,11 @@ function _scale(::Type{T}, v::V, exp, neg) where {T, V <: UInt128}
500
546
end
501
547
502
548
const BIGEXP10 = [1 / exp10 (BigInt (e)) for e = 309 : 327 ]
503
- const BIGFLOAT = BigFloat[]
504
- if VERSION > v " 1.5"
505
- const BIGFLOATEXP10 = [exp10 (BigFloat (i; precision= 64 )) for i = 1 : 308 ]
506
- else
507
- const BIGFLOATEXP10 = [exp10 (BigFloat (i)) for i = 1 : 308 ]
508
- end
549
+ const BIGFLOATS = BigFloat[]
550
+ const BIGFLOATEXP10 = [exp10 (BigFloat (i; precision= 256 )) for i = 1 : 308 ]
509
551
510
552
function _scale (:: Type{T} , v:: V , exp, neg) where {T, V <: BigInt }
511
- x = access_threaded (BigFloat, BIGFLOAT )
553
+ x = access_threaded (BigFloat, BIGFLOATS )
512
554
513
555
ccall ((:mpfr_set_z , :libmpfr ), Int32,
514
556
(Ref{BigFloat}, Ref{BigInt}, Int32),
@@ -531,7 +573,11 @@ function _scale(::Type{T}, v::V, exp, neg) where {T, V <: BigInt}
531
573
x, x, y, MPFR. ROUNDING_MODE[])
532
574
else
533
575
# v * exp10(V(exp))
534
- y = BIGFLOATEXP10[exp]
576
+ if exp <= 308
577
+ y = BIGFLOATEXP10[exp]
578
+ else
579
+ y = exp10 (BigFloat (exp; precision= 256 ))
580
+ end
535
581
ccall ((:mpfr_mul , :libmpfr ), Int32,
536
582
(Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32),
537
583
x, x, y, MPFR. ROUNDING_MODE[])
0 commit comments