/*
    zeta.c, an attempt to compute the Riemann Zeta function
    Copyright (C) 2000 Robert P. Munafo

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

    If you got ries.c from the website www.mrob.com, the
    GNU General Public License may be retrieved from the following URL:

    http://www.mrob.com/ries/COPYING.txt

    to compile (under Linux):

    gcc -o zeta -lm zeta.c 

    If given an argument T, it prints a graph of the Zeta function on the
    critical line (z = 0.5 + i t) from t = T-100 to t = T.

    If not given an argument, it looks for a file called "cdda.wav"
    (which should be 44100-Hz, PCM, 16-bit, stereo) and overwrites it
    with the Zeta function, samples spaced 0.01 apart on the complex
    plane. Thus, each second of audio will cover 441.0 units on the
    imaginary axis.

    This program does not "know" enough about the WAV file format to
    be able to create a WAV file from scratch -- however I suspect that
    would be fairly easy to do.
    I have a description of the AIFF format in "AIFF-1.3.pdf" (in
    another project) but do not currently have a WAV format description.

 */

/*
 Revision history:

 20000302 Initial version
 20000303 Post query to sci.math
 20000304 Add notes from a reply from sci.math

 20050804 Lots of fixes to make it compatible with official ISO
implementation of __complex__. Make little-endian-ness of WAV format
explicit (so it works on a PowerPC). Run it on a 2-GHz G5, the speed
is truly impressive! Increase the number of terms from 1000 to 10000.

 */

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

typedef signed char    s8;
typedef unsigned char  u8;
typedef short          s16;
typedef unsigned short u16;
typedef long           s32;
typedef unsigned long  u32;
typedef long long int  s64;

// Size of wav file header
#define WAV_START 44

#define LOUDNESS 1000.0

/* FREQUENCY is the number of units on the imaginary axis that will be
   heard each second. Imagine that the Zeta function has been carved
   into a phonograph record and we're playing the record. This constant
   tells how many units on the real axis the needle moves each second. */
#define FREQUENCY 441.0

/* INTERVAL is the spacing between values of the argument we pass to the
   Zeta function to produce the audio samples. Note the constant 44100 --
   we assume that our output sound file is 44100 samples per second. */
#define INTERVAL (FREQUENCY / 44100.0)

/* Function prototypes */
extern "C" {
  int main(int, char **);
  double logaddexp(double a, double b);
  double cabssqr(__complex__ double z);
  double cabs(__complex__ double z);
  __complex__ double clog(__complex__ double z);
  __complex__ double cexp(__complex__ double z);
  __complex__ double cpow(double b, __complex__ double e);
 
  void cprint(__complex__ double z);
  __complex__ double gamma1(__complex__ double z);

  __complex__ double zetacrit1(double zi, s32 terms);
}


double g_pi = 3.1415926535897932;


/* logaddexp returns ln(exp(a) + exp(b)), but works even when both exp's
   would overflow and/or underflow */
double logaddexp(double a, double b)
{
  double l, s;

  if (a > b) {
    l = a; s = b;
  } else {
    l = b; s = a;
  }

  /* log(a+b) = log(a) + log(1 + b/a)
              = log(a) + log(1 + exp(log(b)-log(a)))
              = l + log(1 + exp(s-l)) */
  if (s-l < -50.0) {
    /* exp(s-l) reduces to 0 */
    return l;
  } else {
    return (l + log(1.0 + exp(s-l)));
  }
}

double cabssqr(__complex__ double z)
{
  return((__real__ z)*(__real__ z) + (__imag__ z)*(__imag__ z));
}

double cabs(__complex__ double z)
{
  return(sqrt(cabssqr(z)));
}

__complex__ double clog(__complex__ double z)
{
  return(log(cabs(z)) + (atan2(__imag__ z, __real__ z) * 1.0fi));
}

/* complex exponential function */
__complex__ double cexp(__complex__ double z)
{
  __complex__ double rv;

  rv = exp(__real__ z) * (cos(__imag__ z) + (1.0fi * sin(__imag__ z)));
  return(rv);
}


/* real base to a complex power */
__complex__ double cpow(double b, __complex__ double e)
{
  __complex__ double rv;

  //printf("base %10.6lg  log(base) %10.6lg\n", b, log(b));
  //printf("  exponent  %10.6lg + %10.6lg i\n", (__real__ e), (__imag__ e));

  rv = cexp(log(b) * e);
  return(rv);
}

void cprint(__complex__ double z)
{
  printf("%14.10lg + %14.10lg i", __real__ z, __imag__ z);
}

/* This is an attempt to apply analytic extension to the Stirling series
   for the Gamma function. I haven't verified its accuracy yet. I am using
   the Stirling series (generalized to analytic complex functions) and
   the identity:

   gamma(z+1) = z gamma(z)
   gamma(z) = gamma(z+1) / z

   and the (similar but not identical) identity given on the Weisstein page:

   ln(gamma(x + iy + 1)) = ln(x^2 + y^2) + i arctan(y/x) + ln(gamma(x + iy))
   ln(gamma(x + iy)) = ln(gamma(x + iy + 1)) - ln(x^2 + y^2) - i arctan(y/x)

   this second identity would be identical to the first, if only they had
   "ln(sqrt(x^2 + y^2))" instead of "ln(x^2 + y^2)". I suspect it's an error.

   The whole thing is computed in logs and the answer computed with exp()
   at the end. Thus, it is readily converted to a log(gamma(z)) function.
*/
__complex__ double gamma1(__complex__ double z)
{
  __complex__ double scale;
  __complex__ double l;
  __complex__ double zp;

  /* need to test for negative integers and zero */

  scale = 0.0i + 0.0;
  /* On each identity transformation, set its limit to "0.0" to disable,
     around 5.0-10.0 to enable */
  while(cabs(z) < 0.0) { /* strange identity, disabled because it's suspect */
    scale = scale + log(cabssqr(z)) + (1.0fi * atan2(__imag__ z, __real__ z));
    z = z + 1.0;
  }
  while (((__real__ z) <= 0.0) || (cabs(z) < 8.0)) {
    scale = scale + clog(z);
    z = z + 1.0;
  }
  l = 0.0 - scale;

  /* The series I use for gamma is actually the Stirling series adapted for
     computing factorials, that is, x! = gamma(x+1). Rather than change the
     series back to a real gamma series, I just subtract 1 from z. Works just
     fine (-: */
  z = z - 1.0;

  l = l + 0.5 * log(2.0 * g_pi);
  l = l + (z + 0.5) * clog(z);
  l = l - z;
  l = l + 1.0 / (12.0 * z);
  //printf("gamma 1 "); cprint(cexp(l)); printf("\n");
  zp = z * z * z; /* z^3 */
  l = l - 1.0 / (360.0 * zp);
  //printf("gamma 2 "); cprint(cexp(l)); printf("\n");
  zp = zp * z * z; /* z^5 */
  l = l + 1.0 / (1260.0 * zp);
  //printf("gamma 3 "); cprint(cexp(l)); printf("\n");
  zp = zp * z * z; /* z^7 */
  l = l - 1.0 / (1680.0 * zp);
  //printf("gamma 4 "); cprint(cexp(l)); printf("\n");
  zp = zp * z * z; /* z^9 */
  l = l + 1.0 / (1188.0 * zp);
  printf("gamma 5 "); cprint(cexp(l)); printf("\n");

  return(cexp(l));
}

/* Here are some formulas to try, but some are probably only valid for
   real values, and most probably won't be useful:

   zeta(z) = 1/(1-2^(-z)) PRODUCT(p prime) 1/(1-p^(-z))
   don't know if "p prime" includes p=2 -- if not, the first term could
   be moved inside the product

   zeta(1-z) = (2/(2 pi)^z) cos(pi z/2) gamma(z) zeta(z)
   zeta(1/2 - it) = (2/(2 pi)^(1/2 + it)) cos(pi (1/2 + it)/2)
                    gamma(1/2 + it) zeta(1/2 + it)

   The "main term" of the Riemann-Seigel Z function:
   Z(t) ~= 2 SIGMA[k=1..v(t)] cos(theta(t) - t ln k) / sqrt(k)
   theta(t) = arg gamma(1/4 + it/2) - (t ln pi)/2
   v(t) = floor(sqrt(t/(2 pi)))
*/

/*
   This is a possibly-futile attempt to compute the Riemann-Siegel z(t)
   function, that is, |zeta(1/2 + i t)|

   For complex s with real(s) > 0, we have:

     Z(s) = 1/(1-2^(1-s)) SIGMA[n=1..inf] -1^(n-1) n^(-s)
          = 1/(1-2^(1-s)) {  SIGMA[odd n>0] n^(-s)
                           - SIGMA[even n>0] n^(-s)  }

   The terms converge at a rate of 1/sqrt(n), which is very slowly.
   They're smaller if you add consecutive terms, but it still converges
   just as slowly.
 */
__complex__ double zetacrit1(double zi, s32 terms)
{
  s32 i;
  __complex__ double z;
  __complex__ double m; /* the first multiplier, 1/(1-2^(1-s)) */
  __complex__ double t; /* one term in the sum */
  double n;
  __complex__ double sum;
  __complex__ double t1;
  __complex__ double rv;

  /* Compute the value of z on the critical line */
  z = 0.5 + (zi * 1.0fi);

  /* compute the first multiplier */
  m = 1.0 / (1.0 - cpow(2.0, (1.0 - z)));

  if (0) {
    printf("m = %10.6lg + %10.6lg i\n",
           (__real__ m), (__imag__ m));
  }

  /* Compute first term */
  n = 1.0;
  t = cpow(n, 0.0 - z);
  sum = t;

  /* Add some terms */
  i = 0;
  while (i < terms) {
    /* a negative term */
    n += 1.0;
    t1 = cpow(n, 0.0 - z);

    /* and a positive term */
    n += 1.0;
    t = cpow(n, 0.0 - z);
    t = t - t1;
    sum += t;
    if(0) {
      printf("n=%7lg  2 terms = %10.6lg + %10.6lg i", n,
             (__real__ t), (__imag__ t));
      printf("    sum = %10.6lg + %10.6lg i\n",
             (__real__ sum), (__imag__ sum));
    }
    i++;
  }

  //  printf("    sum = %10.6lg + %10.6lg i\n", (__real__ sum), (__imag__ sum));
  rv = m * sum;
  //  printf("    rv = %10.6lg + %10.6lg i\n", (__real__ rv), (__imag__ rv));
  return(rv);
}



/*
 From: Raymond Manzoni <raymman@club-internet.fr>
 Date: 03/04/2000
       sci.math
  Hi,
  
   Thank you for the take-off sound, 
   For which country? I don't know!
  
   Here are some of the papers that could interest you:
  
 A nice historical introduction may be found in 
 "Computational Number Theory at CWI in 1970-1994"
 but the paper (Riele-Lune.ps) doesn't seem available online any longer. 

 At CECM you'll find :
 "Computational strategies for the Riemann zeta function" at    
 http://www.cecm.sfu.ca/preprints/1998pp.html#98:118

 consult too Andrew Odlyzko's "Papers on Zeros of the Riemann Zeta
 Function and Related Topics" at
 http://www.research.att.com/~amo/doc/zeta.html (for example "Fast
 algorithms for multiple evaluations of the Riemann zeta function")

 and Andrew's "Tables of zeros of the Riemann zeta function" at
 http://www.research.att.com/~amo/zeta_tables/index.html

 For more papers concerning Riemann zeta look at :
 http://www.archetipo.web66.com/zeta/briefintro.htm


 Or forget all that and use simply the quick pari/gp for evaluation of
 zeta in the critical strip (available with C source at
 http://www.parigp-home.de/ )

 ftp://megrez.math.u-bordeaux.fr/pub/pari/unix/

  
 Hope it helped a little,

                 Raymond Manzoni

   I added too a resume of my earlier attempts using Euler Mac-Laurin : 

     I used the following formula  for that purpose:
 
     zeta(x) ~ sum(1/k^x, k=1..N)
               + 1/((x-1)*N^(x-1))
               -1/(2*N^x)
               +x/(12*N^(x+1))
               -x*(x+1)*(x+2)/(720*N^(x+3))
               +x*(x+1)*(x+2)*(x+3)*(x+4)/(30240*N^(x+5))
               +...

     The greater the N the better the result (but even N=10 is good) 
  
     The accuracy is rather fine when |Im(x)| < 2*PI*N and becomes
     very bad after that !
  
     You may find more terms using the more general Euler Mac-Laurin
     formula : This allows to explore better the negative values of x
     (in fact you will get the exact values for x=-1,-2, until -2
     times the number of Bernoulli terms!). With more terms the
     precision in the imaginary plane is better under 2*PI*N and worse
     over it!
  
    f(x) ~  sum(f(k), k=1..N)
             + int(f(u),u=N..oo)
             -f(N)/2
             -B1/2!*f'(N)
             +B2/4!*f'''(N)
             -...
     where Bn means the nth Bernoulli number

     This formula was very fruitful for me even on my little TI 92.
     (N=100 was enough)

 */



int main(int nargs, char ** argv)
{
  char ** av;
  int   args;
  char * arg;
  double zi;
  s16   got_zi;  /* got z? */

  got_zi = 0;

  if (nargs > 1) {
    av = argv;
    args = nargs;

    av++; args--;  /* skip our program name/path */

    while(args) {
      arg = *av; av++; args--;
      
      /* Test for number */
      if (((arg[0] >= '0') && (arg[0] <= '9'))
          || (arg[0] == '-') || (arg[0] == '.')) {
        sscanf(arg, "%lf", &zi);
        got_zi = 1;
      }
    }
  }

  if (0) {
    printf("zetacrit1(0.5 + b i) = %14.10lg\n", cabs(zetacrit1(zi, 100)));
    printf("                1000   %14.10lg\n", cabs(zetacrit1(zi, 1000)));
    printf("               10000   %14.10lg\n", cabs(zetacrit1(zi, 10000)));
    printf("              100000   %14.10lg\n", cabs(zetacrit1(zi, 100000)));
    printf("             1000000   %14.10lg\n", cabs(zetacrit1(zi, 1000000)));
    exit(0);
  }

  if(0) {
    gamma1(0.25 + 0.0i);
    gamma1(3.0 + 0.0i); /* expect 2.0 */
    gamma1(5.0 + 0.0i); /* expect 24.0 */
    gamma1(0.25 + (1.0fi * zi));
    gamma1(0.25 + (1.0fi * (zi + 1.0)));
    gamma1(0.25 + (1.0fi * (zi + 2.0)));
    exit(0);
  }

  if(got_zi) {
    __complex__ double z;
    double t, z1, z2r, z2i;
    s32 terms;

    terms = 500;
    t = zi - 100.0;
    if (t < 0.0) {
      t = 0.0;
    }
    for(; t<zi; t += 0.1) {
      s16 i, i1, i2i, i2r, lim;

      z1 = __real__ zetacrit1(t, terms);
      z = zetacrit1(t, 10 * terms);
      //      printf(" %10.6lg + %10.6lg i\n", (__real__ z), (__imag__ z));
      z2r = __real__ z; z2i = __imag__ z;

      i1 = ((int)(z1 * 10.0)) + 40;
      i2i = ((int)(z2i * 10.0)) + 40;
      i2r = ((int)(z2r * 10.0)) + 40;

      printf("%6.4lg", t);
      lim = i1;
      if (i2i > lim) {
        lim = i2i;
      }
      if (i2r > lim) {
        lim = i2r;
      }
      if (40 > lim) {
        lim = 40;
      }
      if(lim > 80) {
        lim = 80;
      }
      for(i=0; i<=lim; i++) {
        if (i == i2i) {
          printf("i");
        } else if (i == i2r) {
          printf("r");
        } else if (i == i1) {
          printf(".");
        } else if (i == 40) {
          printf("|");
        } else {
          printf(" ");
        }
      }
      printf("\n");
    }
  } else {
    __complex__ double z;
    double t, z1r, z1i;
    s32  nsamp, ctr, cs;
    s16  samp[2];
    s8   dbuf[4];
    FILE * wav;

    wav = fopen("cdda.wav", "r+");
    fseek(wav, 0, SEEK_END);
    nsamp = ftell(wav);
    nsamp -= WAV_START;
    nsamp /= 4;
    nsamp--;

    printf("Length: %ld samples\n", nsamp);
    printf("Evaluating Zeta from 0.0 to %10f\n", ((double) nsamp) * INTERVAL);

    fseek(wav, WAV_START, SEEK_SET);

    t = 0.0; ctr = 0; cs = 0;
    while(nsamp > 0) {
      z = zetacrit1(t, 10000);
      z1r = __real__ z; z1i = __imag__ z;
      samp[0] = ((int) (z1i * LOUDNESS));
      samp[1] = ((int) (z1r * LOUDNESS));
      dbuf[0] = samp[0] & 0xff; dbuf[1] = samp[0] >> 8;
      dbuf[2] = samp[1] & 0xff; dbuf[3] = samp[1] >> 8;
      fwrite((void *) dbuf, 2, 2, wav);
      ctr++;
      if (ctr >= 441) {
        printf("."); fflush(stdout);
        ctr = 0;
        cs++;
        if (cs >= 68) {
          printf("\n%10.6lg", t); fflush(stdout);
          cs = 0;
        }
      }
      nsamp--;
      t += INTERVAL;
    }

    fclose(wav);
  }

  return(0);
}