Besselj1.cin

// BesselJ[1,z] implementation:

z_type BesselJ1o(z_type z){ int n; z_type c,s,t; s=1.; c=1.; t=-z*z/4.; for(n=1;n<30;n++) {c/=0.+n*n; c*=t; s+=c/(1.+n);} return s*z/2.;}

z_type BesselJ1B(z_type z){ int n; z_type a,f,c,s,t,u,v; f=M_PI/4.+z; c=cos(f); s=sin(f); a=sqrt(2./M_PI/z); t=1./(z*z); u= (( (( (( (( ( 11587123363874217382000061956546875./302231454903657293676544.*t - 531595142508493798628972203125./1180591620717411303424.)*t + 61394155884765868567528125./9223372036854775808.)*t     -1149690375852815671875./9007199254740992.)*t + 232376754295310625./70368744177664.) * t          - 33424574007825./274877906944. ) * t             + 14783093325./2147483648. ) * t               - 2837835./4194304.      ) * t                  + 4725./32768. ) * t                  - 15./128.          ) *t -1. ; v= (((((((( ( - 38190914185478633427818266171875./9444732965739290427392.*t +3918391713821821611515765625./73786976294838206464.)*t -64152722972587114490625./72057594037927936.)*t  + 11100458801337530625./562949953421312.)*t - 1327867167401775./2199023255552. )*t        + 468131288625./17179869184.   )*t - 66891825./33554432.     )*t              + 72765./262144.        )*t - 105./1024.         )*t              + 3./8.             )/z; return (c*u+s*v)*a; }

z_type BesselJ1(z_type z){ DB x,y; x=Re(z)/12.6; y=Im(z)/18.5; if(x*x+y*y<1.) return BesselJ1o(z); if(x<0) return - BesselJ1B(-z); return  BesselJ1B(z); }

Keywords
BesselJ1, C++