#include <stdio.h>
#include <math.h> 
#include <string.h>
#define MAXLINE 1000

double pi = 3.14159265358979323846;  /* Pi */

double a=0.;
double b=0.;

void arctan(double x, int n);
void errf(double x, int n);  
double confrac(double x, double limit, void (*func)(double, int));


void arctan(double x, int n)
/*For arctan(x) we have a_0 = 0.  For n >=1, a_n = 2n-1
b_1 = x.  For n >= 2, b_n = (n-1)^2*x^2.                      */

{
  if (0 == n) {
    a=0.;
  }
  if (1==n) {
    a=1.;
    b=x;
  }
  if (n>1) {
    a=2.*n-1.;
    b=(n-1)*(n-1)*x*x;
  }
   printf("arctan: n = %d   a = %f   b = %f \n", n, a, b);
   return;
}

void errf(double x, int n) {
  if (0 == n) {
    a = 1.0;
  }
  if (1 == n) {
    a = x;
    b = -1.0*exp(-1.*x*x)/sqrt(pi);
  }
  if (1< n) {
    a = x*(2. - n%2);
    b = n-1;
  }
  return;
}

double confrac(double x, double limit, void (*func)(double, int))
  /* This function calculates a continued fraction.  Which fraction it
calculates is determined by the name "func" which it is told to calculate.

Suppose f(x) is an infinite continued fraction.
Write f as  f = a_0 + b_1(a_1 + b_2/(a_2 + ...
I.e. let f_n := a_0 + b_1/(a_1 + b_2 /(a_2 + ... a_(n-1) + b_n/a_n )) ...))
So f = Limit as n -> oo of f_n.
Write f_n = P_n/Q_n.
We take P_(-1) = 1 and Q_(-1) = 0.  We also take P_0 = a_0 and Q_0 = 1.
Then for n >= 1 we use this recursion relationship:
P_(n+1) = a_(n+1)*P_n + b_(n+1)*P_(n-1)
Q_(n+1) = a_(n+1)*Q_n + b_(n+1)*Q_(n-1) */

{
  double pold, qold, p, q, pnew, qnew;
  double out, newout;
  int n=0;
  double bignum, smallnum;
  bignum=1e40;
  smallnum=1/bignum;
  pold=1.;
  qold=0.;

  func(x, 0);
  p =a;
  q=1.;
  newout=a;

  do {
    out=newout;
    n++;
    func(x, n);
    pnew=a*p+b*pold;
    qnew=a*q+b*qold;
    newout=pnew/qnew;
/*    printf("n = %d   a = %f   b = %f \n", n, a, b);

    printf("pold = %f  qold = %f \n", pold, qold);
    printf("p = %f  q = %f \n", p, q);
    printf("pnew = %f  qnew = %f \n", pnew, qnew);
    printf("out = %f   newout = %f \n", out, newout); */
    pold=p;
    qold=q;
    p=pnew;
    q=qnew;

    /* if necessary rescale to prevent overflows or underflows */
    if((fabs(p)>bignum) ||(fabs(q)>bignum)) {
      p=p*smallnum;
      q=q*smallnum;
      pold=pold*smallnum;
      qold=qold*smallnum;
    }
    else if ((fabs(p)<smallnum) || (fabs(q)<smallnum)) {
      p=p*bignum;
      q=q*bignum;
      pold=pold*bignum;
      qold=qold*bignum;
    }
  }
  while(fabs(out-newout)>limit);
  printf("Number of iterations required = %d    \n", n);
  return out;
}

int main(void)
{

  double dval, x;
  double limit = 1e-9;
  double result;

  int num;
  int instruc;
  char line[MAXLINE], comstring[MAXLINE];
  

  for(;;)  {
    printf(" Enter command.  (Use h for help.) \n");
    fgets(line, MAXLINE, stdin);
    num = sscanf(line, "%s%lf", comstring, &dval);
    if (0 >= num)
      instruc='h';
    else
      instruc = comstring[0];
    if (2 == num) {
          x = dval;
    }

    if ('a' == instruc) {
      if (num == 1) {
        printf(" Enter number to find the arctan of. \n");
        fgets(line, MAXLINE, stdin);
        num = sscanf(line, "%lf", &x);
      }
      if (fabs(x)<1)
        result = confrac(x, limit, arctan);
      else
        result = pi*0.5 - confrac(1/x, limit, arctan);
      printf("arctan(%f) = %f\n", x, result);
    }

    if ('e' == instruc)
      exit(0);
    if ('h' == instruc) {
      printf("Help:\n");
      printf("a:  Calculate an arctan.\n");
      printf("e:  End the program.\n");
      printf("h:  Print this list.\n");
      printf("l:  Change the convergence limit.\n");
      printf("r: Calculate an error function.\n");
    }
    if ('l' == instruc) {
    printf("The convergence limit is %e", limit);
    printf(" To change it, enter a new value.\n");
    fgets(line, MAXLINE, stdin);
    num = sscanf(line, "%lf", &dval);
    if (1 == num)
     limit=dval;
    }
    if ('r' == instruc) {
      if (num==1) {
        printf(" Enter number to find the error function of. \n");
        fgets(line, MAXLINE, stdin);
        num = sscanf(line, "%lf", &x);
      }
      if (0 <= x) {
        result = confrac(x, limit, errf);
        printf("erf(%f) = %f\n", x, result);
      }
      else
        printf("x should be positive \n");

    }

  }
}