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++