How do you do arbitrary precision calculations in Julia?

I come from a Python background and an important point using Python’s Decimal module for arbitrary precision calculations is knowing when you do and don’t have to specify numbers of type Decimal in, for example, probability calculations.

How often do you have to specify arbitrary precision types for numbers in a Julia calculation?

I learn best by examples and the example I want to solve follows.

Calculate the probability of at least one success when randomly sampling without replacement a bucket of 26 ABC to Z blocks 26!, 3(26)!, and 5(26)! times.

Given:

Success on one trial = randomly drawing the blocks out in the correct order: A, B, C, to Z.

The number of random trials is n = 26! or 3n or 5n.

The probability of success on one random trial, p = 1/n

The probability of failure on one random trial is f = 1 – p.

Calculate (using type specification as little as necessary):

The confidence level, CL, of at least one success in n trials using CL = 1 – f^n

The CL for 3n trials using CL = 1 - f^(3n)

The CL for 5n trials using CL = 1 - f^(5n)

2 Answers
2

Julia works very well at flowing the "big" number type specification through the calculations.
For example, in the following code I only specify one number as big"26", everything else is automatic.

I have just started using Julia, but have been doing arbitrary precision calculations in various ways for years. Julia provides, hands-down, the most pleasant experience with this kind of thing I have ever had.

# Set precision to 150 bits which should be adequate precision.
setprecision(150)

# Note that we only have to specify a "big" number here.
n = factorial(big"26")
println("n = factorial(big"26") = ", n)
println("Note that we never have to use "big" again in the following code.")
println("typeof(n) = ", typeof(n), "n")

# p is the probability of success on 1 trial.
p = 1/n
println("p = 1/n = ", p)
# Note we did not have to specify the type of p.
println("typeof(p) = ", typeof(p), "n")

# f is the probability of failure on 1 trial.
f = 1 - p
println("f = 1 - p = ", f)
println("typeof(f) = ", typeof(f), "n")

# CL is the probability of at least 1 success in n trials.
# CL stands for confidence level.
CL = 1 - f^n
println("The 63% CL for n trials = 1 - f^n = ", CL)
println("typeof(CL) = ", typeof(CL), "n")

# Here is the 95% conf. level using 3n random trials.
CL95 = 1 - f^(3n)
println("The 95% CL for 3n trials = ", CL95)
println("typeof(CL95) = ", typeof(CL95), "n")

# Here is the 99% conf. level using 5n random trials.
CL99 = 1 - f^(5n)
println("The 99% CL for 5n trials = ", CL99)
println("typeof(CL99) = ", typeof(CL99), "n")

""" ============================= Output ==============================
n = factorial(big"26") = 403291461126605635584000000
Note that we never have to use "big" again in the following code.
typeof(n) = BigInt

p = 1/n = 2.4795962632247974600749435458479566174226555415e-27
typeof(p) = BigFloat

f = 1 - p = 9.9999999999999999999999999752040373677520254001e-01
typeof(f) = BigFloat

The 63% CL for n trials = 1 - f^n = 6.3212055882855767839219205578358958187929158048e-01
typeof(CL) = BigFloat

The 95% CL for 3n trials = 9.5021293163213605701567013782477488392169554992e-01
typeof(CL95) = BigFloat

The 99% CL for 5n trials = 9.9326205300091453290223898909666750856240017783e-01
typeof(CL99) = BigFloat
"""


If you want a Precision Decimal Floating Point like the Decimal in Python, I have written one myself. You can set the decimal precision to any positive integer value. Contact me if you want a copy of the module.

julia> factorial(BigInt(26))
403291461126605635584000000

julia> println("How many decimal digits do we need?")
How many decimal digits do we need?

julia> length(string( factorial(BigInt(26)) ))
27

julia> using PDFPs

julia> PDFP_setDefaultPrecision(40)
40

julia> n = PDFP( factorial(BigInt(26)) )
PDFP(0, 26, [4, 0, 3, 2, 9, 1, 4, 6, 1, 1, 2, 6, 6, 0, 5, 6, 3, 5, 5, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

julia> disp(a::PDFP) = println(PDFP_toShortCommonString(a))
disp (generic function with 1 method)

julia> p = 1/n
PDFP(0, -27, [2, 4, 7, 9, 5, 9, 6, 2, 6, 3, 2, 2, 4, 7, 9, 7, 4, 6, 0, 0, 7, 4, 9, 4, 3, 5, 4, 5, 8, 4, 7, 9, 5, 6, 6, 1, 7, 4, 2, 2])

julia> disp(p)
2.47960E-27

julia> println("Let's call failure q as per convention in probability questions")
Let's call failure q as per normal in probability questions

julia> q = 1 - p
PDFP(0, -1, [9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 7, 5, 2, 0, 4, 0, 3, 7, 3, 6, 7, 7, 5, 2])

julia> disp(q)
1.00000

julia> PDFP_toFortranString(q)
"9.999999999999999999999999975204037367752E-1"

julia> at_least_one_success_in_n_trials = 1 - q^n
PDFP(0, -1, [6, 3, 2, 1, 2, 0, 5, 5, 8, 8, 2, 8, 5, 5, 8, 0, 5, 5, 6, 0, 2, 0, 5, 7, 9, 5, 3, 4, 5, 6, 9, 0, 2, 1, 9, 8, 7, 6, 4, 9])

julia> disp(at_least_one_success_in_n_trials)
0.632121

julia> at_least_one_success_in_3n_trials = 1 - q^(3*n)
PDFP(0, -1, [9, 5, 0, 2, 1, 2, 9, 3, 1, 6, 3, 2, 1, 3, 6, 2, 1, 0, 1, 6, 5, 3, 1, 1, 5, 3, 8, 4, 6, 6, 4, 3, 9, 0, 4, 8, 2, 9, 1, 0])

julia> disp(at_least_one_success_in_3n_trials)
0.950213

julia> at_least_one_success_in_5n_trials = 1 - q^(5*n)
PDFP(0, -1, [9, 9, 3, 2, 6, 2, 0, 5, 3, 0, 0, 0, 9, 1, 4, 5, 6, 7, 4, 4, 6, 4, 3, 7, 4, 3, 4, 3, 4, 4, 7, 4, 8, 6, 8, 9, 6, 2, 8, 1])

julia> disp(at_least_one_success_in_5n_trials)
0.993262


Thanks for contributing an answer to Stack Overflow!

But avoid …

To learn more, see our tips on writing great answers.

Required, but never shown

Required, but never shown

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy