// Calculate 'e' and 'pi'...

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

double arctan(double X1) {
  const double        R1 =  0x1.9310cfe85307cp+3;  /*  1.2595802263029547e+1  EMAS HEX = 41C98867F42983DF */
  const double        R2 = -0x1.58beca04f1805p+6;  /* -8.6186317517509520e+1  EMAS HEX = C2562FB2813C6014 */
  const double        R3 = -0x1.46d547fed8a3dp+0;  /* -1.2766919133361079e+0  EMAS HEX = C1146D547FED8A3D */
  const double        R4 = -0x1.57bd961f06c89p-4;  /* -8.3921038065840512e-2  EMAS HEX = C0157BD961F06C89 */
  const double        S1 =  0x1.b189e39236635p+4;  /*  2.7096164294378656e+1  EMAS HEX = 421B189E39236635 */
  const double        S2 =  0x1.a3b86f7830dc0p+2;  /*  6.5581320451487386e+0  EMAS HEX = 4168EE1BDE0C3700 */
  const double        S3 =  0x1.1273f9e5eff20p+1;  /*  2.1441643116703661e+0  EMAS HEX = 41224E7F3CBDFE41 */
  const double        S4 =  0x1.44831dafbf542p+0;  /*  1.2676256708212610e+0  EMAS HEX = 41144831DAFBF542 */
  const double       RT3 =  0x1.bb67ae8584caap+0;  /*  1.7320508075688772e+0  EMAS HEX = 411BB67AE8584CAA */
  const double     PIBY6 =  0x1.0c152382d7365p-1;  /*  5.2359877559829887e-1  EMAS HEX = 40860A91C16B9B2C */
  const double   PIBY2M1 =  0x1.243f6a8885a30p-1;  /*  5.7079632679489661e-1  EMAS HEX = 40921FB54442D184 */
  const double     RT3M1 =  0x1.76cf5d0b09955p-1;  /*  7.3205080756887728e-1  EMAS HEX = 40BB67AE8584CAA7 */
  const double TANPIBY12 =  0x1.126145e9ecd56p-2;  /*  2.6794919243112271e-1  EMAS HEX = 404498517A7B3559 */
  const double       ONE =  0x1.0000000000000p+0;  /*  1.0000000000000000e+0  EMAS HEX = 4110000000000000 */
  double XX1,XSQ,CONSTANT=0.0;
  int SIGN=0,INV=0;

  if (X1<0) {SIGN=1; XX1=-X1;} else XX1=X1;
  if (XX1>ONE) {XX1=1.0/XX1; INV=1;}
  if (XX1>TANPIBY12) {XX1=(RT3M1*XX1-1.0+XX1)/(XX1+RT3); CONSTANT=PIBY6;}

  XSQ=XX1*XX1;
  XX1 = XX1*(R1/(XSQ+S1+R2/(XSQ+S2+R3/(XSQ+S3+R4/(XSQ+S4)))));  // <--- taylor series expansion
  XX1 = XX1 + CONSTANT;

  if (INV) XX1=1.0-XX1+PIBY2M1;
  if (SIGN) XX1=-XX1;
  return XX1;
}

//double pi(int a, int n, int b, int limit) {
//  if (numerator == limit) return ...;
//  return ...;  // need a good expansion series here.
//}

double e(double numerator, int limit) {
  if (numerator == limit) return 0;
  return (numerator==0?2:numerator) + ((numerator==0?1:numerator) / e(numerator+1, limit));
}

int main(int argc, char **argv) {
  double r, prev_r = 0.0;
  int i = 2;
  for (;;) {
    r = e(0, i);
    if (r == prev_r) break; else fprintf(stderr, " e = %0.15lf (%d)\n", r, i-1);
    prev_r = r;
    i = i + 1;
  }
  fprintf(stderr, "Algorithm converged after %d calls\n\n", i-2);
  fprintf(stderr, " e = %0.15lf\n", e(0,i-1));

  //  prev_r = 0.0; i = 2;
  //  for (;;) {
  //    r = pi(1, i);
  //    if (r == prev_r) break; else fprintf(stderr, "pi = %0.12lf (%d)\n", 4*r, i-1);
  //    prev_r = r;
  //    i = i + 1;
  //  }
  //  fprintf(stderr, "Algorithm converged after %d calls at pi = %0.12lf\n", i-3, 4*r);

  fprintf(stderr, "Pi = %0.15lf\n", arctan(1)*4);
  
  exit(EXIT_SUCCESS);
}