BBP Pi mod By ShuffleSource

/*
    This program implements the BBP algorithm to generate a few hexadecimal
    digits beginning immediately after a given position id, or in other words
    beginning at position id + 1.  On most systems using IEEE 64-bit floating-
    point arithmetic, this code works correctly so long as d is less than
    approximately 1.18 x 10^7.  If 80-bit arithmetic can be employed, this limit
    is significantly higher.  Whatever arithmetic is used, results for a given
    position id can be checked by repeating with id-1 or id+1, and verifying
    that the hex digits perfectly overlap with an offset of one, except possibly
    for a few trailing digits.  The resulting fractions are typically accurate
    to at least 11 decimal digits, and to at least 9 hex digits.
*/

/*  David H. Bailey     2006-09-08 */
/* !!!MODIFIED BY BALÁZS KÓTI     2013-03-14!!! */

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

int main(int argc, char *argv[])
{
  if (argc >= 2)
  {
    //paraméterek
    int szamjegyek = atoi(argv[1]);
    int kezdo; //hanyadik rész?
    int kiir;
    switch (argc)
    {
      case 2:
        kezdo = 1;
        kiir = 1;
      case 3:
        kezdo = atoi(argv[2]);
        kiir = 1;
        break;
      case 4:
        kezdo = atoi(argv[2]);
        kiir = atoi(argv[3]);
        break;
      default:
        printf("hiba!\n");
        return 0;
        break;
    }

    //csak 10-el osztható darab számot lehet generálni
    if (szamjegyek%10 != 0)
    {
      printf("10-el osztható darab számjegyet lehet csak generálni!\n");
      return 0;
    }

    //fájlkezelés
    char fileName[FILENAME_MAX];
    sprintf(fileName, "PiGen_%d_%d.txt", szamjegyek, kezdo);
    FILE *file;
    file = fopen(fileName ,"a+");

    //számítás
    double pid, s1, s2, s3, s4;
    double sorozat (int m, int n);
    void ihex (double x, int m, char c[]);
    #define NHX 16
    char chx[NHX];
    float toredek,folyamat;
    int kiirszamlalo = 0;
    //
    printf("====PiGen by Shuffle=====\n");
    printf("Kezdo : %d, számjegyek: %d\n", kezdo, szamjegyek);
    printf("Mentett fájl: %s\n", fileName );
    if (kiir == 0)
      printf("Kiíratás kikapcsolva!\n");
    for(toredek=0;toredek<szamjegyek/10;toredek++)
    {
      int toredekMul = ((kezdo - 1) * szamjegyek) + (toredek * 10);
      s1 = sorozat (1, toredekMul);
      s2 = sorozat (4, toredekMul);
      s3 = sorozat (5, toredekMul);
      s4 = sorozat (6, toredekMul);
      pid = 4. * s1 - 2. * s2 - s3 - s4;
      pid = pid - (int) pid + 1.;
      ihex (pid, NHX, chx);
      fprintf(file, "%10.10s", chx );
      //folyamat kiírása
      if (kiir == 1 && kiirszamlalo == 20)
      {
        if (system("cls")) system("clear");
        printf("====PiGen by Shuffle=====\n");
        printf("Kezdo : %d, számjegyek: %d\n", kezdo, szamjegyek);
        printf("Mentett fájl: %s\n", fileName );
        folyamat = (toredek/(szamjegyek/10))*100;
        printf("Gerenálás...%.2f\n",folyamat);
        printf("Aktuális csonk: %10.10s\n",chx);
        kiirszamlalo = 0;
      }

      kiirszamlalo++;
    }
    printf("KÉSZ!\n");
    fclose(file);
    /* printf (" position = %i\n fraction = %.15f \n hex digits =  %10.10s\n",
    toredek, pid, chx);*/
  }
  else
    printf("Nincs megadva paraméter! Paraméterek: <számjegyek száma> <hanyadik rész> <kiíratás ki/be (0/1)> \n");
  return 0;
}

void ihex (double x, int nhx, char chx[])

/*  This returns, in chx, the first nhx hex digits of the fraction of x. */

{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
    y = 16. * (y - floor (y));
    chx[i] = hx[(int) y];
  }
}

double sorozat (int m, int id)

/*  This routine evaluates the sorozat  sum_k 16^(id-k)/(8*k+m)
    using the modular exponentiation technique. */

{
  int k;
  double ak, eps, p, s, t;
  double expm (double x, double y);
#define eps 1e-17

  s = 0.;

/*  Sum the sorozat up to id. */

  for (k = 0; k < id; k++){
    ak = 8 * k + m;
    p = id - k;
    t = expm (p, ak);
    s = s + t / ak;
    s = s - (int) s;
  }

/*  Compute a few terms where k >= id. */

  for (k = id; k <= id + 100; k++){
    ak = 8 * k + m;
    t = pow (16., (double) (id - k)) / ak;
    if (t < eps) break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}

double expm (double p, double ak)

/*  expm = 16^p mod ak.  This routine uses the left-to-right binary
    exponentiation scheme. */

{
  int i, j;
  double p1, pt, r;
#define ntp 25
  static double tp[ntp];
  static int tp1 = 0;

/*  If this is the first call to expm, fill the power of two table tp. */

  if (tp1 == 0) {
    tp1 = 1;
    tp[0] = 1.;

    for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
  }

  if (ak == 1.) return 0.;

/*  Find the greatest power of two less than or equal to p. */

  for (i = 0; i < ntp; i++) if (tp[i] > p) break;

  pt = tp[i-1];
  p1 = p;
  r = 1.;

/*  Perform binary exponentiation algorithm modulo ak. */

  for (j = 1; j <= i; j++){
    if (p1 >= pt){
      r = 16. * r;
      r = r - (int) (r / ak) * ak;
      p1 = p1 - pt;
    }
    pt = 0.5 * pt;
    if (pt >= 1.){
      r = r * r;
      r = r - (int) (r / ak) * ak;
    }
  }

  return r;
}