# And Julia keeps amazing me: high precision computation

**URL:** <https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740>\
**Category:** Community\
**Tags:** precision, bigfloat\
**Created:** [December 17, 2023, 5:21pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740 "2023-12-17T17:21:15Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 17, 2023, 5:21pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/1 "2023-12-17T17:21:15Z")

</div>

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](https://github.com/JuliaSpace/SatelliteAnalysis.jl/blob/main/src/frozen_orbits.jl)) and to obtain the frozen eccentricity we need to compute this function:

 ![Captura de Tela 2023-12-17 às 14.16.56](https://global.discourse-cdn.com/julialang/original/3X/3/8/38b4c1696c5b7059231b790094ac011a05186109.png)

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:

```julia
        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!

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [December 17, 2023, 5:49pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/2 "2023-12-17T17:49:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![devel-chm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devel-chm/32/3572_2.png) [@devel-chm](https://discourse.julialang.org/u/devel-chm)\
**Post date:** [December 17, 2023, 6:04pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/3 "2023-12-17T18:04:05Z")

</div>

You can use [Stirling’s approximation](https://en.wikipedia.org/wiki/Stirling%27s_approximation) to the factorial  
to calculate the number of decimal places needed.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 17, 2023, 6:12pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/4 "2023-12-17T18:12:46Z")

</div>

> 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](https://en.wikipedia.org/wiki/Stirling%27s_approximation) to the factorial  
> to calculate the number of decimal places needed.

Thanks for the info!

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [December 18, 2023, 1:50pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/5 "2023-12-18T13:50:03Z")

</div>

You can also use [Arblib.jl](https://github.com/kalmarek/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 😉

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 18, 2023, 2:51pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/6 "2023-12-18T14:51:03Z")

</div>

> [@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](https://discourse.julialang.org/t/special-functions-associated-legendre-type-3/2163/5)
- [Sinc function - and sinc neural networks - and function approximation for the former - #4 by stevengj](https://discourse.julialang.org/t/sinc-function-and-sinc-neural-networks-and-function-approximation-for-the-former/46483/4)

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 18, 2023, 3:26pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/7 "2023-12-18T15:26:39Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 18, 2023, 3:29pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/8 "2023-12-18T15:29:49Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 18, 2023, 4:15pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/9 "2023-12-18T16:15:30Z")

</div>

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

```julia
# 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 🙂

> <https://github.com/JuliaSpace/SatelliteAnalysis.jl/blob/32ad240d229e2869e0594d80bf127190ecace50f/src/frozen_orbits.jl#L205>

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [December 18, 2023, 5:24pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/10 "2023-12-18T17:24:41Z")

</div>

Looking at [https://github.com/JuliaSpace/SatelliteAnalysis.jl/blob/32ad240d229e2869e0594d80bf127190ecace50f/src/frozen\_orbits.jl#L224-L228](https://github.com/JuliaSpace/SatelliteAnalysis.jl/blob/32ad240d229e2869e0594d80bf127190ecace50f/src/frozen_orbits.jl#L224-L228)

```julia
    # 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`.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 18, 2023, 5:30pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/11 "2023-12-18T17:30:03Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [December 18, 2023, 5:35pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/12 "2023-12-18T17:35:30Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 18, 2023, 5:37pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/13 "2023-12-18T17:37:11Z")

</div>

Thanks! I will try!

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [December 18, 2023, 5:41pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/14 "2023-12-18T17:41:52Z")

</div>

> [@](#):
>
> ```julia
> # 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,

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

```

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [December 18, 2023, 8:10pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/15 "2023-12-18T20:10:41Z")

</div>

Awesome! Thank you very much @danielwe !

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [December 18, 2023, 8:24pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/16 "2023-12-18T20:24:01Z")

</div>

On the topic of “why binomial”:

```julia-repl
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
[...]

```

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [December 18, 2023, 11:51pm UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/17 "2023-12-18T23:51:19Z")

</div>

> [@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)`.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [December 22, 2023, 10:46am UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/18 "2023-12-22T10:46:10Z")

</div>

3 posts were split to a new topic: [Latex parens sizing](https://discourse.julialang.org/t/latex-parens-sizing/107940)

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [December 19, 2023, 1:12am UTC](https://discourse.julialang.org/t/and-julia-keeps-amazing-me-high-precision-computation/107740/19 "2023-12-19T01:12:06Z")

</div>

@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

```julia
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

```julia
ab = big(a)
ib = big(i)
num = zero(BigFloat)
den = zero(BigFloat)

```

by

```julia
num = zero(ab)
den = zero(ab)

```

- the `for p in 1:p_max` loop should be replaced with

```julia
    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)

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

```