#include <complex.h>
#include <iostream.h>

complex<double> HKQ(complex<double>z)
{
    complex<double>Hxy,zp,w,q,a,e;
    double eps,PI,u,v,x,y,tanp,tanm;
    double xph,xmh,alp,alm;
    int i;

    eps = 1e-8;
    PI = M_PI;
    x = real(z);
    y = imag(z);
    v = x+x;
    xph = x + 0.5;
    xmh = x - 0.5;
    alp = log(y*y + xph*xph);
    alm = log(y*y + xmh*xmh);
    u = (xmh*alm - xph*alp)/PI;
    if (fabs(y) > 1e-8) {
        tanp = atan(xph/y);
        tanm = atan(xmh/y);
        v = 2.0*(xph*tanp-xmh*tanm-0.5*y*(alp-alm))/PI;
        u -= (2.0*y*(tanp-tanm)/PI);
    }
    else if (fabs(x) > 0.5)
        v = 1.0;
    zp = PI*complex<double>(0.0,1.0);
// Newton's method starts here
    w = complex<double>(u,v);
    for (i=0;i<100;i++) {
        e = exp(-PI*w);
        q = sqrt(1.0+e);
        a = 2.0*(q-0.5*log((1.0+q)/(1.0-q))+(z-0.5)*zp)/(q*PI);
        w += a;
        if (abs(a) < eps) break;
    }
    if (i > 99) {
        cout << "Convergence failure!" << endl;
        return complex<double>(0.0,0.0);
    }
    cout << w << endl;
    a = complex<double>(0.0,1.0) / q;
    Hxy = complex<double>(imag(a),real(a));

    return Hxy;
}