//  cpsi.cpp -- Complex psi function.
//      Algorithms and coefficient values from "Computation of Special
//      Functions", Zhang and Jin, John Wiley and Sons, 1996.
//
//  (C) 2003, C. Bond. All rights reserved.
//
#include <iostream.h>
#include <complex.h>

complex<double> cpsi(complex<double> z)
{
    double x,y,x0,x1,y1,psr,psi;
    double q0,q2,rr,ri,th,tn,tm,ct2;
    complex<double>ps;

    int n,k;
    static double a[] = {
        -0.8333333333333e-01,
         0.83333333333333333e-02,
        -0.39682539682539683e-02,
         0.41666666666666667e-02,
        -0.75757575757575758e-02,
         0.21092796092796093e-01,
        -0.83333333333333333e-01,
         0.4432598039215686};

    x = real(z);
    y = imag(z);
    x1 = 0.0;
    n = 0;
    if ((y == 0.0) && (x == (int)x) && (x <= 0.0)) {
        ps = complex<double> (1e308,0);
    }
    else {
        if (x < 0.0) {
            x1 = x;
            y1 = y;
            x = -x;
            y = -y;
        }
        x0 = x;
        if (x < 8.0) {
            n = 8 - (int)x;
            x0 = x + n;
        }
        if ((x0 == 0.0) && (y != 0.0))
            th = 0.5*M_PI;
        if (x0 != 0.0)
            th = atan(y/x0);
        q2 = x0*x0+y*y;
        q0 = sqrt(q2);
        psr = log(q0)-0.5*x0/q2;
        psi = th+0.5*y/q2;
        for (k=1;k<=8;k++) {
            psr += (a[k-1]*pow(q2,-k)*cos(2.0*k*th));
            psi -= (a[k-1]*pow(q2,-k)*sin(2.0*k*th));
        }
        if (x < 8.0) {
            rr = 0.0;
            ri = 0.0;
            for (k=1;k<=n;k++) {
                rr += ((x0-k)/(pow(x0-k,2.0)+y*y));
                ri += (y/(pow(x0-k,2.0)+y*y));
            }
            psr -= rr;
            psi += ri;
        }
        if (x1 < 0.0) {
            tn = tan(M_PI*x);
            tm = tanh(M_PI*y);
            ct2 = tn*tn+tm*tm;
            psr = psr+x/(x*x+y*y)+M_PI*(tn-tn*tm*tm)/ct2;
            psi = psi-y/(x*x+y*y)-M_PI*tm*(1.0+tn*tn)/ct2;
            x = x1;
            y = y1;
        }
        ps = complex<double>(psr,psi);
    }
    return ps;
}