#include <iostream>
#include <iomanip>
#include <cmath>
#include "nr.h"
using namespace std;
extern int idum;
void NR::vegas(Vec_I_DP ®n, DP fxn(Vec_I_DP &, const DP), const int init,
const int ncall, const int itmx, const int nprn, DP &tgral, DP &sd,
DP &chi2a)
{
const int NDMX=50, MXDIM=10;
const DP ALPH=1.5, TINY=1.0e-30;
static int i,it,j,k,mds,nd,ndo,ng,npg;
static DP calls,dv2g,dxg,f,f2,f2b,fb,rc,ti;
static DP tsi,wgt,xjac,xn,xnd,xo,schi,si,swgt;
static Vec_INT ia(MXDIM),kg(MXDIM);
static Vec_DP dt(MXDIM),dx(MXDIM),r(NDMX),x(MXDIM),xin(NDMX);
static Mat_DP d(NDMX,MXDIM),di(NDMX,MXDIM),xi(MXDIM,NDMX);
int ndim=regn.size()/2;
if (init <= 0) {
mds=ndo=1;
for (j=0;j<ndim;j++) xi[j][0]=1.0;
}
if (init <= 1) si=swgt=schi=0.0;
if (init <= 2) {
nd=NDMX;
ng=1;
if (mds != 0) {
ng=int(pow(ncall/2.0+0.25,1.0/ndim));
mds=1;
if ((2*ng-NDMX) >= 0) {
mds = -1;
npg=ng/NDMX+1;
nd=ng/npg;
ng=npg*nd;
}
}
for (k=1,i=0;i<ndim;i++) k *= ng;
npg=MAX(int(ncall/k),2);
calls=DP(npg)*DP(k);
dxg=1.0/ng;
for (dv2g=1,i=0;i<ndim;i++) dv2g *= dxg;
dv2g=SQR(calls*dv2g)/npg/npg/(npg-1.0);
xnd=nd;
dxg *= xnd;
xjac=1.0/calls;
for (j=0;j<ndim;j++) {
dx[j]=regn[j+ndim]-regn[j];
xjac *= dx[j];
}
if (nd != ndo) {
for (i=0;i<MAX(nd,ndo);i++) r[i]=1.0;
for (j=0;j<ndim;j++)
rebin(ndo/xnd,nd,r,xin,xi,j);
ndo=nd;
}
if (nprn >= 0) {
cout << " Input parameters for vegas";
cout << " ndim= " << setw(4) << ndim;
cout << " ncall= " << setw(8) << calls << endl;
cout << setw(34) << " it=" << setw(5) << it;
cout << " itmx=" << setw(5) << itmx << endl;
cout << setw(34) << " nprn=" << setw(5) << nprn;
cout << " ALPH=" << setw(9) << ALPH << endl;
cout << setw(34) << " mds=" << setw(5) << mds;
cout << " nd=" << setw(5) << nd << endl;
for (j=0;j<ndim;j++) {
cout << setw(30) << " x1[" << setw(2) << j;
cout << "]= " << setw(11) << regn[j] << " xu[";
cout << setw(2) << j << "]= ";
cout << setw(11) << regn[j+ndim] << endl;
}
}
}
for (it=0;it<itmx;it++) {
ti=tsi=0.0;
for (j=0;j<ndim;j++) {
kg[j]=1;
for (i=0;i<nd;i++) d[i][j]=di[i][j]=0.0;
}
for (;;) {
fb=f2b=0.0;
for (k=0;k<npg;k++) {
wgt=xjac;
for (j=0;j<ndim;j++) {
xn=(kg[j]-ran2(idum))*dxg+1.0;
ia[j]=MAX(MIN(int(xn),NDMX),1);
if (ia[j] > 1) {
xo=xi[j][ia[j]-1]-xi[j][ia[j]-2];
rc=xi[j][ia[j]-2]+(xn-ia[j])*xo;
} else {
xo=xi[j][ia[j]-1];
rc=(xn-ia[j])*xo;
}
x[j]=regn[j]+rc*dx[j];
wgt *= xo*xnd;
}
f=wgt*fxn(x,wgt);
f2=f*f;
fb += f;
f2b += f2;
for (j=0;j<ndim;j++) {
di[ia[j]-1][j] += f;
if (mds >= 0) d[ia[j]-1][j] += f2;
}
}
f2b=sqrt(f2b*npg);
f2b=(f2b-fb)*(f2b+fb);
if (f2b <= 0.0) f2b=TINY;
ti += fb;
tsi += f2b;
if (mds < 0) {
for (j=0;j<ndim;j++) d[ia[j]-1][j] += f2b;
}
for (k=ndim-1;k>=0;k--) {
kg[k] %= ng;
if (++kg[k] != 1) break;
}
if (k < 0) break;
}
tsi *= dv2g;
wgt=1.0/tsi;
si += wgt*ti;
schi += wgt*ti*ti;
swgt += wgt;
tgral=si/swgt;
chi2a=(schi-si*tgral)/(it+0.0001);
if (chi2a < 0.0) chi2a = 0.0;
sd=sqrt(1.0/swgt);
tsi=sqrt(tsi);
if (nprn >= 0) {
cout << " iteration no. " << setw(3) << (it+1);
cout << " : integral = " << setw(14) << ti;
cout << " +/- " << setw(9) << tsi << endl;
cout << " all iterations: " << " integral =";
cout << setw(14) << tgral << "+-" << setw(9) << sd;
cout << " chi**2/IT n =" << setw(9) << chi2a << endl;
if (nprn != 0) {
for (j=0;j<ndim;j++) {
cout << " DATA FOR axis " << setw(2) << j << endl;
cout << " X delta i X delta i";
cout << " X deltai" << endl;
for (i=nprn/2;i<nd;i += nprn+2) {
cout << setw(8) << xi[j][i] << setw(12) << di[i][j];
cout << setw(12) << xi[j][i+1] << setw(12) << di[i+1][j];
cout << setw(12) << xi[j][i+2] << setw(12) << di[i+2][j];
cout << endl;
}
}
}
}
for (j=0;j<ndim;j++) {
xo=d[0][j];
xn=d[1][j];
d[0][j]=(xo+xn)/2.0;
dt[j]=d[0][j];
for (i=2;i<nd;i++) {
rc=xo+xn;
xo=xn;
xn=d[i][j];
d[i-1][j] = (rc+xn)/3.0;
dt[j] += d[i-1][j];
}
d[nd-1][j]=(xo+xn)/2.0;
dt[j] += d[nd-1][j];
}
for (j=0;j<ndim;j++) {
rc=0.0;
for (i=0;i<nd;i++) {
if (d[i][j] < TINY) d[i][j]=TINY;
r[i]=pow((1.0-d[i][j]/dt[j])/
(log(dt[j])-log(d[i][j])),ALPH);
rc += r[i];
}
rebin(rc/xnd,nd,r,xin,xi,j);
}
}
}
C++数值算法（第二版）的例程
需积分: 10 47 浏览量
2010-02-10
20:23:13
上传
评论 1
收藏 377KB RAR 举报
C++数值算法（第二版）的例程 （627个子文件） 
vegas.cpp 5KB xfourfs.cpp 5KB xqrupdt.cpp 5KB fourfs.cpp 4KB svdcmp.cpp 4KB bessjy.cpp 3KB xqrdcmp.cpp 3KB broydn.cpp 3KB xfour1.cpp 3KB xludcmp.cpp 3KB stifbs.cpp 3KB xjacobi.cpp 3KB hqr.cpp 3KB bsstep.cpp 3KB xmrqmin.cpp 3KB xhuffman.cpp 3KB indexx.cpp 3KB bessik.cpp 3KB xamebsa.cpp 3KB xlfit.cpp 2KB xarcode.cpp 2KB xsvbksb.cpp 2KB xsvdcmp.cpp 2KB machar.cpp 2KB xgaussj.cpp 2KB mpops.cpp 2KB stiff.cpp 2KB miser.cpp 2KB xspear.cpp 2KB xqromo.cpp 2KB xmachar.cpp 2KB xfitexy.cpp 2KB xkendl2.cpp 2KB xlinbcg.cpp 2KB xqrsolv.cpp 2KB xrlft3.cpp 2KB xran.cpp 2KB dbrent.cpp 2KB xpccheb.cpp 2KB xcntab2.cpp 2KB xtridag.cpp 2KB xmrqcof.cpp 2KB xsprstm.cpp 2KB dftcor.cpp 2KB xlubksb.cpp 2KB mgfas.cpp 2KB xmnewt.cpp 2KB xtqli.cpp 2KB arcode.cpp 2KB dfpmin.cpp 2KB xcrank.cpp 2KB rlft3.cpp 2KB linbcg.cpp 2KB xsprspm.cpp 2KB xpolint.cpp 2KB gaujac.cpp 2KB jacobi.cpp 2KB simplx.cpp 2KB xfourn.cpp 2KB fitexy.cpp 2KB solvde.cpp 2KB xtred2.cpp 2KB xeigsrt.cpp 2KB xsplint.cpp 2KB difeq.cpp 2KB xbsstep.cpp 2KB xodeint.cpp 2KB xrealft.cpp 2KB xpolcof.cpp 2KB xvegas.cpp 2KB xtwofft.cpp 2KB newt.cpp 2KB xpolcoe.cpp 2KB xflmoon.cpp 2KB xcntab1.cpp 2KB amebsa.cpp 2KB xgaujac.cpp 2KB period.cpp 2KB xcosft2.cpp 2KB xorthog.cpp 2KB xcosft1.cpp 2KB xfixrts.cpp 2KB sfroid.cpp 2KB xdecchk.cpp 2KB fasper.cpp 2KB rj.cpp 2KB anneal.cpp 2KB xtptest.cpp 2KB xsinft.cpp 2KB lnsrch.cpp 2KB xsimplx.cpp 2KB xsavgol.cpp 2KB xfrprmn.cpp 2KB xdfpmin.cpp 2KB xzroots.cpp 2KB xdbrent.cpp 2KB xmoment.cpp 2KB xmmid.cpp 2KB brent.cpp 2KB xbessjy.cpp 2KB共 627 条
- 1
- 2
- 3
- 4
- 5
- 6
- 7
资源评论
xxxyyy333
- 粉丝: 27
- 资源: 1
最新资源
- AIS2024 valid
- 最入门的爬虫代码 python.docx
- 爬虫零基础入门-爬取天气预报.pdf
- 最通俗易懂的 MongoDB 非结构化文档存储数据库教程.zip
- 以mongodb为数据库的订单物流小项目.zip
- 腾讯云-mongodb数据库， 项目部署.zip
- 腾讯 APIJSON 的 MongoDB 数据库插件.zip
- 理解非关系型数据库和关系型数据库的区别.zip
- 操作简单的Mongodb网页web管理工具,基于Spring Boot2.0支持mongodb集群.zip
- tms-mongodb-web,提供访问mongodb数据的REST API和可灵活扩展的mongodb web 客户端.zip
资源上传下载、课程学习等过程中有任何疑问或建议，欢迎提出宝贵意见哦~我们会及时处理！
点击此处反馈
