And Julia keeps amazing me: high precision computation

[Ronis_BR]

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!

[Benny]

@Ronis_BR

setting setprecision accordingly

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.

[devel-chm]

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

[Ronis_BR]

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!

[abulak]

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

[stevengj]

@Ronis_BR

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.

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:

• Special functions , associated Legendre type 3 - #5 by stevengj
• Sinc function - and sinc neural networks - and function approximation for the former - #4 by stevengj

[Ronis_BR]

@stevengj

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.

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.

[Ronis_BR]

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.

[Ronis_BR]

@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

JuliaSpace/SatelliteAnalysis.jl/blob/32ad240d229e2869e0594d80bf127190ecace50f/src/frozen_orbits.jl#L205

[danielwe]

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.

[Ronis_BR]

Is there any advantage to use binomial in this case?

[danielwe]

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?

[Ronis_BR]

Thanks! I will try!

[danielwe]

    #          4 ⋅ (l - t + 1) ⋅ (p - t + 1) ⋅ (l - p - t + 1)
    #  k_t = - ─────────────────────────────────────────────── ⋅ k_{t - 1} ,
    #            t ⋅ (2l - 2t + 2) ⋅ (2l - 2t + 1) ⋅ sin(i)^2

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

[Ronis_BR]

Awesome! Thank you very much @danielwe !

[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
[...]

[Benny]

@stevengj

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.

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).

[abulak]

@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

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])