NAME

mathr.h - Real numbers math library

SYNOPSIS

Macros

#define PI 3.1415926535897932384626433832795028841971693994
#define EULERS 2.7182818284590452353602874713526624977572470937
#define INFP 0x7FF0000000000000ull
#define INFN 0xFFF0000000000000ull
#define NAN 0x7FFFFFFFFFFFFFFFull
#define sgn(x) (x > 0 ? 1 : x < 0 ? -1 : 0)
#define abs(x) (x > 0 ? x : -x)

Functions

double fabs (double x)
Returns the absolute value of x. double ceil (double x)
Ceiling function. double floor (double x)
Floor function. double round (double x)
Round function. double trunc (double x)
Truncate function. double exp (double x)
Returns the exponential of x. double sqrt (double x)
Square root function. double cbrt (double x)
Cube root function. double log (double x)
Natural logarithm function (base e) double log10 (double x)
Base 10 logarithm function. double cos (double x)
Cosine function. double sin (double x)
Sine function. double tan (double x)
Tangent function. double sec (double x)
Secant function. double csc (double x)
Cosecant function. double cot (double x)
Cotangent function. double exs (double x)
Exsecant function. double exc (double x)
Excosecant function. double crd (double x)
Chord function. double acos (double x)
Inverse cosine function. double asin (double x)
Inverse sine function. double atan (double x)
Inverse tangent function. double asec (double x)
Inverse secant function. double acsc (double x)
Inverse cosecant function. double acot (double x)
Inverse cotangent function. double aexs (double x)
Inverse exsecant function. double aexc (double x)
Inverse excosecant function. double acrd (double x)
Inverse chord function. double cosh (double x)
Hyperbolic cosine function. double sinh (double x)
Hyperbolic sine function. double tanh (double x)
Hyperbolic tangent function. double sech (double x)
Hyperbolic secant function. double csch (double x)
Hyperbolic cosecant function. double coth (double x)
Hyperbolic cotangent function. double acosh (double x)
Inverse hyperbolic cosine function. double asinh (double x)
Inverse hyperbolic sine function. double atanh (double x)
Inverse hyperbolic tangent function. double asech (double x)
Inverse hyperbolic secant function. double acsch (double x)
Inverse hyperbolic cosecant function. double acoth (double x)
Inverse hyperbolic cotangent function. double ver (double x)
Versed sine function. double vcs (double x)
Versed cosine function. double cvs (double x)
Coversed sine function. double cvc (double x)
Coversed cosine function. double hv (double x)
Haversed sine function. double hvc (double x)
Haversed cosine function. double hcv (double x)
Hacoversed sine function. double hcc (double x)
Hacoversed cosine function. double aver (double x)
Inverse versed sine function. double avcs (double x)
Inverse versed sine. double acvs (double x)
Inverse coversed sine function. double acvc (double x)
Inverse versed cosine. double ahv (double x)
Inverse haversed sine. double ahvc (double x)
Inverse haversed cosine. double ahcv (double x)
Inverse hacoversed sine. double ahcc (double x)
Inverse hacoversed cosine. double pow (double x, double y)
Expontation function. double fmod (double x, double y)
Return x mod y in exact arithmetic. double atan2 (double y, double x)
Inverse tangent function. double hypot (double x, double y)
hypot double log2p (double x, double y)
double log1p (double x)
double expm1 (double x)
double scalbn (double x, int n)
double copysign (double x, double y)
Returns a value with the magnitude of x and with the sign bit of y. int rempio2 (double x, double *y)
unsigned int log2i (unsigned int x)

DESCRIPTION

Method


    acos(x)  = pi/2 - asin(x)


    acos(-x) = pi/2 + asin(x)



    For |x|<=0.5


    acos(x)  = pi/2 - (x + x*x^2*R(x^2))        (see asin.c)



    For x>0.5


    acos(x)  = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))


             = 2asin(sqrt((1-x)/2))


             = 2s + 2s*z*R(z)   ...z=(1-x)/2, s=sqrt(z)


             = 2f + (2c + 2s*z*R(z))


    where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term


    for f so that f+c ~ sqrt(z).



    For x<-0.5


    acos(x)  = pi - 2asin(sqrt((1-|x|)/2))


             = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)

Special cases


    if x is NaN, return NaN


    if |x|>1, return NaN

Definition at line 92 of file acos.c.

double acosh (double x)

Inverse hyperbolic cosine function.

Based on


    acosh(x) = log [ x + sqrt(x*x-1) ]

we have


    acosh(x) = log(x)+ln2, if x is large; else


    acosh(x) = log(2x-1/(sqrt(x*x-1)+x)) if x>2; else


    acosh(x) = log1p(t+sqrt(2.0*t+t*t)); where t=x-1

Special cases


    acosh(x) is NaN if x<1


    acosh(NaN) is NaN

Definition at line 69 of file acosh.c.

double acot (double x)

Inverse cotangent function.

Method


    arccot(x) = arctan(1/x)

Definition at line 45 of file acot.c.

double acoth (double x)

Inverse hyperbolic cotangent function.

Method


                          x + 1


    acoth(x) = 0.5 * ln( ------- )


                          x - 1


    when x in [-1, 1]


    acoth(x) = NaN

Definition at line 49 of file acoth.c.

double acrd (double x)

Inverse chord function.

Method


    arccrd(x) = 2*arcsin(x/2)

Definition at line 45 of file acrd.c.

double acsc (double x)

Inverse cosecant function.

Method


    arccsc(x) = arcsin(1/x)

Definition at line 45 of file acsc.c.

double acsch (double x)

Inverse hyperbolic cosecant function.

Method


                    1+sqrt(1+x*x)


    acsch(x) = ln( --------------- )


                          x


    when x is 0


    acsch(x) = NaN

Definition at line 49 of file acsch.c.

double acvc (double x)

Inverse versed cosine.

Method


    acvc(x) = asin(1+x)

Definition at line 45 of file acvc.c.

double acvs (double x)

Inverse coversed sine function.

Definition at line 40 of file acvs.c.

double aexc (double x)

Inverse excosecant function.

Method


    aexcsc(x) = arccsc(x+1)


              = arcsin(1/(x+1))

Definition at line 46 of file aexc.c.

double aexs (double x)

Inverse exsecant function.

Method


    aexsec(x) = arcsec(x+1)


              = arccos(1/(x+1))


              = arctan(sqrt(x^2+2*X))

Definition at line 47 of file aexs.c.

Method


    arcsec(x) = arccos(1/x)

Definition at line 45 of file asec.c.

double asech (double x)

Inverse hyperbolic secant function.

Method


                    1+sqrt(1-x*x)


    asech(x) = ln( --------------- )


                          x


    when x <= 0


    asech(x) = NaN



    when x > 1


    asech(x) = NaN

Definition at line 52 of file asech.c.

double asin (double x)

Inverse sine function.

Method

Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
we approximate asin(x) on [0,0.5] by


     asin(x) = x + x*x^2*R(x^2)
where


     R(x^2) is a rational approximation of (asin(x)-x)/x^3
and its remez error is bounded by


     |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)

For x in [0.5,1]


     asin(x) = pi/2-2*asin(sqrt((1-x)/2))
Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
then for x>0.98


     asin(x) = pi/2 - 2*(s+s*z*R(z))


             = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)

For x<=0.98, let pio4_hi = pio2_hi/2, then


     f = hi part of s;


     c = sqrt(z) - f = (z-f*f)/(s+f)     ...f+c=sqrt(z)
and


     asin(x) = pi/2 - 2*(s+s*z*R(z))


             = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)


             = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))

Special cases


    if x is NaN, return NaN


    if |x|>1, return NaN

Definition at line 100 of file asin.c.

double asinh (double x)

Inverse hyperbolic sine function.

Method


    Based on


    asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]



    we have


    asinh(x) = x  if  1+x*x=1,


             = sign(x)*(log(x)+ln2)) for large |x|, else


             = sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else


             = sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))

Definition at line 68 of file asinh.c.

double atan (double x)

Inverse tangent function.

Method


    1. Reduce x to positive by atan(x) = -atan(-x).


    2. According to the integer k=4t+0.25 chopped, t=x, the argument


       is further reduced to one of the following intervals and the


       arctangent of t is evaluated by the corresponding formula:



    [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)


    [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )


    [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )


    [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )


    [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )

Constants


    The hexadecimal values are the intended ones for the following


    constants. The decimal values may be used, provided that the


    compiler will convert from decimal to binary accurately enough


    to produce the hexadecimal values shown.

Definition at line 103 of file atan.c.

double atan2 (double y, double x)

Inverse tangent function.

Parameters:

y,x

Method


    1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).


    2. Reduce x to positive by (if x and y are unexceptional):


       ARG (x+iy) = arctan(y/x)           ... if x > 0,


       ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,

Special cases


    ATAN2((anything), NaN ) is NaN;


    ATAN2(NAN , (anything) ) is NaN;


    ATAN2(+-0, +(anything but NaN)) is +-0  ;


    ATAN2(+-0, -(anything but NaN)) is +-pi ;


    ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;


    ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;


    ATAN2(+-(anything but INF and NaN), -INF) is +-pi;


    ATAN2(+-INF,+INF ) is +-pi/4 ;


    ATAN2(+-INF,-INF ) is +-3pi/4;


    ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;

Constants


    The hexadecimal values are the intended ones for the following


    constants. The decimal values may be used, provided that the


    compiler will convert from decimal to binary accurately enough


    to produce the hexadecimal values shown.

Definition at line 87 of file atan2.c.

double atanh (double x)

Inverse hyperbolic tangent function.

Method



1.Reduced x to positive by atanh(-x) = -atanh(x)
2.For x>=0.5


              1               2x                         x


  atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)


              2             1 - x                      1 - x


  For x<0.5


  atanh(x) = 0.5*log1p(2x+2x*x/(1-x))

Special cases


    atanh(x) is NaN if |x| > 1


    atanh(NaN) is that NaN


    atanh(+-1) is +-INF

Definition at line 72 of file atanh.c.

double avcs (double x)

Inverse versed sine.

avcs(x) = acos(1+x)

Definition at line 44 of file avcs.c.

Returns:

x rounded toward -inf to integral value

Method


    Bit twiddling

Exception


    Inexact flag raised if x not equal to ceil(x).

Definition at line 63 of file ceil.c.

double copysign (double x, double y)

Returns a value with the magnitude of x and with the sign bit of y.

Definition at line 47 of file csign.c.

double cos (double x)

Cosine function.

Parameters:

Returns:

Cosine function of x

Kernel function:



__kernel_sin        ... sine function on [-pi/4,pi/4]
__kernel_cos        ... cose function on [-pi/4,pi/4]
__ieee754_rem_pio2  ... argument reduction routine

Method:



Let S,C and T denote the sin, cos and tan respectively on
[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
in [-pi/4 , +pi/4], and let n = k mod 4.
We have


     n        sin(x)      cos(x)        tan(x)
----------------------------------------------------------


     0          S           C             T


     1          C          -S           -1/T


     2         -S          -C             T


     3         -C           S           -1/T
----------------------------------------------------------

Special cases:



Let trig be any of sin, cos, or tan.
trig(+-INF)  is NaN
trig(NaN)    is that NaN

Accuracy:



TRIG(x) returns trig(x) nearly rounded

Definition at line 87 of file cos.c.

double cosh (double x)

Hyperbolic cosine function. Mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2

Method



1. Replace x by |x| (cosh(x) = cosh(-x))
2.


                                              [ exp(x) - 1 ]^2


   0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------


                                                2*exp(x)


                                         exp(x) +  1/exp(x)


   ln2/2    <= x <= 22     :  cosh(x) := -------------------


                                                 2


   22       <= x <= lnovft :  cosh(x) := exp(x)/2


   lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)


   ln2ovft  <  x           :  cosh(x) := huge*huge (overflow)

Special cases:



cosh(x) is |x| if x is +INF, -INF, or NaN
only cosh(0)=1 is exact for finite x

Definition at line 83 of file cosh.c.

double cot (double x)

Cotangent function.

Parameters:

cot(x) = 1/tan(x)


       = cos(x)/sin(x)


       = sin(2*x)/(cos(2*x)-1)

Definition at line 47 of file cot.c.

double coth (double x)

Hyperbolic cotangent function.

coth(x) = cosh(x)/sinh(x)

Definition at line 44 of file coth.c.

double crd (double x)

Chord function.

crd(x) = 2*sin(x/2)

Definition at line 44 of file crd.c.

double csc (double x)

Cosecant function.

csc = sin(1/x)


    = -2*sin(x)/(cos(2*x) - 1)

Definition at line 45 of file csc.c.

double csch (double x)

Hyperbolic cosecant function.

csch(x) = 1/sinh(x)

Definition at line 44 of file csch.c.

double cvc (double x)

Coversed cosine function.

cvc(x) = 1+sin(x)

Definition at line 44 of file cvc.c.

double cvs (double x)

Coversed sine function.

cvs(x) = 1-sin(x)

Definition at line 44 of file cvs.c.

double exc (double x)

Excosecant function.

excsc(x) = csc(x)-1


         = (1-sin(x))/sin(x)


         = cvs(x)/sin(x)


         = cvs(x)*csc(x)

Definition at line 47 of file exc.c.

double exp (double x)

Returns the exponential of x.

Method



1. Argument reduction:


   Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.


   Given x, find r and integer k such that


           x = k*ln2 + r,  |r| <= 0.5*ln2.


   Here r will be represented as r = hi-lo for better


   accuracy.
2. Approximation of exp(r) by a special rational function on


   the interval [0,0.34658]:


   Write


       R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...


   We use a special Remes algorithm on [0,0.34658] to generate


   a polynomial of degree 5 to approximate R. The maximum error


   of this polynomial approximation is bounded by 2**-59. In


   other words,


       R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5


       (where z=r*r, and the values of P1 to P5 are listed below)


   and


       |                  5          |     -59


       | 2.0+P1*z+...+P5*z   -  R(z) | <= 2


       |                             |


   The computation of exp(r) thus becomes


                      2*r


       exp(r) = 1 + -------


                     R - r


                          r*R1(r)


              = 1 + r + ----------- (for better accuracy)


                         2 - R1(r)


   where


                2       4                     10


       R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
3. Scale back to obtain exp(x):


   From step 1, we have


   exp(x) = 2^k * exp(r)

Special cases:



exp(INF) is INF, exp(NaN) is NaN;
exp(-INF) is 0, and
for finite argument, only exp(0)=1 is exact.

Accuracy:



according to an error analysis, the error is always less than
1 ulp (unit in the last place).

Misc. info:



For IEEE double


    if x >  7.09782712893383973096e+02 then exp(x) overflow


    if x < -7.45133219101941108420e+02 then exp(x) underflow

Constants:

The hexadecimal values are the intended ones for the following
constants. The decimal values may be used, provided that the
compiler will convert from decimal to binary accurately enough
to produce the hexadecimal values shown.

Definition at line 138 of file exp.c.

double expm1 (double x)

Definition at line 153 of file expm1.c.

double exs (double x)

Exsecant function.

exsec(x) = sec(x)-1


         = (1-cos(x))/cos(x)


         = ver(x)/cos(x)


         = ver(x)*sec(x)


         = 2*sin(x/2)*sin(x/2)*sec(x)

Definition at line 48 of file exs.c.

double fabs (double x)

Returns the absolute value of x.

Definition at line 51 of file fabs.c.

double floor (double x)

Floor function.

Returns:

x rounded toward -inf to integral value

Method:


    Bit twiddling

Exception:


    Inexact flag raised if x not equal to floor(x)

Definition at line 62 of file floor.c.

double fmod (double x, double y)

Return x mod y in exact arithmetic. Method: Shift and subtract

Definition at line 58 of file fmod.c.

double hcc (double x)

Hacoversed cosine function.

hcc(x) = (1+sin(x))/2

Definition at line 44 of file hcc.c.

double hcv (double x)

Hacoversed sine function.

hcv(x) = (1-sin(x))/2

Definition at line 44 of file hcv.c.

double hv (double x)

Haversed sine function.

hv(x) = (1-cos(x))/2

Definition at line 44 of file hv.c.

double hvc (double x)

Haversed cosine function.

hvc(x) = (1+cos(x))/2

Definition at line 44 of file hvc.c.

double hypot (double x, double y)

hypot

Method


 If (assume round-to-nearest) z=x*x+y*y


 has error less than sqrt(2)/2 ulp, than


 sqrt(z) has error less than 1 ulp (exercise).



 So, compute sqrt(x*x+y*y) with some care as


 follows to get the error below 1 ulp:



 Assume x>y>0;


 (if possible, set rounding to round-to-nearest)


 1. if x > 2y  use


    x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y


 where x1 = x with lower 32 bits cleared, x2 = x-x1; else


 2. if x <= 2y use


    t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))


 where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,


 y1= y with lower 32 bits chopped, y2 = y-y1.



 NOTE: scaling may be necessary if some argument is too


       large or too tiny

Special cases:


 hypot(x,y) is INF if x or y is +INF or -INF; else


 hypot(x,y) is NAN if x or y is NAN.

Accuracy:


    hypot(x,y) returns sqrt(x^2+y^2) with error less


    than 1 ulps (units in the last place)

Definition at line 81 of file hypot.c.

double log (double x)

Natural logarithm function (base e)

Method


  1. Argument Reduction: find k and f such that


        x = 2^k * (1+f),


    where  sqrt(2)/2 < 1+f < sqrt(2) .



  2. Approximation of log(1+f).


 Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)


     = 2s + 2/3 s**3 + 2/5 s**5 + .....,


         = 2s + s*R


     We use a special Remes algorithm on [0,0.1716] to generate


    a polynomial of degree 14 to approximate R The maximum error


 of this polynomial approximation is bounded by 2**-58.45. In


 other words,


            2      4      6      8      10      12      14


     R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s


    (the values of Lg1 to Lg7 are listed in the program)


 and


     |      2          14          |     -58.45


     | Lg1*s +...+Lg7*s    -  R(z) | <= 2


     |                             |


 Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.


 In order to guarantee error in log below 1ulp, we compute log


 by


    log(1+f) = f - s*(f - R)    (if f is not too large)


    log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)



 3. Finally,  log(x) = k*ln2 + log(1+f).


            = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))


    Here ln2 is split into two floating point number:


        ln2_hi + ln2_lo,


    where n*ln2_hi is always exact for |n| < 2000.

Special cases:


 log(x) is NaN with signal if x < 0 (including -INF) ;


 log(+INF) is +INF; log(0) is -INF with signal;


 log(NaN) is that NaN with no signal.

Accuracy:


 according to an error analysis, the error is always less than


 1 ulp (unit in the last place).

Constants:
The hexadecimal values are the intended ones for the following
constants. The decimal values may be used, provided that the
compiler will convert from decimal to binary accurately enough
to produce the hexadecimal values shown.

Definition at line 109 of file log.c.

double log10 (double x)

Base 10 logarithm function.

Method:



Let log10_2hi = leading 40 bits of log10(2) and


    log10_2lo = log10(2) - log10_2hi,


    ivln10    = 1/log(10) rounded.
Then


    n = ilogb(x),


    if(n<0)  n = n+1;


    x = scalbn(x,-n);


    log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
Note 1:
To guarantee log10(10**n)=n, where 10**n is normal, the rounding
mode must set to Round-to-Nearest.
Note 2:
[1/log(10)] rounded to 53 bits has error .198 ulps;
log10 is monotonic at all binary break points.

Special cases:



log10(x) is NaN with signal if x < 0;
log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
log10(NaN) is that NaN with no signal;
log10(10**N) = N  for N=0,1,...,22.

Constants:
The hexadecimal values are the intended ones for the following constants.
The decimal values may be used, provided that the compiler will convert
from decimal to binary accurately enough to produce the hexadecimal values
shown.

Definition at line 93 of file log10.c.

Method:  Let x =  2   * (1+f)


 1. Compute and return log2(x) in two pieces:


    log2(x) = w1 + w2,


    where w1 has 53-24 = 29 bit trailing zeros.


 2. Perform y*log2(x) = n+y' by simulating muti-precision


    arithmetic, where |y'|<=0.5.


 3. Return x**y = 2**n*exp(y'*log2)

Special cases:


 1.  (anything) ** 0  is 1


 2.  (anything) ** 1  is itself


 3.  (anything) ** NAN is NAN


 4.  NAN ** (anything except 0) is NAN


 5.  +-(|x| > 1) **  +INF is +INF


 6.  +-(|x| > 1) **  -INF is +0


 7.  +-(|x| < 1) **  +INF is +0


 8.  +-(|x| < 1) **  -INF is +INF


 9.  +-1         ** +-INF is NAN


 10. +0 ** (+anything except 0, NAN)               is +0


 11. -0 ** (+anything except 0, NAN, odd integer)  is +0


 12. +0 ** (-anything except 0, NAN)               is +INF


 13. -0 ** (-anything except 0, NAN, odd integer)  is +INF


 14. -0 ** (odd integer) = -( +0 ** (odd integer) )


 15. +INF ** (+anything except 0,NAN) is +INF


 16. +INF ** (-anything except 0,NAN) is +0


 17. -INF ** (anything)  = -0 ** (-anything)


 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)


 19. (-anything except 0 and inf) ** (non-integer) is NAN

Accuracy:


 pow(x,y) returns x**y nearly rounded. In particular


        pow(integer,integer)


 always returns the correct integer provided it is


 representable.

Constants :
The hexadecimal values are the intended ones for the following
constants. The decimal values may be used, provided that the
compiler will convert from decimal to binary accurately enough
to produce the hexadecimal values shown.

Definition at line 138 of file pow.c.

sec(x) = 1/cos(x)


       = 1/sqrt(cos(x)*cos(x))

Definition at line 45 of file sec.c.

double sech (double x)

Hyperbolic secant function.

sech(x) = 1/cosh(x)

Definition at line 44 of file sech.c.

double sin (double x)

Sine function.

Returns:

Sine function of x

Kernel function:



__kernel_sin        ... sine function on [-pi/4,pi/4]
__kernel_cos        ... cose function on [-pi/4,pi/4]
__ieee754_rem_pio2  ... argument reduction routine

Method:



Let S,C and T denote the sin, cos and tan respectively on
[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
in [-pi/4 , +pi/4], and let n = k mod 4.
We have


     n        sin(x)      cos(x)        tan(x)
----------------------------------------------------------


     0          S           C             T


     1          C          -S           -1/T


     2         -S          -C             T


     3         -C           S           -1/T
----------------------------------------------------------

Special cases:



Let trig be any of sin, cos, or tan.
trig(+-INF)  is NaN
trig(NaN)    is that NaN

Accuracy:



TRIG(x) returns trig(x) nearly rounded

Definition at line 86 of file sin.c.

double sinh (double x)

Hyperbolic sine function.

Method
mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2


 1. Replace x by |x| (sinh(-x) = -sinh(x)).


 2.


                                        E + E/(E+1)


     0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)


                                2



     22       <= x <= lnovft :  sinh(x) := exp(x)/2


     lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)


     ln2ovft  <  x      :  sinh(x) := x*shuge (overflow)

Special cases


    sinh(x) is |x| if x is +INF, -INF, or NaN.


    only sinh(0)=0 is exact for finite x.

Definition at line 77 of file sinh.c.

double sqrt (double x)

Square root function.

Returns:

Correctly rounded sqrt

Method:


  Bit by bit method using integer arithmetic. (Slow, but portable)


  1. Normalization


     Scale x to y in [1,4) with even powers of 2:


     find an integer k such that  1 <= (y=x*2^(2k)) < 4, then


     sqrt(x) = 2^k * sqrt(y)


  2. Bit by bit computation


     Let q  = sqrt(y) truncated to i bit after binary point (q = 1),


          i                          0


                                    i+1         2


         s  = 2*q , and y  =  2   * ( y - q  ).     (1)


          i      i            i                 i



     To compute q    from q , one checks whether


        i+1       i





  -(i+1) 2

(q + 2 ) <= y. (2)
i
-(i+1)
If (2) is false, then q = q ; otherwise q = q + 2 .
i+1 i i+1 i



     With some algebric manipulation, it is not difficult to see


     that (2) is equivalent to


                            -(i+1)


        s  +  2       <= y          (3)


         i                i



     The advantage of (3) is that s  and y  can be computed by


                  i      i


     the following recurrence formula:


         if (3) is false



         s     =  s  ,  y    = y   ;            (4)


          i+1      i         i+1    i



         otherwise,


                        -i                     -(i+1)


         s    =  s  + 2  ,  y    = y  -  s  - 2         (5)


          i+1      i          i+1    i     i



     One may easily use induction to prove (4) and (5).


     Note. Since the left hand side of (3) contain only i+2 bits,


         it does not necessary to do a full (53-bit) comparison


         in (3).


  3. Final rounding


     After generating the 53 bits result, we compute one more bit.


     Together with the remainder, we can decide whether the


     result is exact, bigger than 1/2ulp, or less than 1/2ulp


     (it will never equal to 1/2ulp).


     The rounding mode can be detected by checking whether


     huge + tiny is equal to huge, and whether huge - tiny is


     equal to huge for some floating point number 'huge' and 'tiny'.

Special cases:


  sqrt(+-0) = +-0    ... exact


  sqrt(inf) = inf


  sqrt(-ve) = NaN    ... with invalid signal


  sqrt(NaN) = NaN    ... with invalid signal for signaling NaN

Definition at line 119 of file sqrt.c.

double tan (double x)

Tangent function.

Returns:

Tangent function of x

Kernel function:



__kernel_sin        ... sine function on [-pi/4,pi/4]
__kernel_cos        ... cose function on [-pi/4,pi/4]
__ieee754_rem_pio2  ... argument reduction routine

Method:



Let S,C and T denote the sin, cos and tan respectively on
[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
in [-pi/4 , +pi/4], and let n = k mod 4.
We have


     n        sin(x)      cos(x)        tan(x)
----------------------------------------------------------


     0          S           C             T


     1          C          -S           -1/T


     2         -S          -C             T


     3         -C           S           -1/T
----------------------------------------------------------

Special cases:



Let trig be any of sin, cos, or tan.
trig(+-INF)  is NaN
trig(NaN)    is that NaN

Accuracy:



TRIG(x) returns trig(x) nearly rounded

Definition at line 87 of file tan.c.

double tanh (double x)

Hyperbolic tangent function.

Returns:

The Hyperbolic Tangent of x

Method


                                   x    -x


                                  e  - e


    0. tanh(x) is defined to be -----------


                                   x    -x


                                  e  + e



    1. reduce x to non-negative by tanh(-x) = -tanh(x)



    2.  0      <= x <= 2**-55 : tanh(x) = x*(one+x)


                                           -t


        2**-55 <  x <=  1     : tanh(x) = -----; t = expm1(-2x)


                                          t + 2


                                                  2


        1      <= x <=  22.0  : tanh(x) = 1-  ----- ; t=expm1(2x)


                                                t + 2


        22.0   <  x <= INF    : tanh(x) = 1

Special cases


    tanh(NaN) is NaN


    only tanh(0)=0 is exact for finite argument.

Definition at line 81 of file tanh.c.

double trunc (double x)

Truncate function.

when x > 0
trunc(0)   = floor(x)

when x < 0
trunc(x)   = ceil(x)

Special case
trunc(0)   = 0
trunc(NaN) = NaN

Definition at line 52 of file trunc.c.

double vcs (double x)

Versed cosine function.

vercos = 1+cos(x)


       = 2*cos(x/2)*cos(x/2)

Definition at line 45 of file vcs.c.

double ver (double x)

Versed sine function.

versin = 1-cos(x)


       = 2*sin(x/2)*sin(x/2)

Definition at line 45 of file ver.c.

HOMEPAGE

https://amath.innolan.net/

AUTHORS

The library is based on fdlib by Sun Microsystems. The original library can be downloaded at netlib.org:
http://www.netlib.org/fdlibm/

or from mirror site hensa.ac.uk:
http://www.hensa.ac.uk/