On: Cosine Implementation in C

[nhatcher]

I recently had to implement some Bessel functions numerically. It's fantastic to see that someone at Sun Microsystems implemented all this procedures in the early 90's (I'm sure even long before that). You will find this implementations in R, Julia, musl, you name it!

I like it specially the use the "GET_HIGH_WORD" and friends:

https://git.musl-libc.org/cgit/musl/tree/src/math/j0.c

It's a lot of fun to figure out what they meant by things like:

> if (ix >= 0x7ff00000) return 1/(x*x);

[adgjlsfhk1 → nhatcher]

you should check out Bessels.jl which has really high performance implications.

[nhatcher → adgjlsfhk1]

https://github.com/JuliaMath/Bessels.jl/blob/master/src/bess...

Thanks! I love it, so easy to understand and follow.

My favourite work on the subject is Fredrik Johansson's:

https://github.com/fredrik-johansson/arb

Whenever I feel down and without energy I just read something in there

[adgjlsfhk1 → nhatcher]

yeah, Fredrik's papers are great. A lot of the reason Bessels.jl is so good is because ARB makes it easy for us to have test against a known source of truth (which make things like minimax polynomials and picking optimal tradeoffs a ton easier).

[nhatcher → adgjlsfhk1]

You work with Julia? Nice! I actually use Julia to generate my test cases:

  using Nemo

  CC = AcbField(100)

  values = [1, 2, 3, -2, 5, 30, 2e-8]

  for value in values
      y_acb = besseli(CC(1), CC(value))
      real64 = convert(Float64, real(y_acb))
      im64 = convert(Float64, real(y_acb))
      println("(", value, ", ", real64, "),")
  end