/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sasinh.h" #include #include float /*FUNCTION*/ sasinh( float x) { float sasinh_v; /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2001-05-25 SASINH Krogh Minor change for making .f90 version. *>> 1996-04-01 SASINH Krogh Converted to Fortran 77 from SFTRAN. *>> 1995-10-24 SASINH Krogh Removed blanks in numbers for C conversion. *>> 1994-10-19 SASINH Krogh Changes to use M77CON *>> 1994-04-19 SASINH CLL Edited type stmts to make DP & SP similar. *>> 1992-03-13 SASINH FTK Removed implicit statements. *>> 1986-03-18 SASINH Lawson Initial code. *--S replaces "?": ?ACOSH,?ACSCH,?ACTNH,?ASECH,?ASINH,?ATANH *--& ?HYPC, ?HYPS, ?HYPH, ?HYPER, ?ERM1 * * This program unit has six entry points for computing * the six inverse hyperbolic functions. * Using the SFTRAN3 preprocessor. * Uses R1MACH(3) to obtain machine precision. * All critical constants are given correctly rounded * to 25 decimal places. Accuracy of the subprogram * adapts to the machine precision up to 25 decimal * places precision. * * ----------------------------------------------------------- * *-- Begin mask code changes * The single precision code was tested in May 1983, on a * Univac 1100 using the JPL program ST4 that compares a * single precision function with a double precision function. * Each function was tested on at least 6000 points. Max * relative errors are expressed as a multiple of R1MACH(3), * which for the Univac 1100 is 2**(-27) = 0.745E-8 . * * FUNCTION RANGE MAX REL ERROR * -------- ----- ------------- * SASINH ALL X 2.4 * * SACOSH 1.0 .LE. X .LE. 1.03 LARGE * 1.03 .LE. X .LE. 1.21 3.4 * 1.21 .LE. X 2.3 * * SATANH ABS(X) .LE. 0.44 3.2 * 0.44 .LE. ABS(X) .LE. 0.92 5.1 * 0.92 .LE. ABS(X) .LE. 1.0 LARGE * * SACSCH ABS(X) .NE. 0 3.4 * * SASECH 0.0 .LT. X .LE. 0.24 2.0 * 0.24 .LE. X .LE. 0.68 3.0 * 0.68 .LE. X .LE. 0.88 10.0 * 0.88 .LE. X .LE. 1.0 LARGE * * SACTNH 1.0 .LE. ABS(X) .LE. 1.16 LARGE * 1.16 .LE. ABS(X) .LE. 2.20 7.1 * 2.20 LE. ABS(X) 3.7 * * ----------------------------------------------------------- * * The D.P. and S.P. codes have each been tested using * identities such as X - SINH( ASINH(X) ) = 0. The * D.P. codes have about the same accuracy relative to * D1MACH(3) as do the S.P. codes relative to R1MACH(3). * *-- End mask code changes * Error Messages * 1) Arg .LT. 1 ( IN SACOSH ) * 2) Arg .LT. or .EQ. to 1 ( IN SATANH ) * 3) ABS(Arg) .LT. or .EQ. to 1 ( IN SACTNH ) * 4) Arg .EQ. 0 ( IN SACSCH ) * 5) Arg .LT. or .EQ. 0 or Arg .GT. 1 ( IN SASECH ) * * ----------------------------------------------------------- * * C.L.Lawson and Stella Chan,JPL,Feb 23,1983. * */ /* Begin computation for SASINH(x) * Defined for all x.The value u ranges from negative infinity to * positive infinity as x ranges from negative infinity to positive i */ if (x == 0.e0) { sasinh_v = 0.e0; } else { sasinh_v = signf( shyps( fabsf( x ) ), x ); } return( sasinh_v ); } /* end of function */ float /*FUNCTION*/ sacosh( float x) { float sacosh_v; /* Defined (float -valued) for all X greater than or equal to 1. * We compute the nonnegative value U ranging from zero to infinity * as X ranges from 1 to infinity.The second value would be -U. */ if (x < 1.e0) { sacosh_v = 0.e0; serm1( "SACOSH", 1, 0, "ARG.LT.1", "ARG", x, '.' ); } else { sacosh_v = shypc( x ); } return( sacosh_v ); } /* end of function */ float /*FUNCTION*/ satanh( float x) { float satanh_v; /* Defined for -1 < x < +1. The value U ranges from negative * infinity to positive infinity. */ if (fabsf( x ) >= 1.e0) { satanh_v = 0.e0; serm1( "SATANH", 2, 0, "ABS(ARG).GE.1", "ARG", x, '.' ); } else if (x == 0.e0) { satanh_v = 0.e0; } else { satanh_v = signf( shyph( fabsf( x ) ), x ); } return( satanh_v ); } /* end of function */ float /*FUNCTION*/ sactnh( float x) { float sactnh_v; /* Defined for ABS(X) > 1. The value U ranges from zero to negative * infinity as X ranges from negative infinity to -1,and from * positive infinity to zero as X ranges from 1 to positive * infinity. */ if (fabsf( x ) <= 1.e0) { sactnh_v = 0.e0; serm1( "SACTNH", 3, 0, "ABS(ARG).LE.1", "ARG", x, '.' ); } else { sactnh_v = signf( shyph( 1.e0/fabsf( x ) ), x ); } return( sactnh_v ); } /* end of function */ float /*FUNCTION*/ sacsch( float x) { float sacsch_v; /* Defined for all X not equal to 0. The value varies from * zero to negative infinity as X ranges from negative * infinity to zero,jumps from negative infinity to positive * infinity at zero, and varies from positive infinity to * zero as X ranges from zero to positive infinity. */ if (x == 0.e0) { sacsch_v = 0.e0; serm1( "SACSCH", 4, 0, "ARG.EQ.0", "ARG", x, '.' ); } else { sacsch_v = signf( shyps( 1.e0/fabsf( x ) ), x ); } return( sacsch_v ); } /* end of function */ float /*FUNCTION*/ sasech( float x) { float sasech_v; /* Defined (float valued) for X greater than zero and X less than * or equal to 1.E0. We compute the nonnegative value that varies * from infinity to zero as X varies from zero to 1.E0. The other * value would be the negative of this. */ if (x <= 0.e0 || x > 1.e0) { sasech_v = 0.e0; serm1( "SASECH", 5, 0, "ARG.GT.1 OR ARG.LE.0", "ARG", x, '.' ); } else { sasech_v = shypc( 1.e0/x ); } return( sasech_v ); } /* end of function */ float /*FUNCTION*/ shyps( float p) { float cu, shyps_v; static float cut2 = 1.176e0; static float hicut = 1.e16; static float loge2 = 0.6931471805599453094172321e0; /* Here P > 0. * We break the p range into three intervals separated by CUT2 * and HICUT. In the middle range between CUT2 and HICUT we use * SHYPS=LOG(P+SQRT(1+P**2)) * CUT2 is set to 1.176 to keep the argument of LOG() greater * than e = 2.718 to avoid amplification of relative error by LOG(). * We set HICUT = 10.**16 assuming that all machines on which * this code is to run have overflow limit greater than 10.**32 * and precision not greater than 32 decimal places. Thus we * assume that HICUT**2 would not overflow and HICUT**2 plus or * minus 1.E0 would evaluate to HICUT**2. * The idea of using SHYPS = LOG(P) + LOG(2) for P .GT. HICUT is * copied from an ASINH subroutine due to WAYNE FULLERTON. * */ if (p > hicut) { shyps_v = logf( p ) + loge2; } else { cu = sqrtf( 1.e0 + SQ(p) ); if (p < cut2) { shyps_v = shyper( p, cu ); } else { shyps_v = logf( p + cu ); } } return( shyps_v ); } /* end of function */ float /*FUNCTION*/ shypc( float p) { float shypc_v, su; static float cut2 = 1.176e0; static float hicut = 1.e16; static float loge2 = 0.6931471805599453094172321e0; /* Here P .GE. 1. * See comments in function SHYPS * The mid-range formula for SACOSH is: * U=LOG(P+SQRT(P**2-1)) * */ if (p > hicut) { shypc_v = logf( p ) + loge2; } else { su = sqrtf( (p - 1.e0)*(p + 1.e0) ); if (su < cut2) { shypc_v = shyper( su, p ); } else { shypc_v = logf( su + p ); } } return( shypc_v ); } /* end of function */ float /*FUNCTION*/ shyph( float q) { float dif, qsq, shyph_v, term; static float cut1 = 0.463e0; /* Here Q satisfies 0 < Q < 1. * When Q exceeds CUT1 = 0.463 then (1+Q)/(1-Q) > e = 2.718, * and this is large enough to be used as an argument * to LOG() without amplification of relative error. * When Q is less than CUT1 we transform Q to COSH(U) * and SINH(U) and use procedure (LOW RANGE). * * The simple formulas for CU and SU would be * CU = 1.E0 / SQRT(1.E0 - Q**2) * and SU = Q / SQRT(1.E0 - Q**2) * The more complicated formulas used here have less * accumulation of round-of error since * 0.0 .LT. TERM .LT. 0.1283 * */ if (q < cut1) { qsq = q*q; dif = 1.e0 - qsq; term = qsq/(dif + sqrtf( dif )); shyph_v = shyper( q + q*term, 1.e0 + term ); } else { shyph_v = 0.5e0*logf( (1.e0 + q)/(1.e0 - q) ); } return( shyph_v ); } /* end of function */ /* do (LOW RANGE) ==> DYPER(SU, CU) * */ float /*FUNCTION*/ shyper( float su, float cu) { long int _l0, i, ii, _i, _r; static long int nmax; float r, rsq, shyper_v, u; static float c[13-(0)+1], cv[8], sv[8], v[8]; static float b1 = 0.52344e0; static float b2 = 4.4306e0; static LOGICAL32 first = TRUE; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Cv = &cv[0] - 1; float *const Sv = &sv[0] - 1; float *const V = &v[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ c[0] = 0.5e0; c[1] = -.1666666666666666666666667e00; c[2] = .7500000000000000000000000e-01; c[3] = -.4464285714285714285714286e-01; c[4] = .3038194444444444444444444e-01; c[5] = -.2237215909090909090909091e-01; c[6] = .1735276442307692307692308e-01; c[7] = -.1396484375000000000000000e-01; c[8] = .1155180089613970588235294e-01; c[9] = -.9761609529194078947368421e-02; c[10] = .8390335809616815476190476e-02; c[11] = -.7312525873598845108695652e-02; c[12] = .6447210311889648437500000e-02; c[13] = -.5740037670841923466435185e-02; V[1] = .1250000000000000000000000e00; Sv[1] = .1253257752411154569820575e00; Cv[1] = .1007822677825710859846950e01; V[2] = .2500000000000000000000000e00; Sv[2] = .2526123168081683079141252e00; Cv[2] = .1031413099879573176159295e01; V[3] = .3750000000000000000000000e00; Sv[3] = .3838510679136145687542957e00; Cv[3] = .1071140346704586767299498e01; V[4] = .5000000000000000000000000e00; Sv[4] = .5210953054937473616224256e00; Cv[4] = .1127625965206380785226225e01; V[5] = .6250000000000000000000000e00; Sv[5] = .6664922644566160822726066e00; Cv[5] = .1201753692975606324229229e01; V[6] = .7500000000000000000000000e00; Sv[6] = .8223167319358299807036616e00; Cv[6] = .1294683284676844687841708e01; V[7] = .8750000000000000000000000e00; Sv[7] = .9910066371442947560531743e00; Cv[7] = .1407868656822803158638471e01; V[8] = .1000000000000000000000000e01; Sv[8] = .1175201193643801456882382e01; Cv[8] = .1543080634815243778477906e01; _aini = 0; } /* This procedure controls computation of U, given * SU=SINH(U) and CU=COSH(U).These will satisfy * 0 .LE. SU .LE. 1.176 * AND 0 .LE. CU .LE. 1.5437 * so the result, U , will satisfy * 0 .LE. U .LE. 1.00052 * * * Here the argument R satisfies 0 .LE. R .LE. SINH(0.125) = 0.125326 * * ------------------------------------------------------------------ * * On the first call to this procedure the value NMAX * will be set and saved. * FOR MACHINE PRECISION BETWEEN 4.43062 AND 6.45985 USE NMAX = 2 * FOR MACHINE PRECISION BETWEEN 6.45985 AND 8.43090 USE NMAX = 3 * FOR MACHINE PRECISION BETWEEN 8.43090 AND 10.36773 USE NMAX = 4 * FOR MACHINE PRECISION BETWEEN 10.36773 AND 12.28199 USE NMAX = 5 * FOR MACHINE PRECISION BETWEEN 12.28199 AND 14.18024 USE NMAX = 6 * FOR MACHINE PRECISION BETWEEN 14.18024 AND 16.06654 USE NMAX = 7 * FOR MACHINE PRECISION BETWEEN 16.06654 AND 17.94359 USE NMAX = 8 * FOR MACHINE PRECISION BETWEEN 17.94359 AND 19.81325 USE NMAX = 9 * FOR MACHINE PRECISION BETWEEN 19.81325 AND 21.67688 USE NMAX = 10 * FOR MACHINE PRECISION BETWEEN 21.67688 AND 23.53550 USE NMAX = 11 * FOR MACHINE PRECISION BETWEEN 23.53550 AND 25.38987 USE NMAX = 12 * */ /* Coeffs for SASINH series. C(0) is half of its true value. * This series will be used only for arguments, R , in the * range 0 .LE. R .LE. SINH(0.125) = 0.125326 . * */ /* SV(I) = SINH(V(I)) * CV(I) = COSH(V(I)) * */ if (first) { first = FALSE; /* The following linear expression approximate the NMAX * selection criteria described above. The expression * never sets NMAX too small but will set NMAX too * large, by one, in boundary cases. */ nmax = 2 + (long)( b1*(-log10f( FLT_EPSILON/FLT_RADIX ) - b2) ); nmax = max( 2, min( nmax, 12 ) ); } if (su <= Sv[1]) { r = su; } else { ii = min( 7, (long)( su/Sv[1] ) ); if (su < Sv[ii]) ii -= 1; /* Here we transform to an argument, R , in the range * 0 .LE. R .LE. SINH(0.125). The simple formula would * be : R = SU * CV(II) - SV(II) * CU * However we transform this to extract the factor * (SU - SV(II)) for better control of round-off error. * */ r = (su - Sv[ii])*(Cv[ii] - (Sv[ii] + su)*Sv[ii]/(Cv[ii] + cu)); } rsq = r*r; u = c[nmax]; for (i = nmax - 1; i >= 0; i--) { u = rsq*u + c[i]; } shyper_v = r*(u + 0.5e0); if (su > Sv[1]) shyper_v += V[ii]; return( shyper_v ); } /* end of function */