including Jacobi polynomials and Wigner functions

[mmikhasenko]

I have the functions coded in julia in the PartialWaveFunctions package.

jacobi_pols
https://github.com/mmikhasenko/PartialWaveFunctions.jl/blob/master/src/wignerd.jl#L10

wignerd, wignerD
https://github.com/mmikhasenko/PartialWaveFunctions.jl/blob/master/src/wignerd.jl#L60

Does it fit the SpecialFunctions.jl package? Do you accept PR?

What is about clebsch_gordan coefficients?

[stevengj]

Why aren't you calculating it using a Clenshaw algorithm using the 3-term recurrence? That's the typical way to evaluate orthogonal polynomials, and seems a lot more efficient than evaluating lots of exp and log functions.

[stevengj]

See also #175, #124, #24, #136.

I feel like evaluating orthogonal polynomials is kind of a different category of problem than evaluating transcendental functions, and perhaps there should be a separate package (SpecialPolynomials.jl?). I think that there should be a package for these, in part because a lot of people want them, and also because they are easy to implement badly, but combining them with SpecialFunctions.jl may complicate maintenance more than it is worth…

I think @dlfivefifty may have been working on something along these lines in https://github.com/JuliaApproximation/OrthogonalPolynomialsQuasi.jl and similar?

[dlfivefifty]

Yes, and it works pretty well:

julia> a,b = 0.1,0.2; x = 0.3;

julia> @time Jacobi(a,b)[x,1:100_000_000] # Evaluate P_n^(a,b)(x) for n=0:99_999_999
  1.066723 seconds (36 allocations: 762.942 MiB, 15.73% gc time)
100000000-element Array{Float64,1}:
  1.0
  0.295
 -0.40136250000000007
 -0.4040010625000001
  0.09190393023437506
  0.3755560868932033
  0.12997481357103505
 -0.25062955253876074
  ⋮
 -6.12822549513539e-5
 -8.111091463677328e-5
  1.2615705799800099e-5
  8.868033726769704e-5
  4.059249642093418e-5
 -6.432483865011064e-5
 -7.918739901210108e-5

There is some discussion of pulling the evaluation code into a separate, more lightweight package that doesn't use the ContinuumArrays.jl language (CC @guilgautier):

JuliaApproximation/OrthogonalPolynomialsQuasi.jl#54 (comment)

This would still use InfiniteArrays.jl and LazyArrays.jl as it allows for fast recurrence evaluations without having to rewrite the forward recurrence for each polynomial.

[mmikhasenko]

Why aren't you calculating it using a Clenshaw algorithm using the 3-term recurrence? That's the typical way to evaluate orthogonal polynomials, and seems a lot more efficient than evaluating lots of exp and log functions.

That is interesting, I was not aware of the Clenshaw algorithm. When it is beneficial?
I would be curious to benchmark.

[mmikhasenko]

@dlfivefifty do you also have wigner d-functions codded with quasi arrays?

[mmikhasenko]

I was also wondering what is the best way to allow for the call of the polynomial functions with symbolic input.
As far as I understand SpecialFunctions.jl just work.
@dlfivefifty, would your implementation work with SymPy.Sym(x) input?

[dlfivefifty]

@dlfivefifty do you also have wigner d-functions codded with quasi arrays?

No but @MikaelSlevinsky might have something

That is interesting, I was not aware of the Clenshaw algorithm. When it is beneficial?

Anytime you need to generate all OPs in a sequence: even cos(n*acos(x)) is much slower than the recurrence T_{n+1} = 2x*T_n - T_{n-1} which involves no transcendentals.

@dlfivefifty, would your implementation work with SymPy.Sym(x) input?

Probably

[stevengj]

Even for evaluating a single Jacobi polynomial of degree n, both the Clenshaw algorithm and @mmikhasenko's implementation are O(n), but the Clenshaw algorithm should be faster because it involves no transcendental functions.

The Clenshaw algorithm also works better for a broader range of input types, e.g. symbolic inputs, because it only involves * and + operations.

[mmikhasenko]

Anytime you need to generate all OPs in a sequence: even cos(n*acos(x)) is much slower than the recurrence T_{n+1} = 2x*T_n - T_{n-1} which involves no transcendentals.

I do not use acos, but I use log(1+cos), and exp. Getting rid of these might indeed lead to speed up.
Then, my plan would be to check if the current implementation of Clenshaw is faster than mine and switch to it.

[MikaelSlevinsky]

As for Gaunt (Clebsch--Gordan) coefficients, you may find this implementation interesting https://github.com/JuliaApproximation/FastTransforms.jl/blob/master/src/gaunt.jl. There may be others in other packages.

[mmikhasenko]

Thanks! I have the coefficients coded in Julia, it works fine.
My post is more about convergence in the ecosystem to a finite set of well documented packages.

[mmikhasenko]

The package description suggests that it is c++ wrapper. That is not the right place for my code.

[MikaelSlevinsky]

In my opinion, JuliaApproximation is the right place: ApproxFun.jl is the umbrella package (that @dlfivefifty is refactoring via QuasiArrays.jl-ContinuumArrays.jl-OrthogonalPolynomialsQuasi.jl), but there's also FastGaussQuadrature.jl and FastTransforms.jl that are standalone and reasonably well used/documented.

[mmikhasenko]

@dlfivefifty , do you think the Wigner d fit to OrthogonalPolynomiasQuasi.jl ?

From here https://en.m.wikipedia.org/wiki/Wigner_D-matrix

[dlfivefifty]

https://github.com/JuliaApproximation/SphericalHarmonics.jl

might make more sense. Though unfortunately I don't have time at the moment to flesh that out.

[mmikhasenko]

you can pin me once you get there

[mmikhasenko]

@dlfivefifty hey,
did you have any chance to look at the wigner d?

now, as we have got Symbolics.jl, I think, it should work very nicely with your way of computing polynomials.
I am very interested to explore that.

[dlfivefifty]

No I haven't. I don't think I have a need for Wigner d at the moment so no plans on implementing it, but happy to help if you are interested

[mmikhasenko]

ok, thanks.
For numerical computation, I have PartialWaveFunctions.jl.

SymPy works very well for me with analytic expressions.
I hope to use pure Julia one day, but it is not clear to me yet how and when I would be able to do similar things.
E.g. getting :(sqrt(5)) instead of a Float.

I believe that the recursive way of computing will be the way to go, but the details will change.