/*
  Author:  Pate Williams (c) 1997

  "Exercise 3.1. The prime zeta-function.
  The function P(s) = sum(p, p ^ - s),
  summed over all primes, is called the
  prime zeta-function. It is indispensable
  for the evaluation of expressions
  containing all primes such as the twin
  prime constant, occurring in (3.1), or
  the infinite product (3.6). Values of
  P(s) may be computed from values of the
  Riemann zeta-function zeta(s) in the
  following way: taking logarithms of
  both sides of (2.7) gives
  ln(zeta(s)) = sum(k = 1 to inf, P(ks) / k).
  Inversion of this formula gives
  P(s) = sum(k = 1 to inf, mu(k) ln zeta(ks) / k).
  Here mu(k) is Moebius function and the
  values of the zeta function are to be
  computed as in exercise 2.2 on p. 45. For
  large values of ks the following shortcut
  may be convienient:
  ln zeta(ks) = ln(1 + 2 ^ - ks + 3 ^ - ks + ...)
              = 2 ^ - ks + 3 ^ - ks + ...
  Write a function P(s) calling the function
  zeta(s) from exercise 2.2 and giving the
  value of P(s) s >= 2. Incorporate P(s) in
  a test program and check your programming
  of P(s) by re-computing ln zeta(s) from
  P(s) by the formula above and comparing the
  result with the result of a direct
  computation of ln zeta(s)." -Hans Riesel-
  See "Prime Numbers and Computer Methods for
  Factorization by Hans Riesel second edition
  page 65.
*/

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

#define K 64
#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 mu(long n)
{
  long e, i, k = 0, p;
  long prime[5] = {2, 3, 5, 7, 11};

  /* factor n by trial division */
  if (n == 1) return 1;
  for (i = 0; i < 5; i++) {
    p = prime[i];
    if (n % p == 0) {
      k++;
      e = 0;
      do {
        e++;
        n /= p;
      } while (n % p == 0);
      if (e > 1) return 0;
      if (n == 1) return (k & 1) ? - 1 : + 1;
    }
  }
  return 0;
}

double P(double s)
{
  double sum = 0;
  long k;

  for (k = 1; k <= K; k++)
    sum += mu(k) * log(zeta(k * s)) / k;
  return sum;
}

double log_zeta(double s)
{
  double sum = 0;
  long k;

  for (k = 1; k <= K; k++)
    sum += P(k * s) / k;
  return sum;
}

int main(void)
{
  double r, s, t;

  for (s = 2.0; s < 10.0; s += 0.5) {
    r = log(zeta(s));
    t = log_zeta(s);
    printf("%2.1lf %lf %lf %lf %lf\n", s, P(s),
           r, t, fabs(r - t));
  }
  return 0;
}