# Timing Clebsch-Gordan-coefficients calls in Julia (beating c++ again)

**URL:** <https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736>\
**Category:** Numerics\
**Created:** [December 27, 2019, 7:54am UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736 "2019-12-27T07:54:27Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [December 27, 2019, 7:54am UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/1 "2019-12-27T07:54:27Z")

</div>

I have compared 4 different implementations of the [Clebsch-Gordan coefficients](https://en.wikipedia.org/wiki/Clebsch%E2%80%93Gordan_coefficients)

1. [freshly coded](https://github.com/mmikhasenko/ThreeBodyDecay.jl/blob/master/src/angular_functions.jl#L71) in my [`ThreeBodyDecay.jl`](https://github.com/mmikhasenko/ThreeBodyDecay.jl) package
2. from `GSL.jl` call (c++ implementation)
3. from `SymPy.jl` (for the sake of completeness)
4. julia package [`WignerSymbols.jl`](https://github.com/Jutho/WignerSymbols.jl)

The timing for 1000 calls of random coefficient with j₁, j₂ \< 16 follows.

```julia
ThreeBodyDecay.jl 187 μs
GSL 743 μs
WignerSymbols.jl 764 μs
SymPy 1.61 s

```

Details are [here](https://github.com/mmikhasenko/ThreeBodyDecay.jl/blob/master/test/clebsches_native.jl#L79).

- I was glad to use the GSL until I discovered that the evaluation of the CG coefficients in some scattering amplitudes takes 50% of the total time. A native Julia implementation (perhaps, case-specific) turned to be 4 times (!) faster.
- It seems that there is a room for optimization in the `WignerSymbols.jl`. I am happy to help if authors can point out what the right strategy could be. Having a dedicated, tested package is obviously preferable to the customary implementation.
- I have a similar implementation of the Wigner-d functions. AFAIK, they do not exist in julia yet (?).

I am wondering if creating a new package (e.g. HEP-specific) is a good idea.

---

<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 27, 2019, 2:42pm UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/2 "2019-12-27T14:42:38Z")

</div>

It looks like most of your 1000 calls are trivial cases where `ClGd` returns `0.0`; maybe it would be better to benchmark only the nontrivial cases?

It looks like you can get some additional speedup by eliminating your `f_logfact(n)` and just doing `logfact[n+1]` directly with an `@inbounds` directive (assuming that you’ve previously checked the arguments to ensure in-bounds accesses). (Even without the `@inbounds`, the `n < 0` check in `f_logfact` is redundant with bounds-checking.) (To get some slight additional speedup, you could replace `logfact[1+div(x,2)]` with `logfact2[1+x]`, where `logfact2` is a double-length array in which each entry is repeated twice.) Also, isn’t `res = 0` type-unstable?

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [December 27, 2019, 7:47pm UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/3 "2019-12-27T19:47:07Z")

</div>

Hi Steven, these are great suggestions!  
The amount of zeros is ≈ 17%, and 14% are equal to 1.0. I tried to exclude obviosly-trivial Clebsches in the `rand_clebsch()` function (+fixed a typo).

Just checked the other suggestions:

- type stability + `logfact2` gives 1-2 μs gain
- `@inbounds` removes extra 2-4 μs for the same tests.

Not so much, but still cool. Thanks a lot for pointing out!

Did you notice if I am doing something wrong with `SymPy` calls? My colleagues running tensorflow used to calculate both the Clebsch-Gordan coefficientss and the Wigner d using the `sympy.physics.wigner`. It is hard to believe that the call time which they are working with is x1000. I have not benchmarked it thoroughly, tho.

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [May 9, 2020, 2:57pm UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/4 "2020-05-09T14:57:58Z")

</div>

update: the package `PartialWaveFunctions.jl` has been registered!

> **[GitHub - mmikhasenko/PartialWaveFunctions.jl: Partial Wave Analysis specific...](https://github.com/mmikhasenko/PartialWaveFunctions.jl)**
>
> Partial Wave Analysis specific special funtions. Contribute to mmikhasenko/PartialWaveFunctions.jl development by creating an account on GitHub.

timing is here

> <https://github.com/mmikhasenko/PartialWaveFunctions.jl/blob/master/test/timing_different_packages.jl#L108>

---

<div class="post-metadata">

**Author:** ![xzackli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xzackli/32/38301_2.png) [@xzackli](https://discourse.julialang.org/u/xzackli)\
**Post date:** [May 9, 2020, 7:40pm UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/5 "2020-05-09T19:40:17Z")

</div>

WignerSymbols.jl could use a bit of performance tuning, but I would mention that the prime factorization method in that package is able to compute symbols to very high quantum numbers _exactly_, in the form of roots of `Rational{BigInt}`.

This is somewhat orthogonal to the methods you’ve benchmarked, but I thought it might be helpful to mention – I’ve recently implemented (yet another) method in [WignerFamilies.jl](https://github.com/xzackli/WignerFamilies.jl). In my field (cosmology) the Wigner 3j-symbols are not usually used via random access, but instead in summations over all nontrivial symbols for some fixed quantum numbers. The Wigner 3j-symbols are subject to convenient recurrence relations which can provide all symbols in a family with relatively few floating-point operations per symbol.

For example, there are 2001 nontrivial 3j symbols varying j with fixed j\_2 = 1000, j\_3 = 1000, m\_2=2, m\_3=-2, m\_1 = -m\_2 -m\_3.

f(j) = \begin{pmatrix} j & j\_2 & j\_3 \\ 0 & 2 & -2 \end{pmatrix}

```julia
using WignerFamilies, BenchmarkTools
@btime wigner3j_f(1000,1000,2,-2)
  22.369 μs (2 allocations: 15.84 KiB)
2001-element WignerSymbolVector{Float64,Int64,Array{Float64,1}} with indices WignerFamilies.URange(0,2000):
  0.022355091700494854
  4.468784506164827e-5
 -0.011177416041038234
  ⋮
 -0.001741253798177275
  0.0003163901836069519
  0.0025009714219007464

```

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [May 9, 2020, 8:02pm UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/6 "2020-05-09T20:02:43Z")

</div>

> [@xzackli](#):
>
> I’ve recently implemented (yet another) method in [WignerFamilies.jl](https://github.com/xzackli/WignerFamilies.jl)

Hi, thanks for the comment. Good to know about `WignerFamilies.jl`, I should add a link to the package in my README. Indeed, my clebsches would not work above j=50 without extending precalculated factorials (well, it is never needed in HEP).

---

<div class="post-metadata">

**Author:** ![oschulz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oschulz/32/2998_2.png) [@oschulz](https://discourse.julialang.org/u/oschulz)\
**Post date:** [May 12, 2020, 7:50am UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/7 "2020-05-12T07:50:29Z")

</div>

Hi @misha_mikhasenko, your package will be very welcome here, if you’re interested: [https://github.com/JuliaHEP](https://github.com/JuliaHEP)

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [May 12, 2020, 8:45am UTC](https://discourse.julialang.org/t/timing-clebsch-gordan-coefficients-calls-in-julia-beating-c-again/32736/8 "2020-05-12T08:45:34Z")

</div>

Hi Oliver,  
yes, makes sense, It should be moved there. Thanks