/*
* Euler - a numerical lab
*
* platform : neutral
*
* file : special.h -- special math functions
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <ctype.h>
#include "output.h"
#include "spread.h"
static double polevl(double x, double coef[], int N )
{
double ans;
int i;
double *p;
p = coef;
ans = *p++;
i = N;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}
static double p1evl(double x, double coef[], int N )
{
double ans;
double *p;
int i;
p = coef;
ans = x + *p++;
i = N-1;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}
static double P[] = {
1.60119522476751861407E-4,
1.19135147006586384913E-3,
1.04213797561761569935E-2,
4.76367800457137231464E-2,
2.07448227648435975150E-1,
4.94214826801497100753E-1,
9.99999999999999996796E-1
};
static double Q[] = {
-2.31581873324120129819E-5,
5.39605580493303397842E-4,
-4.45641913851797240494E-3,
1.18139785222060435552E-2,
3.58236398605498653373E-2,
-2.34591795718243348568E-1,
7.14304917030273074085E-2,
1.00000000000000000320E0
};
#define MAXGAM 171.624376956302725
static double LOGPI = 1.14472988584940017414;
static double STIR[5] = {
7.87311395793093628397E-4,
-2.29549961613378126380E-4,
-2.68132617805781232825E-3,
3.47222221605458667310E-3,
8.33333333333482257126E-2,
};
#define MAXSTIR 143.01608
static double SQTPI = 2.50662827463100050242E0;
static int sgngam = 0;
//extern int sgngam;
extern double MAXLOG,MAXNUM;
#define PI M_PI
/* Gamma function computed by Stirling's formula.
* The polynomial STIR is valid for 33 <= x <= 172.
*/
static double stirf (double x)
{ double y, w, v;
w = 1.0/x;
w = 1.0 + w * polevl( w, STIR, 4 );
y = exp(x);
if( x > MAXSTIR )
{ v = pow( x, 0.5 * x - 0.25 );
y = v * (v / y);
}
else
{ y = pow( x, x - 0.5 ) / y;
}
y = SQTPI * y * w;
return( y );
}
static double gamm (double x)
{
double p, q, z;
int i;
sgngam = 1;
q = fabs(x);
if( q > 33.0 )
{ if( x < 0.0 )
{ p = floor(q);
if( p == q )
goto goverf;
i = (int)p;
if( (i & 1) == 0 )
sgngam = -1;
z = q - p;
if( z > 0.5 )
{
p += 1.0;
z = q - p;
}
z = q * sin( PI * z );
if( z == 0.0 )
{
goverf:
output("Overflow in gamma\n");
error=1; return 0;
}
z = fabs(z);
z = PI/(z * stirf(q) );
}
else
{
z = stirf(x);
}
return( sgngam * z );
}
z = 1.0;
while( x >= 3.0 )
{
x -= 1.0;
z *= x;
}
while( x < 0.0 )
{
if( x > -1.E-9 )
goto small;
z /= x;
x += 1.0;
}
while( x < 2.0 )
{
if( x < 1.e-9 )
goto small;
z /= x;
x += 1.0;
}
if( (x == 2.0) || (x == 3.0) )
return(z);
x -= 2.0;
p = polevl( x, P, 6 );
q = polevl( x, Q, 7 );
return( z * p / q );
small:
if( x == 0.0 )
{
output("Wrong argument for gamma.\n");
error=1; return 0;
}
else
return( z/((1.0 + 0.5772156649015329 * x) * x) );
}
static double A[] = {
8.11614167470508450300E-4,
-5.95061904284301438324E-4,
7.93650340457716943945E-4,
-2.77777777730099687205E-3,
8.33333333333331927722E-2
};
static double B[] = {
-1.37825152569120859100E3,
-3.88016315134637840924E4,
-3.31612992738871184744E5,
-1.16237097492762307383E6,
-1.72173700820839662146E6,
-8.53555664245765465627E5
};
static double C[] = {
/* 1.00000000000000000000E0, */
-3.51815701436523470549E2,
-1.70642106651881159223E4,
-2.20528590553854454839E5,
-1.13933444367982507207E6,
-2.53252307177582951285E6,
-2.01889141433532773231E6
};
/* log( sqrt( 2*pi ) ) */
static double LS2PI = 0.91893853320467274178;
#define MAXLGM 2.556348e305
/* Logarithm of gamma function */
static double gammln (double x)
{
double p, q, w, z;
int i;
sgngam = 1;
if( x < -34.0 )
{
q = -x;
w = gammln(q); /* note this modifies sgngam! */
p = floor(q);
if( p == q )
goto loverf;
i = (int)p;
if( (i & 1) == 0 )
sgngam = -1;
else
sgngam = 1;
z = q - p;
if( z > 0.5 )
{
p += 1.0;
z = p - q;
}
z = q * sin( PI * z );
if( z == 0.0 )
goto loverf;
/* z = log(PI) - log( z ) - w;*/
z = LOGPI - log( z ) - w;
return( z );
}
if( x < 13.0 )
{
z = 1.0;
while( x >= 3.0 )
{
x -= 1.0;
z *= x;
}
while( x < 2.0 )
{
if( x == 0.0 )
goto loverf;
z /= x;
x += 1.0;
}
if( z < 0.0 )
{
sgngam = -1;
z = -z;
}
else
sgngam = 1;
if( x == 2.0 )
return( log(z) );
x -= 2.0;
p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
return( log(z) + p );
}
if( x > MAXLGM )
{
loverf:
output("Overflow in loggamma\n");
error=1; return 0;
}
q = ( x - 0.5 ) * log(x) - x + LS2PI;
if( x > 1.0e8 )
return( q );
p = 1.0/(x*x);
if( x >= 1000.0 )
q += (( 7.9365079365079365079365e-4 * p
- 2.7777777777777777777778e-3) *p
+ 0.0833333333333333333333) / x;
else
q += polevl( p, A, 4 ) / x;
return( q );
}
void mgammaln (header *hd)
{ spread1(gammln,0,hd);
test_error("gammaln");
}
void mgamma (header *hd)
{ spread1(gamm,0,hd);
test_error("gamma");
}
#define ITMAX 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30
static double gser (double a, double x)
{ int n;
double sum,del,ap;
if (x <= 0.0) return 0;
else
{ ap=a;
del=sum=1.0/a;
for (n=1;n<=ITMAX;n++) {
++ap;
del *= x / ap;
sum += del;
if (fabs(del) < fabs(sum)*EPS) {
return sum*exp(-x+a *log(x));
}
}
output("Iteration fails in gamma\n");
error=1;
return 0;
}
}
static double gcf (double a, double x)
{
int i;
double an,b,c,d,del,h;
b=x+1.0-a;
c=1.0/FPMIN;
d=1.0/b;
h=d;
for (i=1;i<=ITMAX;i++) {
an = -i*(i-a);
b += 2.0;
d=an*d+b;
if (fabs(d) < FPMIN) d=FPMIN;
c=b+an/c;
if (fabs(c) < FPMIN) c=FPMIN;
d=1.0/d;
del=d*c;
h *= del;
if (fabs(del-1.0) < EPS) break;
}
if (i > ITMAX)
{ output("a too large, ITMAX too small in gcf\n");
error=1; return 0;
}
return exp(-x+a*log(x))*h;
}
static void gammp (double *a, double *x, double *res)
{
if (*x < 0.0 || *a <= 0.0)
{ output("Invalid arguments in routine gammp\n");
error=1; return;
}
if (*x < (*a+1.0)) {
*res=gser(*a,*x);
return;
} else {
*res=gamm(*a)-gcf(*a,*x);
return;
}
}
void mgammai (header *hd)
{ spreadf2(gammp,0,0,hd);
test_error("gamma");
}
static double gauss (double z)
{ double res,a=0.5;
double x=z*z/2;
gammp(&a,&x,&res);
res/=gamm(0.5)*2;
if (z<0)
return 0.5-res;
else
return 0.5+res;
}
void mgauss (header *hd)
{ spread1(gauss,0,hd);
test_error("normaldis");
}
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define RNMX (1.0-EPS)
static long randseed=1234512345;
#define IDUM2 123456789
static long idum2=IDUM2;
static long iy=0;
static long iv[NTAB];
static double ran2 (void)
{
int j;
long k;
double temp;
if (randseed <= 0) {
if (-(randseed) < 1) randseed=1;
else randseed = -(randseed);
idum2=(randseed);
for (j=NTAB+7;j>=0;j--) {
k=(randseed)/IQ1;
randseed=IA1*(randseed-k*IQ1)-k*IR1;
if (randseed < 0) randseed += IM1;
if (j < NTAB) iv[j] = randseed;
}
iy=iv[0];
}
k=(randseed)/IQ1;
randseed=IA1*(randseed-k*IQ1)-k*IR1;
if (randseed < 0) randseed += IM1;
k=idum2/IQ2;
idum2=IA2*(idum2-k*IQ2)-k*IR2;
if (idum2 < 0) idum2 += IM2;
j=iy/NDIV;
iy=iv[j]-idum2;
iv[j] = randseed;
if (iy < 1) iy += IMM1;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
void mrandom (header *hd)
{ header *st=hd,*result;
double *m;
int r,c;
long k,n;
hd=getvalue(hd); if (error) return;
if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2
|| *(m=matrixof(hd))<0 || *m>=INT_MAX
|| *(m+1)<0 || *(m+1)>INT_MAX)
wrong_arg_in("random");
r=(int)*m;
c=(int)*(m+1);
result=new_matrix(r,c,""); if (error) return;
m=matrixof(result);
n=(long)c*r;
for (k=0; k<n; k++) *m++=(double)ran2();
moveresult(st,result);
}
void mshuffle (header *hd)
{ header *st=hd,*result;
double *m,*mr,x;
int i,j,n;
hd=getvalue(hd); if (error) return;
if (hd->type!=s_matrix || dimsof(hd)->r!=1)
wrong_arg_in("shuffle");
n=dimsof(hd)->c;
result=new_matrix(1,n,"");
m=matrixof(hd); mr=matrixof(result);
for (i=0; i<n; i++) *mr++=*m++;
mr=matrixof(result);
for (i=n-1; i>0; i--)
{ j=(int)floor(ran2()*(i+1));
if (i!=j)
{ x=*(mr+i); *(mr+i)=*(mr+j); *(mr+j)=x;
}
}
moveresult(st,result);
}
static double gasdev (void)
{ static int iset=0;
static double gset;
double fac,rsq,v1,v2;
if (iset == 0)
{ do
{ v1=2.0*ran2()-1.0;
v2=2.0*ran2()-1.0;
rsq=v1*v1+v2*v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac=sqrt(-2.0*log(rsq)/rsq);
gset=v1*fac;
iset=1;
return v2*fac;
}
else
{ iset=0;
return gset;
}
}
void mnormalnew (header *hd)
{ header *st=hd,*result;
double *m;
int r,c;
long k,n;
hd=getvalue(hd); if (error) return;
if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2
|| *(m=matrixof(hd))<0 || *m>=INT_MAX
|| *(m+1)<0 || *(m+1)>INT_MAX)
wrong_arg_in("normal");
r=(int)*m;
c=(int)*(m+1);
result=new_matrix(r,c,""); if (error) return;
m=matrixof(result);
n=(long)c*r;
for (k=0; k<n; k++) *m++=(double)gasdev();
moveresult(st,result);
}
void mseed (header *hd)
{ header *st=hd,*result;
hd=getvalue(hd); if (error) return;
if (hd->type!=s_real) wrong_arg_in("seed");
result=new_real(*realof(hd),"");
randseed=-labs((long)(*realof(hd)*LONG_MAX));
moveresult(st,result);
}
/************ BETA Functions *******************/
#undef ITMAX
#define ITMAX 100
#undef FPMIN
#define FPMIN 1e-30
static double betacf (double a, double b, double x)
{
int m,m2;
double aa,c,d,del,h,qab,qam,qap;
qab=a+b;
qap=a+1.0;
qam=a-1.0;
c=1.0;
d=1.0-qab*x/qap;
if (fabs(d) < FPMIN) d=FPMIN;
d=1.0/d;
h=d;
for (m=1;m<=ITMAX;m++) {
m2=2*m;
aa=m*(b-m)*x/((qam+m2)*(a+m2));
d=1.0+aa*d;
if (fabs(d) < FPMIN) d=FPMIN;
c=1.0+aa/c;
if (fabs(c) < FPMIN) c=FPMIN;
d=1.0/d;
h *= d*c;
aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
d=1.0+aa*d;
if (fabs(d) < FPMIN) d=FPMIN;
c=1.0+aa/c;
if (fabs(c) < FPMIN) c=FPMIN;
d=1.0/d;
del=d*c;
h *= del;
if (fabs(del-1.0) < EPS) break;
}
if (m > ITMAX)
{ output( "a or b too big, or ITMAX too small in betai\n" );
error=1;
}
return h;
}
static double betai (double x, double a, double b)
{ double bt;
if (x < 0.0 || x > 1.0)
{ output( "x not in [0,1] in betai\n" );
error=1; return 0;
}
if (x == 0.0 || x == 1.0) bt=0.0;
else
bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
if (x < (a+1.0)/(a+b+2.0))
return bt*betacf(a,b,x)/a;
else
return 1.0-bt*betacf(b,a,1.0-x)/b;
}
void mbetai (header *hd)
{ header *result,*st=hd,*hda,*hdb;
double x;
hda=nextof(hd);
hdb=nextof(hda);
hd=getvalue(hd); if (error) return;
hda=getvalue(hda); if (error) return;
hdb=getvalue(hdb); if (error) return;
if (hd->type!=s_real || hda->type!=s_real ||
hdb->type!=s_real)
wrong_arg_in("betai");
x=betai(*realof(hd),*realof(hda),*realof(hdb));
if (error) return;
result=new_real(x,""); if (error) return;
moveresult(st,result);
}
static double chebev (double a, double b, double c[], int m, double x)
{
double d=0.0,dd=0.0,sv,y,y2;
int j;
if ((x-a)*(x-b) > 0.0)
{ output("x not in range in routine bessel\n");
error = 1; return 0;
}
y2=2.0*(y=(2.0*x-a-b)/(b-a));
for (j=m-1;j>=1;j--) {
sv=d;
d=y2*d-dd+c[j];
dd=sv;
}
return y*d-dd+0.5*c[0];
}
#define NUSE1 7
#define NUSE2 8
static void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi)
{ double xx;
static double c1[] = {
-1.142022680371172e0,6.516511267076e-3,
3.08709017308e-4,-3.470626964e-6,6.943764e-9,
3.6780e-11,-1.36e-13};
static double c2[] = {
1.843740587300906e0,-0.076852840844786e0,
1.271927136655e-3,-4.971736704e-6,-3.3126120e-8,
2.42310e-10,-1.70e-13,-1.0e-15};
xx=8.0*x*x-1.0;
*gam1=chebev(-1.0,1.0,c1,NUSE1,xx); if (error) return;
*gam2=chebev(-1.0,1.0,c2,NUSE2,xx); if (error) return;
*gampl= *gam2-x*(*gam1);
*gammi= *gam2+x*(*gam1);
}
#undef NUSE1
#undef NUSE2
static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
(imaxarg1) : (imaxarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#undef EPS
#undef FPMIN
#define EPS 1.0e-16
#define FPMIN 1.0e-30
#define MAXIT 10000
#define XMIN 2.0
static void bessjy (double xnu, double x, double *rj, double *ry,
double *rjp, double *ryp)
{
int i,isign,l,nl;
double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl,
rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1,
temp,w,x2,xi,xi2,xmu,xmu2;
if (x <= 0.0 || xnu < 0.0)
{
output("bad arguments in bessel\n");
error = 1; return;
}
nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5)));
xmu=xnu-nl;
xmu2=xmu*xmu;
xi=1.0/x;
xi2=2.0*xi;
w=xi2/M_PI;
isign=1;
h=xnu*xi;
if (h < FPMIN) h=FPMIN;
b=xi2*xnu;
d=0.0;
c=h;
for (i=1;i<=MAXIT;i++) {
b += xi2;
d=b-d;
if (fabs(d) < FPMIN) d=FPMIN;
c=b-1.0/c;
if (fabs(c) < FPMIN) c=FPMIN;
d=1.0/d;
del=c*d;
h=del*h;
if (d < 0.0) isign = -isign;
if (fabs(del-1.0) < EPS) break;
}
if (i > MAXIT)
{ output("x too large in bessel\n");
error=1; return;
}
rjl=isign*FPMIN;
rjpl=h*rjl;
rjl1=rjl;
rjp1=rjpl;
fact=xnu*xi;
for (l=nl;l>=1;l--) {
rjtemp=fact*rjl+rjpl;
fact -= xi;
rjpl=fact*rjtemp-rjl;
rjl=rjtemp;
}
if (rjl == 0.0) rjl=EPS;
f=rjpl/rjl;
if (x < XMIN) {
x2=0.5*x;
pimu=M_PI*xmu;
fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
d = -log(x2);
e=xmu*d;
fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
beschb(xmu,&gam1,&gam2,&gampl,&gammi);
ff=2.0/M_PI*fact*(gam1*cosh(e)+gam2*fact2*d);
e=exp(e);
p=e/(gampl*M_PI);
q=1.0/(e*M_PI*gammi);
pimu2=0.5*pimu;
fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2);
r=M_PI*pimu2*fact3*fact3;
c=1.0;
d = -x2*x2;
sum=ff+r*q;
sum1=p;
for (i=1;i<=MAXIT;i++) {
ff=(i*ff+p+q)/(i*i-xmu2);
c *= (d/i);
p /= (i-xmu);
q /= (i+xmu);
del=c*(ff+r*q);
sum += del;
del1=c*p-i*del;
sum1 += del1;
if (fabs(del) < (1.0+fabs(sum))*EPS) break;
}
if (i > MAXIT)
{ output("bessrl series failed to converge\n");
error=1; return;
}
rymu = -sum;
ry1 = -sum1*xi2;
rymup=xmu*xi*rymu-ry1;
rjmu=w/(rymup-f*rymu);
} else {
a=0.25-xmu2;
p = -0.5*xi;
q=1.0;
br=2.0*x;
bi=2.0;
fact=a*xi/(p*p+q*q);
cr=br+q*fact;
ci=bi+p*fact;
den=br*br+bi*bi;
dr=br/den;
di = -bi/den;
dlr=cr*dr-ci*di;
dli=cr*di+ci*dr;
temp=p*dlr-q*dli;
q=p*dli+q*dlr;
p=temp;
for (i=2;i<=MAXIT;i++) {
a += 2*(i-1);
bi += 2.0;
dr=a*dr+br;
di=a*di+bi;
if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN;
fact=a/(cr*cr+ci*ci);
cr=br+cr*fact;
ci=bi-ci*fact;
if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN;
den=dr*dr+di*di;
dr /= den;
di /= -den;
dlr=cr*dr-ci*di;
dli=cr*di+ci*dr;
temp=p*dlr-q*dli;
q=p*dli+q*dlr;
p=temp;
if (fabs(dlr-1.0)+fabs(dli) < EPS) break;
}
if (i > MAXIT)
{ output("cf2 failed in bessjy\n");
error=1; return;
}
gam=(p-f)/q;
rjmu=sqrt(w/((p-f)*gam+q));
rjmu=SIGN(rjmu,rjl);
rymu=rjmu*gam;
rymup=rymu*(p+q/gam);
ry1=xmu*xi*rymu-rymup;
}
fact=rjmu/rjl;
*rj=rjl1*fact;
*rjp=rjp1*fact;
for (i=1;i<=nl;i++) {
rytemp=(xmu+i)*xi2*ry1-rymu;
rymu=ry1;
ry1=rytemp;
}
*ry=rymu;
*ryp=xnu*xi*rymu-ry1;
}
#undef EPS
#undef FPMIN
#undef MAXIT
#undef XMIN
#undef PI
static void besselj (double *x, double *k, double *r)
{ double h1,h2;
bessjy(*k,*x,r,&h1,&h2,&h2);
}
void mbesselj (header *hd)
{ spreadf2(besselj,0,0,hd);
test_error("besselj");
}
static void bessely (double *x, double *k, double *r)
{ double h1,h2;
bessjy(*k,*x,&h1,r,&h2,&h2);
}
void mbessely (header *hd)
{ spreadf2(bessely,0,0,hd);
test_error("bessely");
}
void mbesselall (header *hdx)
{ header *st=hdx,*hd=nextof(hdx);
double bj,by,bdj,bdy;
hd=getvalue(hd); if (error) return;
hdx=getvalue(hdx); if (error) return;
if (hd->type!=s_real || hdx->type!=s_real)
wrong_arg_in("besselall");
bessjy(*realof(hd),*realof(hdx),&bj,&by,&bdj,&bdy);
if (error) return;
newram=(char *)st;
new_real(bj,""); if (error) return;
new_real(by,""); if (error) return;
new_real(bdj,""); if (error) return;
new_real(bdy,"");
}
#define EPS 1.0e-16
#define FPMIN 1.0e-30
#define MAXIT 10000
#define XMIN 2.0
#define PI M_PI
static void bessik(double xnu, double x, double *ri, double *rk, double *rip, double *rkp)
{
void beschb(double x, double *gam1, double *gam2, double *gampl,
double *gammi);
void nrerror(char error_text[]);
int i,l,nl;
double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,gam1,gam2,
gammi,gampl,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,ripl,
ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2;
if (x <= 0.0 || xnu < 0.0)
{
output("bad arguments in bessik\n");
error = 1; return;
}
nl=(int)(xnu+0.5);
xmu=xnu-nl;
xmu2=xmu*xmu;
xi=1.0/x;
xi2=2.0*xi;
h=xnu*xi;
if (h < FPMIN) h=FPMIN;
b=xi2*xnu;
d=0.0;
c=h;
for (i=1;i<=MAXIT;i++) {
b += xi2;
d=1.0/(b+d);
c=b+1.0/c;
del=c*d;
h=del*h;
if (fabs(del-1.0) < EPS) break;
}
if (i > MAXIT)
{ output("x too large in bessik\n");
error=1; return;
}
ril=FPMIN;
ripl=h*ril;
ril1=ril;
rip1=ripl;
fact=xnu*xi;
for (l=nl;l>=1;l--) {
ritemp=fact*ril+ripl;
fact -= xi;
ripl=fact*ritemp+ril;
ril=ritemp;
}
f=ripl/ril;
if (x < XMIN) {
x2=0.5*x;
pimu=PI*xmu;
fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
d = -log(x2);
e=xmu*d;
fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
beschb(xmu,&gam1,&gam2,&gampl,&gammi);
ff=fact*(gam1*cosh(e)+gam2*fact2*d);
sum=ff;
e=exp(e);
p=0.5*e/gampl;
q=0.5/(e*gammi);
c=1.0;
d=x2*x2;
sum1=p;
for (i=1;i<=MAXIT;i++) {
ff=(i*ff+p+q)/(i*i-xmu2);
c *= (d/i);
p /= (i-xmu);
q /= (i+xmu);
del=c*ff;
sum += del;
del1=c*(p-i*ff);
sum1 += del1;
if (fabs(del) < fabs(sum)*EPS) break;
}
if (i > MAXIT)
{
output("bessk series failed to converge\n");
error = 1; return;
}
rkmu=sum;
rk1=sum1*xi2;
} else {
b=2.0*(1.0+x);
d=1.0/b;
h=delh=d;
q1=0.0;
q2=1.0;
a1=0.25-xmu2;
q=c=a1;
a = -a1;
s=1.0+q*delh;
for (i=2;i<=MAXIT;i++) {
a -= 2*(i-1);
c = -a*c/i;
qnew=(q1-b*q2)/a;
q1=q2;
q2=qnew;
q += c*qnew;
b += 2.0;
d=1.0/(b+a*d);
delh=(b*d-1.0)*delh;
h += delh;
dels=q*delh;
s += dels;
if (fabs(dels/s) < EPS) break;
}
if (i > MAXIT)
{
output("bessik: failure to converge in cf2\n");
error = 1; return;
}
h=a1*h;
rkmu=sqrt(PI/(2.0*x))*exp(-x)/s;
rk1=rkmu*(xmu+x+0.5-h)*xi;
}
rkmup=xmu*xi*rkmu-rk1;
rimu=xi/(f*rkmu-rkmup);
*ri=(rimu*ril1)/ril;
*rip=(rimu*rip1)/ril;
for (i=1;i<=nl;i++) {
rktemp=(xmu+i)*xi2*rk1+rkmu;
rkmu=rk1;
rk1=rktemp;
}
*rk=rkmu;
*rkp=xnu*xi*rkmu-rk1;
}
#undef EPS
#undef FPMIN
#undef MAXIT
#undef XMIN
#undef PI
static void besseli (double *x, double *k, double *r)
{ double h1,h2;
bessik(*k,*x,r,&h1,&h2,&h2);
}
void mbesseli (header *hd)
{ spreadf2(besseli,0,0,hd);
test_error("besseli");
}
static void besselk (double *x, double *k, double *r)
{ double h1,h2;
bessik(*k,*x,&h1,r,&h2,&h2);
}
void mbesselk (header *hd)
{ spreadf2(besselk,0,0,hd);
test_error("besselk");
}
void mbesselmodall (header *hdx)
{ header *st=hdx,*hd=nextof(hdx);
double bj,by,bdj,bdy;
hd=getvalue(hd); if (error) return;
hdx=getvalue(hdx); if (error) return;
if (hd->type!=s_real || hdx->type!=s_real)
wrong_arg_in("besselmodall");
bessik(*realof(hd),*realof(hdx),&bj,&by,&bdj,&bdy);
if (error) return;
newram=(char *)st;
new_real(bj,""); if (error) return;
new_real(by,""); if (error) return;
new_real(bdj,""); if (error) return;
new_real(bdy,"");
}
syntax highlighted by Code2HTML, v. 0.9.1