RE: Bessel series

[Jacques Gélinas <jacquesg00@hotmail.com>]

> If you need the value at 0

Actually, I need values near 0, and the besselj series provides an elegant solution !
This seems to work:

H(p,x) ={
  if( abs(x) < 1. << - (bitprecision(1.)/2),
      eval(Pol(besselj(p,'x+O('x^3)))),
      gamma(p+1)*besselj(p,x)/(x/2)^p );
}


Thanks,

Jacques Gélinas







De : Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>
Envoyé : 17 octobre 2017 05:49
À : pari-users@pari.math.u-bordeaux.fr
Objet : Re: Bessel series
    
On Tue, Oct 17, 2017 at 12:17:17AM +0000, Jacques Gélinas wrote:
> How can I compute the entire function  H_p(x) = Gamma(p+1) J_p(x) / (x/2)^p 
> where J_p denotes the Bessel-J function, and p >  -1 ?
> For example, H_{-1/2} = cos, H_{1/2} = sinc.
> 
> # besselj(-1/2,0)
>   *** besselj: domain error in gpow(0,n): n <= 0

If you need the value at 0, you need to use the power series expanssion
taking into account the documentation:

  If x converts to a power series,  the initial factor
  (x/2)^nu/Gamma(nu+1)  is omitted  (since it cannot be represented in
  PARI when nu is not integral).

So 
? besselj(-1/2,O(x))
%11 = 1+O(x^2)

Maybe you want this:
H(p,x) = if(x,gamma(p+1)*besselj(p,x)/(x/2)^p,polcoeff(besselj(p,O('x)),0))

Cheers,
Bill.