/usr/share/yacas/scripts/specfunc.rep/bessel.ys is in yacas 1.3.6-2+b1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 | /// coded by Jonathan Leto
// When x is <= 1, the series is monotonely decreasing from the
// start, so we don't have to worry about loss of precision from the
// series definition.
// When {n} is an integer, this is fast.
// When {n} is not, it is pretty slow due to Gamma()
Function("BesselNsmall",{n,x,modified})
[
Local(term,result,k);
Local(prec,eps,tmp);
prec:=Builtin'Precision'Get();
Builtin'Precision'Set(Ceil(1.2*prec)); // this is a guess
eps:=5*10^(-prec);
term:=1;
k:=0;
result:=0;
While( Abs(term) >= eps )[
term:=x^(2*k+n);
// The only difference between BesselJ and BesselI
// is an alternating term
If( k%2=1 And modified=0 , term:=term*-1 );
term:=N(term/(2^(2*k+n)* k! * Gamma(k+n+1) ));
//Echo({"term is ",term});
result:=result+term;
k:=k+1;
];
Builtin'Precision'Set(prec);
// This should not round, only truncate
// some outputs will be off by one in the last digit
RoundTo(result,prec);
];
// Seems to get about 8 digits precision for most real numbers
// Only about 2 digits precision for complex
// This is just a temporary implementation, I would not want to
// expose users to it until it is much more robust
// I am still looking for a good arbitrary precision algorithm.
Function("BesselJN0",{x})
[
Local(ax,z,xx,y,result,res1,res2);
Local(c1,c2,c3,c4);
// Coefficients of the rational polynomials to
// approx J_0 for x < 8
c1:={57568490574.0,-13362590354.0,651619640.7,
-11214424.18,77392.33017,-184.9052456};
c2:={57568490411.0,1029532985.0,9494680.718,
59272.64853,267.8532712};
// Coefficients of the rational polynomials to
// approx J_0 for x >= 8
c3:={-0.001098628627,0.00002734510407,-0.000002073370639,
0.0000002093887211};
c4:={-0.01562499995,0.0001430488765,-0.000006911147651,
0.0000007621095161,0.0000000934935152};
ax:=Abs(x);
If( ax < 8.0,[
y:=x^2;
res1:=c1[1]+y*(c1[2]+y*c1[3]+y*(c1[4]+y*(c1[5]+y*(c1[6]))));
res2:=c1[1]+y*(c2[2]+y*c2[3]+y*(c2[4]+y*(c2[5]+y*1.0)));
result:=res1/res2;
],[
z:=8/ax;
y:=z^2;
xx:=ax-0.785398164;
res1:=1.0+y*(c3[1]+y*(c3[2]+y*(c3[3]+y*c4[4])));
res2:=c4[1]+y*(c4[2]+y*(c4[3]+y*(c4[4]-y*c4[5])));
result:=Sqrt(2/(Pi*x))*(Cos(xx)*res1-z*Sin(xx)*res2);
] );
];
|