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!

33 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 :wink:

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 :slight_smile:

10 Likes

Looking at https://github.com/JuliaSpace/SatelliteAnalysis.jl/blob/32ad240d229e2869e0594d80bf127190ecace50f/src/frozen_orbits.jl#L224-L228

    # 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, BigFloats and Arbs :slight_smile:

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