numerical_recipes.h

一些通信编程中常用到的数值解决方法的头文件

{
	int i,imax,j,k;
	float big,dum,sum,temp;
	float *vv;

	vv=vector(1,n);
	*d = 1.0;
	for (i=1;i<=n;i++)
	{
		big=0.0;
		for (j=1;j<=n;j++)
		{
			if ((temp=fabs(a[i][j])) > big)
				big=temp;
		}
		if (big == 0.0) 
			nrerror("Singular matrix in routine ludcmp");
		vv[i] = 1.0/big;
	}
	for (j=1;j<=n;j++)
	{
		for (i=1;i<j;i++) 
		{
			sum=a[i][j];
			for (k=1;k<i;k++) 
				sum -= a[i][k]*a[k][j];
			a[i][j] = sum;
		}
		big=0.0;
		for (i=j;i<=n;i++)
		{
			sum=a[i][j];
			for (k=1;k<j;k++)
				sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*fabs(sum)) >= big)
			{
				big = dum;
				imax = i;
			}
		}
		if (j != imax) 
		{
			for (k=1;k<=n;k++)
			{
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			*d = -(*d);
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (a[j][j] == 0.0)
			a[j][j] = TINY;
		if (j != n)
		{
			dum=1.0/(a[j][j]);
			for (i=j+1;i<=n;i++) 
				a[i][j] *= dum;
		}
	}
	free_vector(vv,1,n);
}

/** ******lubksb************* **/
void lubksb(float **a, int n, int *indx, float b[])
{
	int i,ii=0,ip,j;
	float sum;

	for (i=1;i<=n;i++) 
	{
		ip=indx[i];
		sum=b[ip];
		b[ip]=b[i];
		if (ii)
			for (j=ii;j<=i-1;j++) 
				sum -= a[i][j]*b[j];
		else if (sum) ii=i;
		b[i]=sum;
	}
	for (i=n;i>=1;i--)
	{
		sum=b[i];
		for (j=i+1;j<=n;j++)
			sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];
	}
}

/** ******maint function to compute Ax=b********* **/
void inverse(float ** a,int n, float *b)
{
	int *index;
	float d;
    void ludcmp(float **a, int n, int *indx, float *d);
    void lubksb(float **a, int n, int *indx, float b[]);
    
	index=ivector(1,n);
	ludcmp(a,n,index,&d);
    lubksb(a,n,index,b); 

}

/* *************************** */
/* *** bessel function***** ** */
/* *************************** */

/** ***zero order function*** **/
float bessj0(float x)
{
	float ax,z;
	float xx,y,ans,ans1,ans2;

	if ((ax = fabs(x)) < 8.0)
	{
		y=x*x;
		ans1 = 57568490574.0+y*(-13362590354.0+y*(651619640.7
			+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
		ans2 = 57568490411.0+y*(1029532985.0+y*(9494680.718
			+y*(59272.64853+y*(267.8532712+y*1.0))));
		ans = ans1/ans2;
	} 
	else 
	{
		z = 8.0/ax;
		y = z*z;
		xx = ax-0.785398164;
		ans1 = 1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			+y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			-y*0.934935152e-7)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
	}
	return ans;
}

/** ***first order function*** **/
float bessj1(float x)
{
	float ax,z;
	float xx,y,ans,ans1,ans2;

	if ((ax = fabs(x)) < 8.0) 
	{
		y = x*x;
		ans1 = x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
			+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
		ans2 = 144725228442.0+y*(2300535178.0+y*(18583304.74
			+y*(99447.43394+y*(376.9991397+y*1.0))));
		ans = ans1/ans2;
	} 
	else
	{
		z = 8.0/ax;
		y = z*z;
		xx = ax-2.356194491;
		ans1 = 1.0+y*(0.183105e-2+y*(-0.3516396496e-4
			+y*(0.2457520174e-5+y*(-0.240337019e-6))));
		ans2 = 0.04687499995+y*(-0.2002690873e-3
			+y*(0.8449199096e-5+y*(-0.88228987e-6
			+y*0.105787412e-6)));
		ans = sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
		if (x < 0.0) ans = -ans;
	}
	return ans;
}
/** ******n-th order bessel function**** **/
float bessj(int n, float x)
{
	float bessj0(float x);
	float bessj1(float x);
	void nrerror(char error_text[]);
	int j,jsum,m;
	float ax,bj,bjm,bjp,sum,tox,ans;

	if (n < 2) nrerror("Index n less than 2 in bessj");
	ax = fabs(x);
	if (ax == 0.0)
		return 0.0;
	else if (ax > (float) n) 
	{
		tox=2.0/ax;
		bjm=bessj0(ax);
		bj=bessj1(ax);
		for (j=1;j<n;j++) {
			bjp=j*tox*bj-bjm;
			bjm=bj;
			bj=bjp;
		}
		ans = bj;
	} 
	else
	{
		tox = 2.0/ax;
		m = 2*((n+(int) sqrt(ACC*n))/2);
		jsum=0;
		bjp=ans=sum=0.0;
		bj = 1.0;
		for (j=m;j>0;j--)
		{
			bjm=j*tox*bj-bjp;
			bjp=bj;
			bj=bjm;
			if (fabs(bj) > BIGNO)
			{
				bj *= BIGNI;
				bjp *= BIGNI;
				ans *= BIGNI;
				sum *= BIGNI;
			}
			if (jsum)
				sum += bj;
			jsum=!jsum;
			if (j == n) 
				ans=bjp;
		}
		sum = 2.0*sum-bj;
		ans /= sum;
	}
	return x < 0.0 && (n & 1) ? -ans : ans;
}

/* *********** */
/* *** SVD *** */
/* *********** */

void svdcmp(float **a, int m, int n, float w[], float **v)
{
	float pythag(float a, float b);
	int flag,i,its,j,jj,k,l,nm;
	float anorm,c,f,g,h,s,scale,x,y,z,*rv1;

	rv1=vector(1,n);
	g=scale=anorm=0.0;
	for (i=1;i<=n;i++)
	{
		l=i+1;
		rv1[i]=scale*g;
		g=s=scale=0.0;
		if (i <= m) 
		{
			for (k=i;k<=m;k++) scale += fabs(a[k][i]);
			if (scale)
			{
				for (k=i;k<=m;k++)
				{
					a[k][i] /= scale;
					s += a[k][i]*a[k][i];
				}
				f=a[i][i];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][i]=f-g;
				for (j=l;j<=n;j++)
				{
					for (s=0.0,k=i;k<=m;k++)
						s += a[k][i]*a[k][j];
					f=s/h;
					for (k=i;k<=m;k++)
						a[k][j] += f*a[k][i];
				}
				for (k=i;k<=m;k++) 
					a[k][i] *= scale;
			}
		}
		w[i] = scale *g;
		g=s=scale=0.0;
		if (i <= m && i != n)
		{
			for (k=l;k<=n;k++)
				scale += fabs(a[i][k]);
			if (scale) 
			{
				for (k=l;k<=n;k++)
				{
					a[i][k] /= scale;
					s += a[i][k]*a[i][k];
				}
				f=a[i][l];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][l]=f-g;
				for (k=l;k<=n;k++) 
					rv1[k]=a[i][k]/h;
				for (j=l;j<=m;j++) 
				{
					for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
					for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
				}
				for (k=l;k<=n;k++) 
					a[i][k] *= scale;
			}
		}
		anorm = FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
	}
	for (i=n;i>=1;i--)
	{
		if (i < n)
		{
			if (g)
			{
				for (j=l;j<=n;j++)
					v[j][i]=(a[i][j]/a[i][l])/g;
				for (j=l;j<=n;j++) 
				{
					for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
					for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
				}
			}
			for (j=l;j<=n;j++) 
				v[i][j] = v[j][i] = 0.0;
		}
		v[i][i] = 1.0;
		g = rv1[i];
		l = i;
	}
	for (i=IMIN(m,n);i>=1;i--) 
	{
		l=i+1;
		g=w[i];
		for (j=l;j<=n;j++) 
			a[i][j]=0.0;
		if (g)
		{
			g=1.0/g;
			for (j=l;j<=n;j++)
			{
				for (s=0.0,k=l;k<=m;k++)
					s += a[k][i]*a[k][j];
				f=(s/a[i][i])*g;
				for (k=i;k<=m;k++) 
					a[k][j] += f*a[k][i];
			}
			for (j=i;j<=m;j++) 
				a[j][i] *= g;
		} 
		else
			for (j=i;j<=m;j++) 
				a[j][i] = 0.0;
		++a[i][i];
	}
	for (k=n;k>=1;k--)
	{
		for (its=1;its<=30;its++)
		{
			flag=1;
			for (l=k;l>=1;l--)
			{
				nm=l-1;
				if ((float)(fabs(rv1[l])+anorm) == anorm)
				{
					flag=0;
					break;
				}
				if ((float)(fabs(w[nm])+anorm) == anorm) 
					break;
			}
			if (flag) 
			{
				c=0.0;
				s=1.0;
				for (i=l;i<=k;i++)
				{
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					if ((float)(fabs(f)+anorm) == anorm) 
						break;
					g = w[i];
					h = pythag(f,g);
					w[i] = h;
					h = 1.0/h;
					c = g*h;
					s = -f*h;
					for (j=1;j<=m;j++)
					{
						y=a[j][nm];
						z=a[j][i];
						a[j][nm]=y*c+z*s;
						a[j][i]=z*c-y*s;
					}
				}
			}
			z = w[k];
			if (l == k) 
			{
				if (z < 0.0) 
				{
					w[k] = -z;
					for (j=1;j<=n;j++)
						v[j][k] = -v[j][k];
				}
				break;
			}
			if (its == 30) 
				nrerror("no convergence in 30 svdcmp iterations");
			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 = pythag(f,1.0);
			f = ((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
			c=s=1.0;
			for (j=l;j<=nm;j++)
			{
				i = j+1;
				g = rv1[i];
				y = w[i];
				h = s*g;
				g = c*g;
				z = pythag(f,h);
				rv1[j] = z;
				c = f/z;
				s = h/z;
				f = x*c+g*s;
				g = g*c-x*s;
				h = y*s;
				y *= c;
				for (jj=1;jj<=n;jj++)
				{
					x = v[jj][j];
					z = v[jj][i];
					v[jj][j] = x*c+z*s;
					v[jj][i] = z*c-x*s;
				}
				z = pythag(f,h);
				w[j] = z;
				if (z) 
				{
					z = 1.0/z;
					c = f*z;
					s = h*z;
				}
				f = c*g+s*y;
				x = c*y-s*g;
				for (jj=1;jj<=m;jj++)
				{
					y = a[jj][j];
					z = a[jj][i];
					a[jj][j] = y*c+z*s;
					a[jj][i] = z*c-y*s;
				}
			}
			rv1[l] = 0.0;
			rv1[k] = f;
			w[k] = x;
		}
	}
	free_vector(rv1,1,n);
}

/* ******************************************************* */
/* *** Random generator(Uniform distribution in (0,1)) *** */
/* ******************************************************* */

/** ****** different from NR********* **/
float randu(float *r)
{ 
	int m;
    float s,u,v,p;
    s = 65536.0;
    u = 2053.0;
    v = 13849.0;
    m = (int)(*r/s); 
    *r = *r-m*s;
    *r = u*(*r)+v; 
    m = (int)(*r/s);
    *r = *r-m*s; 
    p = *r/s;
    return p;
}

/* ********************************************* */
/* *** Random generator(Normal distribution) *** */
/* ********************************************* */

/** ****** different from NR******** **/
float randn(float mean,float std,float *r)
{ 
    int i,m;
    float s,w,v,t;
    s = 65536.0; 
    w = 2053.0;
    v = 13849.0;
    t = 0.0;
    for (i=1; i<=12; i++)
	{
	*r = (*r)*w+v; m=(int)(*r/s);
        *r = *r-m*s; t=t+(*r)/s;
	}
    t = mean+std*(t-6.0);
    return t;
}