/*
  Author:  Pate Williams (c) 1997

  "Exercise 2.2 Computing zeta(s).
  zeta(s) may be computed by aid of the
  Euler-Maclaurin sum formula:
  zeta(s) = sum(n = 1 to N - 1, n ^ - s)
          + N ^ (1 - s) / (s - 1)
          + N ^ - s / 2
          + s N ^ (- s - 1) / 12
          - s(s + 1)(s + 2) N ^ (- s - 3) / 720
          + s(s + 1)...(s + 4) N ^ (- s - 5) / 30240
          - ...
  This is the so-called semiconvergent asymptotic
  series. the truncation error made in breaking
  off the summation is always less than the
  immediately following term. Write a function
  zeta(s) performing the summation of the series.
  In single precision arithmetic (about 8
  decimal digits), a suitable value if N is 10.
  Break off the series immediately before the
  last term written out above. Check your values
  against known exact values of
  zeta(2) = pi ^ 2 / 6 , zeta(4) = pi ^ 4 / 90,
  zeta(6) = pi ^ 6 / 945, zeta(8) = pi ^ 8 / 9450
  and zeta(10) = pi ^ 10 / 93555." -Hans Riesel-
  See "Primes Numbers and Computer Methods for
  Factorization" by Hans Riesel second edition
  pages 45-46.
*/

#include <math.h>
#include <stdio.h>

#define N 64

double zeta(double s)
{
  double sum = 0, x = - s;
  long n;

  for (n = 1; n < N; n++) sum += pow(n, x);
  sum += pow(N, 1 - s) / (s - 1)
      + 0.5 * pow(N, - s) + s * pow(N, - s - 1) / 12
      - s * (s + 1) * (s + 2) * pow(N, - s - 3) / 720;
  return sum;
}

int main(void)
{
  double pi = M_PI;
  double zeta2 = pi * pi / 6;
  double zeta4 = pow(pi, 4) / 90;
  double zeta6 = pow(pi, 6) / 945;
  double zeta8 = pow(pi, 8) / 9450;
  double zeta10 = pow(pi, 10) / 93555.0;
  double s[5] = {2, 4, 6, 8, 10};
  double z, zetai[5];
  long i;

  zetai[0] = zeta2;
  zetai[1] = zeta4;
  zetai[2] = zeta6;
  zetai[3] = zeta8;
  zetai[4] = zeta10;
  for (i = 0; i < 5; i++) {
    z = zeta(s[i]);
    printf("%2.0lf %lf %lf %lf\n", s[i], z, zetai[i],
           fabs(z - zetai[i]));
  }
  return 0;
}