# And Julia keeps amazing me: high precision computation

Hi!

This application is probably very common for some people, but it is the very first time I really needed high precision computation (`BigInt` and `BigFloat`).

I am implementing an algorithm for designing frozen orbits in SatelliteAnalysis.jl (https://github.com/JuliaSpace/SatelliteAnalysis.jl/blob/main/src/frozen_orbits.jl) and to obtain the frozen eccentricity we need to compute this function:

The problem is when you try using very high degrees, leading to 300! and so on. I was thinking how I could possible handle this problem, like interactively computing the terms using a nice βguessβ to always divide things to reach lower numbers.

However, it turns out that just using `BigInt` and `BigFloat`, setting `setprecision` accordingly, everything just works out-of-the-box without requiring too much thinking, Hence, I could write:

``````        k_t = factorial(2lb - 2tb) / (
factorial(tb) *
factorial(lb - tb) *
factorial(pb - tb) *
factorial(lb - pb - tb) *
big(2)^(2lb - 2tb)
)
``````

which is exactly the same equation as I needed.

Conclusion: I implemented everything is very small time, which was awesome!

After all these years, I still find amazing things that Julia does extremely well. Thanks, all the devs!

31 Likes

Got any tips for estimating how much precision is needed, given an expression with factorials and inputs? That is, assuming you do the straightforward factorials instead of cancelling terms and reordering operations to restrain the magnitude.

1 Like

You can use Stirlingβs approximation to the factorial
to calculate the number of decimal places needed.

3 Likes

Got any tips for estimating how much precision is needed, given an expression with factorials and inputs?

In my case, I just use a huge number (1024 I think). I know it is correct because the difference between the values computed for an increasing degree was always decreases.

That is, assuming you do the straightforward factorials instead of cancelling terms and reordering operations to restrain the magnitude.

After the implementation, I tried to be more βcleverβ. However, due to the lookup table that is probably used in factorials, the algorithm was way slower than just naively computing the factorials.

You can use Stirlingβs approximation to the factorial
to calculate the number of decimal places needed.

Thanks for the info!

2 Likes

You can also use Arblib.jl to see how much precision you really need (youβll get a certified error bound);

This will require a little bit of coding, as youβll need to implement `factorial` using `Arblib.fac!`, `binomial` using `Arblib.bin!` (PRs most wellcome ;). Let me know if you need help with that

6 Likes

In series expressions like this, you should almost never call the `factorial` function, nor should you compute the powers for each term separately. Instead, compute each term iteratively from the previous term, starting with the smallest power and the smallest factorial. Not only is this computationally cheaper, but it also allows you to circumvent spurious over/underflow and avoid having to use expensive arbitrary-precision arithmetic.

This is something that many people forget. See, for example:

9 Likes

I fully agree! If I obtain those terms and keep everything I have calculated in an accumulator, I will be able to reduce the computational cost significantly. However, and this is what amazed me, the current version is taking some ms to run for 360 degrees. This algorithm is computing the frozen eccentricity and it will be run once or twice per analysis. Hence, performance is not a problem. Furthermore, the precision down to `Float64` (I am converting the result at the end) is equal to what we have in our references.

4 Likes

Ah, by the way, I tried to compute those factorials in each term of the summation using a βcleverβ approach, avoiding calculating the same thing multiple times. However, considering this approach that I do not save the factors for the next summation term, the code was almost 100x slower than calling `factorial`. I have no idea why but I guess `factorial` is using some lookup table.

1 Like

@stevengj On second thoughts, it is better to optimize the code as other people might have different use cases. I followed your advice and implemented those factors interactively. It indeed decreased the computational workload by a lot:

``````# Before
julia> @btime frozen_orbit(7130.982e3, 98.410 |> deg2rad; gravity_model=egm96)
1.728 ms (84088 allocations: 2.25 MiB)
(0.0011641853028456078, 1.5707963267948966)

# After
julia> @btime frozen_orbit(7130.982e3, 98.410 |> deg2rad; gravity_model = egm96)
1.103 ms (64248 allocations: 1.57 MiB)
(0.0011641853028456078, 1.5707963267948966)
``````

I mentioned your name in the comments to thanks for the suggestion. I hope you do not mind

10 Likes
``````    # where the initialization is:
#
#                    2l!
#  k_0 =  βββββββββββββββββββββββββββ β sin(i)^(l) β (-1)^{p - k} .
#         l! β p! β (l - p)! β 2^(2l)
``````

note that the factorials here are equivalent to `binomial(2l, l) * binomial(l, p)`, so you donβt need `factorial` to get `k_0` either. Also, you might try to combine `1 / 2^(2l)` and `sin(i)^l` into `(sin(i) / 4)^l`.

4 Likes

Is there any advantage to use `binomial` in this case?

I havenβt checked, but Iβd be surprised if it isnβt computationally much more efficient, and it doesnβt overflow as long as the result is representable, so maybe you can completely avoid using BigInt?

1 Like

Thanks! I will try!

Looks like the first parenthesized factors in the numerator and denominator cancel and return 1/2. Hence,

``````    #          2 β (p - t + 1) β (l - p - t + 1)
#  k_t = - βββββββββββββββββββββββββββββββββ β k_{t - 1} ,
#            t β (2l - 2t + 1) β sin(i)^2
``````
1 Like

Awesome! Thank you very much @danielwe !

On the topic of βwhy binomialβ:

``````julia> using BenchmarkTools

julia> naivebinomial(n, k) = factorial(n) Γ· (factorial(n - k) * factorial(k));

julia> @btime binomial(20, 7);
2.672 ns (0 allocations: 0 bytes)

julia> @btime binomial(21, 7);
2.673 ns (0 allocations: 0 bytes)

julia> @btime naivebinomial(20, 7);
19.002 ns (0 allocations: 0 bytes)

julia> @btime naivebinomial(21, 7);
ERROR: OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead
[...]
``````
5 Likes

Is it at all feasible for a different `factorial` function to construct a lazy product that combines with others to cancel and combine terms, and materializing it would compute the terms while minimizing the magnitude like this? I imagine bounds with actual values would be needed to combine terms properly, like switching iteration order `1*2*3` vs `(4-3)*(4-2)*(4-1)`.

1 Like

3 posts were split to a new topic: Latex parens sizing

@Ronis_BR I decided to jump the wagon of optimizing the code with the aim of having the same code version working for both `Float64`, `BigFloat`s and `Arb`s

What is needed for Arbs to play βlike floatsβ here

``````function Base.Int(x::Arb)
Arblib.isexact(x) || throw(InexactError(:Int, Int, x))
return Int(midpoint(x))
end

Base.div(x::Arb, n) = Arb(div(Int(x), n); prec=precision(x))

function Base.sincos(x::Arb)
s, c = zero(x), zero(x)
Arblib.sin_cos!(s,c,x)
return (s,c)
end

function Base.binomial(n::Arb, k::Arb)
z = zero(n)
Arblib.bin!(z, n, unsigned(Int(k)))
return z
end
``````

then `frozen_orbit` needs only basic changes:

• change the signature to `frozen_orbit(ab::Number, ib::Number; ...)`
• replace
``````ab  = big(a)
ib  = big(i)
num = zero(BigFloat)
den = zero(BigFloat)
``````

by

``````num = zero(ab)
den = zero(ab)
``````
• the `for p in 1:p_max` loop should be replaced with
``````    for _p in 1:p_max
p = oftype(ab, _p)

``````
• drop the Float64 conversion from `e = Float64(2 * num / den)` and return just `e`

Finally replace `_F_and_βF_l0p` by (sorry for getting rid of your beautiful comments)

``````function _F_and_βF_l0p(lb::Number, pb::Number, ib::Number)
F  = zero(ib)
βF = zero(ib)

kb = div(lb, 2)

sin_i, cos_i = sincos(ib)
sinΒ²_i = sin_i * sin_i

fact = isinteger((pb - kb)/ 2) == 0 ? one(lb) : -one(lb)
k_t  = binomial(2lb, lb) * binomial(lb, pb) * sin_i^lb * fact / 2^(2lb)

F  = k_t
βF = k_t * lb / sin_i

for tb in 1:min(Int(kb), Int(pb))
k_t *= -2 * (pb - tb + 1) * (lb - pb - tb + 1) / (tb * (2lb - 2tb + 1) * sinΒ²_i)
F   += k_t
βF  += k_t * (lb - 2tb) / sin_i
end

return F, βF * cos_i
end
``````

(btw, youβre still using `l` in the computation of `k_t`, looks like an unintentional oversight).

I didnβt run the `frozen_orbit`, but I can confirm that the latter function agrees with the original and runs for all types:

``````julia> let (l,p,i) = (16, 8, 55.0)
for f in (identity, big, Arb)
@info _F_and_βF_l0p(f(l), f(p), f(i))
end
end
[ Info: (-0.03602523798271656, 0.22704216792630078)
[ Info: (-0.0360252379832088852388914537148436296233180397111850594594684606756466552166546, 0.2270421679267800856096976709245996971231420231252597114888955736847305713261805)
[ Info: ([-0.0360252379832088852388914537148436296233180397111850594594684606756467 +/- 5.69e-71], [0.22704216792678008560969767092459969712314202312525971148889557368473057 +/- 4.46e-72])
``````
2 Likes