📄 numerical_recipes.h
字号:
/* ******************************************************************* */
/* Numerical recipes C header file,numerical_recipes.h */
/* */
/* This file contains a set of mathimatical functions that are commonly*/
/* used,including the cholesky decomposition for both real or complex */
/* matrix,the real matrix sigular value decomposition algorithm, */
/* the generation of a random variable which is uniformly distributed */
/* or normally distributed, the bessel functions of o,1,n-th order,the */
/* allocation and disallocation of a vector or a matrix for both real */
/* and complex type, etc. */
/* */
/* Original Source: "Numerical Recipes in C", 2nd Edition by */
/* William H.Press,Saul A. Teukolsky, et. al. */
/* Revised by : Xu Xiaojun, Zeng Meng */
/* National Communications Lab,UESTC, 25 March 2005 */
/* */
/* ******************************************************************* */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include "complex.h"
#define NR_END 1
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
#define FMAX(a,b) (maxarg1 = (a),maxarg2 = (b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
#define IMIN(a,b) (iminarg1 = (a),iminarg2 = (a),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define FREE_ARG char*
#define SWAP(a,b) tempr = (a);(a) = (b);(b) = tempr
#define ITMAX 100
#define FPMIN 1.0e-30
#define EPS 3.0e-7
#define TINY 1.0e-20
#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10
#define Coa 0.3275911
#define CoA 0.254829592
#define CoB -0.284496736
#define CoC 1.421413741
#define CoD -1.453152027
#define CoE 1.061405429
#define Pi 3.1415926
static float sqrarg;
static float maxarg1,maxarg2;
static int iminarg1,iminarg2;
void nrerror(char error_text[]);
int *ivector(long nl,long nh);
float *vector(long nl, long nh);
complex *cvector(long nl, long nh);
void free_vector_int(int *vector,long nl, long nh);
void free_vector(float *v, long nl, long nh);
void free_vector_cmpx(complex *v, long nl, long nh);
float **matrix(long nrl, long nrh, long ncl, long nch);
complex **cmatrix(long nrl, long nrh, long ncl, long nch);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
void free_matrix_cmpx(complex **m, long nrl, long nrh, long ncl, long nch);
float pythag(float a, float b);
/** ******Functions used in erfcomp()************ **/
float erff(float x);
complex erfcomp(complex z);
float gammp(float a, float x);
void gcf(float *gammcf, float a, float x, float *gln);
void gser(float *gamser, float a, float x, float *gln);
float gammln(float xx);
/** ******Bessel function*********** **/
float bessj0(float x);
float bessj1(float x);
float bessj(int n, float x);
/** ******Matrix computation******** **/
/** cholesky decomposition of real matrix **/
void choldc(float **a, int n, float p[]);
/** cholesky decomposition of complex matrix **/
void choldc_cmpx(complex **a, int n, complex p[]);
/** LU decomposition **/
void ludcmp(float **a, int n, int *indx, float *d);
void lubksb(float **a, int n, int *indx, float b[]);
void inverse(float ** a,int n, float *b);
/** SVD decomposition **/
void svdcmp(float **a, int m, int n, float w[], float **v);
/** ******Random variable production**** */
/** randu() for uniform, randn() for normal **/
float randu(float *r);
float randn(float mean,float std,float *r);
/* *********************** */
/* *** Error messaging *** */
/* *********************** */
void nrerror(char error_text[])
/** Numerical Recipes standard error handler **/
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
/* ********************************* */
/* *** Dynamic vector allocation *** */
/* ********************************* */
int *ivector(long nl, long nh)
/** allocate a int vector with subscript range v[nl..nh] **/
{
int *v;
v = (int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v)
nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
/** allocate a float vector with subscript range v[nl..nh] **/
float *vector(long nl, long nh)
{
float *v;
v = (float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v)
nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
/* **************************************** */
/* ***complex Dynamic vector allocation *** */
/* **************************************** */
/** allocate a complex vector with subscript range v[nl..nh] **/
complex *cvector(long nl, long nh)
{
complex *v;
v = (complex *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(complex)));
if (!v)
nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
/* ************************************ */
/* *** Dynamic vector desallocation *** */
/* ************************************ */
void free_vector_int(int *v, long nl, long nh)
/** free a int vector allocated with vector() **/
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_vector(float *v, long nl, long nh)
/** free a float vector allocated with vector() **/
{
free((FREE_ARG) (v+nl-NR_END));
}
/* ******************************************** */
/* *** Complex dynamic vector disallocation *** */
/* ******************************************** */
void free_vector_cmpx(complex *v, long nl, long nh)
{
free((FREE_ARG) (v+nl-NR_END));
}
/* ********************************* */
/* *** Dynamic matrix allocation *** */
/* ********************************* */
float **matrix(long nrl, long nrh, long ncl, long nch)
/** allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] **/
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/** allocate pointers to rows **/
m = (float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m)
nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/** allocate rows and set pointers to them **/
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl])
nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/** return pointer to array of pointers to rows **/
return m;
}
/* ***************************************** */
/* *** complex Dynamic matrix allocation *** */
/* ***************************************** */
complex **cmatrix(long nrl, long nrh, long ncl, long nch)
/** allocate a complex matrix with subscript range m[nrl..nrh][ncl..nch] **/
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
complex **m;
/** allocate pointers to rows **/
m = (complex **) malloc((size_t)((nrow+NR_END)*sizeof(complex *)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/** allocate rows and set pointers to them **/
m[nrl]=(complex *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(complex *)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/** return pointer to array of pointers to rows **/
return m;
}
/* ************************************ */
/* *** Dynamic matrix disallocation *** */
/* ************************************ */
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/** free a float matrix allocated by matrix() **/
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
/* ******************************************** */
/* *** Complex dynamic matrix disallocation *** */
/* ******************************************** */
void free_matrix_cmpx(complex **m, long nrl, long nrh, long ncl, long nch)
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
/* **************************************** */
/* *** Cholesky decomposition algorithm *** */
/* **************************************** */
void choldc(float **a, int n, float p[])
{
void nrerror(char error_text[]);
int i,j,k;
float sum;
for (i=1;i<=n;i++)
{
for (j=i;j<=n;j++)
{
for (sum=a[i][j],k=i-1;k>=1;k--)
sum -= a[i][k]*a[j][k];
if (i == j)
{
if (sum <= 0.0)
nrerror("choldc failed");
p[i]=sqrt(sum);
}
else
a[j][i]=sum/p[i];
}
}
}
/* ************************************************ */
/* *** complex Cholesky decomposition algorithm *** */
/* ************************************************ */
void choldc_cmpx(complex **a, int n, complex p[])
{
void nrerror(char error_text[]);
int i,j,k;
complex sum;
for (i=1;i<=n;i++)
{
for (j=i;j<=n;j++)
{
for (sum=a[i][j],k=i-1;k>=1;k--)
sum -= conj(a[i][k])*a[j][k];
if (i == j)
{
if (norm(sum) <= 0.0)
nrerror("choldc failed");
p[i]=sqrt(sum);
}
else
a[j][i]=sum/p[i];
}
}
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
a[i][j] = conj(a[i][j]);
}
/* ****************************** */
/* *** Computes sqrt(a^2+b^2) *** */
/* ****************************** */
float pythag(float a, float b)
{
float absa,absb;
absa=fabs(a);
absb=fabs(b);
if (absa > absb)
return absa*sqrt(1.0+SQR(absb/absa));
else
return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
}
/* **************************** */
/* *** Error Function erf() *** */
/* **************************** */
/* ******gammln**************** */
float gammln(float xx)
{
float x,y,tmp,ser;
static float cof[6]={76.18009172947146,-86.50532032941677,
24.01409824083091,-1.231739572450155,
0.1208650973866179e-2,-0.5395239384953e-5};
int j;
y=x=xx;
tmp = x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.000000000190015;
for (j=0;j<=5;j++) ser += cof[j]/++y;
return -tmp+log(2.5066282746310005*ser/x);
}
/** ******gser function********* **/
void gser(float *gamser, float a, float x, float *gln)
{
float gammln(float xx);
void nrerror(char error_text[]);
int n;
float sum,del,ap;
*gln = gammln(a);
if (x <= 0.0)
{
if (x < 0.0)
nrerror("x less than 0 in routine gser");
*gamser = 0.0;
return;
}
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) {
*gamser=sum*exp(-x+a*log(x)-(*gln));
return;
}
}
nrerror("a too large, ITMAX too small in routine gser");
return;
}
}
/** ******gcf funtion************ **/
void gcf(float *gammcf, float a, float x, float *gln)
{
float gammln(float xx);
void nrerror(char error_text[]);
int i;
float an,b,c,d,del,h;
*gln = gammln(a);
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)
nrerror("a too large, ITMAX too small in gcf");
*gammcf=exp(-x+a*log(x)-(*gln))*h;
}
/** ******gamma funtion********* **/
float gammp(float a, float x)
{
void gcf(float *gammcf, float a, float x, float *gln);
void gser(float *gamser, float a, float x, float *gln);
void nrerror(char error_text[]);
float gamser,gammcf,gln;
if (x < 0.0 || a <= 0.0)
nrerror("Invalid arguments in routine gammp");
if (x < (a+1.0))
{
gser(&gamser,a,x,&gln);
return gamser;
}
else
{
gcf(&gammcf,a,x,&gln);
return 1.0-gammcf;
}
}
/** ******erf function main part******** **/
float erff(float x)
{
float gammp(float a, float x);
return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}
/** ******complex error function******** **/
complex erfcomp(complex z)
{
float x,y,w,p,q,abserr,fn,gn;
int n = 1;
complex E=complex(0.0,0.0);
complex S=complex(0.0,0.0);
complex Sakt=complex(0.0,0.0);
x = z.getReal();
y = z.getIm();
w = 1.0/(1.0+Coa*fabs(x));
p = CoE*w+CoD;
p = p*w+CoC;
p = p*w+CoB;
p = p*w+CoA;
p = p*w;
q = exp(-x*x);
E = (1-q*p)*(x/fabs(x));
abserr = 100.0;
while(abserr>1e-12)
{
fn = 2*x-2*x*cosh(n*y)*cos(2*x*y)+n*sinh(n*y)*sin(2*x*y);
gn = 2*x*cosh(n*y)*sin(2*x*y)+n*sinh(n*y)*cos(2*x*y);
Sakt.setValue((2/Pi)*exp(-x*x)*exp(-0.25*n*n)/(n*n+4*x*x)*fn,(2/Pi)
*exp(-x*x)*exp(-0.25*n*n)/(n*n+4*x*x)*gn);
abserr = sqrt(norm(Sakt));
n++;
S += Sakt;
}
if(x == 0.0&&y == 0.0)
E.setValue(0.0,0.0);
if(x == 0.0&&y != 0.0)
E.setValue(0.0,y/Pi);
if(x != 0.0&&y != 0.0)
E.setValue(erff(x)+exp(-x*x)*(1-cos(2*x*y))/(2*Pi*x),exp(-x*x)*sin(2*x*y)/(2*Pi*x));
E += S;
return E;
}
/* **************************** */
/* *** sovle AX=b********** *** */
/* **************************** */
/** ******ludcmp************** **/
void ludcmp(float **a, int n, int *indx, float *d)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -