Zeta.cpp

#include "Headers.h"
#include "Zeta.h"


/* coefficients for Maclaurin summation in hzeta()
 * B_{2j}/(2j)!
 */
static double hzeta_c[15] = {
  1.00000000000000000000000000000,
  0.083333333333333333333333333333,
 -0.00138888888888888888888888888889,
  0.000033068783068783068783068783069,
 -8.2671957671957671957671957672e-07,
  2.0876756987868098979210090321e-08,
 -5.2841901386874931848476822022e-10,
  1.3382536530684678832826980975e-11,
 -3.3896802963225828668301953912e-13,
  8.5860620562778445641359054504e-15,
 -2.1748686985580618730415164239e-16,
  5.5090028283602295152026526089e-18,
 -1.3954464685812523340707686264e-19,
  3.5347070396294674716932299778e-21,
 -8.9535174270375468504026113181e-23
};


int HZeta(const double s, const double q, MathResult * result)
{
    /* CHECK_POINTER(result) */

    if(s <= 1.0 || q <= 0.0)
    {
        DOMAIN_ERROR(result);
    }
    else
    {
        const double max_bits = 54.0;
        const double ln_term0 = -s * log(q);  

        if(ln_term0 < MATH_LOG_DBL_MIN + 1.0)
        {
            UNDERFLOW_ERROR(result);
        }
        else if(ln_term0 > MATH_LOG_DBL_MAX - 1.0)
        {
            OVERFLOW_ERROR (result);
        }
        else if((s > max_bits && q < 1.0) || (s > 0.5*max_bits && q < 0.25))
        {
            result->val = pow(q, -s);
            result->err = 2.0 * MATH_DBL_EPSILON * fabs(result->val);
            return MATH_SUCCESS;
        }
        else if(s > 0.5*max_bits && q < 1.0)
        {
            const double p1 = pow(q, -s);
            const double p2 = pow(q/(1.0+q), s);
            const double p3 = pow(q/(2.0+q), s);
            result->val = p1 * (1.0 + p2 + p3);
            result->err = MATH_DBL_EPSILON * (0.5*s + 2.0) * fabs(result->val);
            return MATH_SUCCESS;
        }
        else
        {
            /* Euler-Maclaurin summation formula 
            * [Moshier, p. 400, with several typo corrections]
            */
            const int jmax = 12;
            const int kmax = 10;
            int j, k;
            const double pmax  = pow(kmax + q, -s);
            double scp = s;
            double pcp = pmax / (kmax + q);
            double ans = pmax*((kmax+q)/(s-1.0) + 0.5);

            for(k=0; k<kmax; k++)
            {
                ans += pow(k + q, -s);
            }

            for(j=0; j<=jmax; j++)
            {
                double delta = hzeta_c[j+1] * scp * pcp;
                ans += delta;
                if(fabs(delta/ans) < 0.5*MATH_DBL_EPSILON) break;
                scp *= (s+2*j+1)*(s+2*j+2);
                pcp /= (kmax + q)*(kmax + q);
            }

            result->val = ans;
            result->err = 2.0 * (jmax + 1.0) * MATH_DBL_EPSILON * fabs(ans);
            return MATH_SUCCESS;
        }
    }
}