/***********************************************************************/
/*                                                                     */
/*                                          O                          */
/*                             OOOOOOOOOOOOO                           */
/*                            OOOOOOOOOOOOO                            */
/*                           O  OO     OO                              */
/*                              OO     OO                              */
/*                              OO     OO                              */
/*                              OO     OO                              */
/*                             OO       OO                             */
/*                            OO         OO                            */
/*                                                                     */
/*                                                                     */
/*                                                                     */
/*                                                                     */
/* Program:                                                            */
/* --------                                                            */
/*                                                                     */
/* sspi.c                                                              */
/*                                                                     */
/*                                                                     */
/* Purpose:                                                            */
/* --------                                                            */
/*                                                                     */
/* This program calculates Pi to many decimal digits.                  */
/*                                                                     */
/*                                                                     */
/* Version 1.08:                                                       */
/* -------------                                                       */
/*                                                                     */
/* First version: June 2003                                            */
/* Latest update: 2007-08-13                                           */
/*                                                                     */
/*                                                                     */
/* Author of the program:                                              */
/* ----------------------                                              */
/*                                                                     */
/* Stefan Spaennare, Lund, Sweden                                      */
/*                                                                     */
/* E-mail: stefan@spaennare.se                                         */
/*                                                                     */
/* Home page:    http://www.spaennare.se/index.html                    */
/* Program page: http://www.spaennare.se/ssprog.html                   */
/*                                                                     */
/* If you have problems with the program please let me know!           */
/*                                                                     */
/*                                                                     */
/* References:                                                         */
/* -----------                                                         */
/*                                                                     */
/* 1. Some of the algorithms are found on the web-page "Mathematical   */
/* constants and computation" by Xavier Gourdon and Pascal Sebah:      */
/*                                                                     */
/* http://numbers.computation.free.fr/Constants/constants.html         */
/*                                                                     */
/* 2. The web-page "Fast Algorithms and the FEE Method",               */
/* by Ekatherina A. Karatsuba:                                         */
/*                                                                     */
/* http://www.ccas.ru/personal/karatsuba/algen.htm                     */
/*                                                                     */
/* 3. The web-page "The Karatsuba Method 'Divide and Conquer'",        */
/* by Ekatherina A. Karatsuba:                                         */
/*                                                                     */
/* http://www.ccas.ru/personal/karatsuba/divcen.htm                    */
/*                                                                     */
/* 4. The article "Fast evaluation of transcendental functions",       */
/* Problems of Information Transmission, vol. 27, (1991), p. 339-360,  */
/* by Ekatherina A. Karatsuba.                                         */
/*                                                                     */
/* 5. The article "Fast multiprecision evaluation of series of         */
/* rational numbers", by Bruno Haible and Thomas Papanikolaou, 1997.   */
/*                                                                     */
/* 6. The program now use a very fast FFT (fftsg_h.c) by Takuya Ooura: */  
/*                                                                     */
/* http://momonga.t.u-tokyo.ac.jp/~ooura/                              */
/*                                                                     */
/* 7. "CLN - Class Library for Numbers" by Bruno Haible et. al:        */
/*                                                                     */
/* http://www.ginac.de/CLN/                                            */
/*                                                                     */
/* 8. The documentation "Programs to Calculate some Mathematical       */
/* Constants to Large Precision", by Stefan Spaennare.                 */
/*                                                                     */
/*                                                                     */
/* Program history:                                                    */
/* ----------------                                                    */
/*                                                                     */
/* 2007-08-13:                                                         */
/* The code has been somewhat prettified and shorter (in main).        */
/* The program has the same speed as before.                           */
/*                                                                     */
/* 2004-01-09:                                                         */
/* The name of the program has been changed to sspi.c.                 */
/*                                                                     */
/* 2003-10-09:                                                         */
/* Very minor changes.                                                 */
/*                                                                     */
/* 2003-08-27:                                                         */
/* Minor improvements to some precedures. Updated benchmarks.          */
/*                                                                     */
/* 2003-07-18:                                                         */
/* The multiplication version variable (mulversion) is now set         */
/* (to 1 or 2) in the file "mulver.txt".                               */
/*                                                                     */
/* 2003-06-30:                                                         */
/* Improved FFT-multiplication (lmul). Inversions, divisions and       */
/* square roots up to 25 % faster. Inversions and divisions with       */
/* 4-th order convergence.                                             */
/*                                                                     */
/* 2003-06-24:                                                         */
/* Minor changes to this text header.                                  */
/*                                                                     */
/* First version June 2003.                                            */
/*                                                                     */
/*                                                                     */
/* Copy-right:                                                         */
/* -----------                                                         */
/*                                                                     */
/* This program is free to copy if this information follows the        */
/* program. You are also allowed to modify the program for your        */
/* own calculations if you give a reference to the author.             */
/* See also notice below!                                              */
/*                                                                     */
/*                                                                     */
/* IMPORTANT!                                                          */
/* ----------                                                          */
/*                                                                     */
/* Please don't use the result from this program for important         */
/* purposes without doing an independent check of the result!          */
/* See also WARNING and notice below!                                  */
/*                                                                     */
/*                                                                     */
/* Algorithms:                                                         */
/* -----------                                                         */
/*                                                                     */
/*                                                                     */
/* The program uses (see also the references above):                   */
/*                                                                     */
/* * Chudnovsky or Ramanujan fast series for Pi.                       */
/*                                                                     */
/* * FFT for multiplications.                                          */
/*                                                                     */
/* * 4-th order methods for inversions and divisions.                  */
/*                                                                     */
/* * 4-th order methods for square and cube roots.                     */
/*                                                                     */
/*                                                                     */
/* Chudnovsky series for Pi is given by:                               */ 
/* -------------------------------------                               */
/*                                                                     */
/*                                                                     */
/*   1      12    inf (-1)^n (6n)! (A+Bn)                              */
/*   -- = ------- SUM -------------------                              */
/*   pi   C^(3/2) n=0 (3n)! (n!)^3 C^(3n)                              */
/*                                                                     */
/* where A=13591409, B=545140134, C=640320. Each term in the series    */
/* adds about 14 digits to Pi.                                         */
/*                                                                     */
/*                                                                     */
/* Ramanujan series for Pi is given by:                                */ 
/* ------------------------------------                                */
/*                                                                     */
/*                                                                     */
/*   1    2^(3/2) inf (4n)! (1103+26390n)                              */
/*   -- = ------- SUM -------------------                              */
/*   pi    9801   n=0 (n!)^4 (396)^(4n)                                */
/*                                                                     */
/* Each term in the series adds about 8 digits to Pi.                  */
/*                                                                     */
/*                                                                     */
/* Accuracy:                                                           */
/* ---------                                                           */
/*                                                                     */
/* The program calculates about 0.995*2^(power of 2) decimal digits.   */
/*                                                                     */
/* All printed digits are probably correct. The result for Pi has      */
/* been verified to 0.995*2^24 digits (for mulversion = 1).            */
/*                                                                     */
/*                                                                     */
/* FFT max error                                                       */
/* -------------                                                       */
/*                                                                     */
/* Set the variable mulversion=2 if "FFT max error" > 0.25.            */
/* This makes the program two times slower and requires more memory.   */
/*                                                                     */
/* The variable mulversion is present in the file "mulver.txt"         */
/* containing the single value 1 or 2. Default value is 1.             */
/*                                                                     */
/*                                                                     */
/* WARNING!                                                            */
/* --------                                                            */
/*                                                                     */
/* If you want to make calculations with much more than 2^24           */
/* (16 million) digits the convolutions in the FFT multiplications     */
/* can be saturated which will give erroneous results or no            */
/* convergence at all. If this occur then set mulversion=2 which       */
/* is 2 times slower than mulversion=1 but use the base 100            */
/* instead of 10000. This means that the convolution will be           */
/* saturated for a much larger number of digits (> 100 millions).      */
/*                                                                     */
/*                                                                     */
/* Compilation:                                                        */
/* ------------                                                        */
/*                                                                     */
/* The program works in Unix, Linux and Windows (DOS using DJGPP).     */
/* It is recommended that you optimize the program (i.e. use gcc -O3   */
/* option or similar option) when you compile:                         */
/*                                                                     */
/*>gcc -O3 -o sspi sspi.c fftsg_h.c -lm                                */
/*                                                                     */
/* The time wraps around after 4295 seconds on some computers.         */
/*                                                                     */
/*                                                                     */
/* Output:                                                             */
/* -------                                                             */
/*                                                                     */
/* Output is to the file "pi.txt". This file name can easily be        */
/* changed.                                                            */
/*                                                                     */
/*                                                                     */
/* Timing, benchmark and memory:                                       */
/* -----------------------------                                       */
/*                                                                     */
/* The CPU-times are given for a 1400 MHz Intel Celeron computer       */
/* (100 MHz memory buss) and 512 Mbyte memory. The cache memory is     */
/* 256 kbyte at full CPU speed. The operating system was Red Hat 9     */
/* Linux with the gcc 3.22 C-compiler. The CPU load was on average     */
/* 1.00 (i.e. no other programs running).                              */
/*                                                                     */
/* The program was compiled using the line:                            */
/*                                                                     */
/*>gcc -O3 -o sspi sspi.c fftsg_h.c -lm                                */
/*                                                                     */
/* The times have been measured for "power of 2" equals 17 and 20 with */
/* the "time" function in Linux. TC below means "Time Complexity".     */
/* The variable mulversion = 1.                                        */
/*                                                                     */
/* Intel Celeron 1400 MHz                                              */
/* ----------------------                                              ************/
/*                                                                                */
/* Algorithm           2^17       2^20        2^24         TC       Memory/(2^20) */
/*--------------------------------------------------------------------------------*/
/* Pi, Chudnovsky      3.78 s    57.26  s  1533.77 s  O(n*log(n)^3)  13.4 Mbyte   */
/* Pi, Ramanujan       5.81 s    72.26  s  2449.46 s  O(n*log(n)^3)  13.4 Mbyte   */
/*--------------------------------------------------------------------------------*/
/*                                                                                */
/**********************************************************************************/

/**********************************************************************************/
/*                                                                                */
/* Notice                                                                         */
/* ======                                                                         */
/*                                                                                */
/* I make no warranties that this program is (1) free of errors, (2) consistent   */
/* with any standard merchantability, or (3) meeting the requirements of a        */
/* particular application. This software shall not, partly or as a whole,         */
/* participate in a process, whose outcome can result in injury to a person or    */
/* loss of property. It is solely designed for analytical work. Permission to     */
/* use, copy and distribute is hereby granted without fee, providing that the     */
/* header above including this notice appears in all copies.                      */
/*                                                                                */
/*                                                            Stefan Spaennare    */
/*                                                                                */
/**********************************************************************************/


/**********************************************************************************/
/*                                                                                */
/* Comments by Stefan Spaennare January 2010                                      */
/* =========================================                                      */
/*                                                                                */
/* Many of the procedures in the common block of C-code of the program now        */
/* have got "how to use" comments.                                                */
/*                                                                                */
/* These invloves lfloat floating point functions for huge operands and bigint    */
/* for huge integer operands. The package also involves some procedures for       */
/* reading date and time from the computer clock and write to a C-program.        */ 
/*                                                                                */
/**********************************************************************************/


/**********************************************************************************/
/*                                                                                */
/* Initial includes and declarations and the struct for bigint etc.               */
/*                                                                                */
/**********************************************************************************/

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


double timetic=(double)(CLOCKS_PER_SEC);

double pi=3.14159265358979323;

char version_string[80]="Program sspi2.c, version 1.00, 2010-01-31 (based on sspi.c version 1.08)";

char file_name[20]="pi.txt";         /* Output filename */

char mv_file_name[20]="mulver.txt";  /* Multiplication version file name */

int mulversion;                      /* Multiplication version (1 or 2) */

double bas=10000.0;
double basi=0.0001;

double fftmaxerror;

typedef struct biginttag {
   short *coef;
   int size,maxsize,exp0,sgn0;
} bigint;

bigint tmp,tmp2,tmp3,tmp4;


/**********************************************************************************/
/*                                                                                */
/* mul: Auxiliary procedure for FFT multiplication used by lmul and bigmul.       */
/*                                                                                */
/**********************************************************************************/


void mul(d,n)
double *d;
int n;
{
   int nn,nn3,nn2,j;
   double x1,x2,y1,y2;
   
   nn=n+1;
   nn2=2*n+2;
   nn3=nn2+1;

   x1=d[0]; y1=d[1];
   x2=y1; y1=0.0; y2=0.0;
   d[0]=x1*x2-y1*y2;
   d[1]=x1*y2+y1*x2;

   for (j=3; j<=nn; j=j+2) {
      x1=0.5*(d[j-1]+d[nn2-j-1]);
      x2=0.5*(d[j-1]-d[nn2-j-1]);
      y1=0.5*(d[j]+d[nn3-j-1]);
      y2=0.5*(d[j]-d[nn3-j-1]);
      
      d[j-1]=x1*y1-y2*(-x2);
      d[j]=x1*(-x2)+y2*y1;
      d[nn2-j-1]=x1*y1-(-y2)*x2;
      d[nn3-j-1]=x1*x2+(-y2)*y1;
   } /* for j */

} /* mul */


/**********************************************************************************/
/*                                                                                */
/* lzero: Sets an lfloat operand a to zero and cleares the decimal array.         */
/*                                                                                */
/**********************************************************************************/

void lzero(a,nn)
short *a;
int *nn;
{
   int i,l,il,dl,ib;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   for (i=l;i>=ib; i--) {
      a[i]=0;
   } /* for i */
   
   a[l+1]=1;

} /* lzero */


/**********************************************************************************/
/*                                                                                */
/* linttol: Sets the integer part of an lfloat operand to an integer value        */
/* < 10000 and cleares the decimal array. The first word                          */
/* of the integer part of the lfloat operand will be affected.                    */
/*                                                                                */
/**********************************************************************************/

void linttol(a,man0,exp0,nn)
short *a;
int man0,exp0;
int *nn;
{
   int i,l,il,dl,ib;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   for (i=l;i>=ib; i--) {
      a[i]=0;
   } /* for i */

   a[l-il+exp0+1]=abs(man0);
  
   if (man0 >= 0) { 
      a[l+1]=1;
   }
   else {
      a[l+1]=-1;
  } /* if */

} /* linttol */


/**********************************************************************************/
/*                                                                                */
/* linttol2: Sets the integer part of an lfloat operand to an integer value       */
/* < 2000000000 and cleares the decimal array. The three first words              */
/* of the integer part of the lfloat operand will be affected.                    */
/*                                                                                */
/**********************************************************************************/
      
void linttol2(a,man0,exp0,nn)
short *a;
int man0,exp0;
int *nn;
{
   int i,l,il,dl,ib;
   int tt1,tt2,t1,t2,b;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   for (i=l;i>=ib; i--) {
      a[i]=0;
   } /* for i */

   b=abs(man0);

   if (b < 10000) {
      a[l-il+exp0+1]=b;
   }
   else {
      tt1=b % 100000000;
      tt2=b / 100000000;
      t1=tt1 % 10000;
      t2=tt1 / 10000;
      a[l-il+exp0+1]=t1;
      a[l-il+exp0+2]=t2;
      a[l-il+exp0+3]=tt2;
   } /* if */
  
   if (man0 >= 0) { 
      a[l+1]=1;
   }
   else {
      a[l+1]=-1;
  } /* if */

} /* linttol2 */


/**********************************************************************************/
/*                                                                                */
/* ppp: Prints an integer word to standard output.                                */
/*                                                                                */
/**********************************************************************************/

void ppp(t)
int t; 
{
  if (t > 999) {
     printf("%4d",t);
  }
  else if (t > 99) {
     printf("0%3d",t);
  } 
  else if (t > 9) {
     printf("00%2d",t);
  }
  else {
     printf("000%1d",t);
  } /* if */

} /* ppp */


/**********************************************************************************/
/*                                                                                */
/* lprint: Prints a lfloat operand to standard output.                            */
/*                                                                                */
/**********************************************************************************/

void lprint(a,nn)
short *a;
int *nn;
{

   int i,j,m,l,il,dl,t;

   l=nn[0]; il=nn[1]; dl=nn[3];

   if (a[l+1]==-1) {
      printf("-");
   } /* if */

   i=l;
   while ((a[i]==0) && (i>(l-il+1))) {
      i--;
   } /* while */

   m=0;

   if (i==(l-il+1)) {

      printf("%4d.",a[i]);

      m++;

      for (j=(l-il); j>=(l-il-dl+1); j--) {
            t=a[j];
            ppp(t);
            m++;
            if ((m % 18)==0) {
              printf("\n");
            } /* if */
      } /* for j */ 
   }
   else {

      printf("%d",a[i]);
      m++;

      for (j=(i-1); j>=(l-il+1); j--) {
            t=a[j];
            ppp(t);
            m++;
            if ((m % 18)==0) {
               printf("\n");
            } /* if */
      } /* for j */ 
     
      printf(".");

      for (j=(l-il); j>=(l-il-dl+1); j--) {
            t=a[j];
            ppp(t);
            m++;
            if ((m % 18)==0) {
               printf("\n");
            } /* if */
      } /* for j */ 
   } /* else */

   printf("\n\n");

} /* lprint */


/**********************************************************************************/
/*                                                                                */
/* fppp: Prints an integer word to a file.                                        */
/*                                                                                */
/**********************************************************************************/

void fppp(t,outfile)
int t; 
FILE *outfile;
{
  if (t > 999) {
     fprintf(outfile,"%4d",t);
  }
  else if (t > 99) {
     fprintf(outfile,"0%3d",t);
  } 
  else if (t > 9) {
     fprintf(outfile,"00%2d",t);
  }
  else {
     fprintf(outfile,"000%1d",t);
  } /* if */

} /* fppp */


/**********************************************************************************/
/*                                                                                */
/* lfprint: Prints a lfloat operand to a file.                                    */
/*                                                                                */
/**********************************************************************************/

void lfprint(a,nn,outfile)
short *a;
int *nn;
FILE *outfile;
{

   int i,j,m,l,il,dl,t;

   l=nn[0]; il=nn[1]; dl=nn[3];

   if (a[l+1]==-1) {
      fprintf(outfile,"-");
   } /* if */

   i=l;
   while ((a[i]==0) && (i>(l-il+1))) {
      i--;
   } /* while */

   m=0;

   if (i==(l-il+1)) {

      fprintf(outfile,"%4d.",a[i]);

      m++;

      for (j=(l-il); j>=(l-il-dl+1); j--) {
            t=a[j];
            fppp(t,outfile);
            m++;
            if ((m % 18)==0) {
              fprintf(outfile,"\n");
            } /* if */
      } /* for j */ 
   }
   else {

      fprintf(outfile,"%d",a[i]);
      m++;

      for (j=(i-1); j>=(l-il+1); j--) {
            t=a[j];
            fppp(t,outfile);
            m++;
            if ((m % 18)==0) {
               fprintf(outfile,"\n");
            } /* if */
      } /* for j */ 
     
      fprintf(outfile,".");

      for (j=(l-il); j>=(l-il-dl+1); j--) {
            t=a[j];
            fppp(t,outfile);
            m++;
            if ((m % 18)==0) {
               fprintf(outfile,"\n");
            } /* if */
      } /* for j */ 
   } /* else */

   fprintf(outfile,"\n");

} /* lfprint */


/**********************************************************************************/
/*                                                                                */
/* lass: a = b; for lfloat operands.                                              */
/*                                                                                */
/**********************************************************************************/

void lass(a,b,nn)
short *a;
short *b;
int *nn;
{
  int i,l,il,dl,ib;

  l=nn[0]; il=nn[1]; dl=nn[2];
  
  ib=l-il-dl+1;

  for (i=l+1; i>=ib; i--) {
     a[i]=b[i];
  } /* for i */

} /* lass */


/**********************************************************************************/
/*                                                                                */
/* add: Help function for ladd for lfloat operands. Performs the actual addition  */
/* without taking into account the sign of the operands (a=bb+cc). a, bb and cc   */
/* can be the same variable.                                                      */
/*                                                                                */
/**********************************************************************************/

void add(a,bb,cc,nn)
short *a;
short *bb;
short *cc;
int *nn;
{
  int i,ib,il,dl,l;
  int q,r,basint;

  short *b,*c;

  l=nn[0]; il=nn[1]; dl=nn[2];

  ib=l-il-dl+1;

  b = (short *)calloc(l+2,sizeof(short));
  c = (short *)calloc(l+2,sizeof(short));

  lass(b,bb,nn);
  lass(c,cc,nn);

  basint=10000;
  r=0;

  for (i=ib; i<=l; i++) {
     q=b[i]+c[i]+r;
     if (q >= basint) {
        a[i]=q-basint;
	r=1;
     }
     else {
        a[i]=q;
	r=0;
     } /* if */
  } /* for i */

  free(b);
  free(c);

} /* add */

      
/**********************************************************************************/
/*                                                                                */
/* sub: Help function for lsub for lfloat operands. Performs the actual           */
/* subtraction without taking into account the sign of the operands (a=bb-cc).    */
/* a, bb and cc can be the same variable.                                         */
/*                                                                                */
/**********************************************************************************/

void sub(a,bb,cc,nn)
short *a;
short *bb;
short *cc;
int *nn;
{
  int i,ib,il,dl,l;
  int q,r,basint;

  short *b,*c;

  l=nn[0]; il=nn[1]; dl=nn[2];

  ib=l-il-dl+1;

  b = (short *)calloc(l+2,sizeof(short));
  c = (short *)calloc(l+2,sizeof(short));

  lass(b,bb,nn);
  lass(c,cc,nn);

  basint=10000;
  r=0;

  for (i=ib; i<=l; i++) {
     q=b[i]-c[i]-r;
     if (q < 0) {
        a[i]=basint+q;
	r=1;
     }
     else {
        a[i]=q;
	r=0;
     } /* if */
  } /* for i */

  free(b);
  free(c);

} /* sub */


/**********************************************************************************/
/*                                                                                */
/* labscomp: Compares the absolute values of the a and b lfloats operands.        */
/* Returns 1 if a >= b and -1 if a < b.                                           */
/*                                                                                */
/**********************************************************************************/

int labscomp(a,b,nn)
short *a;
short *b;
int *nn;
{
  int i,l,il,dl,ib;

  l=nn[0]; il=nn[1]; dl=nn[2];

  ib=l-il-dl+1;
  i=l;

  while ((a[i]==b[i]) && (i>=ib)) {
     i--;
  } /* while */

  if (a[i]>=b[i]) {
     return(1);
  }
  else {
     return(-1);
  } /* if */

} /* lasbcomp */


/**********************************************************************************/
/*                                                                                */
/* ladd: a = b + c; for lfloat operands. a, b and c can be the same variable.     */
/*                                                                                */
/**********************************************************************************/

void ladd(a,b,c,nn)
short *a;
short *b;
short *c;
int *nn;
{
   int b2,c2,co,l;

   l=nn[0];

   co=labscomp(b,c,nn);
  
   b2=b[l+1];
   c2=c[l+1]; 

      if ((b2>0) && (c2>0)) {
         add(a,b,c,nn);
         a[l+1]=1;
      } /* if */
      if ((b2>0) && (c2<0) && (co==1)) {
         sub(a,b,c,nn);
         a[l+1]=1;
      } /* if */
      if ((b2>0) && (c2<0) && (co==-1)) {
         sub(a,c,b,nn);
         a[l+1]=-1;
      } /* if */
      if ((b2<0) && (c2>0) && (co==1)) {
         sub(a,b,c,nn);
         a[l+1]=-1;
      } /* if */
      if ((b2<0) && (c2>0) && (co==-1)) {
         sub(a,c,b,nn);
         a[l+1]=1;
      } /* if */
      if ((b2<0) && (c2<0)) {
         add(a,b,c,nn);
         a[l+1]=-1;
      } /* if */

} /* ladd */


/**********************************************************************************/
/*                                                                                */
/* lsub: a = b - c for; lfloat operands. a, b and c can be the same variable.     */
/*                                                                                */
/**********************************************************************************/

void lsub(a,b,c,nn)
short *a;
short *b;
short *c;
int *nn;
{
   int b2,c2,co,l;

   l=nn[0];

   co=labscomp(b,c,nn);
  
   b2=b[l+1];
   c2=c[l+1]; 

      if ((b2>0) && (c2>0) && (co==1)) {
         sub(a,b,c,nn);
         a[l+1]=1;
      } /* if */
      if ((b2>0) && (c2>0) && (co==-1)) {
         sub(a,c,b,nn);
         a[l+1]=-1;
      } /* if */
      if ((b2>0) && (c2<0)) {
         add(a,b,c,nn);
         a[l+1]=1;
      } /* if */
      if ((b2<0) && (c2>0)) {
         add(a,b,c,nn);
         a[l+1]=-1;
      } /* if */
      if ((b2<0) && (c2<0) && (co==1)) {
         sub(a,c,b,nn);
         a[l+1]=-1;
      } /* if */
      if ((b2<0) && (c2<0) && (co==-1)) {
         sub(a,b,c,nn);
         a[l+1]=1;
      } /* if */

} /* lsub */


/**********************************************************************************/
/*                                                                                */
/* lcmul: a = bb * c*word^exp0; Here a and bb are lfloat operands and             */
/* c is an integer 32-bit variable.                                               */
/*                                                                                */
/**********************************************************************************/


void lcmul(a,bb,c,exp0,nn)
short *a;
short *bb;
int c;
int exp0;
int *nn;
{
   int i,j,k,l,il,dl,ib;
   double q,r,cr;

   short *b;

   l=nn[0]; il=nn[1]; dl=nn[2];
   ib=l-il-dl+1;

   b = (short *)calloc(l+2,sizeof(short));

   lass(b,bb,nn);

   k=ib+exp0;

   for (i=ib; i<k; i++) {
      a[i]=0;
   } /* for i */

   cr=(double)(abs(c));
   r=0.0;

   for (j=ib; j<=l; j++) {
      q=(double)(b[j])*cr+r+0.5;
      r=(double)((int)(basi*q));
      if (k>=ib) {
         a[k]=(int)(q-bas*r);
      } /* if */
      k++;
   } /* for j */

   a[k]=(int)(r+0.5);

   k++;
   for (i=k; i<=l; i++) {
      a[i]=0;
   } /* for i */

   if (c>=0) {
      a[l+1]=b[l+1];
   }
   else {
      a[l+1]=-b[l+1];
   } /* if */

   free(b);

} /* lcmul */


/**********************************************************************************/
/*                                                                                */
/* lcmul: a = bb * c*word^exp0; Here a and bb are lfloat operands and             */
/* c is a double floating point variable.                                         */
/*                                                                                */
/**********************************************************************************/

void lcmul2(a,bb,c,exp0,nn)
short *a;
short *bb;
double c;
int exp0;
int *nn;
{
   int i,j,k,l,il,dl,ib;
   double q,r,cr;

   short *b;

   l=nn[0]; il=nn[1]; dl=nn[2];
   ib=l-il-dl+1;

   b = (short *)calloc(l+2,sizeof(short));

   lass(b,bb,nn);

   k=ib+exp0;

   for (i=ib; i<k; i++) {
      a[i]=0;
   } /* for i */

   cr=fabs(c);
   r=0.0;

   for (j=ib; j<=l; j++) {
      q=(double)(b[j])*cr+r+0.5;
      r=floor(basi*q);
      if (k>=ib) {
         a[k]=(int)(q-bas*r);
      } /* if */
      k++;
   } /* for j */

   a[k]=(int)(r+0.5);

   k++;
   for (i=k; i<=l; i++) {
      a[i]=0;
   } /* for i */

   if (c>=0.0) {
      a[l+1]=b[l+1];
   }
   else {
      a[l+1]=-b[l+1];
   } /* if */

   free(b);

} /* lcmul2 */



/**********************************************************************************/
/*                                                                                */
/* ltoreal: man0*10^exp0 ~ a where man0 and exp0 are double floating point        */
/* varibles and a a lfloat operand. man0 is thus highly truncated to fit in       */
/* a double variable.                                                             */
/*                                                                                */
/**********************************************************************************/


void ltoreal(a,man0,exp0,nn)
short *a;
double *man0,*exp0;
int *nn;
{
   char s0[20];
   char s1[10];
   char s2[10];
   char s3[10];
   char s4[10];

   char st2[10];
   char st3[10];
   char st4[10];

   int i,l,il,dl,ib,length,le;

   double man,ex;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   i=l;

   while ((a[i]==0) && (i>=ib)) {
      i--;
   } /* while */

   if ((i==ib) && (a[i]==0)) {
      man=0.0;
      ex=0.0;
   }
   else {
      sprintf(s1,"%d",a[i]);
      length=strlen(s1);
      ex=(double)(4*(i-(l-il+1))+length);

      sprintf(s2,"%d",a[i-1]);
      sprintf(s3,"%d",a[i-2]);
      sprintf(s4,"%d",a[i-3]);
 
      strcpy(st2,"");
      le=strlen(s2);
      for (i=1; i<=(4-le); i++) {
         strcat(st2,"0");
      } /* for i */
      strcat(st2,s2);

      strcpy(st3,"");
      le=strlen(s3);
      for (i=1; i<=(4-le); i++) {
         strcat(st3,"0");
      } /* for i */
      strcat(st3,s3);

      strcpy(st4,"");
      le=strlen(s4);
      for (i=1; i<=(4-le);i++) {
         strcat(st4,"0");
      } /* for i */
      strcat(st4,s4);

      sprintf(s0,"%s%s%s%s%s","0.",s1,st2,st3,st4);
      man=atof(s0);
   } /* if */
   
   *man0=(double)(a[l+1])*man;
   *exp0=ex;

} /* ltoreal */


/**********************************************************************************/
/*                                                                                */
/* lrealtol: a = man0*10^exp0. Here a is a lfloat operand and man0 and exp0       */
/* are double floating point variables. The rest of the a array is set zero.      */
/*                                                                                */
/**********************************************************************************/

void lrealtol(a,man0,exp0,nn)
short *a;
double man0,exp0;
int *nn;
{
    char s0[18];
    char s1[16];
 
    char z1[7];
    char z2[7];
    char z3[7];
    char z4[7];

    double man,ex,exp2,man2;

    int i,k,l,il,dl,ib,l1,l2,ll;

    l=nn[0]; il=nn[1]; dl=nn[2];

    ib=l-il-dl+1;

    for (i=l; i>=ib; i--) {
       a[i]=0;
    } /* for i */   

    man2=man0;
    ex=exp0;
    man=fabs(man2);

    while (man>=1.0) {
       man=man/10.0;
       ex=ex+1.0;
    } /* while */

    sprintf(s0,"%17.16f",man);
 
    for (i=2; i<=17; i++) {
       s1[i-2]=s0[i];
    } /* for i */


    exp2=fabs(ex/(double)(4));
    l1=(int)(ex/(double)(4));
    l2=(int)(4*(exp2-(double)(abs(l1))));
    
    if (ex<0) {
       l2=4-l2;
       l1=l1-1;
    } /* if */

    if (l2==0) {
       l2=4;
       l1=l1-1; 
    } /* if */

    for (i=0; i<=l2-1; i++) {
       z1[i]=s1[i];
    } /* for i */

    z1[i]='\0';

    for (i=0; i<=4-1; i++) {
       z2[i]=s1[l2+i];
    } /* for i */

    z2[i]='\0';

    k=0;

    for (i=0; i<=4-1; i++) {
       if ((l2+i+4)<=16) {
          z3[i]=s1[l2+i+4];
          k++;
       } /* if */
    } /* for i */
 
    z3[k]='\0'; 

    for (i=1; i<=4-k; i++) {
       strcat(z3,"0");
    } /* for i */

    k=0;

    for (i=0; i<=4-1; i++) {
       if ((l2+i+2*4)<=16) {
          z4[i]=s1[l2+i+2*4];
          k++;
       } /* if */
    } /* for i */
 
    z4[k]='\0'; 

    for (i=1; i<=4-k; i++) {
       strcat(z4,"0");
    } /* for i */
    
    ll=l-il+l1+1;
 
    a[ll]=atoi(z1);
    a[ll-1]=atoi(z2);
    a[ll-2]=atoi(z3);
    a[ll-3]=atoi(z4);

    if (man2>=0) {
       a[l+1]=1;
    }
    else {
       a[l+1]=-1;
    } /* if */

} /* lrealtol */


/**********************************************************************************/
/*                                                                                */
/* lcompeps: Checks if the lfloat variable a is zero down to the nn[3] decimal    */
/* word calculated down from the "decimal point" defined by nn[0]-nn[1].          */
/* The function returns 1 if this is true and -1 if it is not.                    */
/*                                                                                */
/**********************************************************************************/

int lcompeps(a,nn)
short *a;
int *nn;
{
   int i,l,il,dl,k,ib,n,ll;

   l=nn[0]; il=nn[1]; dl=nn[2]; n=nn[3];

   ib=l-il-dl+1;

   k=0;

   ll=l-il-n;

   for (i=l; i>=ll; i--) {
      if (a[i]==0) {
         k++;
      } /* if */
   } /* for i */

   if (k==(il+n+1)) {
      return(-1);
   }
   else {
      return(1);
   } /* if */

}  /* lcompeps */



/**********************************************************************************/
/*                                                                                */
/* lgetdec: Returns the number of decimals (in base 10) the lfloat variable a     */
/* is zero below the decmal point.                                                */
/*                                                                                */
/**********************************************************************************/

int lgetdec(a,nn)
short *a;
int *nn;
{
   int i,l,il,dl,k,ib;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   k=0;

   i=l-il;

   while ((a[i]==0) && (i>ib)) {
      k++;
      i--;
   } /* while */

   return(4*k);

}  /* lgetdec */


/**********************************************************************************/
/*                                                                                */
/* lcheckdec: Checks how many words (from the most significant) that is identical */
/* of lfloat operand a.                                                           */
/*                                                                                */
/**********************************************************************************/

int lcheckdec(a,nn)
short *a;
int *nn;
{
   int i,l,il,dl,k,ib,ss;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   k=0;

   i=l-il;
   ss=a[i];

   while ((a[i]==ss) && (i>ib)) {
      k++;
      i--;
   } /* while */

   if ((i==ib) && (ss==0)) {
      k=0;
   } /* if */

   return(k);

}  /* lcheckdec */



/**********************************************************************************/
/*                                                                                */
/* lmul: a = b * c; for lfloat operands. a, b and c can be the same variable.     */
/* The result is truncated to the same size as b and c.                           */
/*                                                                                */
/**********************************************************************************/


void lmul(a,b,c,nn)
short *a;
short *b;
short *c;
int *nn;
{


   int i,k,l,il,dl,ib,n,n2;
   int bi,b1,c1;
   int bn,cn,kkk,bnib,cnib;

   double isize,rr,xx;
   double br,bri,maxerror,temp,qqq;

   double *q;


if (mulversion == 1) {

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   i=l;

   while ((b[i] == 0) && (i > ib)) {
      i--;
   } /* while */

   bnib=i;
   bn=bnib-ib;

   i=l;

   while ((c[i] == 0) && (i > ib)) {
      i--;
   } /* while */

   cnib=i;
   cn=cnib-ib;

   kkk=bn+cn;

   if (kkk > dl) {

   n=1;

   while (n < kkk) {
      n=2*n;
   } /* while */

   n2=2*n;

   q=(double *)malloc((n2+1)*sizeof(double));

   isize=1.0/sqrt((double)(n));

   for (i=0; i<=n2; i++) {
      q[i]=0.0;
   } /* for i */
 
   k=0;

   for (i=ib; i<=bnib; i++) {
      q[k]=(double)(b[i])*isize;
      k=k+2;
   } /* for i */

   k=0;

   for (i=ib; i<=cnib; i++) {
      q[k+1]=(double)(c[i])*isize;
      k=k+2;
   } /* for i */


   cdft(n2,1,q);

   mul(q,n);

   cdft(n2,-1,q);

   maxerror=0.0;

   rr=0.0;

   for (i=0; i<=n2; i=i+2) {
      qqq=q[i];
      temp=fabs(qqq-floor(qqq+0.5));
      if (temp > maxerror) {
         maxerror=temp;
      } /* if */
      xx=qqq+rr;
      rr=floor((xx+0.5)*basi);
      q[i]=xx-bas*rr; 
   } /* for i */

   if (maxerror > fftmaxerror) {
      fftmaxerror=maxerror;
   } /* if */

   k=2*dl;

   for (i=ib; i<=l; i++) {
      a[i]=0;
      if (k <= n2) {
         a[i]=(int)(q[k]+0.5);
      } /* if */
      k=k+2;
   } /* for i */

   free(q);

   }
   else {

   for (i=ib; i<=l; i++) {
      a[i]=0;
   } /* for i */

   } /* if */

   a[l+1]=b[l+1]*c[l+1]; 

}
else {

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   i=l;

   while ((b[i] == 0) && (i > ib)) {
      i--;
   } /* while */

   bnib=i;
   bn=bnib-ib;

   i=l;

   while ((c[i] == 0) && (i > ib)) {
      i--;
   } /* while */

   cnib=i;
   cn=cnib-ib;

   kkk=bn+cn;

   if (kkk > dl) {

   n=1;

   while (n < kkk) {
      n=2*n;
   } /* while */

   n=2*n;
   n2=2*n;

   q=(double *)malloc((n2+1)*sizeof(double));

   br=100.0;
   bri=0.01;
   bi=100;

   isize=1.0/sqrt((double)(n));

   for (i=0; i<=n2; i++) {
      q[i]=0.0;
   } /* for i */
 
   k=0;

   for (i=ib; i<=bnib; i++) {
      b1=b[i] % bi;
      q[k]=(double)(b1)*isize;
      b1=b[i] / bi;
      q[k+2]=(double)(b1)*isize;
      k=k+4;
   } /* for i */

   k=0;

   for (i=ib; i<=cnib; i++) {
      c1=c[i] % bi;
      q[k+1]=(double)(c1)*isize;
      c1=c[i] / bi;
      q[k+3]=(double)(c1)*isize;
      k=k+4;
   } /* for i */


   cdft(n2,1,q);

   mul(q,n);

   cdft(n2,-1,q);

   maxerror=0.0;

   rr=0.0;

   for (i=0; i<=n2; i=i+2) {
      qqq=q[i];
      temp=fabs(qqq-floor(qqq+0.5));
      if (temp > maxerror) {
         maxerror=temp;
      } /* if */
      xx=qqq+rr;
      rr=floor((xx+0.5)*bri);
      q[i]=xx-br*rr; 
   } /* for i */

   if (maxerror > fftmaxerror) {
      fftmaxerror=maxerror;
   } /* if */

   k=4*dl;

   for (i=ib; i<=l; i++) {
      a[i]=0;
      if (k <= n2) {
         a[i]=bi*(int)(q[k+2]+0.5)+(int)(q[k]+0.5);
      } /* if */
      k=k+4;
   } /* for i */

   free(q);

   }
   else {

   for (i=ib; i<=l; i++) {
      a[i]=0;
   } /* for i */

   } /* if */

   a[l+1]=b[l+1]*c[l+1]; 

} /* if */


} /* lmul */


/**********************************************************************************/
/*                                                                                */
/* linv: a = 1 / b; b != 0; for lfloat operands.                                  */
/*                                                                                */
/**********************************************************************************/

void linv(a,b,nn)
short *a;
short *b;
int *nn;
{

   double man,ex;
   
   int s,l,il,dl,dd,kk,bb,bbb,i;

   short *x1,*x2,*t1,*t2;

   int ss[4];
 
   l=nn[0]; il=nn[1]; dl=nn[2];
   dd=il+dl;

   x1=(short *)malloc((l+2)*sizeof(short));
   x2=(short *)malloc((l+2)*sizeof(short));
   t1=(short *)malloc((l+2)*sizeof(short));
   t2=(short *)malloc((l+2)*sizeof(short));

   lzero(x1,nn);
   lzero(x2,nn);
   lzero(t1,nn);
   lzero(t2,nn);
 
   s=b[l+1];
   b[l+1]=1;

   ltoreal(b,&man,&ex,nn);

   man=1.0/man;
   ex=-ex;

   lrealtol(x1,man,ex,nn);

   bbb=lcheckdec(b,nn);

   if (bbb < 16) {
      bb=32;
   }
   else {
      kk=(int)(log((double)(bbb+2))/log(2.0));
      bb=2;
      for (i=1; i<=kk; i++) {
         bb=2*bb;
      } /* for i */
   } /* if */

   kk=bb;

   ss[0]=l; ss[1]=il;

   do {
      if (kk > dd) {
         kk=dd;
      } /* if */
      ss[2]=kk-il;

      lmul(t1,x1,b,ss);
      linttol(t2,1,0,ss);
      lsub(x2,t2,t1,ss);

      lmul(t1,x2,x2,ss);
      lmul(a,x2,t1,ss);
      ladd(t2,t1,a,ss);
      ladd(t1,x2,t2,ss);
      lmul(t2,t1,x1,ss);
      ladd(t1,x1,t2,ss);
      lass(x1,t1,ss);

      kk=4*kk;
   } while (lcompeps(x2,nn)==1); 

   lass(a,x1,nn);

   a[l+1]=s;

   free(x1);
   free(x2);
   free(t1);
   free(t2);

} /* linv */ 

 
/**********************************************************************************/
/*                                                                                */
/* ldivv: a = b / c; c != 0; for lfloat operands. a, b and c can be the same      */
/* variable.                                                                      */
/*                                                                                */
/**********************************************************************************/
 
void ldivv(a,b,c,nn)
short *a;
short *b;
short *c;
int *nn;
{

   linv(a,c,nn);

   lmul(a,a,b,nn);

} /* ldivv */    

 
/**********************************************************************************/
/*                                                                                */
/* lsqrtinv: a = 1 / sqrt(b); b != 0; for lfloat operands.                        */
/*                                                                                */
/**********************************************************************************/

void lsqrtinv(a,b,nn)
short *a;
short *b;
int *nn;
{

   double re0,re1,re2,s10,man,ex,eps;

   int s,l,il,dl,dd,kk,bb,bbb,i;

   short *x1,*x2,*t1,*t2;
    
   int ss[4]; 
 
   l=nn[0]; il=nn[1]; dl=nn[2];
   dd=il+dl;

   x1=(short *)malloc((l+2)*sizeof(short));
   x2=(short *)malloc((l+2)*sizeof(short));
   t1=(short *)malloc((l+2)*sizeof(short));
   t2=(short *)malloc((l+2)*sizeof(short));

   lzero(x1,nn);
   lzero(x2,nn);
   lzero(t1,nn);
   lzero(t2,nn);

   eps=0.001; 
   s=b[l+1];
   b[l+1]=1;

   ltoreal(b,&man,&ex,nn);

   re0=1.0/sqrt(man);
   re1=-ex/2.0;
   re2=fabs(re1);

   s10=sqrt(10.0);

   if ((re2-(double)((int)(re2)))>eps) {
      re0=re0/s10;
      re1=re1+0.5;
   } /* if */ 

   man=re0;
   ex=re1;

   lrealtol(x1,man,ex,nn);

   bbb=lcheckdec(b,nn);

   if (bbb < 16) {
      bb=32;
   }
   else {
      kk=(int)(log((double)(bbb+2))/log(2.0));
      bb=2;
      for (i=1; i<=kk; i++) {
         bb=2*bb;
      } /* for i */
   } /* if */

   kk=bb;

   ss[0]=l; ss[1]=il; 

   do {
      if (kk > dd) {
         kk=dd;
      } /* if */ 
      ss[2]=kk-il;

      lmul(t1,x1,b,ss);
      lmul(t1,x1,t1,ss);
      linttol(t2,1,0,ss);
      lsub(x2,t2,t1,ss);

      lmul(t1,x2,x2,ss);
      lmul(a,t1,x2,ss);
      lcmul(t2,a,5,0,ss);
      lcmul(a,t1,6,0,ss);
      ladd(t1,a,t2,ss);
      lcmul(t2,x2,8,0,ss);
      ladd(a,t1,t2,ss);
      lmul(t1,a,x1,ss);
      lcmul(a,t1,625,-1,ss);
      ladd(t2,x1,a,ss);
      lass(x1,t2,ss);

      kk=4*kk;
   } while (lcompeps(x2,nn)==1); 

   lass(a,x1,nn);

   free(x1);
   free(x2);
   free(t1);
   free(t2);

} /* lsqrtinv */


/**********************************************************************************/
/*                                                                                */
/* lsqrtinv: a = sqrt(b); b != 0; for lfloat operands.                            */
/*                                                                                */
/**********************************************************************************/

void lsqrt(a,b,nn)
short *a;
short *b;
int *nn;
{

   lsqrtinv(a,b,nn);

   lmul(a,a,b,nn);

} /* lsqrt */


/**********************************************************************************/
/*                                                                                */
/* lsqrtsqrtinv: a = 1 / b^(1/4); b != 0; for lfloat operands.                    */
/*                                                                                */
/**********************************************************************************/


void lsqrtsqrtinv(a,b,nn)
short *a;
short *b;
int *nn;
{

   double re0,re1,re2,s10,man,ex,eps;

   int s,l,il,dl,bb,dd,kk,bbb,i;

   short *x1,*x2,*t1,*t2;

   int ss[4];
 
   l=nn[0]; il=nn[1]; dl=nn[2];
   dd=il+dl;

   x1=(short *)malloc((l+2)*sizeof(short));
   x2=(short *)malloc((l+2)*sizeof(short));
   t1=(short *)malloc((l+2)*sizeof(short));
   t2=(short *)malloc((l+2)*sizeof(short));

   lzero(x1,nn);
   lzero(x2,nn);
   lzero(t1,nn);
   lzero(t2,nn);
 
   eps=0.001;
   s=b[l+1];
   b[l+1]=1;

   ltoreal(b,&man,&ex,nn);

   re0=1.0/sqrt(sqrt(man));
   re1=-ex/4.0;
   re2=fabs(re1);

   s10=sqrt(sqrt(10.0));

   while ((re2-(double)((int)(re2)))>eps) {
      re0=re0/s10;
      re1=re1+0.25;
      re2=fabs(re1);
   } /* while */

   man=re0;
   ex=re1;

   lrealtol(x1,man,ex,nn);

   bbb=lcheckdec(b,nn);

   if (bbb < 16) {
      bb=32;
   }
   else {
      kk=(int)(log((double)(bbb+2))/log(2.0));
      bb=2;
      for (i=1; i<=kk; i++) {
         bb=2*bb;
      } /* for i */
   } /* if */

   kk=bb;

   ss[0]=l; ss[1]=il; 

   do {
      if (kk > dd) {
         kk=dd;
      } /* if */
      ss[2]=kk-il;

      lmul(t1,x1,x1,ss);
      lmul(t1,t1,t1,ss);
      lmul(t1,t1,b,ss);
      linttol(t2,1,0,ss);
      lsub(x2,t2,t1,ss);

      lmul(t1,x2,x2,ss);
      lmul(a,t1,x2,ss);
      lcmul(t2,a,15,0,ss);
      lcmul(a,t1,20,0,ss);
      ladd(t1,a,t2,ss);
      lcmul(t2,x2,32,0,ss);
      ladd(a,t1,t2,ss);
      lmul(t1,a,x1,ss);
      lcmul(a,t1,781250,-2,ss);
      ladd(t2,x1,a,ss);
      lass(x1,t2,ss);

      kk=4*kk;
   } while (lcompeps(x2,nn)==1); 

   lass(a,x1,nn);

   free(x1);
   free(x2);
   free(t1);
   free(t2);

} /* lsqrtsqrtinv */


/**********************************************************************************/
/*                                                                                */
/* lsqrtsqrt: a = b^(1/4); b != 0; for lfloat operands.                           */
/*                                                                                */
/**********************************************************************************/

void lsqrtsqrt(a,b,nn)
short *a;
short *b;
int *nn;
{

  int l,il,dl,dd;

   short *t1;

   l=nn[0]; il=nn[1]; dl=nn[2];
   dd=il+dl;

   lsqrtsqrtinv(a,b,nn);

   t1=(short *)malloc((l+2)*sizeof(short));

   lmul(t1,a,a,nn);
   lmul(t1,t1,a,nn);
   lmul(a,t1,b,nn);

   free(t1);

} /* lsqrtsqrt */



/**********************************************************************************/
/*                                                                                */
/* linv3: a = 1 / 3; for lfloat operand a. A special function not used for        */
/* the moment.                                                                    */
/*                                                                                */
/**********************************************************************************/

void linv3(a,nn)
short *a;
int *nn;
{
   int i,l,il,dl,ib,ll;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ll=l-il+1;
   ib=l-il-dl+1;

   for (i=l;i>=ll; i--) {
      a[i]=0;
   } /* for i */

   for (i=ll-1;i>=ib; i--) {
      a[i]=3333;
   } /* for i */

   a[l+1]=1;

} /* linv3 */


/**********************************************************************************/
/*                                                                                */
/* linv162: a = 1 / 162; for lfloat operand a. A special function used by         */
/* lcubertinv.                                                                    */
/*                                                                                */
/**********************************************************************************/



void linv162(a,nn)
short *a;
int *nn;
{
   int i,l,il,dl,ib,ll,kk;

   char cc[12];
   char ss[12];

   cc[0]='7';
   cc[1]='2';
   cc[2]='8';
   cc[3]='3';
   cc[4]='9';
   cc[5]='5';
   cc[6]='0';
   cc[7]='6';
   cc[8]='1';
   cc[9]='\0';

   l=nn[0]; il=nn[1]; dl=nn[2];

   ll=l-il+1;
   ib=l-il-dl+1;

   for (i=l;i>=ll; i--) {
      a[i]=0;
   } /* for i */

   a[ll-1]=61;

   kk=0;

   for (i=ll-2;i>=ib; i--) {

      ss[0]=cc[(kk % 9)];
      ss[1]=cc[((kk+1) % 9)];
      ss[2]=cc[((kk+2) % 9)];
      ss[3]=cc[((kk+3) % 9)];
      ss[4]='\0';

      a[i]=atoi(ss);

      kk=kk+4;

   } /* for i */

   a[l+1]=1;

} /* linv162 */


/**********************************************************************************/
/*                                                                                */
/* lcubertinv: a = 1 / b^(1/3); b != 0; for lfloat operands.                      */
/*                                                                                */
/**********************************************************************************/

void lcubertinv(a,b,nn)
short *a;
short *b;
int *nn;
{

   double re0,re1,re2,s10,man,ex,eps;

   int s,l,il,dl,bb,dd,kk,bbb,i;

   short *x1,*x2,*t1,*t2,*t3;

   int ss[4];
 
   l=nn[0]; il=nn[1]; dl=nn[2];
   dd=il+dl;

   x1=(short *)malloc((l+2)*sizeof(short));
   x2=(short *)malloc((l+2)*sizeof(short));
   t1=(short *)malloc((l+2)*sizeof(short));
   t2=(short *)malloc((l+2)*sizeof(short));
   t3=(short *)malloc((l+2)*sizeof(short));

   lzero(x1,nn);
   lzero(x2,nn);
   lzero(t1,nn);
   lzero(t2,nn);

   linv162(t3,nn);

   eps=0.001;
   s=b[l+1];
   b[l+1]=1;

   ltoreal(b,&man,&ex,nn);

   re0=1.0/exp(log(man)/3.0);
   re1=-ex/3.0;
   re2=fabs(re1);

   s10=exp(log(10.0)/3.0);

   while ((re2-(double)((int)(re2)))>eps) {
      re0=re0/s10;
      re1=re1+1.0/3.0;
      re2=fabs(re1);
   } /* while */

   man=re0;
   ex=re1;

   lrealtol(x1,man,ex,nn);

   bbb=lcheckdec(b,nn);

   if (bbb < 16) {
      bb=32;
   }
   else {
      kk=(int)(log((double)(bbb+2))/log(2.0));
      bb=2;
      for (i=1; i<=kk; i++) {
         bb=2*bb;
      } /* for i */
   } /* if */

   kk=bb;

   ss[0]=l; ss[1]=il; 

   do {
      if (kk > dd) {
         kk=dd;
      } /* if */
      ss[2]=kk-il;

      lmul(t1,x1,x1,ss);
      lmul(t1,t1,x1,ss);
      lmul(t1,t1,b,ss);
      linttol(t2,1,0,ss);
      lsub(x2,t2,t1,ss);

      lmul(t1,x2,x2,ss);
      lmul(a,t1,x2,ss);
      lcmul(t2,a,28,0,ss);
      lcmul(a,t1,36,0,ss);
      ladd(t1,a,t2,ss);
      lcmul(t2,x2,54,0,ss);
      ladd(a,t1,t2,ss);
      lmul(t1,a,x1,ss);
      lmul(a,t3,t1,ss);
      ladd(t2,x1,a,ss);
      lass(x1,t2,ss);

      kk=4*kk;
   } while (lcompeps(x2,nn)==1); 

   lass(a,x1,nn);

   free(x1);
   free(x2);
   free(t1);
   free(t2);
   free(t3);

} /* lcubertinv */


/**********************************************************************************/
/*                                                                                */
/* lcubert: a = b^(1/3); b != 0; for lfloat operands.                             */
/*                                                                                */
/**********************************************************************************/

void lcubert(a,b,nn)
short *a;
short *b;
int *nn;
{

  int l,il,dl,dd;

   short *t1;

   l=nn[0]; il=nn[1]; dl=nn[2];
   dd=il+dl;

   lcubertinv(a,b,nn);

   t1=(short *)malloc((l+2)*sizeof(short));

   lmul(t1,a,a,nn);
   lmul(a,t1,b,nn);

   free(t1);

} /* lcubert */


void biginit(a,maxsize)
bigint *a;
int maxsize;
{
   a->coef=(short *)malloc((maxsize+1)*sizeof(short));
   a->size=1;
   a->maxsize=maxsize;
   a->exp0=0;
   a->sgn0=1;

   a->coef[1]=0;
} /* biginit */


void biginit2(a,maxsize)
bigint *a;
int maxsize;
{
   int i;

   a->coef=(short *)malloc((maxsize+1)*sizeof(short));
   a->size=1;
   a->maxsize=maxsize;
   a->exp0=0;
   a->sgn0=1;

   for (i=1; i<=maxsize; i++) {
      a->coef[i]=0;
   } /* for i */

} /* biginit2 */


void bigfree(a)
bigint *a;
{
  free(a->coef);
} /* bigfree */


void inttobig(a,bb)
bigint *a;
int bb;
{
   int tt1,tt2,t1,t2,b;

   b=abs(bb);

   if (b < 10000) {
      a->coef[1]=b;
      a->size=1;
   }
   else {
      tt1=b % 100000000;
      tt2=b / 100000000;
      t1=tt1 % 10000;
      t2=tt1 / 10000;
      a->coef[1]=t1;
      a->coef[2]=t2;
      a->coef[3]=tt2;
      a->size=1;
      if (t2 > 0) {
         a->size=2;
      } /* if */
      if (tt2 > 0) {
         a->size=3;
      } /* if */
   } /* if */

   a->exp0=0;

   if (bb >= 0) {
      a->sgn0=1;
   }
   else {
      a->sgn0=-1;
   } /* if */

} /* inttobig */


void bigass(a,b)
bigint *a;
bigint *b;
{
   int i,n;

   n=b->size;

   for (i=1; i<=n; i++) {
      a->coef[i]=b->coef[i];
   } /* for i */

   a->size=n;
   a->exp0=b->exp0;
   a->sgn0=b->sgn0;

} /* bigass */


int bigcomp2(b,c)
bigint *b;
bigint *c;
{

  int i,n,bn,cn,bb,cc,maxsize,maxs;
  int bexp,cexp,expmax;
  int bk,ck,bi,ci,flag;

  bn=b->size;
  cn=c->size;

  maxsize=b->maxsize;
  maxs=maxsize-2;

  bexp=b->exp0;
  cexp=c->exp0;

  if ((bn+bexp) > (cn+cexp)) {
     expmax=bexp;
  }
  else {
     expmax=cexp;
  } /* if */

  if (bn > cn) {
     n=bn;
  }
  else {
     n=cn;
  } /* if */

  bi=expmax-bexp;
  ci=expmax-cexp;

  i=n;
  flag=0;

  while ((flag == 0) && (i >= 1)) {

     bk=i+bi;
     if ((bk <= bn) && (bk >= 1)) {
        bb=b->coef[bk];
     }
     else {
        bb=0;
     } /* if */

     ck=i+ci;
     if ((ck <= cn) && (ck >= 1)) {
        cc=c->coef[ck];
     } 
     else {
        cc=0;
     } /* if */

     if (bb != cc) {
        flag=1;
     } /* if */

     i--;

  } /* while*/

  if (bb >= cc) {
     return(1);
  }
  else {
     return(-1);
  } /* if */


} /* bigcomp2 */
      

void bigadd2(a,bbb,ccc)
bigint *a;
bigint *bbb;
bigint *ccc;
{

  int i,n,bn,cn,bb,cc,maxsize,maxs;
  int q,r,basint,bexp,cexp,expmax;
  int bk,ck,bi,ci;

  bigint b,c;

  bn=bbb->size;
  cn=ccc->size;

  maxsize=bbb->maxsize;
  maxs=maxsize-2;

  biginit(&b,maxsize);
  biginit(&c,maxsize);

  bigass(&b,bbb);
  bigass(&c,ccc);

  bexp=bbb->exp0;
  cexp=ccc->exp0;

  if ((bn+bexp) > (cn+cexp)) {
     expmax=bexp;
  }
  else {
     expmax=cexp;
  } /* if */

  if (bn > cn) {
     n=bn;
  }
  else {
     n=cn;
  } /* if */

  bi=expmax-bexp;
  ci=expmax-cexp;

  basint=10000;
  r=0;

  for (i=1; i<=n; i++) {
     bk=i+bi;
     if ((bk <= bn) && (bk >= 1)) {
        bb=b.coef[bk];
     }
     else {
        bb=0;
     } /* if */

     ck=i+ci;
     if ((ck <= cn) && (ck >= 1)) {
        cc=c.coef[ck];
     } 
     else {
        cc=0;
     } /* if */

     q=bb+cc+r;
     if (q >= basint) {
        a->coef[i]=q-basint;
	r=1;
     }
     else {
        a->coef[i]=q;
	r=0;
     } /* if */
  } /* for i */


  if (n >= maxs) {

     if (r > 0) {

        for (i=1; i < n; i++) {
           a->coef[i]=a->coef[i+1];
        } /* for i */

        a->coef[n]=r;
        a->size=n;
        a->exp0=expmax+1;
     }
     else {
        a->size=n;
        a->exp0=expmax;
     } /* if */

  }
  else {
  
     if (r > 0) {
        a->coef[n+1]=r;
        a->size=n+1;
        a->exp0=expmax;
     }
     else {
        a->size=n;
        a->exp0=expmax;
     } /* if */

  } /* if */

  bigfree(&b);
  bigfree(&c);

} /* bigadd2 */
      

void bigsub2(a,bbb,ccc)
bigint *a;
bigint *bbb;
bigint *ccc;
{

  int i,n,bn,cn,bb,cc,maxsize,maxs;
  int q,r,basint,bexp,cexp,expmax;
  int bk,ck,bi,ci;

  bigint b,c;

  bn=bbb->size;
  cn=ccc->size;

  maxsize=bbb->maxsize;
  maxs=maxsize-2;

  biginit(&b,maxsize);
  biginit(&c,maxsize);

  bigass(&b,bbb);
  bigass(&c,ccc);

  bexp=bbb->exp0;
  cexp=ccc->exp0;

  if ((bn+bexp) > (cn+cexp)) {
     expmax=bexp;
  }
  else {
     expmax=cexp;
  } /* if */

  if (bn > cn) {
     n=bn;
  }
  else {
     n=cn;
  } /* if */

  bi=expmax-bexp;
  ci=expmax-cexp;

  basint=10000;
  r=0;

  for (i=1; i<=n; i++) {
     bk=i+bi;
     if ((bk <= bn) && (bk >= 1)) {
        bb=b.coef[bk];
     }
     else {
        bb=0;
     } /* if */

     ck=i+ci;
     if ((ck <= cn) && (ck >= 1)) {
        cc=c.coef[ck];
     } 
     else {
        cc=0;
     } /* if */

     q=bb-cc-r;
     if (q < 0) {
        a->coef[i]=basint+q;
	r=1;
     }
     else {
        a->coef[i]=q;
	r=0;
     } /* if */

  } /* for i */

  i=n;

  while ((a->coef[i] == 0) && (i > 1)) {
     i--;
  } /* while */

  a->size=i;
  a->exp0=expmax;

  bigfree(&b);
  bigfree(&c);

} /* bigsub2 */


void bigadd(a,b,c)
bigint *a;
bigint *b;
bigint *c;
{
   int b2,c2,co;

   co=bigcomp2(b,c);
  
   b2=b->sgn0;
   c2=c->sgn0; 

      if ((b2>0) && (c2>0)) {
         bigadd2(a,b,c);
	 a->sgn0=1;
      } /* if */
      if ((b2>0) && (c2<0) && (co==1)) {
         bigsub2(a,b,c);
	 a->sgn0=1;
      } /* if */
      if ((b2>0) && (c2<0) && (co==-1)) {
         bigsub2(a,c,b);
         a->sgn0=-1;
      } /* if */
      if ((b2<0) && (c2>0) && (co==1)) {
         bigsub2(a,b,c);
         a->sgn0=-1;
      } /* if */
      if ((b2<0) && (c2>0) && (co==-1)) {
         bigsub2(a,c,b);
	 a->sgn0=1;
      } /* if */
      if ((b2<0) && (c2<0)) {
         bigadd2(a,b,c);
	 a->sgn0=-1;
      } /* if */

} /* bigadd */


void bigsub(a,b,c,nn)
bigint *a;
bigint *b;
bigint *c;
int *nn;
{
   int b2,c2,co;

   co=bigcomp2(b,c);
  
   b2=b->sgn0;
   c2=c->sgn0; 

      if ((b2>0) && (c2>0) && (co==1)) {
         bigsub2(a,b,c);
         a->sgn0=1;
      } /* if */
      if ((b2>0) && (c2>0) && (co==-1)) {
         bigsub2(a,c,b);
         a->sgn0=-1;
      } /* if */
      if ((b2>0) && (c2<0)) {
         bigadd2(a,b,c);
         a->sgn0=1;
      } /* if */
      if ((b2<0) && (c2>0)) {
         bigadd2(a,b,c);
         a->sgn0=-1;
      } /* if */
      if ((b2<0) && (c2<0) && (co==1)) {
         bigsub2(a,c,b);
         a->sgn0=-1;
      } /* if */
      if ((b2<0) && (c2<0) && (co==-1)) {
         bigsub2(a,b,c);
         a->sgn0=1;
      } /* if */

} /* bigsub */


void bigmul(a,b,c)
bigint *a;
bigint *b;
bigint *c;
{

   double *q;

   int i,k,n,n2;
   int bi,b1,c1;
   int bn,cn,kkk,maxsize;
   int bexp,cexp,exptot,gg,maxs;

   double isize,rr,xx;
   double br,bri,maxerror,qqq,temp;


if (mulversion == 1) {

   bn=b->size;
   cn=c->size;

   bexp=b->exp0;
   cexp=c->exp0;
   exptot=bexp+cexp;

   maxsize=b->maxsize;
   maxs=maxsize-2;

   kkk=bn+cn;

   n=1;

   while (n < kkk) {
      n=2*n;
   } /* while */

   n2=2*n;

   q=(double *)malloc((n2+1)*sizeof(double));

   for (i=0; i<=n2; i++) {
      q[i]=0.0;
   } /* for i */

   isize=1.0/sqrt((double)(n));
 
   k=0;

   for (i=1; i<=bn; i++) {
      q[k]=(double)(b->coef[i])*isize;
      k=k+2;
   } /* for i */

   k=0;

   for (i=1; i<=cn; i++) {
      q[k+1]=(double)(c->coef[i])*isize;
      k=k+2;
   } /* for i */

   cdft(n2,1,q);

   mul(q,n);

   cdft(n2,-1,q);

   maxerror=0.0;

   rr=0.0;

   for (i=0; i<=n2; i=i+2) {
      qqq=q[i];
      temp=fabs(qqq-floor(qqq+0.5));
      if (temp > maxerror) {
         maxerror=temp;
      } /* if */
      xx=qqq+rr;
      rr=floor((xx+0.5)*basi);
      q[i]=xx-bas*rr; 
   } /* for i */

   if (maxerror > fftmaxerror) {
      fftmaxerror=maxerror;
   } /* if */

   gg=1;
   if (kkk > maxs) {
   /*
      printf("* %8d %8d\n",kkk,maxs);
      */
      gg=kkk-maxs+1;
   } /* if */

   k=2*gg-2;

   for (i=gg; i<=kkk; i++) {
      a->coef[i-gg+1]=(int)(q[k]+0.5);
      k=k+2;
   } /* for i */


   if (a->coef[kkk-gg+1] > 0) {
      a->size=kkk-gg+1;
      a->exp0=exptot+gg-1;
   }
   else {
      a->size=kkk-gg;
      a->exp0=exptot+gg-1;
   } /* if */

   a->sgn0=(b->sgn0)*(c->sgn0);


   free(q);


} 
else {

   br=100.0;
   bri=0.01;
   bi=100;

   bn=b->size;
   cn=c->size;

   bexp=b->exp0;
   cexp=c->exp0;
   exptot=bexp+cexp;

   maxsize=b->maxsize;
   maxs=maxsize-2;

   kkk=bn+cn;

   n=1;

   while (n < kkk) {
      n=2*n;
   } /* while */

   n=2*n;
   n2=2*n;

   q=(double *)malloc((n2+1)*sizeof(double));

   for (i=0; i<=n2; i++) {
      q[i]=0.0;
   } /* for i */

   isize=1.0/sqrt((double)(n));
 
   k=0;

   for (i=1; i<=bn; i++) {
      b1=b->coef[i] % bi;
      q[k]=(double)(b1)*isize;
      b1=b->coef[i] / bi;
      q[k+2]=(double)(b1)*isize;
      k=k+4;
   } /* for i */

   k=0;

   for (i=1; i<=cn; i++) {
      c1=c->coef[i] % bi;
      q[k+1]=(double)(c1)*isize;
      c1=c->coef[i] / bi;
      q[k+3]=(double)(c1)*isize;
      k=k+4;
   } /* for i */

   cdft(n2,1,q);

   mul(q,n);

   cdft(n2,-1,q);

   maxerror=0.0;

   rr=0.0;

   for (i=0; i<=n2; i=i+2) {
      qqq=q[i];
      temp=fabs(qqq-floor(qqq+0.5));
      if (temp > maxerror) {
         maxerror=temp;
      } /* if */
      xx=qqq+rr;
      rr=floor((xx+0.5)*bri);
      q[i]=xx-br*rr; 
   } /* for i */

   if (maxerror > fftmaxerror) {
      fftmaxerror=maxerror;
   } /* if */

   gg=1;
   if (kkk > maxs) {
   /*
      printf("* %8d %8d\n",kkk,maxs);
      */
      gg=kkk-maxs+1;
   } /* if */

   k=4*gg-4;

   for (i=gg; i<=kkk; i++) {
      a->coef[i-gg+1]=bi*(int)(q[k+2]+0.5)+(int)(q[k]+0.5);
      k=k+4;
   } /* for i */

   if (a->coef[kkk-gg+1] > 0) {
      a->size=kkk-gg+1;
      a->exp0=exptot+gg-1;
   }
   else {
      a->size=kkk-gg;
      a->exp0=exptot+gg-1;
   } /* if */

   a->sgn0=(b->sgn0)*(c->sgn0);

   free(q);

} /* if */

} /* bigmul */


void big2l(a_big,b_big,a_l,b_l,nn)
bigint *a_big;
bigint *b_big;
short *a_l;
short *b_l;
int *nn;
{
   int i,ii,l,il,dl,ib,k,kk,jj,maxexp;

   l=nn[0]; il=nn[1]; dl=nn[2];

   ib=l-il-dl+1;

   if (b_big->exp0 > a_big->exp0) {
      maxexp=b_big->exp0;
   } 
   else {
      maxexp=a_big->exp0;
   } /* if */

   if (b_big->size > a_big->size) {
      ii=b_big->size;
   }
   else {
      ii=a_big->size;
   } /* if */

   k=l-il;

   for (i=ii; i >= 1; i--) {
      if (k >= 1) {
         jj=i+maxexp-b_big->exp0;
         if (jj <= b_big->size) {
            b_l[k]=b_big->coef[jj];
	 } /* if */
         kk=i+maxexp-a_big->exp0;
         if (kk <= a_big->size) {
            a_l[k]=a_big->coef[kk];
	 } /* if */
      } /* if */
      k--;
   } /* for i */

} /* big2l */


double julian_day(year,month,day)
int year,month,day;
{
   double f,g,a,jdd;

   if (month>=3) {
      f=(double)(year);
      g=(double)(month);
   }
   else {
      f=(double)(year-1);
      g=(double)(month+12);
   } /* if */

   a=2.0-floor(f/100.0)+floor(f/400.0);

   jdd=floor(365.25*f)+floor(30.6001*(g+1.0))+
       (double)(day)+a+1720994.5;

   return(jdd);

} /* julian_day */


double get_float_days(ttime)
time_t *ttime;
{
   int year,month,day,hour,min,sec;

   double fdays,fhours;

   struct tm *tt;

   tt=localtime(ttime);

   year=tt->tm_year+1900;
   month=tt->tm_mon+1;
   day=tt->tm_mday;

   hour=tt->tm_hour;
   min=tt->tm_min;
   sec=tt->tm_sec;
  
   fhours=(double)(hour)+(double)(min)/60.0+(double)(sec)/3600.0;

   fdays=julian_day(year,month,day)+fhours/24.0;

   return(fdays);

} /* get_float_days */


void day_and_time(ttime)
time_t *ttime;
{

   int year,month,day,hour,min,sec;

   struct tm *tt;

   tt=localtime(ttime);

   year=tt->tm_year+1900;
   month=tt->tm_mon+1;
   day=tt->tm_mday;

   hour=tt->tm_hour;
   min=tt->tm_min;
   sec=tt->tm_sec;

   printf("%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,min,sec);

} /* day_and_time */


void fday_and_time(outfile,ttime)
FILE *outfile;
time_t *ttime;
{

   int year,month,day,hour,min,sec;

   struct tm *tt;

   tt=localtime(ttime);

   year=tt->tm_year+1900;
   month=tt->tm_mon+1;
   day=tt->tm_mday;

   hour=tt->tm_hour;
   min=tt->tm_min;
   sec=tt->tm_sec;

   fprintf(outfile,"%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,min,sec);

} /* fday_and_time */


void sum_pi_chudnovsky(rP,rQ,rT,n1,n2,maxsize)
bigint *rP;
bigint *rQ;
bigint *rT;
int n1;
int n2;
int maxsize;
{

   int nm,n;

   bigint LP,LQ,LT;
   bigint RP,RQ,RT;

   if ((n2-n1) == 1) {

      n=n1;

      if (n == 0) {

         inttobig(rP,1);
         inttobig(rQ,1);
         inttobig(rT,13591409);

      } 
      else {


         inttobig(&tmp,(6*n-5));
         inttobig(&tmp2,(2*n-1));
	 bigmul(&tmp,&tmp,&tmp2);
         inttobig(&tmp2,(6*n-1));
	 bigmul(rP,&tmp,&tmp2);
	 rP->sgn0=-1;

         inttobig(&tmp,545140134);
         inttobig(&tmp2,n);
         bigmul(&tmp,&tmp,&tmp2);
         inttobig(&tmp2,13591409);
         bigadd(&tmp3,&tmp,&tmp2);

	 bigmul(rT,rP,&tmp3);

	 inttobig(&tmp,26680);
         inttobig(&tmp2,640320);
	 bigmul(&tmp,&tmp2,&tmp);
         inttobig(&tmp2,640320);
	 bigmul(&tmp,&tmp2,&tmp);
         inttobig(&tmp2,n);
	 bigmul(&tmp,&tmp2,&tmp);
         inttobig(&tmp2,n);
	 bigmul(&tmp,&tmp2,&tmp);
         inttobig(&tmp2,n);
	 bigmul(rQ,&tmp2,&tmp);

     } /* if */


   }
   else {

      biginit(&LP,maxsize);
      biginit(&LQ,maxsize);
      biginit(&LT,maxsize);

      biginit(&RP,maxsize);
      biginit(&RQ,maxsize);
      biginit(&RT,maxsize);

      nm=(n1+n2)/2;

      sum_pi_chudnovsky(&LP,&LQ,&LT,n1,nm,maxsize);

      sum_pi_chudnovsky(&RP,&RQ,&RT,nm,n2,maxsize);


      bigmul(rP,&LP,&RP);
      bigmul(rQ,&LQ,&RQ);

      bigmul(&tmp,&RQ,&LT);

      bigmul(&tmp2,&LP,&RT);

      bigadd(rT,&tmp,&tmp2);

      bigfree(&LP);
      bigfree(&LQ);
      bigfree(&LT);
      
      bigfree(&RP);
      bigfree(&RQ);
      bigfree(&RT);


   } /* if */

} /* sum_pi_chudnovsky */


void sum_pi_ramanujan(rP,rQ,rT,n1,n2,maxsize)
bigint *rP;
bigint *rQ;
bigint *rT;
int n1;
int n2;
int maxsize;
{

   int nm,n;

   bigint LP,LQ,LT;
   bigint RP,RQ,RT;

   if ((n2-n1) == 1) {

      n=n1;

      if (n == 0) {

         inttobig(rP,1);
         inttobig(rQ,1);
         inttobig(rT,1103);

      } 
      else {


         inttobig(&tmp,(4*n-3));
         inttobig(&tmp2,(4*n-2));
	 bigmul(&tmp,&tmp,&tmp2);
         inttobig(&tmp2,(4*n-1));
	 bigmul(&tmp,&tmp,&tmp2);
         inttobig(&tmp2,(4*n));
	 bigmul(rP,&tmp,&tmp2);

         inttobig(&tmp,26390);
         inttobig(&tmp2,n);
         bigmul(&tmp,&tmp,&tmp2);
         inttobig(&tmp2,1103);
         bigadd(&tmp3,&tmp,&tmp2);

	 bigmul(rT,rP,&tmp3);

	 inttobig(&tmp,156816);
	 bigmul(&tmp3,&tmp,&tmp);
         inttobig(&tmp2,n);
	 bigmul(&tmp2,&tmp2,&tmp2);
	 bigmul(&tmp2,&tmp2,&tmp2);
	 bigmul(rQ,&tmp2,&tmp3);

     } /* if */


   }
   else {

      biginit(&LP,maxsize);
      biginit(&LQ,maxsize);
      biginit(&LT,maxsize);

      biginit(&RP,maxsize);
      biginit(&RQ,maxsize);
      biginit(&RT,maxsize);

      nm=(n1+n2)/2;

      sum_pi_ramanujan(&LP,&LQ,&LT,n1,nm,maxsize);

      sum_pi_ramanujan(&RP,&RQ,&RT,nm,n2,maxsize);


      bigmul(rP,&LP,&RP);
      bigmul(rQ,&LQ,&RQ);

      bigmul(&tmp,&RQ,&LT);

      bigmul(&tmp2,&LP,&RT);

      bigadd(rT,&tmp,&tmp2);

      bigfree(&LP);
      bigfree(&LQ);
      bigfree(&LT);
      
      bigfree(&RP);
      bigfree(&RQ);
      bigfree(&RT);


   } /* if */

} /* sum_pi_ramanujan */


int main(argc,argv) 
int argc;
char *argv[];
{

   FILE *outfile,*mvfile;

   int bigsize,lsize,maxmem;
   int i,n,length,ilength,pot,pot0;
   int mem1,mem2,decimals,decimals2;
   int n1,n2,method;

   double elapseddays,frac,temp;

   int nn[5];

   clock_t time1,timediff1;
   clock_t time2,timediff2;

   time_t start_time,stop_time;

   short *a,*b,*c,*d,*ab,*pi_result; 

   bigint rP;
   bigint rQ;
   bigint rT;


   if ((mvfile = fopen(mv_file_name,"r")) == NULL) {
      printf("\n");
      printf("Warning, file %s not found, mulversion now set to 1!\n",mv_file_name);
      printf("\n");
      mulversion=1;
   } 
   else {
      fscanf(mvfile,"%d",&mulversion);
      fclose(mvfile);

      if ((mulversion != 1) && (mulversion != 2)) {
         printf("\n");
         printf("Warning, wrong value of mulversion in %s, mulversion now set to 1!\n",mv_file_name);
         printf("\n");
	 mulversion=1;
      } /* if */

   } /* if */


   if (argc != 3) {
      printf("\n");
      printf("Usage: %s methodi power-of-2i\n",argv[0]);
      printf("\n");
      printf("%s\n",version_string);
      printf("\n");
      printf("method = 1 => Pi using Chudnovsky fast series.\n");
      printf("method = 2 => Pi using Ramanujan fast series.\n");
      printf("\n");
      printf("Pi is calculated with about 2^(power-of-2) decimal digits.\n");
      printf("\n");
      printf("The variable mulversion is at present set to %1d.\n",mulversion);
      printf("It is possible to change the value (to %1d) by editing the file \"%s\".\n",
             3-mulversion,mv_file_name);
      printf("Set the variable to 2 for calculations with more than 2^24 (16 million) digits,\n");
      printf("or more important at least when \"FFT max error\" > 0.25.\n");
      printf("\n");
      exit(0);
   } /* if */

   method=atoi(argv[1]);
   pot0=atoi(argv[2]);

   if ((method < 1) || (method > 2)) {
      printf("\n");
      printf("Wrong method number.\n");
      printf("\n");
      exit(0);
   } /* if */

   if (pot0 < 8) {
      pot0=8;
   } /* if */


if (method == 1) {   


   printf("\n");
   printf("Pi, Chudnovsky\n");
   printf("============== \n");
   printf("\n");
   printf("%s\n",version_string);
   printf("\n");

   pot=pot0;

   n=1;

   for (i=1; i<=pot; i++) {
      n=2*n;
   } /* for i */

   decimals=n;
   decimals2=decimals;

   frac=0.995;
   decimals=decimals/4;
   decimals=(int)(frac*(double)(decimals-8))-2;
   decimals=4*decimals;

   n1=0;
   n2=decimals2/14+1;

   bigsize=n/4;

   length=n/4;
   ilength=8;
   frac=0.995;

   /* nn:    Short array with data for all lfloat variables. 0 <= word <= 9999.         */
   /* nn[0]: Total length of lfloat variable in words.                                  */
   /* nn[1]: Length of integer part of lfloat variable in words.                        */
   /* nn[2]: Length of decimal part of lfloat variable in words.                        */
   /* nn[3]: Fraction of the decimals that for sure will be correct during calulations. */
   /*        Here fraction is 0.995.                                                    */

   nn[0]=length;
   nn[1]=ilength;
   nn[2]=length-ilength;
   nn[3]=(int)(frac*(double)(length-ilength))-2;

   lsize=nn[0];

   printf("mulversion = %10d\n",mulversion);
   printf("Power of 2 = %10d\n",pot);
   printf("n          = %10d\n",n);
   printf("FFT        = %10d\n",2*mulversion*bigsize);
   printf("maxsize    = %10d\n",bigsize);
   printf("decimals   = %10d  (%d)\n",decimals,decimals2);
   printf("n1         = %10d\n",n1);
   printf("n2         = %10d\n",n2);
   printf("\n");

   printf("Number of calculated decimals: %10d\n",decimals);
   printf("\n");

   temp=(double)(bigsize);
   mem1=(int)(2.4*1048576.0*log(temp)*temp/(log(1048576.0/8.0)*1048576.0/8.0)+0.5);

   mem2=8*2*2*mulversion*bigsize;
   maxmem=mem1+mem2;
   maxmem=maxmem/1024;

   printf("Max approx memory usage: %10d kbyte;   mulversion: %3d\n",maxmem,mulversion);

   printf("\n");

   time(&start_time);
   printf("Start time: ");
   day_and_time(&start_time);

   time2=clock();

   printf("\n");
   printf("\n");

   biginit2(&rP,bigsize);
   biginit2(&rQ,bigsize);
   biginit2(&rT,bigsize);

   inttobig(&rP,0);
   inttobig(&rQ,0);
   inttobig(&rT,0);

   biginit2(&tmp,bigsize);
   biginit2(&tmp2,bigsize);
   biginit2(&tmp3,bigsize);
   biginit2(&tmp4,bigsize);

   inttobig(&tmp,0);
   inttobig(&tmp2,0);
   inttobig(&tmp3,0);
   inttobig(&tmp4,0);

   printf("Calculating Pi ...\n");

   fftmaxerror=0.0;

   time1=clock();

   sum_pi_chudnovsky(&rP,&rQ,&rT,n1,n2,bigsize);

   timediff1=clock()-time1;

   printf("\n");

   printf("Time Pi, fast series:         %10.2f\n",(double)(timediff1)/timetic);

   bigfree(&rP);

   bigfree(&tmp3);
   bigfree(&tmp4);

   a = (short *)calloc(lsize+2,sizeof(short));
   b = (short *)calloc(lsize+2,sizeof(short));
   c = (short *)calloc(lsize+2,sizeof(short));
   d = (short *)calloc(lsize+2,sizeof(short));
   ab = (short *)calloc(lsize+2,sizeof(short));
   pi_result = (short *)calloc(lsize+2,sizeof(short));

   inttobig(&tmp,12);

   bigmul(&tmp,&tmp,&rT);

   bigass(&tmp2,&rQ);

   bigfree(&rQ);
   bigfree(&rT);

   lzero(a,nn);
   lzero(b,nn);
   lzero(ab,nn);

   big2l(&tmp,&tmp2,a,b,nn);

   bigfree(&tmp);
   bigfree(&tmp2);

   time1=clock();
       
   ldivv(ab,b,a,nn);

   timediff1=clock()-time1;

   printf("Time div:                     %10.2f\n",(double)(timediff1)/timetic);

   linttol2(d,640320,0,nn);

   time1=clock();

   lsqrt(c,d,nn);

   timediff1=clock()-time1;
   printf("Time sqrt:                    %10.2f\n",(double)(timediff1)/timetic);

   lcmul(d,c,640320,0,nn);

   lmul(pi_result,d,ab,nn);

   free(a);
   free(b);
   free(c);
   free(d);
   free(ab);

   printf("\n");

   printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror);

   timediff2=clock()-time2;
   printf("\n");
   printf("Time TOTAL Pi:                %10.2f\n",(double)(timediff2)/timetic);

   printf("\n");

   time(&stop_time); 
   printf("Stop time:  ");
   day_and_time(&stop_time);

   printf("\n");
   printf("\n");

   elapseddays=get_float_days(&stop_time)-get_float_days(&start_time);

   printf("Elapsed clock time: %18.8lf days = %18.6lf hours = %18.2lf seconds\n",
          elapseddays,24.0*elapseddays,24.0*3600.0*elapseddays);

   printf("\n");

   printf("Calculation finished, output to %s.\n",file_name);

   printf("\n");

   outfile=fopen(file_name,"w");

   lfprint(pi_result,nn,outfile);  

   fclose(outfile);


} /* if */


if (method == 2) {


   printf("\n");
   printf("Pi, Ramanujan\n");
   printf("============= \n");
   printf("\n");
   printf("%s\n",version_string);
   printf("\n");

   pot=pot0;

   n=1;

   for (i=1; i<=pot; i++) {
      n=2*n;
   } /* for i */

   decimals=n;
   decimals2=decimals;

   frac=0.995;
   decimals=decimals/4;
   decimals=(int)(frac*(double)(decimals-8))-2;
   decimals=4*decimals;

   n1=0;
   n2=decimals2/8+1;

   bigsize=n/4;

   length=n/4;
   ilength=8;
   frac=0.995;

   /* nn:    Short array with data for all lfloat variables. 0 <= word <= 9999.         */
   /* nn[0]: Total length of lfloat variable in words.                                  */
   /* nn[1]: Length of integer part of lfloat variable in words.                        */
   /* nn[2]: Length of decimal part of lfloat variable in words.                        */
   /* nn[3]: Fraction of the decimals that for sure will be correct during calulations. */
   /*        Here fraction is 0.995.                                                    */

   nn[0]=length;
   nn[1]=ilength;
   nn[2]=length-ilength;
   nn[3]=(int)(frac*(double)(length-ilength))-2;

   lsize=nn[0];

   printf("mulversion = %10d\n",mulversion);
   printf("Power of 2 = %10d\n",pot);
   printf("n          = %10d\n",n);
   printf("FFT        = %10d\n",2*mulversion*bigsize);
   printf("maxsize    = %10d\n",bigsize);
   printf("decimals   = %10d  (%d)\n",decimals,decimals2);
   printf("n1         = %10d\n",n1);
   printf("n2         = %10d\n",n2);
   printf("\n");

   printf("Number of calculated decimals: %10d\n",decimals);
   printf("\n");

   temp=(double)(bigsize);
   mem1=(int)(2.4*1048576.0*log(temp)*temp/(log(1048576.0/8.0)*1048576.0/8.0)+0.5);

   mem2=8*2*2*mulversion*bigsize;
   maxmem=mem1+mem2;
   maxmem=maxmem/1024;

   printf("Max approx memory usage: %10d kbyte;   mulversion: %3d\n",maxmem,mulversion);

   printf("\n");

   time(&start_time);
   printf("Start time: ");
   day_and_time(&start_time);

   time2=clock();

   printf("\n");
   printf("\n");

   biginit2(&rP,bigsize);
   biginit2(&rQ,bigsize);
   biginit2(&rT,bigsize);

   inttobig(&rP,0);
   inttobig(&rQ,0);
   inttobig(&rT,0);

   biginit2(&tmp,bigsize);
   biginit2(&tmp2,bigsize);
   biginit2(&tmp3,bigsize);
   biginit2(&tmp4,bigsize);

   inttobig(&tmp,0);
   inttobig(&tmp2,0);
   inttobig(&tmp3,0);
   inttobig(&tmp4,0);

   printf("Calculating Pi ...\n");

   fftmaxerror=0.0;

   time1=clock();

   sum_pi_ramanujan(&rP,&rQ,&rT,n1,n2,bigsize);

   timediff1=clock()-time1;

   printf("\n");

   printf("Time Pi, fast series:         %10.2f\n",(double)(timediff1)/timetic);

   bigfree(&rP);

   bigfree(&tmp3);
   bigfree(&tmp4);

   a = (short *)calloc(lsize+2,sizeof(short));
   b = (short *)calloc(lsize+2,sizeof(short));
   c = (short *)calloc(lsize+2,sizeof(short));
   d = (short *)calloc(lsize+2,sizeof(short));
   ab = (short *)calloc(lsize+2,sizeof(short));
   pi_result = (short *)calloc(lsize+2,sizeof(short));

   inttobig(&tmp,9801);

   bigmul(&tmp,&tmp,&rQ);

   bigass(&tmp2,&rT);

   bigfree(&rQ);
   bigfree(&rT);

   lzero(a,nn);
   lzero(b,nn);
   lzero(ab,nn);

   big2l(&tmp,&tmp2,a,b,nn);

   bigfree(&tmp);
   bigfree(&tmp2);

   linttol2(d,2,0,nn);

   time1=clock();

   lsqrt(c,d,nn);

   timediff1=clock()-time1;
   printf("Time sqrt:                    %10.2f\n",(double)(timediff1)/timetic);

   lcmul(d,c,2,0,nn);

   lmul(b,b,d,nn);

   time1=clock();
       
   ldivv(pi_result,a,b,nn);

   timediff1=clock()-time1;

   printf("Time div:                     %10.2f\n",(double)(timediff1)/timetic);

   free(a);
   free(b);
   free(c);
   free(d);
   free(ab);

   printf("\n");

   printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror);

   timediff2=clock()-time2;
   printf("\n");
   printf("Time TOTAL Pi:                %10.2f\n",(double)(timediff2)/timetic);

   printf("\n");

   time(&stop_time); 
   printf("Stop time:  ");
   day_and_time(&stop_time);

   printf("\n");
   printf("\n");

   elapseddays=get_float_days(&stop_time)-get_float_days(&start_time);

   printf("Elapsed clock time: %18.8lf days = %18.6lf hours = %18.2lf seconds\n",
          elapseddays,24.0*elapseddays,24.0*3600.0*elapseddays);

   printf("\n");

   printf("Calculation finished, output to %s.\n",file_name);

   printf("\n");

   outfile=fopen(file_name,"w");

   lfprint(pi_result,nn,outfile);  

   fclose(outfile);


} /* if */


} /* End */