Continued Fraction Expansion

// Continued Fraction Expansion yields seven transcendental 
// values from one function.
// Original Turbo C, modified and tested with Pelles C  vegaseat  8oct2004

#include <stdio.h>   // in/out functions
#include <math.h>   // math functions

#define SIN   1
#define COS   2
#define TAN   3
#define EXP   4
#define SINH  5
#define COSH  6
#define TANH  7
#define PI    3.1415926536  // or use atan(1)*4.0

double trigfunc(int p, int n, double x);

//
// test trigfunc() against compiler functions
//
int main()
{
  double x;

  x = 30 * PI/180;  // convert 30 degrees to radians
  puts("Continued Fraction Expansion can yield seven transcendental values ... \n");
  printf("x = 30 degrees = %11.10f radians \n\n",x);
  printf("compiler sin(x)    = %11.10f \n",sin(x));
  printf("trigfunc(SIN,2,x)  = %11.10f \n",trigfunc(SIN,2,x));
  printf("trigfunc(SIN,3,x)  = %11.10f \n\n",trigfunc(SIN,3,x));
  printf("compiler cos(x)    = %11.10f \n",cos(x));
  printf("trigfunc(COS,3,x)  = %11.10f \n\n",trigfunc(COS,3,x));
  printf("compiler tan(x)    = %11.10f \n",tan(x));
  printf("trigfunc(TAN,3,x)  = %11.10f \n\n",trigfunc(TAN,3,x));
  printf("compiler exp(x)    = %11.10f \n",exp(x));
  printf("trigfunc(EXP,3,x)  = %11.10f \n\n",trigfunc(EXP,3,x));
  printf("compiler sinh(x)   = %11.10f \n",sinh(x));
  printf("trigfunc(SINH,3,x) = %11.10f \n\n",trigfunc(SINH,3,x));
  printf("compiler cosh(x)   = %11.10f \n",cosh(x));
  printf("trigfunc(COSH,3,x) = %11.10f \n\n",trigfunc(COSH,3,x));
  printf("compiler tanh(x)   = %11.10f \n",tanh(x));
  printf("trigfunc(TANH,3,x) = %11.10f \n\n",trigfunc(TANH,3,x));
  getchar();
  return 0;
}

//
// generate 7 transcendental values with one function 
// p = function selector,  #define the following:            
// SIN 1   COS 2   TAN 3   EXP 4   SINH 5   COSH 6   TANH 7  
// n = convergence cutoff, n = 3 is typical, 
// higher n gives higher precision, but slower execution
//
double trigfunc(int p, int n, double x)
{
  int     k; 
  double  r, s, t;

  if (p <= 3)
  {
    r = - x * x;    // trig
  } 
  else
  {     
    r = x * x;      // hyperbolic
  } 
  s = 4 * n + 2;
  for (k = n; k > 0; k--)
  {
    s = 4 * k - 2 + r/s;
  }
  switch (p % 4) 
  {
    case 0 : t = (s + x)/(s - x);         // exp
             break;
    case 1 : t = 2 * x * s/(s * s - r);   // sin, sinh
             break;
    case 2 : t = (s * s + r)/(s * s - r); // cos, cosh
             break;
    case 3 : t = 2 * x * s/(s * s + r);   // tan, tanh
             break;
  }
  return (t);
}