static double f(double z){
    return log((1+z) / (1-z));
}

/* integral of z*f(z) */
static double fa(double z){
    return (z*z-1) / 2 * f(z) + z;
}
/* integral of z^3*f(z) */
static double fb(double z){
    double zz;
    zz = z*z;
    return (zz*zz-1) / 4 * f(z) + (zz*z/3 + z) / 2;
}
/* integral of z^5*f(z) */
static double fc(double z){
    double zz;
    zz = z*z;
    return (zz*zz*zz-1) / 6 * f(z) + (zz*zz*z/5 + zz*z/3 + z) / 3;
}

/*
 *  method of least squares
 *  D= integral of {f(z) - (az+bz^3+cz^5)}^2, (-S<= z <= S)
 *  get the value of (a, b, c) which attains the least D
 *  dD/da = dD/db = dD/dc = 0
 *  if and only if (a b*S^2 c*S^4) A = (fa/(2S^3), fb/(2S^5), fc/(2S^7))
 *  Here,
      (1/3 1/5 1/7)
    A=(1/5 1/7 1/9)
      (1/7 1/9 1/11)
    A^-1 = {cij}
 */
static void calc_coef(double *pa, double *pb, double *pc){
    double S, SSS;
    double c11, c12, c13, c22, c23, c33;
    double a, b, c;
    double p, q, r;
    double det;

    S = (sqrt(2)-1);S *= S;

    c11=1.0/77 - 1.0/81;
    c12=1.0/63 - 1.0/55;
    c13=1.0/45 - 1.0/49;
    c22=1.0/33 - 1.0/49;
    c23=1.0/35 - 1.0/27;
    c33=1.0/21 - 1.0/25;
    det = c11/3 + c12/5 + c13/7;
    c11 /= det;
    c12 /= det;
    c13 /= det;
    c22 /= det;
    c23 /= det;
    c33 /= det;

    SSS = S * S * S;
    p = fa(S) / (2*SSS);
    q = fb(S) / (2*SSS*S*S);
    r = fc(S) / (2*SSS*SSS*S);

    a = c11*p + c12*q + c13*r;
    b = c12*p + c22*q + c23*r; b /=(S*S);
    c = c13*p + c23*q + c33*r; c /=(SSS*S);
//  printf("a-1  =%12.8f\n", a-1);
//  printf("b-1/3=%12.8f\n", b - 1.0/3);
//  printf("c-1/5=%12.8f\n", c - 0.2);
    *pa = a;
    *pb = b;
    *pc = c;
}