📄 d9function.java
字号:
System.out.println("Singular matrix");
System.exit(1);
}
PIVINV=1.0/(A[(ICOL-1)*MA+ICOL]);
A[(ICOL-1)*MA+ICOL]=1.0;
for(l=1;l<=N;l++)
{
A[(ICOL-1)*MA+l]=A[(ICOL-1)*MA+l]*PIVINV;
}
B[ICOL]=B[ICOL]*PIVINV;
for(ll=1;ll<=N;ll++)
{
if(ll!=ICOL)
{
DUM=A[(ll-1)*MA+ICOL];
A[(ll-1)*MA+ICOL]=0.0;
for(l=1;l<=N;l++)
{
A[(ll-1)*MA+l]=A[(ll-1)*MA+l]-A[(ICOL-1)*MA+l]*DUM;
}
B[ll]=B[ll]-B[ICOL]*DUM;
}
}
}
for(l=N;l>=1;l--)
{
if(INDXR[l]!=INDXC[l])
{
for(k=1;k<=N;k++)
{
DUM=A[(k-1)*MA+INDXR[l]];
A[(k-1)*MA+INDXR[l]]=A[(k-1)*MA+INDXC[l]];
A[(k-1)*MA+INDXR[l]]=DUM;
}
}
}
}
double gammq_a, gammq_x;
double gammq(double a, double x)
{
double temp,gamser,gammcf,gln;
gamser=0;gln=0;gammcf=0;
if (x < 0.0 || a <= 0.0)
{
System.out.println("pause in gammq");
System.exit(1);
}
if (x < a + 1.0)
{
gser(gamser, a, x, gln);
gamser=gser_gamser;a=gser_a;x=gser_x;gln=gser_gln;
temp= 1.0 - gamser;
}
else
{
gcf(gammcf, a, x, gln);
gammcf=gcf_gammcf;a=gcf_a;x=gcf_x;gln=gcf_gln;
temp = gammcf;
}
gammq_a=a;
gammq_x=x;
return temp;
}
double gcf_gammcf,gcf_a,gcf_x,gcf_gln;
void gcf(double gammcf, double a, double x, double gln)
{
int itmax,n;
double eps,a0,a1,b0,b1,fac,an,ana,anf,gold,g;
itmax = 100;
eps = 0.0000003;
gln = gammln(a);
gold = 0.0;
a0 = 1.0;
a1 = x;
b0 = 0.0;
b1 = 1.0;
fac = 1.0;
for (n = 1;n<=itmax;n++)
{
an = n;
ana = an - a;
a0 = (a1 + a0 * ana) * fac;
b0 = (b1 + b0 * ana) * fac;
anf = an * fac;
a1 = x * a0 + anf * a1;
b1 = x * b0 + anf * b1;
if (a1 != 0.0)
{
fac = 1.0 / a1;
g = b1 * fac;
if (Math.abs((g - gold) / g) < eps )
{
gammcf = Math.exp(-x + a * Math.log(x) - gln) * g;
break;
}
gold = g;
}
}
gcf_gammcf=gammcf;gcf_a=a;gcf_x=x;gcf_gln=gln;
}
double gser_gamser,gser_a,gser_x,gser_gln;
void gser(double gamser, double a, double x, double gln)
{
int itmax,n;
double ap,sum,del,eps;
itmax = 100;
eps = 0.0000003;
gln = gammln(a);
if (x <= 0.0)
{
if (x < 0.0)
{
System.out.println("pause in gser");
System.exit(1);
}
gamser = 0.0;
gser_gamser=gamser;gser_a=a;gser_x=x;gser_gln=gln;
return;
}
ap = a;
sum = 1.0 / a;
del = sum;
for( n = 1;n<=itmax;n++)
{
ap = ap + 1.0;
del = del * x / ap;
sum = sum + del;
if (Math.abs(del) < (Math.abs(sum) * eps))
{// goto loop;
gamser = sum * Math.exp(-x + a * Math.log(x) - gln);
break;
}
}
gser_gamser=gamser;gser_a=a;gser_x=x;gser_gln=gln;
}
double gammln(double xx)
{
int j;
double temp;
double cof[]=new double[7];
double stp,half,one,fpf,x,tmp,ser;
cof[1] = 76.18009173;
cof[2] = -86.50532033;
cof[3] = 24.01409822;
cof[4] = -1.231739516;
cof[5] = 0.00120858003;
cof[6] = -0.00000536382;
stp = 2.50662827465;
half = 0.5;
one = 1.0;
fpf = 5.5;
x = xx - one;
tmp = x + fpf;
tmp = (x + half) * Math.log(tmp) - tmp;
ser = one;
for (j = 1;j<=6;j++)
{
x = x + one;
ser = ser + cof[j] / x;
}
temp = tmp + Math.log(stp * ser);
return temp;
}
static int iset;
static double gset;
int gasdev_idum;
double gasdev(int idum)
{
double t,v1,v2,fac,r;
if (iset == 0.0)
{
do
{
v1 = 2.0 * ran1(idum) - 1.0;
idum=ran1_idum;
v2 = 2.0 * ran1(idum) - 1.0;
idum=ran1_idum;
r = v1 *v1 + v2 *v2;
}while ((r >= 1.0) || (r == 0.0));
fac = Math.sqrt(-2.0 * Math.log(r) / r);
gset = v1 * fac;
t = v2 * fac;
iset = 1;
gasdev_idum=idum;
return t;
}
else
{
t = gset;
gasdev_idum=idum;
iset = 0;
return t;
}
}
void sort(int n, double ra[])
{
int i,j,l,ir;
double rra;
l = (int)(n / 2) + 1;
ir = n;
do
{
if (l > 1)
{
l = l - 1;
rra = ra[l];
}
else
{
rra = ra[ir];
ra[ir] = ra[1];
ir = ir - 1;
if (ir == 1)
{
ra[1] = rra;
return;
}
}
i = l;
j = l + l;
while (j <= ir)
{
if (j < ir)
{
if (ra[j] < ra[j + 1])
{
j+=1;
}
}
if (rra < ra[j])
{
ra[i] = ra[j];
i = j;
j += j;
}
else
{
j = ir + 1;
}
}
ra[i] = rra;
}while(true);
}
void svdcmp(double a[],int m,int n,double w[],double v[])
{
int nmax=100;
double rv1[]=new double[nmax];
int i,j,l,k,nm=0;
double g,s,scale,f,h,c,x,y,z;
if(m<n)
{
System.out.println("you must augment a with extra zero rows.");
System.exit(1);
}
g=0.0;
scale=0.0;
double anorm=0.0;
l=0;
for(i=1; i<=n;i++)
{
l=i+1;
rv1[i]=scale*g;
g=0.0;
s=0.0;
scale=0.0;
if(i<=m)
{
for(k=i;k<=m;k++)
{
scale=scale+Math.abs(a[(k-1)*n+i]);
}
if(scale!=0.0)
{
for(k=i;k<=m;k++)
{
a[(k-1)*n+i]=a[(k-1)*n+i]/scale;
s=s+a[(k-1)*n+i]*a[(k-1)*n+i];
}
f=a[(i-1)*n+i];
g=-Math.sqrt(s)*sgn(f);
h=f*g-s;
a[(i-1)*n+i]=f-g;
if(i!=n)
{
for(j=l;j<=n;j++)
{
s=0.0;
for(k=i;k<=m;k++)
{
s=s+a[(k-1)*n+i]*a[(k-1)*n+j];
}
f=s/h;
for(k=i;k<=m;k++)
{
a[(k-1)*n+j]=a[(k-1)*n+j]+f*a[(k-1)*n+i];
}
}
}
for(k=i;k<=m;k++)
{
a[(k-1)*n+i]=scale*a[(k-1)*n+i];
}
}
}
w[i]=scale*g;
g=0.0;
s=0.0;
scale=0.0;
if ((i<=m)&&(i!=n))
{
for(k=l;k<=n;k++)
{
scale=scale+Math.abs(a[(i-1)*n+k]);
}
if (scale!=0.0)
{
for(k=l;k<=n;k++)
{
a[(i-1)*n+k]=a[(i-1)*n+k]/scale;
s=s+a[(i-1)*n+k]*a[(i-1)*n+k];
}
f=a[(i-1)*n+l];
g=-Math.sqrt(s)*sgn(f);
h=f*g-s;
a[(i-1)*n+l]=f-g;
for(k=l;k<=n;k++)
{
rv1[k]=a[(i-1)*n+k]/h;
}
if(i!=m)
{
for(j=l;j<=m;j++)
{
s=0.0;
for(k=l;k<=n;k++)
{
s=s+a[(j-1)*n+k]*a[(i-1)*n+k];
}
for(k=l;k<=n;k++)
{
a[(j-1)*n+k]=a[(j-1)*n+k]+s*rv1[k];
}
}
}
for(k=l;k<=n;k++)
{
a[(i-1)*n+k]=scale*a[(i-1)*n+k];
}
}
}
anorm=Math.max(anorm,(Math.abs(w[i])+Math.abs(rv1[i])));
}
for(i=n;i>=1;i--)
{
if(i<n)
{
if(g!=0.0)
{
for(j=l;j<=n;j++)
{
v[(j-1)*n+i]=(a[(i-1)*n+j]/a[(i-1)*n+l])/g;
}
for(j=l;j<=n;j++)
{
s=0.0;
for(k=l;k<=n;k++)
{
s=s+a[(i-1)*n+k]*v[(k-1)*n+j];
}
for(k=l;k<=n;k++)
{
v[(k-1)*n+j]=v[(k-1)*n+j]+s*v[(k-1)*n+i];
}
}
}
for(j=l;j<=n;j++)
{
v[(i-1)*n+j]=0.0;
v[(j-1)*n+i]=0.0;
}
}
v[(i-1)*n+i]=1.0;
g=rv1[i];
l=i;
}
for(i=n;i>=1;i--)
{
l=i+1;
g=w[i];
if (i<n)
{
for(j=l;j<=n;j++)
{
a[(i-1)*n+j]=0.0;
}
}
if (g!=0.0)
{
g=1.0/g;
if (i!=n)
{
for(j=l;j<=n;j++)
{
s=0.0;
for(k=l;k<=m;k++)
{
s=s+a[(k-1)*n+i]*a[(k-1)*n+j];
}
f=(s/a[(i-1)*n+i])*g;
for(k=i;k<=m;k++)
{
a[(k-1)*n+j]=a[(k-1)*n+j]+f*a[(k-1)*n+i];
}
}
}
for(j=i;j<=m;j++)
{
a[(j-1)*n+i]=a[(j-1)*n+i]*g;
}
}
else
{
for(j=i;j<=m;j++)
{
a[(j-1)*n+i]=0.0;
}
}
a[(i-1)*n+i]=a[(i-1)*n+i]+1.0;
}
for(k=n;k>=1;k--)
{
for(int its=1;its<=30;its++)
{
for(l=k;l>=1;l--)
{
nm=l-1;
if((Math.abs(rv1[l])+anorm)==anorm)
{
break;
}
if((Math.abs(w[nm])+anorm)==anorm)
{
break;
}
}
if(Math.abs(rv1[l])+anorm!=anorm)
{
c=0.0;
s=1.0;
for(i=l;i<=k;i++)
{
f=s*rv1[i];
if((Math.abs(f)+anorm)!=anorm)
{
g=w[i];
h=Math.sqrt(f*f+g*g);
w[i]=h;
h=1.0/h;
c=(g*h);
s=-(f*h);
for(j=1;j<=m;j++)
{
y=a[(j-1)*n+nm];
z=a[(j-1)*n+i];
a[(j-1)*n+nm]=(y*c)+(z*s);
a[(j-1)*n+i]=-(y*s)+(z*c);
}
}
}
}
z=w[k];
if(l==k)
{
if(z<0.0)
{
w[k]=-z;
for(j=1;j<=n;j++)
{
v[(j-1)*n+k]=-v[(j-1)*n+k];
}
}
break;
}
if(its==30)
{
System.out.println("no convergence in 30 iterations");
System.exit(1);
}
x=w[l];
nm=k-1;
y=w[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=Math.sqrt(f*f+1.0);
f=((x-z)*(x+z)+h*((y/(f+Math.abs(g)*sgn(f)))-h))/x;
c=1.0;
s=1.0;
for(j=l;j<=nm;j++)
{
i=j+1;
g=rv1[i];
y=w[i];
h=s*g;
g=g*c;
z=Math.sqrt(f*f+h*h);
rv1[j]=z;
c=f/z;
s=h/z;
f=(x*c)+(g*s);
g=-(x*s)+(g*c);
h=y*s;
y=y*c;
for(int pp=1;pp<=n;pp++)
{
x=v[(pp-1)*n+j];
z=v[(pp-1)*n+i];
v[(pp-1)*n+j]=(x*c)+(z*s);
v[(pp-1)*n+i]=-(x*s)+(z*c);
}
z=Math.sqrt(f*f+h*h);
w[j]=z;
if(z!=0.0)
{
z=1.0/z;
c=f*z;
s=h*z;
}
f=(c*g)+(s*y);
x=-(s*g)+(c*y);
for(int qq=1;qq<=m;qq++)
{
y=a[(qq-1)*n+j];
z=a[(qq-1)*n+i];
a[(qq-1)*n+j]=(y*c)+(z*s);
a[(qq-1)*n+i]=-(y*s)+(z*c);
}
}
rv1[l]=0.0;
rv1[k]=f;
w[k]=x;
}
}
}
void svbksb(double u[], double w[], double v[], int m, int n, double b[], double x[])
{
double tmp[]=new double[100];
int i,j,jj;
double s;
for (j = 1; j<=n; j++)
{
s = 0.0;
if (w[j] != 0.0)
{
for (i = 1; i<=m; i++)
{
s = s + u[(i-1)*n+j] * b[i];
}
s = s / w[j];
}
tmp[j] = s;
}
for (j = 1; j<=n; j++)
{
s = 0.0;
for (jj = 1; jj<=n; jj++)
{
s = s + v[(j-1)*n+jj] * tmp[jj];
}
x[j] = s;
}
}
int ran1_idum;
static double r[]=new double[98];
static int ix1,ix2,ix3;
static int iff;
double ran1(int idum)
{
double rm1,rm2,t;
int m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3,j;
m1 = 259200; ia1 = 7141; ic1 = 54773; rm1 = 0.0000038580247;
m2 = 134456; ia2 = 8121; ic2 = 28411; rm2 = 0.0000074373773;
m3 = 243000; ia3= 4561; ic3 = 51349;
if ((idum < 0) || (iff == 0))
{
iff = 1;
ix1 = (ic1 - idum) % m1;
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = ix1 % m2;
ix1 = (ia1 * ix1 + ic1) % m1;
ix3 = ix1 % m3;
for (j = 1; j<=97; j++)
{
ix1 = (ia1 * ix1 + ic1)% m1;
ix2 = (ia2 * ix2 + ic2)% m2;
r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
}
idum = 1;
}
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = (ia2 * ix2 + ic2) % m2;
ix3 = (ia3 * ix3 + ic3) % m3;
j = 1 + (int)((97 * ix3) / m3);
if ((j > 97) || (j < 1))
{
System.out.println( "abnormal exit");
System.exit(1);
}
t = r[j];
r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
ran1_idum=idum;
return t;
}
int sgn(double pa)
{
if (pa>0.0)
{
return 1;
}
else
{
if (pa<0.0)
{
return -1;
}
}
return 0;
}
double plgndr(double l, double m, double x)
{
int ll,i;
double pmm,somx2,fact,temp,pmmp1,pll=0;
if( m < 0 || m > l || Math.abs(x) > 1.0)
{
System.out.println("bad arguments");
System.exit(1);
}
pmm = 1.0;
if (m > 0)
{
somx2 = Math.sqrt((1.0 - x) * (1.0 + x));
fact = 1.0;
for( i = 1; i<=m; i++)
{
pmm = -pmm * fact * somx2;
fact = fact + 2.0;
}
}
if (l == m)
{
temp = pmm;
}
else
{
pmmp1 = x * (2 * m + 1) * pmm;
if (l == m + 1)
{
temp = pmmp1;
}
else
{
for (ll = (int)m + 2; ll<=l; ll++)
{
pll = x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm;
pll = pll / (ll - m);
pmm = pmmp1;
pmmp1 = pll;
}
temp = pll;
}
}
return temp;
}
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -