/*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 <math.h>
#include "fcrt.h"
#include "sasinh.h"
#include <float.h>
#include <stdlib.h>
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 */