Comments (10)
Further investigation shows, that the powers of big integers are unreliable for ~ n < -2^53
:
The first column uses Int
, the second Float64
exponent type. Third column the correct value.
julia> [45:62 prevfloat(1.0).^(-2 .^ (45:62)) prevfloat(1.0).^float(-2 .^ (45:62)) exp.(log(big(prevfloat(1.0))) .* (-2 .^ (45:62)))]
18×4 Matrix{BigFloat}:
45.0 1.00391 1.00391 1.00391
46.0 1.00781 1.00781 1.00784
47.0 1.01562 1.01562 1.01575
48.0 1.03123 1.03123 1.03174
49.0 1.06233 1.06233 1.06449
50.0 1.12352 1.12352 1.13315
51.0 1.23654 1.23654 1.28403
52.0 1.35914 1.35914 1.64872
53.0 1.10235e-07 1.10235e-07 2.71828
54.0 -54.5981 -54.5981 7.38906
55.0 -8942.87 -8942.87 54.5982
56.0 -6.22028e+07 -6.22028e+07 2980.96
57.0 -1.18444e+15 -1.18444e+15 8.88611e+06
58.0 -1.9329e+29 -1.9329e+29 7.8963e+13
59.0 -2.44925e+57 -2.44925e+57 6.23515e+27
60.0 -1.91951e+113 -1.91951e+113 3.88771e+55
61.0 -5.82523e+224 -5.82523e+224 1.51143e+111
62.0 Inf Inf 2.28441e+222
from julia.
The issue here is that inv(x)
returns a rounded Float64 of the inverse — with 52 bits precision. In fact, we start losing 1ulp accuracy at -2^27 and it just goes downhill from there. Seems we need a better strategy here.
Line 1322 in e07c0f1
from julia.
@mbauman that gets compensated for on
Line 1324 in e07c0f1
I think something must be going wrong in the error compensation though. Investigating now.
from julia.
Ah, I see. We still are missing compensation somewhere here:
julia> ulps(b, x) = abs(b^x - Float64(big(b)^x))/eps(b^x)
ulps (generic function with 3 methods)
julia> [ulps(prevfloat(1.0), -2^e) for e in 26:54]
29-element Vector{Float64}:
0.0
1.0
2.0
8.0
32.0
128.0
512.0
2048.0
8192.0
32768.0
131073.0
524302.0
2.097258e6
8.389461e6
3.3561258e7
1.34272348e8
5.37307984e8
2.15098174e9
8.617942835e9
3.4584177964e10
1.3924045658e11
5.6426422088e11
2.316653480914e12
9.762831672521e12
4.3352610026219e13
2.13851005933614e14
1.304153939836051e15
2.0538755963422595e23
8.723923242651639e15
from julia.
When I fixed the first issue (typemin(Int)
), I found the second (inaccuracy) disappeared. Will drop PR soon.
from julia.
@KlausC that suggests something you were doing is wrong. From what I can tell, this appears to be an issue caused by me skipping a term in the error compensation that is very small but can accumulate to a significant amount. At this point, I believe that I have corrected for the case of power of 2 exponents, but not for the case of exponents with bits set.
from julia.
that suggests something you were doing is wrong. From what I can tell, this appears to be an issue caused by
That is possible. I fixed the first issue primarily. As a side-effect the exploding errors for prevfloat(1.0)
powers disappeared. Of course it is possible, that there are more improvements needed.
from julia.
BTW, I repeated the ULPs counting after my fix:
julia> ulps(b, x) = abs(b^x - Float64(big(b)^x))/eps(b^x)
ulps (generic function with 1 method)
julia> [26:63 [ulps(prevfloat(1.0), -2^e) for e in 26:63]]
38×2 Matrix{Float64}:
26.0 0.0
27.0 0.0
28.0 0.0
29.0 0.0
30.0 0.0
31.0 0.0
32.0 0.0
33.0 0.0
34.0 0.0
35.0 0.0
36.0 0.0
37.0 0.0
38.0 0.0
39.0 0.0
40.0 0.0
41.0 0.0
42.0 0.0
43.0 0.0
44.0 0.0
45.0 0.0
46.0 0.0
47.0 0.0
48.0 0.0
49.0 0.0
50.0 0.0
51.0 0.0
52.0 0.0
53.0 1.0
54.0 1.0
55.0 3.0
56.0 12.0
57.0 34.0
58.0 144.0
59.0 646.0
60.0 3255.0
61.0 10322.0
62.0 51893.0
63.0 NaN
from julia.
what about [26:63 [ulps(prevfloat(1.0), 1-2^e) for e in 26:63]]
?
from julia.
Also fine:
julia> [26:63 [ulps(prevfloat(1.0), -2^e+1) for e in 26:63]]
38×2 Matrix{Float64}:
26.0 0.0
27.0 1.0
28.0 0.0
29.0 0.0
30.0 0.0
31.0 0.0
32.0 0.0
33.0 0.0
34.0 0.0
35.0 0.0
36.0 0.0
37.0 1.0
38.0 0.0
39.0 1.0
40.0 0.0
41.0 0.0
42.0 0.0
43.0 0.0
44.0 0.0
45.0 0.0
46.0 1.0
47.0 0.0
48.0 1.0
49.0 1.0
50.0 1.0
51.0 0.0
52.0 1.0
53.0 0.0
54.0 0.0
55.0 2.0
56.0 11.0
57.0 33.0
58.0 144.0
59.0 646.0
60.0 3254.0
61.0 10321.0
62.0 51892.0
63.0 NaN
from julia.
Related Issues (20)
- Autocompletion triggers error in invalid `AbstractDict`
- Foo{T >: Int} gives ERROR: UndefVarError: `T` not defined HOT 7
- Double counting test broken on i686-linux
- JET linting error from `wrap` on `Array` (on 1.11) HOT 3
- Add `strict` mechanism for opting into stricter subsets of the language HOT 10
- `UnionAll`s with trivial bodies ignore type lower bounds? HOT 2
- RFC: Export versioning
- auto-completion doesn't work correctly in pkg mode HOT 2
- verify_ir desn't catch SSAValue from outside range HOT 1
- Allow type variables with `do` syntax HOT 3
- Type instabilities with `mapslices` and `PermutedDimsArray` functions HOT 1
- maximum(abs, ::Array{Complex{Int}}; dims=...) fails with an internal error HOT 4
- Concatenation of string constants isn't concretely evaluated, but is inferred as a `Const` HOT 4
- Test failure in `Sockets` with `IOError: listen: address already in use (EADDRINUSE)`
- Overrides.toml broken on master
- `GlobalRef`s as `Dict` keys in precompilation are invalid HOT 1
- `Array{T}(::T)` should return a zero-dimensional array HOT 7
- GC is too lazy, sometimes HOT 10
- External Pkg always recompiles upon loading HOT 12
- REPLExt for Pkg no longer loads when loading a custom Pkg HOT 2
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from julia.