/****************************************************
* Program to demonstrate the hyperbolic subroutines *
* ------------------------------------------------- *
* Reference: BASIC Scientific Subroutines, Vol. II  *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
*                                                   *
*                      C++ version by J-P Moreau.   *
*                          (www.jpmoreau.fr)        *
****************************************************/
#include <stdio.h>
#include <math.h>

double  A[10], W[10];
int     N,ligne;
double  aa,E,K,XX,Y,Z;


void Coeffs();
void WCoeffs(double X);

//Modified cordic exponential subroutine
double XExp(double X) {
// This subroutine takes an input value and returns Y:=EXP(X)
// X may be any positive or negative real value
double Y;
int I;
// Get coefficients 
  Coeffs();
// Reduce the range of X
  K=int(X);
  X=X-K;
// Determine the weighting coeffs, W(I)
  WCoeffs(X);
// Calculate products
  Y=1.0;
  for (I=1; I<=N; I++)
    if (W[I]>0.0) Y*=A[I];
// Perform residual multiplication
  Y=Y*(1.0+Z*(1.0+Z/2.0*(1.0+Z/3.0*(1.0+Z/4.0))));
// Account for factor EXP(K)
  if (K<0) E=1.0/E;
  if (fabs(K)<1) return Y;
  for (I=1; I<=fabs(K); I++)  Y*=E;
// Restore X
  X+=K;
  return Y;
}

void WCoeffs(double X) {
  int I;
  aa=0.5;
  Z=X;
  for (I=1; I<=N; I++) {
    W[I]=0.0;
    if (Z>aa) W[I]=1.0;
    Z=Z-W[I]*aa;
    aa/=2.0;
  }
}

void Coeffs() {
  N=9;
  E=2.718281828459045;
  A[1]=1.648721270700128;
  A[2]=1.284025416687742;
  A[3]=1.133148453066826;
  A[4]=1.064494458917859;
  A[5]=1.031743407499103;
  A[6]=1.015747708586686;
  A[7]=1.007843097206448;
  A[8]=1.003913889338348;
  A[9]=1.001955033591003;
}

/*--------------------------------------------*
*          Hyperbolic sine Function           *
* ------------------------------------------- *
* This Procedure uses the definition of the   *
* hyperbolic sine and the modified cordic     *
* exponential Function XExp(X) to approximate *
* SINH(X) over the entire range of real X.    *
* -------------------------------------------*/
double SinH(double X) {
//Label: L10;
double Y; int I;
// Is X small enough to cause round off erroe ?
  if (fabs(X)<0.35) goto L10;
// Calculate SINH(X) using exponential definition
  Y=XExp(X);
  Y=(Y-(1.0/Y))/2.0;
  return Y;
L10: // series approximation (for X small)
  Z=1.0; Y=1.0;
  for (I=1; I<9; I++) {
    Z=Z*X*X/((2*I)*(2*I+1));
    Y=Y+Z;
  }
  Y*=X;
  return Y;
}

// hyperbolic cosine Function
double CosH(double X) {
  double Y;
  Y=XExp(X);
  Y=(Y+(1.0/Y))/2.0;
  return Y;
}

// hyperbolic tangent Function
//   TANH(X]:=SINH(X)/COSH(X)
double TanH(double X) {
  double V,Y;
  V=SinH(X);
  Y=CosH(X);
  return (V/Y);
}


void main()  {
  ligne=1;
  printf("\n   X        SINH(X)         COSH(X)         TANH(X) \n");
  printf(" ---------------------------------------------------\n");
  XX=-5.0;
  do {
    printf("%6.1f  ", XX);
    Y=SinH(XX);
    printf("%14.10f  ", Y);
    Y=CosH(XX);
    printf("%14.10f  ", Y);
    Y=TanH(XX);
    printf("%14.10f  ", Y);
    if (ligne==20)  {
      ligne=0;
      getchar();
    }
	else 
	  printf("\n");
    ligne++;
    XX+=0.2;
  } while (XX<=5.2);
  printf("\n\n");
}

// End of file hyper.cpp