Reputation: 1102
I'm computing the digits of pi using the the Chudnovsky algorithm, and this is my code:
qMax = 1000000
L = Vector{BigFloat}(undef, qMax)
L[1] = 13591409
for q in 2:qMax
L[q] = L[q-1] + 545140134
end
X = Vector{BigFloat}(undef, qMax)
X[1] = 1.0
for q in 2:qMax
X[q] = X[q-1] * (-262537412640768000)
end
M = Vector{BigFloat}(undef, qMax)
M[1] = 1.0
for q in 2:qMax
M[q] = M[q-1] * (((12*(q-2) + 2) * (12*(q-2)+6) * (12*(q-2) + 10))/(q-1)^3)
end
result = 0
for q in 1:qMax
global result += (M[q] * L[q]) / X[q]
end
print(426880 * sqrt(10005) * 1/result)
But when I print the result, I get the same number of digits after the 3 for 1, 10, 1000 etc. values of qMax. How to get arbitrary many digits after the 3?
Upvotes: 1
Views: 1514
Reputation: 44828
You can use setprecision
:
julia> BigFloat(π)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198
julia> setprecision(BigFloat, 1024)
1024
julia> BigFloat(π)
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997
julia>
It's also possible to specify the precision in the constructor:
julia> BigFloat(π, precision=2048)
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526357
However, operators don't seem to respect precision set in the constructor:
julia> setprecision(BigFloat, 10) # low precision
10
julia> BigFloat(1) / 3
0.3335 # OK, low precision, as expected
julia> BigFloat(1, precision=100) / 3
0.3335 # Low precision still!
I would expect higher precision to be propagated through the whole expression, but that doesn't seem to be the case...
EDIT: I opened an issue about this behavior.
EDIT 2:
In the loop that updates M
, all calculations are done on q
, which is a regular integer, so when you divide by (q-1)^3
, that will mean regular floating-point division, not involving BigFloat
.
You can use BigFloat
like this:
for q in 2:qMax
# Use `BigFloat`!
q_f = BigFloat(q)
M[q] = M[q-1] * (
(
(12*(q_f-2) + 2) * (12*(q_f-2)+6) * (12*(q_f-2) + 10)
) / (q_f-1)^3
)
end
Code:
setprecision(1024)
# your code...
println(426880 * sqrt(BigFloat(10005)) * 1/result)
println(BigFloat(π))
Output:
~/test $ julia so2.jl
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587042
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997
Upvotes: 3