C++ Vec_IO_DP类代码示例

void NR::rkqs(Vec_IO_DP &y, Vec_IO_DP &dydx, DP &x, const DP htry,
	const DP eps, Vec_I_DP &yscal, DP &hdid, DP &hnext,
	void derivs(const DP, Vec_I_DP &, Vec_O_DP &))
{
	const DP SAFETY=0.9, PGROW=-0.2, PSHRNK=-0.25, ERRCON=1.89e-4;
	int i;
	DP errmax,h,htemp,xnew;

	int n=y.size();
	h=htry;
	Vec_DP yerr(n),ytemp(n);
	for (;;) {
		rkck(y,dydx,x,h,ytemp,yerr,derivs);
		errmax=0.0;
		for (i=0;i<n;i++) errmax=MAX(errmax,fabs(yerr[i]/yscal[i]));
		errmax /= eps;
		if (errmax <= 1.0) break;
		htemp=SAFETY*h*pow(errmax,PSHRNK);
		h=(h >= 0.0 ? MAX(htemp,0.1*h) : MIN(htemp,0.1*h));
		xnew=x+h;
		if (xnew == x) nrerror("stepsize underflow in rkqs");
	}
	if (errmax > ERRCON) hnext=SAFETY*h*pow(errmax,PGROW);
	else hnext=5.0*h;
	x += (hdid=h);
	for (i=0;i<n;i++) y[i]=ytemp[i];
}

void NR::eulsum(DP &sum, const DP term, const int jterm, Vec_IO_DP &wksp)
{
	int j;
	static int nterm;
	DP tmp,dum;

	if (jterm == 0) {
		nterm=1;
		sum=0.5*(wksp[0]=term);
	} else {
		if (nterm+1 > wksp.size()) nrerror("wksp too small in euler");
		tmp=wksp[0];
		wksp[0]=term;
		for (j=1;j<nterm;j++) {
			dum=wksp[j];
			wksp[j]=0.5*(wksp[j-1]+tmp);
			tmp=dum;
		}
		wksp[nterm]=0.5*(wksp[nterm-1]+tmp);
		if (fabs(wksp[nterm]) <= fabs(wksp[nterm-1]))
			sum += (0.5*wksp[nterm++]);
		else
			sum += wksp[nterm];
	}
}

void NR::cosft2(Vec_IO_DP &y, const int isign)
{
	const DP PI=3.141592653589793238;
	int i;
	DP sum,sum1,y1,y2,ytemp,theta,wi=0.0,wi1,wpi,wpr,wr=1.0,wr1,wtemp;

	int n=y.size();
	theta=0.5*PI/n;
	wr1=cos(theta);
	wi1=sin(theta);
	wpr = -2.0*wi1*wi1;
	wpi=sin(2.0*theta);
	if (isign == 1) {
		for (i=0;i<n/2;i++) {
			y1=0.5*(y[i]+y[n-1-i]);
			y2=wi1*(y[i]-y[n-1-i]);
			y[i]=y1+y2;
			y[n-1-i]=y1-y2;
			wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1;
			wi1=wi1*wpr+wtemp*wpi+wi1;
		}
		realft(y,1);
		for (i=2;i<n;i+=2) {
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
			y1=y[i]*wr-y[i+1]*wi;
			y2=y[i+1]*wr+y[i]*wi;
			y[i]=y1;
			y[i+1]=y2;
		}
		sum=0.5*y[1];
		for (i=n-1;i>0;i-=2) {
			sum1=sum;
			sum += y[i];
			y[i]=sum1;
		}
	} else if (isign == -1) {
		ytemp=y[n-1];
		for (i=n-1;i>2;i-=2)
			y[i]=y[i-2]-y[i];
		y[1]=2.0*ytemp;
		for (i=2;i<n;i+=2) {
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
			y1=y[i]*wr+y[i+1]*wi;
			y2=y[i+1]*wr-y[i]*wi;
			y[i]=y1;
			y[i+1]=y2;
		}
		realft(y,-1);
		for (i=0;i<n/2;i++) {
			y1=y[i]+y[n-1-i];
			y2=(0.5/wi1)*(y[i]-y[n-1-i]);
			y[i]=0.5*(y1+y2);
			y[n-1-i]=0.5*(y1-y2);
			wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1;
			wi1=wi1*wpr+wtemp*wpi+wi1;
		}
	}
}

void NR::fixrts(Vec_IO_DP &d)
{
	bool polish=true;
	int i,j;

	int m=d.size();
	Vec_CPLX_DP a(m+1),roots(m);
	a[m]=1.0;
	for (j=0;j<m;j++)
		a[j]= -d[m-1-j];
	zroots(a,roots,polish);
	for (j=0;j<m;j++)
		if (abs(roots[j]) > 1.0)
			roots[j]=1.0/conj(roots[j]);
	a[0]= -roots[0];
	a[1]=1.0;
	for (j=1;j<m;j++) {
		a[j+1]=1.0;
		for (i=j;i>=1;i--)
			a[i]=a[i-1]-roots[j]*a[i];
		a[0]= -roots[j]*a[0];
	}
	for (j=0;j<m;j++)
		d[m-1-j] = -real(a[j]);
}

void NR::kstwo(Vec_IO_DP &data1, Vec_IO_DP &data2, DP &d, DP &prob)
{
	int j1=0,j2=0;
	DP d1,d2,dt,en1,en2,en,fn1=0.0,fn2=0.0;

	int n1=data1.size();
	int n2=data2.size();
	sort(data1);
	sort(data2);
	en1=n1;
	en2=n2;
	d=0.0;
	while (j1 < n1 && j2 < n2) {
		if ((d1=data1[j1]) <= (d2=data2[j2])) fn1=j1++/en1;
		if (d2 <= d1) fn2=j2++/en2;
		if ((dt=fabs(fn2-fn1)) > d) d=dt;
	}
	en=sqrt(en1*en2/(en1+en2));
	prob=probks((en+0.12+0.11/en)*d);
}

void NR::wt1(Vec_IO_DP &a, const int isign,
	void wtstep(Vec_IO_DP &, const int, const int))
{
	int nn;

	int n=a.size();
	if (n < 4) return;
	if (isign >= 0) {
		for (nn=n;nn>=4;nn>>=1) wtstep(a,nn,isign);
	} else {
		for (nn=4;nn<=n;nn<<=1) wtstep(a,nn,isign);

void NR::sort3(Vec_IO_DP &ra, Vec_IO_DP &rb, Vec_IO_DP &rc)
{
	int j;

	int n=ra.size();
	Vec_INT iwksp(n);
	Vec_DP wksp(n);
	indexx(ra,iwksp);
	for (j=0;j<n;j++) wksp[j]=ra[j];
	for (j=0;j<n;j++) ra[j]=wksp[iwksp[j]];
	for (j=0;j<n;j++) wksp[j]=rb[j];
	for (j=0;j<n;j++) rb[j]=wksp[iwksp[j]];
	for (j=0;j<n;j++) wksp[j]=rc[j];
	for (j=0;j<n;j++) rc[j]=wksp[iwksp[j]];
}

void NR::piksrt(Vec_IO_DP &arr)
{
	int i,j;
	DP a;

	int n=arr.size();
	for (j=1;j<n;j++) {
		a=arr[j];
		i=j;
		while (i > 0 && arr[i-1] > a) {
			arr[i]=arr[i-1];
			i--;
		}
		arr[i]=a;
	}
}

void NR::pcshft(const DP a, const DP b, Vec_IO_DP &d)
{
	int k,j;
	DP fac,cnst;

	int n=d.size();
	cnst=2.0/(b-a);
	fac=cnst;
	for (j=1;j<n;j++) {
		d[j] *= fac;
		fac *= cnst;
	}
	cnst=0.5*(a+b);
	for (j=0;j<=n-2;j++)
		for (k=n-2;k>=j;k--)
			d[k] -= cnst*d[k+1];
}

void NR::four1(Vec_IO_DP &data, const int isign)
{
	int n,mmax,m,j,istep,i;
	DP wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;

	int nn=data.size()/2;
	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j-1],data[i-1]);
			SWAP(data[j],data[i]);
		}
		m=nn;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j-1]-wi*data[j];
				tempi=wr*data[j]+wi*data[j-1];
				data[j-1]=data[i-1]-tempr;
				data[j]=data[i]-tempi;
				data[i-1] += tempr;
				data[i] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}

void NR::ksone(Vec_IO_DP &data, DP func(const DP), DP &d, DP &prob)
{
	int j;
	DP dt,en,ff,fn,fo=0.0;

	int n=data.size();
	sort(data);
	en=n;
	d=0.0;
	for (j=0;j<n;j++) {
		fn=(j+1)/en;
		ff=func(data[j]);
		dt=MAX(fabs(fo-ff),fabs(fn-ff));
		if (dt > d) d=dt;
		fo=fn;
	}
	en=sqrt(en);
	prob=probks((en+0.12+0.11/en)*d);
}

void NR::frprmn(Vec_IO_DP &p, const DP ftol, int &iter, DP &fret,
	DP func(Vec_I_DP &), void dfunc(Vec_I_DP &, Vec_O_DP &))
{
	const int ITMAX=200;
	const DP EPS=1.0e-18;
	int j,its;
	DP gg,gam,fp,dgg;

	int n=p.size();
	Vec_DP g(n),h(n),xi(n);
	fp=func(p);
	dfunc(p,xi);
	for (j=0;j<n;j++) {
		g[j] = -xi[j];
		xi[j]=h[j]=g[j];
	}
	for (its=0;its<ITMAX;its++) {
		iter=its;
		linmin(p,xi,fret,func);
		if (2.0*fabs(fret-fp) <= ftol*(fabs(fret)+fabs(fp)+EPS))
			return;
		fp=fret;
		dfunc(p,xi);
		dgg=gg=0.0;
		for (j=0;j<n;j++) {
			gg += g[j]*g[j];
//		  dgg += xi[j]*xi[j];
			dgg += (xi[j]+g[j])*xi[j];
		}
		if (gg == 0.0)
			return;
		gam=dgg/gg;
		for (j=0;j<n;j++) {
			g[j] = -xi[j];
			xi[j]=h[j]=g[j]+gam*h[j];
		}
	}
	nrerror("Too many iterations in frprmn");
}

void NR::pade(Vec_IO_DP &cof, DP &resid)
{
	const DP BIG=1.0e30;
	int j,k;
	DP d,rr,rrold,sum;

	int n=(cof.size()-1)/2;
	Mat_DP q(n,n),qlu(n,n);
	Vec_INT indx(n);
	Vec_DP x(n),y(n),z(n);
	for (j=0;j<n;j++) {
		y[j]=x[j]=cof[n+j+1];
		for (k=0;k<n;k++) {
			q[j][k]=cof[j-k+n];
			qlu[j][k]=q[j][k];
		}
	}
	ludcmp(qlu,indx,d);
	lubksb(qlu,indx,x);
	rr=BIG;
	do {
		rrold=rr;
		for (j=0;j<n;j++) z[j]=x[j];
		mprove(q,qlu,indx,y,x);
		for (rr=0.0,j=0;j<n;j++)
			rr += SQR(z[j]-x[j]);
	} while (rr < rrold);
	resid=sqrt(rrold);
	for (k=0;k<n;k++) {
		for (sum=cof[k+1],j=0;j<=k;j++)
			sum -= z[j]*cof[k-j];
		y[k]=sum;
	}
	for (j=0;j<n;j++) {
		cof[j+1]=y[j];
		cof[j+1+n] = -z[j];
	}
}

DP NR::select(const int k, Vec_IO_DP &arr)
{
	int i,ir,j,l,mid;
	DP a;

	int n=arr.size();
	l=0;
	ir=n-1;
	for (;;) {
		if (ir <= l+1) {
			if (ir == l+1 && arr[ir] < arr[l])
				SWAP(arr[l],arr[ir]);
			return arr[k];
		} else {
			mid=(l+ir) >> 1;
			SWAP(arr[mid],arr[l+1]);
			if (arr[l] > arr[ir])
				SWAP(arr[l],arr[ir]);
			if (arr[l+1] > arr[ir])
				SWAP(arr[l+1],arr[ir]);
			if (arr[l] > arr[l+1])
				SWAP(arr[l],arr[l+1]);
			i=l+1;
			j=ir;
			a=arr[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP(arr[i],arr[j]);
			}
			arr[l+1]=arr[j];
			arr[j]=a;
			if (j >= k) ir=j-1;
			if (j <= k) l=i;
		}
	}
}

void NR::crank(Vec_IO_DP &w, DP &s)
{
	int j=1,ji,jt;
	DP t,rank;

	int n=w.size();
	s=0.0;
	while (j < n) {
		if (w[j] != w[j-1]) {
			w[j-1]=j;
			++j;
		} else {
			for (jt=j+1;jt<=n && w[jt-1]==w[j-1];jt++);
			rank=0.5*(j+jt-1);
			for (ji=j;ji<=(jt-1);ji++)
				w[ji-1]=rank;
			t=jt-j;
			s += (t*t*t-t);
			j=jt;
		}
	}
	if (j == n) w[n-1]=n;
}