/***************************************************
*  Program to demonstrate the inverse hyperbolic   *
*  functions 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 E; int ligne,N;
double A[16], W[16];
double x,y,Z;

void Coeff(int *);
void WCoeff(double);

/*********************************************************
*        Modified cordic logarithm function
* -------------------------------------------------------
* This subroutine takes an input value and returns Y=LN(X)
* X may be any positive real value.
*********************************************************/
double XLn(double X) {
// Labels: L10,L20,L30
  int I,K;
  double aa,X1,Y;
// Get coefficients
  Coeff(&N);
//If X<=0 then an error exists, return
  if (X<=0) return -1e18;
  K=0;
// Save X
  X1=X;
// Reduce the range of X
L10: if (X<E) goto L20;
// Divide out a power of E
    K=K+1;
    X=X/E;
  goto L10;
// Test if X>=1, if so go to next step 
// Otherwise, bring X to >1
L20: if (X>=1.0) goto L30;
  K=K-1;
  X=X*E;
  goto L20;
// Determine the weighting coefficients, W(I)
L30: WCoeff(X);
// Calculate residual factor based on Z
// Want LN(Z), where Z is near unity
  Z=Z-1;
  Z=Z*(1.0-(Z/2.0)*(1.0+(Z/3.0)*(1.0-Z/4.0)));
// Assemble results
  aa=0.5;
  for (I=1; I<=N; I++) {
    Z=Z+W[I]*aa;
    aa=aa/2.0;
  }
// Z is now the mantissa, K the characteristic
  Y=K+Z;
// Restore X
  X=X1;
  return Y;
}

// Weight determination function
void WCoeff(double X) {
  int I;
  Z=X;
  for (I=1; I<=N; I++) {
    W[I]=0.0;
    if (Z>A[I]) W[I]=1.0;
    if (W[I]==1.0) Z/=A[I];
  }
}

// Exponential coefficients function
void Coeff(int *N) {
  *N=15;
  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;
  A[10]=1.000977039492417;
  A[11]=1.000488400478694;
  A[12]=1.000244170429748;
  A[13]=1.000122077763384;
  A[14]=1.000061037018933;
  A[15]=1.000030518043791;
}

/*---------------------------------------------
 *    Inverse hyperbolic sine subroutine
 * -------------------------------------------
 * This subroutine calculates the inverse of
 * the hyperbolic sine using the modified cordic
 * natural logarithm subroutine 500.
 * Input: X - Output: Y = ARCSINH(X)
 * Formula: ARCSINH(X]=LN(X+SQRT(X*X+1))
 ---------------------------------------------*/
double ArcSinH(double X) {
// Label: L10
  double X2,Y;
// Test for zero argument
  if (X!=0.0) goto L10;
  return 0;
L10: //Save X
  X2=X;
  X=fabs(X);
  X=X+sqrt(X*X+1.0);
// Get LN(X)
  Y=XLn(X);
// Insert sign
  Y=(X2/fabs(X2))*Y;
// Restore X
  X=X2;
  return Y;
}

// Inverse hyperbolic cosine subroutine
//    ARCCOSH(X]=LN(X+SQRT(X*X-1))
double ArcCosH(double X) {
// Label: L10
  double X2,Y;
// Test for argument less than or equal to unity
  if (X>1.0) goto L10;
  return 0.0;
L10: //Save X
  X2=X;
  X=fabs(X);
  X=X+sqrt(X-1.0)*sqrt(X+1.0);
  Y=XLn(X);
// Restore X
  X=X2;
  return Y;
}

// Inverse hyperbolic tangent subroutine
//    ARCTANH(X]=1/2*LN(1+x/1-x)
double ArcTanH(double X) {
// Label: L10,L20
  double X2,Y;
// Test for X>= +/- 1 
  if (fabs(X)<=0.999999) goto L10;
// Y is BIG! (here +/- 1E18)
  Y=(X/fabs(X))*1e18;
  return Y;
L10: // Test for zero argument
  if (X!=0.0) goto L20;
  return 0.0;
L20: // Save X
  X2=X;
  X=(1.0+X)/(1.0-X);
  Y=XLn(X);
// Restore X
  X=X2;
  return Y;
}


void main() {
  ligne=1;
  printf("\n");
  printf("   X      ARCSINH(X)      ARCCOSH(X)      ARCTANH(X) \n"); 
  printf(" ----------------------------------------------------\n");
  x=-3.0;
  while (x<=3.2) {
    printf("%6.1f  ", x);
    y=ArcSinH(x);
    printf("%14.10f  ", y);
    y=ArcCosH(x);
    printf("%14.10f  ", y);
    y=ArcTanH(x);
    if (fabs(y)<1e12)
      printf("%14.10f\n", y);
    else
      if (y<0) 
	printf(" -1E18\n");
      else 
	printf("  1E18\n");
    if (ligne==20) {
      ligne=0;
      getchar();
      printf("\n");
    }
    ligne++;
    x+=0.2;
  }
  printf("\n\n");
}

// End of file invhyper.cpp