/*
  Author:  Pate Williams (c) 1997

  "Exercise 3.2. The twin prime constant.
  The twin prime constant C in (3.1) can be evaluated
  by the following:
  ln C / 2 = - sum(j = 2 to inf, (2 ^ j - 2) / j)
           * sum(p >= 3, p ^ - j).
  Now sum(p >= 3, p ^ - s) = P(s) - 2 ^ - s, where
  P(s) is the prime zeta-function studied in
  exercise 3.1 above. Use your function P(s) in a
  computer program which evaluates ln(C / 2) after
  the formula given above. Compare the result with
  the value C calculated from C given in (3.1)!"
  -Hans Riesel-
  See "Prime Numbers and Computer Methods for
  Factorization" by Hans Riesel second edition
  pages 65-66.
*/

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

#define J 32
#define K 64
#define N 64
#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define SIZE 8192

long get_bit(long i, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;

  return (sieve[c] >> (BITS_PER_LONG_1 - b)) & 1;
}

void set_bit(long i, long v, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;
  long mask = 1 << (BITS_PER_LONG_1 - b);

  if (v == 1)
    sieve[c] |= mask;
  else
    sieve[c] &= ~mask;
}

void Sieve(long n, long *sieve)
{
  long c, i, inc;

  set_bit(0, 0, sieve);
  set_bit(1, 0, sieve);
  set_bit(2, 1, sieve);
  for (i = 3; i <= n; i++)
    set_bit(i, i & 1, sieve);
  c = 3;
  do {
    i = c * c, inc = c + c;
    while (i <= n) {
      set_bit(i, 0, sieve);
      i += inc;
    }
    c += 2;
    while (!get_bit(c, sieve)) c++;
  } while (c * c <= n);
}

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;
}

double constant(void)
{
  double sum = 0;
  long j;

  for (j = 2; j <= J; j++)
    sum += (pow(2, j) - 2) * (P(j) - pow(2, - j)) / j;
  return 2.0 * exp(- sum);
}

double constant3_1(void)
{
  double d, n, product = 1;
  long p = 3, sieve[SIZE];

  Sieve(100000l, sieve);
  while (p < 100000l) {
    while (!get_bit(p, sieve)) p++;
    n = (double) p * (double) (p - 2);
    d = (double) (p - 1) * (double) (p - 1);
    product *= n / d;
    p++;
  }
  return 2 * product;
}

int main(void)
{
  double C = 1.320323632;
  double c = constant(), d = constant3_1();

  printf("calculated from ln(C / 2), c = %lf\n", c);
  printf("calculated from (3.1), D = %lf\n", d);
  printf("given C = %lf\n", C);
  printf("| c - C | = %lf\n", fabs(c - C));
  printf("| c - D | = %lf\n", fabs(c - d));
  return 0;
}