/*************************************************** * 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 #include 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=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