BigInt precision issue


#1

I am finding the cube root of a big number. When I raise the answer to the power of 3, I ended up with a different number.

julia> y = 2205316413931134031046440767620541984801091216351222789180535786851451917462804449135087209259828503848304180574549372616172217553002988241140344023060716738565104171296716554122734607654513009667720334889869007276287692856645210293194853
2205316413931134031046440767620541984801091216351222789180535786851451917462804449135087209259828503848304180574549372616172217553002988241140344023060716738565104171296716554122734607654513009667720334889869007276287692856645210293194853

julia> x = BigInt(BigInt(y) ^ (BigInt(1)/BigInt(3)))
13016382529449106065839070830454998857466392684017754632233814825405652260970880

julia> x ^ BigInt(3)
2205316413931134031046440767620541984801091216351222789180535786851451917468010747269793145541333050162118216925368707794887629451553641866690547562191367342837416208560274377321512108411649807242343108436126660405507526307305662185472000

julia> x ^ BigInt(3) == y
false

Let’s try this:

julia> z = 13016382529449106065839070830454998857466392684017754632233814825405652260960637
13016382529449106065839070830454998857466392684017754632233814825405652260960637

julia> z ^ BigInt(3) 
2205316413931134031046440767620541984801091216351222789180535786851451917462804449135087209259828503848304180574549372616172217553002988241140344023060716738565104171296716554122734607654513009667720334889869007276287692856645210293194853

julia> z ^ BigInt(3) == y
true

You may wonder how I obtain the value of z. It came from our old friend Python:

>>> from decimal import *
>>> getcontext().prec = 100
>>> Decimal(2205316413931134031046440767620541984801091216351222789180535786851451917462804449135087209259828503848304180574549372616172217553002988241140344023060716738565104171296716554122734607654513009667720334889869007276287692856645210293194853) ** (Decimal(1)/Decimal(3))
Decimal('13016382529449106065839070830454998857466392684017754632233814825405652260960636.99999999999999999976')

#2

Dividing two BigInts produces a BigFloat. Further operations are floating point and hence not exact.

julia> BigInt(1)/BigInt(3)
3.333333333333333333333333333333333333333333333333333333333333333333333333333348e-01

julia> typeof(BigInt(1)/BigInt(3))
BigFloat

I’m not sure what the issue is with the second example. EDIT: ok I see, no issue, you’re just showing the correct integer answer.


#3

My little investigation.

Default precision:

julia> BigInt((z ^ big"3") ^ (BigInt(1)/BigInt(3))) - z
10243

Same as here:

julia> setprecision(256) do 
       BigInt((z ^ big"3") ^ (big"1"/big"3")) - z 
       end
10243

We could do better:

julia> setprecision(260) do 
       BigInt((z ^ big"3") ^ (big"1"/big"3")) - z 
       end
643

more better:

julia> setprecision(267) do 
       BigInt((z ^ big"3") ^ (big"1"/big"3")) - z 
       end
-5

But this is too much:

julia> setprecision(268) do 
       BigInt((z ^ big"3") ^ (big"1"/big"3")) - z 
       end
ERROR: InexactError: BigInt(BigInt, 1.30163825294491060658390708304549988574663926840177546322338148254056522609606395e+79)
Stacktrace:
 [1] BigInt(::BigFloat) at ./mpfr.jl:272
 [2] (::getfield(Main, Symbol("##75#76")))() at ./REPL[95]:2
 [3] setprecision(::getfield(Main, Symbol("##75#76")), ::Type{BigFloat}, ::Int64) at ./mpfr.jl:904
 [4] setprecision(::Function, ::Int64) at ./mpfr.jl:910
 [5] top-level scope at none:0

But you could do:

julia> setprecision(1000) do 
       w = y ^ (big"1"/big"3")
       # w = BigInt(w)   # this throw Exception! 
       end
       w = BigInt(w)
       w == z
true

#4
julia> using Printf

julia> y = 2205316413931134031046440767620541984801091216351222789180535786851451917462804449135087209259828503848304180574549372616172217553002988241140344023060716738565104171296716554122734607654513009667720334889869007276287692856645210293194853

julia> setprecision(900) do
       round(exp10( log10(y)/3.0 ) )
       end
1.3016382529449106065839070830454998857466392684017754632233814825405652260960637e+79

julia> @printf("%90d",ans)
          13016382529449106065839070830454998857466392684017754632233814825405652260960637