// cplxde.cpp -- complex integration of function by double exponential
//  transform. Based on algorithm in Applied Numerical Methods with
//  Software, by Nakamura.
//
//  (C) 2001, C. Bond. All rights reserved.
//
//  Note: in these routines the value of the product 'n*h' is limited
//  by the allowed range of floating point numbers before overflow
//  occurs. The maximum floating point value for IEEE double
//  precision arithmetic is:   1.79769E+308.
//

#include <float.h>
#include <math.h>
#include <complex.h>

const double pi4 = M_PI_4;

// This version closely follows the method in Nakamura's book, except
// that it is ported to C, converted to complex<double> and has minor
// improvements in efficiency.
//
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);
}