Adding Reply Archived Topic: It will not be bumped to the top of the forum. Topic Find the problem. Membername   Not a member? Sign Up Free (takes 20 seconds) Password   Forgotten your password? Post Font: Verdana Arial Black Georgia Comic Sans Impact   Size: Tiny Normal Big Huge   Color: Default Dark Red Red Orange Brown Green Olive Cyan Blue Dark Blue Indigo Black [ Check Spelling ] [ Post Preview ] [ LiveTags ] FAQ Keyword Search: Post Options Favorites Manager Notify me of new replies to this topic by email Notify me of new replies to this topic by private message Original Post Curtained Posted at 7:14 am on Mar. 22, 2007 Who can find the problem?   Code: #include <float.h> #include <math.h> #include <complex.h> const double pi4 = M_PI_4; complex<double> cdeint(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double> cdeint2(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double> cdeint3(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double>)(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<double> cdeint4(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double>)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); }

Replies Curtained Posted at 7:14 am on Mar. 22, 2007 Who can find the problem?   Code: #include <float.h> #include <math.h> #include <complex.h> const double pi4 = M_PI_4; complex<double> cdeint(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double> cdeint2(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double> cdeint3(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double>)(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<double> cdeint4(complex<double> (*f)(complex<double>),int n,    complex<double> a,complex<double> b) {    complex<double> 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<double>)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<double>)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); } All 9 previous replies displayed.