# Problems with SpecialFunctions.jl
**URL:** https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705
**Category:** General Usage
**Tags:** question, package
**Created:** [January 28, 2023, 6:33pm UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705 "2023-01-28T18:33:01Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@fda-tome](https://discourse.julialang.org/u/fda-tome)
#### Post date: [January 28, 2023, 6:33pm UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/1 "2023-01-28T18:33:01Z")
Hey guys, im having trouble using BigFloat variables on spherical bessel functions given by the SpecialFunctions,jl package. I’m using the following code and receiving the error below. Thanks beforehand 🙂
```julia
function partialwavexp(psiamp, axiconang, order, r, theta, phi)
k = 1000 * 2 * big(pi) / 1.54
kr = k * r
krho = k * sin(axiconang)
kz = k * cos(axiconang)
phi0 = z0 = rho0 = 0
nmax = Integer(ceil(kr + (big(405) / 100) * (kr^(1/3)) + 2))
psi = 0
for n in 0:nmax
for m in -n:n
BSC = bsccalc(n, m, axiconang, order, krho, kz, phi0, z0, rho0)
psi += BSC * sphericalbesselj(n, kr) * Plm(n, m, cos(theta)) * cis(m * phi)
end
end
return psi * psiamp
end
```
And the error:
```julia
ERROR: LoadError: MethodError: no method matching besselj(::BigFloat, ::Complex{BigFloat})
Closest candidates are:
besselj(::T, ::Complex{T}) where T<:AbstractFloat at ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:590
besselj(::Real, ::Complex) at ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:584
besselj(::Real, ::AbstractFloat) at ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:482
...
Stacktrace:
[1] besselj(k::BigFloat, z::Complex{BigFloat})
@ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:590
[2] besselj(nu::BigFloat, x::BigFloat)
@ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:490
[3] sphericalbesselj(nu::BigInt, x::BigFloat)
@ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:701
[4] partialwavexp(psiamp::Int64, axiconang::BigFloat, order::Int64, r::BigFloat, theta::Float64, phi::BigFloat)
@ Main ~/Documentos/BSCSim/src/naive/naive_implem.jl:30
```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [January 28, 2023, 6:39pm UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/2 "2023-01-28T18:39:53Z")
if you need extended precision Bessel functions, they are provided by Arblib.jl.
---
### Author:  [@fda-tome](https://discourse.julialang.org/u/fda-tome)
#### Post date: [January 28, 2023, 6:57pm UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/3 "2023-01-28T18:57:06Z")
Thanks for the advice 🙂 , but i was not able to find spherical functions on the arb packages and documentation. Are there any other wrappers you could indicate please.
---
### Author:  [@heltonmc](https://discourse.julialang.org/u/heltonmc)
#### Post date: [January 28, 2023, 10:05pm UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/4 "2023-01-28T22:05:28Z")
Just use the relation to the normal Bessel function [http://dlmf.nist.gov/10.47.E3](https://dlmf.nist.gov/10.47.E3).
```julia
sphericalbesselj(nu, x::T) where T = sqrt(T(π)/2) * besselj(nu + one(T)/2, x) / sqrt(x)
```
You could do much better for performance though by moving the `sphericalbesselj` outside of the `m` loop and using recurrence. You should probably also move the `sphericalbesselj` completely outside the loop and use recurrence.
This is only stable in the backward direction so you would need to compute `sphericalbesselj(nmax, x)` and `sphericalbesselj(nmax-1, x)` and sum from `for n in nmax:-1:0` then you would only need two bessel function evaluations… anyway that requires a bit more care in your implementation but the speed gains will be immense if the other calculations in your loop are cheap.
Edit: My previous advice if you are evaluating for complex numbers might need more care for stability in recurrence.
Edit2: I also usually use `ArbNumerics.jl` but the default precision is not always enough if `nu` or `x` is large so you might need to increase the precision to get fully accurate results.
---
### Author:  [@fda-tome](https://discourse.julialang.org/u/fda-tome)
#### Post date: [January 29, 2023, 8:05pm UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/5 "2023-01-29T20:05:30Z")
Hello Helton, i cannot thank you enough for helping with this specific problem and giving valuable advice. bessel.jl developed by you and Oscar that also answered my question is a cornerstone on my new research project. I am still a julia newbie as you can see, but somehow in the future i still want to contribute to some JuliaMath repo due to how impressive they are. Again, thank you very much.
---
### Author:  [@heltonmc](https://discourse.julialang.org/u/heltonmc)
#### Post date: [January 30, 2023, 1:45am UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/6 "2023-01-30T01:45:11Z")
Of course! I’m glad you found it useful 🙂
Just to add another recommendation for you to consider. In general, you usually can avoid directly calculating Bessel (or Gamma) functions within the loop by using recurrence. Though, your problem sets up well for a different approach.
You could apply [Miller’s recurrence algorithm](https://en.wikipedia.org/wiki/Miller%27s_recurrence_algorithm). It will require some care for your desired precision and you need to determine how much higher your trial _M_ should be over _N_. I’m not for sure if you need an arbitrary precision routine but if you only need a set amount (say `Float128` or `Double64`) this number can be determined before hand.
The main advantage of this is that you must calculate `sphericalbesselj(0, kr)` at the start (or end of the loop) which has a [closed form expression](https://dlmf.nist.gov/10.49.E2) (`j_0=sin(x)/x`) this is important because you can use that as your normalization condition at the end. So, now instead of calculating `sphericalbesselj` everytime within your loop you just have to calculate `sin(x)/x` at the end of your loop one time and normalize your sum. The advantage of this is you can find fast quadruple or double-double implementations of `sin(x)` in `DoubleFloats.jl` and `QuadMath.jl` if you don’t need an arbitrary precision routine. I’m away traveling currently so can’t code this up but would be happy to try and answer any questions if you want to give it a try! This approach essentially avoids computing the Bessel function in higher precision which will give you a significant speed gain. Of course if performance doesn’t matter then the above suggestion will be much easier 🙂
---
### Author:  [@fda-tome](https://discourse.julialang.org/u/fda-tome)
#### Post date: [January 30, 2023, 2:05am UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/7 "2023-01-30T02:05:37Z")
Heeeey, Miller’s algorithm seems great. I will implement it tomorrow, probably using float128, as a first approach i don’t think i’ll be needing more precision than that. My concern solely is when x approaches 0 in `j_0=sin(x)/x`, but i’ll try to work around that!
Thank you so much again and i hope you have a great traveling experience.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [January 30, 2023, 2:06am UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/8 "2023-01-30T02:06:13Z")
In that case, the taylor series will save you.
---
### Author:  [@fda-tome](https://discourse.julialang.org/u/fda-tome)
#### Post date: [January 30, 2023, 2:08am UTC](https://discourse.julialang.org/t/problems-with-specialfunctions-jl/93705/9 "2023-01-30T02:08:50Z")
Nice, i’ll implement it using the taylor series expansion to some degree. I cannot stress how much you guys helped. Thanks a lot.