Bessel function

An Bessel function of the first kind of order $\nu$ can be computed using Taylor expansion

$$$J_\nu(z) = \sum\limits_{n=0}^{\infty} \frac{(z/2)^\nu}{\Gamma(k+1)\Gamma(k+\nu+1)} (-z^2/4)^{n}$$$

where $\Gamma(n) = (n-1)!$ is the Gamma function. One can compute the accumulated item iteratively as $s_n = -\frac{z^2}{4} s_{n-1}$.

using NiLang, NiLang.AD
using ForwardDiff: Dual

Since we need to use logarithmic numbers to handle the sequential mutiplication. Let's first add patch about the conversion between ULogarithmic and Dual number.

function Base.convert(::Type{Dual{T,V,N}}, x::ULogarithmic) where {T,V,N}
Dual{T,V,N}(exp(x.log))
end

@i function ibesselj(y!::T, ν, z::T; atol=1e-8) where T
if z == 0
if v == 0
out! += 1
end
else
@routine @invcheckoff begin
k ← 0
@ones ULogarithmic{T} lz halfz halfz_power_2 s
@zeros T out_anc
lz *= convert(z)
halfz *= lz / 2
halfz_power_2 *= halfz ^ 2
# s *= (z/2)^ν/ factorial(ν)
s *= halfz ^ ν
for i=1:ν
s /= i
end
out_anc += convert(s)
@from k==0 while s.log > -25 # upto precision e^-25
k += 1
# s *= 1 / k / (k+ν) * (z/2)^2
s *= halfz_power_2 / (@const k*(k+ν))
if k%2 == 0
out_anc += convert(s)
else
out_anc -= convert(s)
end
end
end
y! += out_anc
~@routine
end
end

y, x = 0.0, 1.0
Grad(ibesselj)(Val(1), y, 2, x)
(Val{1}(), GVar(0.0, 1.0), 2, GVar(1.0, 0.21024361588126536))

Here, Grad(ibesselj) is a callable instance of type Grad{typeof(ibesselj)}}. The first parameter Val(1) indicates the first argument is the loss.

To obtain second order gradients, one can Feed dual numbers to this gradient function.

_, hxy, _, hxx = Grad(ibesselj)(Val(1), Dual(y, zero(y)), 2, Dual(x, one(x)))