|
-- Posted by Curtained at 7:14 am on Mar. 22, 2007
Who can find the problem?  | Code: | #include
#include
#include const double pi4 = M_PI_4;
complex cdeint(complex (*f)(complex),int n,
complex a,complex b)
{
complex h,ss,z,exz,hcos,hsin,s,dxdz,x,w; int k; n = n/2;
h = 5.0/n; // 5.0-6.0 is rough limit of K in exp(exp(K)) ss = 0.0;
for (k=-n;k<=n;k++) {
z = h * (complex)k;
exz = exp(z);
hcos = exz+1.0/exz; hsin = exz-1.0/exz;
s = exp(pi4*hsin);
w = s + 1.0/s;
x = (b*s+a/s)/w; // transformed abscissa
if (x != a && x != b) {
dxdz = 2.0*(b-a)*pi4*hcos/(w*w); // transformed weight ss += h * f(x)*dxdz;
}
}
return ss;
} // Similar to above, except that the main loop has half the number
// of iterations. Positive and negative loop variables are folded
// together so that common quantities only need to be computed
// once.
//
complex cdeint2(complex (*f)(complex),int n,
complex a,complex b)
{
complex h,ss,z,exz,hcos,hsin,s,dxdz,x1,x2,w; int k; n = n/2;
h = 5.0/n; // 5.0-6.0 is rough limit of K in exp(exp(K)) ss = 0.5*f((a+b)/2.0);
for (k=-n;k<0;k++) {
z = h*(complex)k;
exz = exp(z); hcos = exz+1.0/exz;
hsin = exz-1.0/exz;
s = exp(pi4*hsin);
w = s + 1.0/s;
dxdz = hcos/(w*w); // transformed weight
x1 = (b*s+a/s)/w; // transformed abscissa x2 = (a*s+b/s)/w; // " "
if (x1 != a && x1 != b)
ss += dxdz * f(x1);
if (x2 != a && x2 != b)
ss += dxdz * f(x2);
}
return 2.0*(b-a)*pi4*h*ss; } // This version adds a set of in-between points to the integration
complex cdeint3(complex (*f)(complex),int n,
complex a,complex b) {
complex h,ss,z,exz,hcos,hsin,s,dxdz,x,w;
int k; n = n/2;
h = 5.0/n; // 5.0-6.0 is rough limit of K in exp(exp(K)) ss = 0.0;
for (k=-n;k<=n;k++) {
z = h * (complex)k; exz = exp(z);
hcos = exz+1.0/exz;
hsin = exz-1.0/exz;
s = exp(pi4*hsin);
w = s + 1.0/s;
x = (b*s+a/s)/w; // transformed abscissa
if (x != a && x != b) { dxdz = hcos/(w*w); // transformed weight
ss += f(x) * dxdz;
}
} for (k=-n;k<=n;k++) {
z = h * ((complex)(k)+0.5); // use in-between points this pass exz = exp(z);
hcos = exz+1.0/exz;
hsin = exz-1.0/exz;
s = exp(pi4*hsin);
w = s + 1.0/s;
x = (b*s+a/s)/w; // transformed abscissa
if (x != a && x != b) {
dxdz = hcos/(w*w); // transformed weight ss += f(x) * dxdz;
}
}
return h*(b-a)*pi4*ss;
}
// Progressive version
complex cdeint4(complex (*f)(complex),int n,
complex a,complex b) {
complex h,ss,ss1,ss2,z,exz,hcos,hsin,s,dxdz,x1,x2,w;
int k,l; n /= 2;
h = 5.0/n; // 5.0-6.0 is rough limit of K in exp(exp(K))
ss = 0.5*f(0.5*(a+b));
ss1 = 0.0
;
ss2 = 0.0;
for (k=-n;k<0;k++) {
z = h*(complex)k;
exz = exp(z);
hcos = exz+1.0/exz;
hsin = exz-1.0/exz;
s = exp(pi4*hsin);
w = s + 1.0/s;
dxdz = hcos/(w*w); // transformed weight x1 = (b*s+a/s)/w; // transformed abscissa
x2 = (a*s+b/s)/w; // " "
if (x1 != a && x1 != b)
ss1 += dxdz * f(x1);
if (x2 != a && x2 != b) ss2 += dxdz * f(x2);
} for (l=0;l<3;l++) { for (k=-n;k<0;k++) {
z = h*((complex)k+0.5);
exz = exp(z);
hcos = exz+1.0/exz;
hsin = exz-1.0/exz; s = exp(pi4*hsin);
w = s + 1.0/s;
dxdz = hcos/(w*w); // transformed weight
x1 = (b*s+a/s)/w; // transformed abscissa
x2 = (a*s+b/s)/w; // " " if (x1 != a && x1 != b)
ss1 += dxdz * f(x1);
if (x2 != a && x2 != b)
ss2 += dxdz * f(x2);
}
n *= 2;
h /= 2;
}
return 2.0*(b-a)*pi4*h*(ss+ss1+ss2);
} |
|