math_next.h

Go to the documentation of this file.

#ifndef _MATH_NEXT_H_
#define _MATH_NEXT_H_

#include <nlibc.h>

/*==================================================================
!! Syntax extensions for mathematical functions on apeNEXT
!!
!! All constants changed to records which are defined
!! and initialized in mathvar.hzt
!!
!! $Id: math_next.h,v 2.7 2005/06/08 14:06:53 pleiter Exp $
!! from : math_next.hzt,v 1.3 2001/11/02 09:28:41 simma Exp
!!
!! ZZ version : Roberto De Pietri: v. 0.4 26/08/2001
!! C conversion : N. Paschedag : 15/02/2002
!!
!! Implemented functions:
!!
!!   ->  sqrt(x)     <- 
!!   ->  isqrt(x)    <- 1/sqrt(x)
!!   ->  exp(x)      <- e^x
!!   ->  exp1m(x)    <- e^x-1
!!   ->  sinh(x)     <- 
!!   ->  cosh(x)     <- 
!!   ->  tanh(x)     <- 
!!   ->  log(x)      <- log(x)   \__ / Using a degree 33 approximations
!!   ->  log1p(x)    <- log(1+x) /   \ of log(1+y) for -0.5<y<0.5
!!   ->  sin_1oct(x) <- sin(x), -pi/4<x<pi/4  
!!   ->  sin_1oct(x) <- cos(x), -pi/4<x<pi/4  
!!   ->  sin(x)      <- 
!!   ->  cos(x)      <- 
!!   ->  tan(x)      <- 
!!   ->  asin_k(x)   <- asin(x), -1/2<x<1/2 
!!   ->  asin(x)     <- 
!!   ->  acos(x)     <- 
!!   ->  atan_k(x)   <- atan(x), -1<x<1 
!!   ->  atan(x)     <- 
!!
!!==================================================================*/


#ifdef _INLINE_MATH_
#  define math_inline inline
#else
#  define math_inline 
#endif

/*======================================================
!! macro to set a double from its hex representation
!!====================================================*/

#define _math_setdouble(U,D) { \
  register union { double d; unsigned i; } u; \
  u.i=U; D=u.d; \
}

/*======================================================
!!  Function inverted square root (Newton iterations)
!!====================================================*/

math_inline double isqrt( double a ) {
    register double ris, tmp_a2;
    register double three_half, one_half;
    register double aR;

    _math_setdouble(0x3FE0000000000000ULL,one_half);    // 0.5
    _math_setdouble(0x3FF8000000000000ULL,three_half);  // 1.5

    aR = a;
    asm("\tlutinvsqrt.L %0 %1" : "=r" (ris) : "r" (aR));
    asm("\tlutcross.H %0 $ZERO" : "=r" (ris));

    tmp_a2 = aR * one_half;
    ris = ris*( three_half - tmp_a2 * ris * ris);
    ris = ris*( three_half - tmp_a2 * ris * ris);
    ris = ris*( three_half - tmp_a2 * ris * ris);
    ris = ris*( three_half - tmp_a2 * ris * ris);
    return ris;
}

/*======================================================
!!  Function square root (Newton iterations)
!!====================================================*/

inline double sqrt( double a ) {
    register double ris, tmp_a2;
    register double three_half, one_half;
    register double aR;

    _math_setdouble(0x3FE0000000000000ULL,one_half);   // 0.5
    _math_setdouble(0x3FF8000000000000ULL,three_half); // 1.5

    aR = a;
    ris = 0.0;
    where ( aR ) {
        asm("\tlutisqrt.L %0 %1" : "=r" (ris) : "r" (aR));
        asm("\tlutcross.H %0 $ZERO" : "=r" (ris));
    }
    tmp_a2 = aR * one_half;
    ris = ris*( three_half - tmp_a2 * ris * ris);
    ris = ris*( three_half - tmp_a2 * ris * ris);
    ris = ris*( three_half - tmp_a2 * ris * ris);
    ris = ris*( three_half - tmp_a2 * ris * ris);
    ris *= aR;
    return ris;
}


/*-----------------------------------------------------------------
!!
!! The invsqrt function compute 1/Sqrt(x) using a table based guess  
!! y0 = lutisqrt(x) such that r = 1 - x y0^2 is |r|<2^(-6.72)
!!
!! 1/Sqrt[x] = y0 ( 1 + r Q(r) )
!!  
!!        (    1          )    1
!! Q(r) = ( --------- - 1 ) x ---
!!        ( Sqrt[1-r]     )    r
!!
!! Q(r) has been approximated by a Remez minimax polynomial:
!!
!! degree = 5
!! [a,b]  = [-0.01104854345604,0.01104854345604]
!! w(x)   = r Sqrt(1-r)
!! N Iter = 10000 
!! Max Error = 6.6*10^-17
!!
!!---------------------------------------------------------------*/


math_inline double isqrt2( double x ) {
    register double xx,y0,q0,q1;
    register double r,r2,Q;
    register double t1,t2,t3;
    register double p0,p1,p2,p3,p4,p5;

    // -- Remez coefficients for 1/sqrt(x) -------------------------------
    _math_setdouble( 0x3fe0000000000176ULL, p0 );
    _math_setdouble( 0x3fd8000000000640ULL, p1 );
    _math_setdouble( 0x3fd3fffffd11f1d0ULL, p2 );
    _math_setdouble( 0x3fd17ffffbfd0ec7ULL, p3 );
    _math_setdouble( 0x3fcf81774c2db32fULL, p4 );
    _math_setdouble( 0x3fcce1962db980e2ULL, p5 );
    
    xx = x;
    
    asm("\tlutisqrt.L %0 %1" : "=r" (y0) : "r" (xx));
    asm("\tlutcross.H %0 $ZERO" : "=r" (y0));
    
    q0 = y0 * y0;
    r  = (1.0) - xx * q0;
    r2 = r * r;
    
    t1 = p5 * r + p4;
    t2 = p3 * r + p2;
    t3 = p1 * r + p0;
    
    t1 = t1 * r2 + t2;
    Q  = t1 * r2 + t3;
    q1 = r  * y0;
    Q  = q1 * Q + y0;

    return Q;
}

math_inline double sqrt2( double x ) { 
    register double xx,y0,q0,q1;
    register double r,r2,Q;
    register double t1,t2,t3;
    register double p0,p1,p2,p3,p4,p5;

    // -- Remez coefficients for 1/sqrt(x) -------------------------------
    _math_setdouble( 0x3fe0000000000176ULL, p0 );
    _math_setdouble( 0x3fd8000000000640ULL, p1 );
    _math_setdouble( 0x3fd3fffffd11f1d0ULL, p2 );
    _math_setdouble( 0x3fd17ffffbfd0ec7ULL, p3 );
    _math_setdouble( 0x3fcf81774c2db32fULL, p4 );
    _math_setdouble( 0x3fcce1962db980e2ULL, p5 );
    
    xx = x;
    y0 = xx;

    where ( xx ) {
        asm("\tlutisqrt.L %0 %1" : "=r" (y0) : "r" (xx));
        asm("\tlutcross.H %0 $ZERO" : "=r" (y0));
    }

    q0 = y0 * y0;
    r  = 1.0 - xx * q0;
    r2 = r * r;
    
    y0 = y0 * xx;
    q1 = r  * y0;
    
    t1 = p5 * r + p4;
    t2 = p3 * r + p2;
    t3 = p1 * r + p0;
    
    t1 = t1 * r2 + t2;
    Q  = t1 * r2 + t3;
    
    Q  = q1 * Q + y0;

    return Q;
}


/*==================================================================
!!==================================================================
!! Exponential and Hyperbolic functions
!!
!! Exp(x) = 2^(1/log[2) b) 
!!        = 2^f 2^(1/Log(2) x -f ) 
!!        = P 2^q'                   where  q' = 1/Log[2] x - f
!!        = P Exp(q)                    and    q  = x - Log[2] f
!!
!! Since 0 < |q'| < 0.5  we will have  0 < |q| < Log[2]/2 ~ 0.347 
!!
!!                   1
!! Exp[q] = 1 + q + --- q^2 + q^3 P(q) )
!!                   2
!!
!! ---------------------------------------------------------------
!!
!!        Exp[q] - 1 - q - q^2/2
!! P(q) = -----------------------
!!                q^3
!! P(q) has been approximated by a Remez minimax polynomial:
!!
!! degree = 8 
!! [a,b]  = [-0.34657359027997,0.34657359027997]
!! w(q)   = q^3/(Exp[q]). 
!! N Iter = 20000 
!! Max Error = 5.3*10^-18  
!!
!!---------------------------------------------------------------*/



/*======================================================
!! Function exponential (CORE) 
!!====================================================*/
math_inline double exp( double x ) { 

    register double PP,QQ,Qh,res;
    register double b,z,f;
    register double q,qlead,qtrail;
    register double q2,q4,q5;
    register double t1,t2,t3;
    register double t1h,t2h,t3h;
    register double log2e,ln2lead,ln2trail,nummag;
    register double p0,p1,p2,p3,p4,p5,p6,p7,p8;
    register double half;

    _math_setdouble(  0x3fe0000000000000 , half );

    _math_setdouble(  0x3ff71547652b82fe, log2e    );
    _math_setdouble(  0x3fe62e42fefa4000, ln2lead  ); 
    _math_setdouble(  0xbd48432a1b0e2634, ln2trail ); 
    _math_setdouble(  0x43300000000003ff, nummag   ); 

    // -- Remez coefficients for exp -------------------------------
    _math_setdouble(  0x3fc5555555555502 , p0 );
    _math_setdouble(  0x3fa555555555320b , p1 );
    _math_setdouble(  0x3f81111111128589 , p2 );
    _math_setdouble(  0x3f56c16c17cc1b0b , p3 );
    _math_setdouble(  0x3f2a01a011cdf3cb , p4 );     
    _math_setdouble(  0x3efa019ab32afe96 , p5 );     
    _math_setdouble(  0x3ec71df5419f16f5 , p6 );     
    _math_setdouble(  0x3e9289f84ba3248e , p7 );     
    _math_setdouble(  0x3e5ad2147cc869b8 , p8 );

    b = log2e * x;
    z = nummag + b;
    f = z - nummag;

    asm("\tlutexp.L %0 %1" : "=r" (PP) : "r" (z));
    asm("\tlutcross.H %0 $ZERO" : "=r" (PP));

    qlead  = x - f * ln2lead;
    qtrail =    - f * ln2trail;
    q  = qlead + qtrail;
    q2 = q  * q;
    q4 = q2 * q2;
    q5 = q  * q4;

    t1h = p7 + q  * p8;
    t2h = p4 + q  * p5;
    t3h = p6 + q  * t1h;
    Qh  = t2h + q2 * t3h;

    t1 = p2   + q  * p3;
    t2 = half + q  * p0;
    t3 = p1   + q  * t1;
    QQ = t2    + q2 * t3;

    QQ = QQ + q5 * Qh;
    QQ = qtrail + q2 * QQ;
    QQ = qlead  + QQ;

    res = PP + PP * QQ;

    return res;
}

/*======================================================
!! Function exponential minus 1 (CORE) 
!!====================================================*/
math_inline double exp1m( double x ) {

    register double PP, Pm1, QQ, Qh, res;
    register double b, z, f;
    register double q, qlead, qtrail;
    register double q2, q4, q5;
    register double t1, t2, t3;
    register double t1h, t2h, t3h;
    register double log2e, ln2lead, ln2trail, nummag;
    register double p0, p1, p2, p3, p4, p5, p6, p7, p8;
    register double half;

    _math_setdouble(  0x3fe0000000000000 , half    );
 
    _math_setdouble(  0x3ff71547652b82fe, log2e    );
    _math_setdouble(  0x3fe62e42fefa4000, ln2lead  ); 
    _math_setdouble(  0xbd48432a1b0e2634, ln2trail ); 
    _math_setdouble(  0x43300000000003ff, nummag   ); 

    // -- Remez coefficients for exp -------------------------------
    _math_setdouble(  0x3fc5555555555502 , p0 );
    _math_setdouble(  0x3fa555555555320b , p1 );
    _math_setdouble(  0x3f81111111128589 , p2 );
    _math_setdouble(  0x3f56c16c17cc1b0b , p3 );
    _math_setdouble(  0x3f2a01a011cdf3cb , p4 );     
    _math_setdouble(  0x3efa019ab32afe96 , p5 );     
    _math_setdouble(  0x3ec71df5419f16f5 , p6 );     
    _math_setdouble(  0x3e9289f84ba3248e , p7 );     
    _math_setdouble(  0x3e5ad2147cc869b8 , p8 );

    b = log2e * x;
    z = nummag + b;
    f = z - nummag;

    asm("\tlutexp.L %0 %1" : "=r" (PP) : "r" (z));
    asm("\tlutcross.H %0 $ZERO" : "=r" (PP));

    Pm1 = PP - 1.0;

    qlead  = x - f * ln2lead;
    qtrail =    - f * ln2trail;
    q  = qlead + qtrail;
    q2 = q  * q;
    q4 = q2 * q2;
    q5 = q  * q4;

    t1h = p7 + q  * p8;
    t2h = p4 + q  * p5;
    t3h = p6 + q  * t1h;
    Qh  = t2h + q2 * t3h;

    t1 = p2   + q  * p3;
    t2 = half + q  * p0;
    t3 = p1   + q  * t1;
    QQ = t2   + q2 * t3;

    QQ = QQ + q5 * Qh;
    QQ = qtrail + q2 * QQ;
    QQ = qlead  + QQ;

    res = Pm1 + PP * QQ;

    return res;
}


/*======================================================
!! Function sinh(x) 
!!====================================================*/
math_inline double sinh( double x ) {

    register double xm, expx, expxm, half, res;

    half = 0.5;
    xm = - x;

    expx  = exp1m(x);
    expxm = exp1m(xm);

    res = half * expx - half * expxm;

    return res;
}

/*======================================================
!! Function cosh(x) 
!!====================================================*/
math_inline double cosh( double x ) {

    register double xm, expx, expxm, half, res;

    half = 0.5;
    xm = - x;

    expx  = exp1m( x );
    expxm = exp1m( xm );

    res = half * expx + half * expxm + 1.0;

    return res;

}

/*======================================================
!! Function tanh(x) 
!!====================================================*/
math_inline double tanh( double xxx ) {

    register double x, x2, exp2x, two, temp;
    register double res;
    two = 1.0 + 1.0;

    x  = xxx;
    x2 = x;
    where ( x < 0.0 ) {
      x2 = -x;
    }

    x2  = two * x2;
    exp2x  = exp1m(x2);

    temp = exp2x + 2.0;
    temp = 1.0 / temp;

    res = exp2x * temp;
    where ( x2 > two ) {
      res = 1.0 - two * temp;
    }

    where ( x < 0.0 ) {
      res = - res;
    }

    return res;
}


/*-------------------------------------------------------
!!
!! Log(x) = e log(2) + log(x0)
!!        = e log(2) + log(1+r) 
!!        = e log(2) + r + r^2 Q(r)
!! 
!! where:  e   = lutloge(x)     x = 2^e x0
!!         x0  = lutlogm(x),    0.5 < x0 < 1.5
!!         r   = x0 - 1  
!!
!! --------------------------------------
!!        Log(1+r) - r
!! Q(r) = ------------- 
!!            r^2
!! has been approximated by a Remez minimax polynomial:
!! ==> p0 + p1 r + p2 r^2 + ..... + p36 r^36
!!
!! degree = 36
!! [a,b]  = [-0.5,0.5]
!! w(r)   = r^2/(Log[1+r]). 
!! N Iter = 20000 
!! Max Error = 1.6*10^-18  
!!-----------------------------------------------------*/



//-----------------------------------------------------------------

math_inline double log( double x ) {
    register double xx, x0, e, e1, r, r2, r4;
    register double P, Pa, Pb, Pc, Pd;
    register double trail;
    register double log2e, ln2lead, ln2trail, nummag;

    register double p0, p1, p2, p3, p4, p5, p6, p7, p8, p9;
    register double p10, p11, p12, p13, p14, p15, p16, p17, p18, p19;
    register double p20, p21, p22, p23, p24, p25, p26, p27, p28, p29;
    register double p30, p31, p32, p33, p34, p35, p36;

    _math_setdouble(  0x3ff71547652b82fe, log2e    );
    _math_setdouble(  0x3fe62e42fefa4000, ln2lead  ); 
    _math_setdouble(  0xbd48432a1b0e2634, ln2trail ); 
    _math_setdouble(  0x43300000000003ff, nummag   ); 

    // -- Remez coefficients for log -------------------------------
    _math_setdouble(  0xbfe0000000000000 , p0  );
    _math_setdouble(  0x3fd5555555555555 , p1  );
    _math_setdouble(  0xbfcfffffffffffd2 , p2  );
    _math_setdouble(  0x3fc9999999999944 , p3  );
    _math_setdouble(  0xbfc555555555839f , p4  );
    _math_setdouble(  0x3fc24924924994cb , p5  );
    _math_setdouble(  0xbfbfffffffd3792a , p6  );
    _math_setdouble(  0x3fbc71c71bfbb905 , p7  );
    _math_setdouble(  0xbfb99999a6067f7a , p8  );
    _math_setdouble(  0x3fb745d1960bf8dd , p9  );
    _math_setdouble(  0xbfb555531a405781 , p10 );
    _math_setdouble(  0x3fb3b1351bb667d8 , p11 );
    _math_setdouble(  0xbfb2496af8fbbbe4 , p12 );
    _math_setdouble(  0x3fb111c6afc8c6e1 , p13 );
    _math_setdouble(  0xbfaff38dd41057ce , p14 );
    _math_setdouble(  0x3fadffb7d7a94696 , p15 );
    _math_setdouble(  0xbfad41ba8a134738 , p16 );
    _math_setdouble(  0x3faccbc9726ece12 , p17 );
    _math_setdouble(  0xbf9ec38f119d3d06 , p18 );
    _math_setdouble(  0x3f785f81962a50fb , p19 );
    _math_setdouble(  0xbfce4be7a8a5d72b , p20 );
    _math_setdouble(  0x3fd9d8d252726e6c , p21 );
    _math_setdouble(  0x3ff5c8aa212fd722 , p22 );
    _math_setdouble(  0xc0026bd5ff44cd4e , p23 );
    _math_setdouble(  0xc01fa5187761a8ad , p24 );
    _math_setdouble(  0x4026d78e0bf83a6c , p25 );
    _math_setdouble(  0x4040ba1a1b262a4f , p26 );
    _math_setdouble(  0xc044570c25ae2da0 , p27 );
    _math_setdouble(  0xc05a8787b4194713 , p28 );
    _math_setdouble(  0x405a009c5d877e88 , p29 );
    _math_setdouble(  0x406e49c81a0007b0 , p30 );
    _math_setdouble(  0xc0666aa1bfe5d0fe , p31 );
    _math_setdouble(  0xc077926b8e4542e3 , p32 );
    _math_setdouble(  0x406767610205d729 , p33 );
    _math_setdouble(  0x40765b8a9462a4c6 , p34 );
    _math_setdouble(  0xc0564d5567c48d46 , p35 );
    _math_setdouble(  0xc06383c8abcc91fa , p36 );
    
    xx = x;

    asm("\tlutlogm.L %0 %1" : "=r" (x0) : "r" (xx));
    asm("\tlutloget.L %0 %1" : "=r" (e1) : "r" (xx));
    asm("\tlutcross.H %0 $ZERO" : "=r" (x0));
    asm("\tlutcross.H %0 $ZERO" : "=r" (e1));

    e = e1 - nummag;
    r = x0 - 1.0;
    r2 = r  * r;
    r4 = r2 * r2;
    
    Pa = p36;
    Pa = r4 * Pa + p32;
    Pa = r4 * Pa + p28;
    Pa = r4 * Pa + p24;
    Pa = r4 * Pa + p20;
    Pa = r4 * Pa + p16;
    Pa = r4 * Pa + p12;
    Pa = r4 * Pa + p8;
    Pa = r4 * Pa + p4;
    Pa = r4 * Pa + p0;
    
    Pb = p33;
    Pb = r4 * Pb + p29;
    Pb = r4 * Pb + p25;
    Pb = r4 * Pb + p21;
    Pb = r4 * Pb + p17;
    Pb = r4 * Pb + p13;
    Pb = r4 * Pb + p9;
    Pb = r4 * Pb + p5;
    Pb = r4 * Pb + p1;
    
    Pc = p34;
    Pc = r4 * Pc + p30;
    Pc = r4 * Pc + p26;
    Pc = r4 * Pc + p22;
    Pc = r4 * Pc + p18;
    Pc = r4 * Pc + p14;
    Pc = r4 * Pc + p10;
    Pc = r4 * Pc + p6;
    Pc = r4 * Pc + p2;
    
    Pd = p35;
    Pd = r4 * Pd + p31;
    Pd = r4 * Pd + p27;
    Pd = r4 * Pd + p23;
    Pd = r4 * Pd + p19;
    Pd = r4 * Pd + p15;
    Pd = r4 * Pd + p11;
    Pd = r4 * Pd + p7;
    Pd = r4 * Pd + p3;
    
    P = Pa + r * Pb + r2 * (Pc + r * Pd);
    
    trail = e * ln2trail;
    trail = r2 * P + trail;
    trail = trail + r;
    P     = trail + e * ln2lead;

    return P;
}


//-----------------------------------------------------------------

math_inline double log1p( double xxx ) {
    register double x,xx,x0,e,e1,r,r2,r4;
    register double P,Pa,Pb,Pc,Pd;
    register double trail;
    register double log2e,ln2lead,ln2trail,nummag  ;

    register double p0,p1,p2,p3,p4,p5,p6,p7,p8,p9;
    register double p10,p11,p12,p13,p14,p15,p16,p17,p18,p19;
    register double p20,p21,p22,p23,p24,p25,p26,p27,p28,p29;
    register double p30,p31,p32,p33,p34,p35,p36;

    _math_setdouble(  0x3ff71547652b82fe, log2e    );
    _math_setdouble(  0x3fe62e42fefa4000, ln2lead  ); 
    _math_setdouble(  0xbd48432a1b0e2634, ln2trail ); 
    _math_setdouble(  0x43300000000003ff, nummag   ); 

    // -- Remez coefficients for log -------------------------------
    _math_setdouble(  0xbfe0000000000000 , p0  );
    _math_setdouble(  0x3fd5555555555555 , p1  );
    _math_setdouble(  0xbfcfffffffffffd2 , p2  );
    _math_setdouble(  0x3fc9999999999944 , p3  );
    _math_setdouble(  0xbfc555555555839f , p4  );
    _math_setdouble(  0x3fc24924924994cb , p5  );
    _math_setdouble(  0xbfbfffffffd3792a , p6  );
    _math_setdouble(  0x3fbc71c71bfbb905 , p7  );
    _math_setdouble(  0xbfb99999a6067f7a , p8  );
    _math_setdouble(  0x3fb745d1960bf8dd , p9  );
    _math_setdouble(  0xbfb555531a405781 , p10 );
    _math_setdouble(  0x3fb3b1351bb667d8 , p11 );
    _math_setdouble(  0xbfb2496af8fbbbe4 , p12 );
    _math_setdouble(  0x3fb111c6afc8c6e1 , p13 );
    _math_setdouble(  0xbfaff38dd41057ce , p14 );
    _math_setdouble(  0x3fadffb7d7a94696 , p15 );
    _math_setdouble(  0xbfad41ba8a134738 , p16 );
    _math_setdouble(  0x3faccbc9726ece12 , p17 );
    _math_setdouble(  0xbf9ec38f119d3d06 , p18 );
    _math_setdouble(  0x3f785f81962a50fb , p19 );
    _math_setdouble(  0xbfce4be7a8a5d72b , p20 );
    _math_setdouble(  0x3fd9d8d252726e6c , p21 );
    _math_setdouble(  0x3ff5c8aa212fd722 , p22 );
    _math_setdouble(  0xc0026bd5ff44cd4e , p23 );
    _math_setdouble(  0xc01fa5187761a8ad , p24 );
    _math_setdouble(  0x4026d78e0bf83a6c , p25 );
    _math_setdouble(  0x4040ba1a1b262a4f , p26 );
    _math_setdouble(  0xc044570c25ae2da0 , p27 );
    _math_setdouble(  0xc05a8787b4194713 , p28 );
    _math_setdouble(  0x405a009c5d877e88 , p29 );
    _math_setdouble(  0x406e49c81a0007b0 , p30 );
    _math_setdouble(  0xc0666aa1bfe5d0fe , p31 );
    _math_setdouble(  0xc077926b8e4542e3 , p32 );
    _math_setdouble(  0x406767610205d729 , p33 );
    _math_setdouble(  0x40765b8a9462a4c6 , p34 );
    _math_setdouble(  0xc0564d5567c48d46 , p35 );
    _math_setdouble(  0xc06383c8abcc91fa , p36 );

    x  = xxx;
    xx = x + 1.0;

    asm("\tlutlogm.L %0 %1" : "=r" (x0) : "r" (xx));
    asm("\tlutloget.L %0 %1" : "=r" (e1) : "r" (xx));
    asm("\tlutcross.H %0 $ZERO" : "=r" (x0));
    asm("\tlutcross.H %0 $ZERO" : "=r" (e1));

    e = e1 - nummag;
    r = x0 - 1.0;
    
    where ( e == 0.0 ) {
        r = x;
    }
    
    r2 = r  * r;
    r4 = r2 * r2;

    Pa = p36;
    Pa = r4 * Pa + p32;
    Pa = r4 * Pa + p28;
    Pa = r4 * Pa + p24;
    Pa = r4 * Pa + p20;
    Pa = r4 * Pa + p16;
    Pa = r4 * Pa + p12;
    Pa = r4 * Pa + p8;
    Pa = r4 * Pa + p4;
    Pa = r4 * Pa + p0;

    Pb = p33;
    Pb = r4 * Pb + p29;
    Pb = r4 * Pb + p25;
    Pb = r4 * Pb + p21;
    Pb = r4 * Pb + p17;
    Pb = r4 * Pb + p13;
    Pb = r4 * Pb + p9;
    Pb = r4 * Pb + p5;
    Pb = r4 * Pb + p1;

    Pc = p34;
    Pc = r4 * Pc + p30;
    Pc = r4 * Pc + p26;
    Pc = r4 * Pc + p22;
    Pc = r4 * Pc + p18;
    Pc = r4 * Pc + p14;
    Pc = r4 * Pc + p10;
    Pc = r4 * Pc + p6;
    Pc = r4 * Pc + p2;

    Pd = p35;
    Pd = r4 * Pd + p31;
    Pd = r4 * Pd + p27;
    Pd = r4 * Pd + p23;
    Pd = r4 * Pd + p19;
    Pd = r4 * Pd + p15;
    Pd = r4 * Pd + p11;
    Pd = r4 * Pd + p7;
    Pd = r4 * Pd + p3;

    P = Pa + r * Pb + r2 * (Pc + r * Pd);

    trail = e * ln2trail;
    trail = r2 * P + trail;
    trail = trail + r;
    P = trail + e * ln2lead;

    return P;
}


/*-----------------------------------------------------------------
!!
!! Kernel function for the computation of Sin[x]
!!
!! Sin[x] = x ( 1 + x^2 Q[x^2] )    ;   -pi/4 < x < pi/4
!!
!!---------------------------------------------------------------*/

math_inline double sin_1oct( double xx ) {
    register double x, x2s, x4s, Qs;
    register double t1s, t2s, t3s;
    register double Sp0, Sp1, Sp2, Sp3, Sp4, Sp5, Sp6;

    // -- Remez coefficients for kernel sin -------------------------------
    _math_setdouble(  0x0 , Sp0 );
    _math_setdouble(  0xbfc5555555555555 , Sp1 );
    _math_setdouble(  0x3f81111111110eb8 , Sp2 );
    _math_setdouble(  0xbf2a01a019f1c947 , Sp3 );
    _math_setdouble(  0x3ec71de384036e7d , Sp4 );
    _math_setdouble(  0xbe5ae60a561eeab5 , Sp5 );
    _math_setdouble(  0x3de5e3c6b7eeb28d , Sp6 );

    x = xx;
    x2s  = x * x;
    x4s = x2s * x2s;
    
    t1s = Sp5 + x2s * Sp6;
    t2s = Sp3 + x2s * Sp4;
    t3s = Sp1 + x2s * Sp2;
    
    t1s = t2s + x4s * t1s;
    Qs  = t3s + x4s * t1s;
    
    Qs  = x2s * Qs;
    Qs  = x + Qs * x;
    
    return Qs;
    
}

/*-----------------------------------------------------------------
!!
!!
!! Kernel function for the computation of Cos[x]
!!
!! Cos[x] = ( 1 + x^2 (1/2 + x^2 Q[x^2] ) )  ;   -pi/4 < x < pi/4
!!
!!---------------------------------------------------------------*/

math_inline double cos_1oct( double x ) { 

    register double x2c, x4c, Qc;
    register double t1c, t2c, t3c;
    register double Cp0, Cp1, Cp2, Cp3, Cp4, Cp5, Cp6;

    // -- Remez coefficients for kernel cos -------------------------------
    _math_setdouble(  0xbfe0000000000000 , Cp0 );
    _math_setdouble(  0x3fa5555555555555 , Cp1 );
    _math_setdouble(  0xbf56c16c16c16910 , Cp2 );
    _math_setdouble(  0x3efa01a019f556d1 , Cp3 );
    _math_setdouble(  0xbe927e4fa28f90c6 , Cp4 );
    _math_setdouble(  0x3e21eeb7c6903ba2 , Cp5 );
    _math_setdouble(  0xbda908b4ef9a7e2e , Cp6 );

    x2c  = x * x;
    x4c = x2c * x2c;

    t1c = Cp5 + x2c * Cp6;
    t2c = Cp3 + x2c * Cp4;
    t3c = Cp1 + x2c * Cp2;

    t1c = t2c + x4c * t1c;
    Qc  = t3c + x4c * t1c;

    Qc  = x2c * Qc + Cp0;
    Qc  = (1.0) + x2c * Qc;

    return Qc;

}


/*-----------------------------------------------------------------
!! sin(x)
!!---------------------------------------------------------------*/

math_inline double sin( double x ) { 

    register double y, siny, cosy;
    register double res;

    register union { 
      double z;
      int    zint;
    } zu;

    register double f;
    register int k;
    register double twooverpi, pihalflead, pihalftrail1, pihalftrail2, magic2;
       
    _math_setdouble(  0x3fe45f306dc9c883 , twooverpi   );
    _math_setdouble(  0x4337ffffffffffff , magic2      );
    _math_setdouble(  0x3ff921fb50000000 , pihalflead  );
    _math_setdouble(  0x3e5110b460000000 , pihalftrail1 );
    _math_setdouble(  0x3c91a62633145c07 , pihalftrail2 );

    zu.z = magic2 + twooverpi * x;
    f = zu.z - magic2;
    y = x - f * pihalflead;
    y = y - f * pihalftrail1;
    y = y - f * pihalftrail2;

    // -----------------------------------------
    // Comute the sin and the cos of y for (-pi/4<y<pi/4) 
    // -----------------------------------------
    
    siny = sin_1oct(y);
    cosy = cos_1oct(y);
    
    /* -----------------------------------------
    !! get the module 4 of _f  The modul is get 
    !! over z that has mantissa _f+0x7fff...f
    !! That means _k = (f%4 + 3)%4
    !! f = 0 -> k= 3 
    !! f = 1 -> k= 0 
    !! f = 2 -> k= 1 
    !! f = 3 -> k= 2 
    !! ---------------------------------------*/
    k = 0x3;
    k = zu.zint & k;

    where ( k == 3 ) {
      res = siny;
    }
    where ( k == 0 ) {
      res = cosy;
    }
    where ( k == 1 ) {
      res = - siny;
    }
    where ( k == 2 ) {
      res = - cosy;
    }
    
    return res;
    
}

/*-----------------------------------------------------------------
!! cos(x)
!!---------------------------------------------------------------*/

math_inline double cos( double x ) {
    register double y, siny, cosy;
    register double res;

    register double f;

    register union { 
      double z;
      int zint;
    } zu;

    register int k;
    register double twooverpi, pihalflead, pihalftrail1, pihalftrail2, magic2;
       
    _math_setdouble(  0x3fe45f306dc9c883 , twooverpi   );
    _math_setdouble(  0x4337ffffffffffff , magic2      );
    _math_setdouble(  0x3ff921fb50000000 , pihalflead  );
    _math_setdouble(  0x3e5110b460000000 , pihalftrail1 );
    _math_setdouble(  0x3c91a62633145c07 , pihalftrail2 );

    zu.z = magic2 + twooverpi * x;
    f = zu.z - magic2;
    y = x - f * pihalflead;
    y = y - f * pihalftrail1;
    y = y - f * pihalftrail2;
    
    // Comute the sin and the cos of (y)
    
    siny = sin_1oct(y);
    cosy = cos_1oct(y);
    
    /* -----------------------------------------
    !! get the module 4 of _f  The modul is get 
    !! over z that has mantissa _f+0x7fff...f
    !! That means _k = (_f%4 + 3)%4
    !! _f = 0 -> _k= 3 
    !! _f = 1 -> _k= 0 
    !! _f = 2 -> _k= 1 
    !! _f = 3 -> _k= 2 
    !! ---------------------------------------*/    
    k = 0x3;
    k = zu.zint & k ;
    
    where ( k == 3 ) {
      res = cosy;
    }
    where ( k == 0 ) {
      res = - siny;
    }
    where ( k == 1 ) {
      res = - cosy;
    }
    where ( k == 2 ) {
      res = siny;
    }
    
    return res;
    
}

/*==================================================================
!!
!! The argument of the function Tan(x) is first reduce to y such that
!!
!! x = k pi/2 + y,   k integer and -pi/4 < y < pi/4  
!!
!! and then:
!!
!! -pi/4 + k pi < x < pi/4 + k pi :   Tan(x) = Tan(y)     
!! -pi/4 + k pi < x < pi/4 + k pi :   Tan(x) = - Cotan(y)
!!
!! Tan(y) and CoTan(y) are computed according to:
!! 
!! -----------------------------------------------------------------
!!
!! P = Tan(y) = y + y^3 P13(y^2) 
!!
!! Where P22(r) is a degree 13 minimax polinomial approximating
!! 
!!        Tan(Sqrt(r))-Sqrt(r)
!! P(r) = ---------------------
!!            Sqrt(r) r 
!!
!! degree = 13
!! [a,b]  = [0,0.61685027506808]
!! w(r)   = r^3/(Tan[r]). 
!! N Iter = 20000 
!! Max Error = 5.3*10^-18  
!!
!! -----------------------------------------------------------------
!!      1
!! Q = --- - CoTan(y) = y + y^3 Q10(y^2) 
!!      y
!!
!! Where Q22(t) is a degree 10 minimax polinomial approximating
!! 
!!        1/Sqrt[r] - CoTan(Sqrt[r]) - Sqrt[r]
!! Q[r] = ------------------------------------
!!                 Sqrt[r] r
!!
!! degree = 10
!! [a,b]  = [0,0.61685027506808]
!! w(r)   = r^3/( 1/Sqrt(r) - CoTan(sqrt(r)) ). 
!! N Iter = 2000 
!! Max Error = 2.36*10^-20  
!!===============================================================*/

math_inline double tan( double x ) {

    register double y, y2, y3, y4, invy;
    register double res, Q, P, Qn, Qp, Pn, Pp;
    
    register double p0, p1, p2, p3, p4, p5, p6, p7, p8, p9;
    register double p10, p11, p12, p13;
    register double q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10;

    register double f;

    register union { 
     double z;
     int    zint;
    } zu;

    register int k;
    register double twooverpi,pihalflead,pihalftrail1,pihalftrail2,magic2;
       
    _math_setdouble(  0x3fe45f306dc9c883 , twooverpi   );
    _math_setdouble(  0x4337ffffffffffff , magic2      );
    _math_setdouble(  0x3ff921fb50000000 , pihalflead  );
    _math_setdouble(  0x3e5110b460000000 , pihalftrail1 );
    _math_setdouble(  0x3c91a62633145c07 , pihalftrail2 );
  
    // -- Remez coefficients for tan -------------------------------
    _math_setdouble(  0x3fd55555555554f2 , p0  );
    _math_setdouble(  0x3fc111111111823f , p1  );
    _math_setdouble(  0x3faba1ba1b472c71 , p2  );
    _math_setdouble(  0x3f9664f49a7087a7 , p3  );
    _math_setdouble(  0x3f8226e1281741ed , p4  );
    _math_setdouble(  0x3f6d6d92e4dc1ec9 , p5  );
    _math_setdouble(  0x3f57d5b9d094e180 , p6  );
    _math_setdouble(  0x3f437f87931093e9 , p7  );
    _math_setdouble(  0x3f2d28899e55dabc , p8  );
    _math_setdouble(  0x3f21df93999dc111 , p9  );
    _math_setdouble(  0xbefb6ea2534187cb , p10 );
    _math_setdouble(  0x3f17656fc1431ef6 , p11 );
    _math_setdouble(  0xbf0778cb8106dd3d , p12 );
    _math_setdouble(  0x3ef5d99c5b37b8fb , p13 );
    // -- Remez coefficients for cotan -------------------------------    
    _math_setdouble(  0x3fd5555555555555 , q0  );
    _math_setdouble(  0x3f96c16c16c16c16 , q1  );
    _math_setdouble(  0x3f61566abc011734 , q2  );
    _math_setdouble(  0x3f2bbd7793321936 , q3  );
    _math_setdouble(  0x3ef66a8f2d1bc68f , q4  );
    _math_setdouble(  0x3ec228059183e28c , q5  );
    _math_setdouble(  0x3e8d6dc6da8b4b97 , q6  );
    _math_setdouble(  0x3e57d86d9f6cba62 , q7  );
    _math_setdouble(  0x3e23722fc10fb082 , q8  );
    _math_setdouble(  0x3ded3f62bcb56407 , q9  );
    _math_setdouble(  0x3dc2113add876256 , q10 );

    zu.z = magic2 + twooverpi * x;
    f = zu.z - magic2;
    y = x - f * pihalflead;
    y = y - f * pihalftrail1;
    y = y - f * pihalftrail2;
    
    y2 = y * y;
    y3 = y * y2;
    y4 = y2 * y2;
    
    invy = y;
    where ( invy == 0.0) {
       invy = 1.0;
    }
    invy = 1.0 / invy;
    
    Pp = p13;
    Pn = p12;
    Pp = y4*Pp + p11;
    Pn = y4*Pn + p10;
    Pp = y4*Pp + p9;
    Pn = y4*Pn + p8;
    Pp = y4*Pp + p7;
    Pn = y4*Pn + p6;
    Pp = y4*Pp + p5;
    Pn = y4*Pn + p4;
    Pp = y4*Pp + p3;
    Pn = y4*Pn + p2;
    Pp = y4*Pp + p1;
    Pn = y4*Pn + p0;
    P = y + y3 *(Pn+y2*Pp);
    
    Qn =  q10;
    Qp =  q9;
    Qn = y4*Qn + q8;
    Qp = y4*Qp + q7;
    Qn = y4*Qn + q6;
    Qp = y4*Qp + q5;
    Qn = y4*Qn + q4;
    Qp = y4*Qp + q3;
    Qn = y4*Qn + q2;
    Qp = y4*Qp + q1;
    Qn = y4*Qn + q0;
    Q = Qn + y2 * Qp;

    /* -----------------------------------------
    !! get the module 2 of f  The modul is get 
    !! over z that has mantissa f+0x7fff...f
    !! That means k = (f%2 + 1)%2
    !! f = 0 -> k= 1 
    !! f = 1 -> k= 0 
    !! ---------------------------------------*/    
    k = 0x1;
    k = zu.zint & k;
    
    where ( k  ==  0 ) {
        P = Q * y  - invy;
    }

    return P;

}


/*===================================================================
!! Kernel function for the computation Atan(x) for |x|<1
!!
!! Res = Atan(x) = x + x^3 P22(x^2) 
!!
!! Where P22(y) is a degree 22 minimax polinomial approximating
!! 
!!        Atan(Sqrt[y])-Sqrt[y]
!! P[y] = ---------------------
!!            Sqrt[y] y
!!
!! degree = 22 
!! [a,b]  = [0,1.0]
!! w(y)   = y^3 / (Atan[Sqrt[y]]). 
!! N Iter = 20000 
!! Max Error = 2.38*10^-20  
!!
!!=================================================================*/

math_inline double atan_k( double xx ) {
    register double x, x2, x3, P;
    register double r2, r4, Pa, Pb, Pc, Pd;
    register double p0, p1, p2, p3, p4, p5, p6, p7, p8, p9;
    register double p10, p11, p12, p13, p14, p15, p16, p17, p18, p19;
    register double p20, p21, p22;

    // -- Remez coefficients for atan -------------------------------
    _math_setdouble(  0xbfd5555555555555 , p0  );
    _math_setdouble(  0x3fc9999999999971 , p1  );
    _math_setdouble(  0xbfc2492492490f58 , p2  );
    _math_setdouble(  0x3fbc71c71c68e470 , p3  );
    _math_setdouble(  0xbfb745d173541a68 , p4  );
    _math_setdouble(  0x3fb3b13affee80cf , p5  );
    _math_setdouble(  0xbfb111100b7ee084 , p6  );
    _math_setdouble(  0x3fae1e0a5cf4626e , p7  );
    _math_setdouble(  0xbfaaf1f6218828be , p8  );
    _math_setdouble(  0x3fa85e4efb1caacf , p9  );
    _math_setdouble(  0xbfa634430a602ac3 , p10 );
    _math_setdouble(  0x3fa446067f9429b9 , p11 );
    _math_setdouble(  0xbfa25977aa870234 , p12 );
    _math_setdouble(  0x3fa0269900f0cc1e , p13 );
    _math_setdouble(  0xbf9ae16c1b259e80 , p14 );
    _math_setdouble(  0x3f946df0d7c219b5 , p15 );
    _math_setdouble(  0xbf8b503895880e16 , p16 );
    _math_setdouble(  0x3f7ee47c1bf074f6 , p17 );
    _math_setdouble(  0xbf6c59076eef60d7 , p18 );
    _math_setdouble(  0x3f54123312f3c886 , p19 );
    _math_setdouble(  0xbf346fc2146a8776 , p20 );
    _math_setdouble(  0x3f0a81635a5fa23b , p21 );
    _math_setdouble(  0xbed063c89c3fdaca , p22 );
    
    x  = xx;
    x2 = x * x;
    x3 = x * x2;
    r2 = x2 * x2;
    r4 = r2 * r2;
    
    Pa = p20;
    Pa = r4 * Pa + p16;
    Pa = r4 * Pa + p12;
    Pa = r4 * Pa + p8;
    Pa = r4 * Pa + p4;
    Pa = r4 * Pa + p0;
    
    Pb = p21;
    Pb = r4 * Pb + p17;
    Pb = r4 * Pb + p13;
    Pb = r4 * Pb + p9;
    Pb = r4 * Pb + p5;
    Pb = r4 * Pb + p1;
    
    Pc =  p22;
    Pc = r4 * Pc + p18;
    Pc = r4 * Pc + p14;
    Pc = r4 * Pc + p10;
    Pc = r4 * Pc + p6;
    Pc = r4 * Pc + p2;
    
    Pd = p19;
    Pd = r4 * Pd + p15;
    Pd = r4 * Pd + p11;
    Pd = r4 * Pd + p7;
    Pd = r4 * Pd + p3;
    
    P = Pa + x2 * Pb + r2 * (Pc + x2 * Pd);
    P = x + x3 *P;

    return P;
}

/* ===================================================================
!! 
!! The computation of Atan(x) is compute according to
!!
!!  1 < x      : Atan(x) =  pi/2 - Atan_k(1/x)
!! -1 < x <  1 : Atan(x) = Atan_k(x)
!!      x < -1 : Atan(x) = -pi/2 - Atan_k(1/x)
!! 
!! Where P22(y) is a degree 22 minimax polinomial approximating
!! (Atan(x)-x)/x^3 in [-1 1]
!!
!! =================================================================*/

math_inline double atan( double x ) { 
    register double xx, absx, y, res;
    register double pihalflead, pihalftrail;

    _math_setdouble(  0x3ff9220000000000 , pihalflead  );
    _math_setdouble(  0xbed2aeef4b9ee59e , pihalftrail );

    xx   = x;
    absx = xx;
    where ( xx < 0.0 ) {
      absx = - absx;
    }

    y = absx;
    where ( absx > 1.0 ) {
      y = 1.0 / absx;
    }
    
    res = atan_k( y );
    
    where ( absx > 1.0 ) {
      res = pihalftrail - res;
      res = pihalflead  + res;
    }
    where ( xx < 0.0 ) {
      res = - res;
    }
 
    return res;
}


/* -----------------------------------------------------------------
!!
!! Main part: Computation of Asin(x) where -0.5 < x < 0.5
!! 
!! Res = Asin(x) = x + x^3 P19(x^2) 
!!
!! Where P19(y) is a degree 19 minimax polinomial approximating
!! 
!!        Asin(Sqrt[y])-Sqrt[y]
!! P[y] = ---------------------
!!            Sqrt[y] y
!!
!! degree = 19
!! [a,b]  = [0,0.5]
!! w(y)   = 1/(Asin[y]). 
!! N Iter = 20000 
!! Max Error = 4.41*10^-19  
!!
!! ---------------------------------------------------------------*/

math_inline double asin_k( double xx ) { 
    register double x, x2, x3, P;
    register double r2, r4, Pa, Pb, Pc, Pd;
    register double p0, p1, p2, p3, p4, p5, p6, p7, p8, p9;
    register double p10, p11, p12, p13, p14, p15, p16, p17, p18, p19;

    // -- Remez coefficients for asin -------------------------------
    _math_setdouble(  0x3fc5555555555555 , p0  );
    _math_setdouble(  0x3fb3333333333ff8 , p1  );
    _math_setdouble(  0x3fa6db6db6c8299f , p2  );
    _math_setdouble(  0x3f9f1c71d2cb0419 , p3  );
    _math_setdouble(  0x3f96e8b839a6eb1f , p4  );
    _math_setdouble(  0x3f91c5217cd11e53 , p5  );
    _math_setdouble(  0x3f8c91ddc9cccb8b , p6  );
    _math_setdouble(  0x3f880fff7f9aa201 , p7  );
    _math_setdouble(  0x3f7ff0699c7911d9 , p8  );
    _math_setdouble(  0x3f97c8c548d85470 , p9  );
    _math_setdouble(  0xbfb43c911218c604 , p10 );
    _math_setdouble(  0x3fd968af9c41272b , p11 );
    _math_setdouble(  0xbff5ecfed8eb21be , p12 );
    _math_setdouble(  0x400e2c596eac2542 , p13 );
    _math_setdouble(  0xc01fb74a405ddd3f , p14 );
    _math_setdouble(  0x4029446d2f84e0c3 , p15 );
    _math_setdouble(  0xc02d7056e867daf6 , p16 );
    _math_setdouble(  0x4027cbd4a391e704 , p17 );
    _math_setdouble(  0xc017e8e8f8f0fb8c , p18 );
    _math_setdouble(  0x3ff6d9ee9188998b , p19 );
    
    x  = xx;
    x2 = x * x;
    x3 = x * x2;
    
    r2 = x2 * x2;
    r4 = r2 * r2;
    
    Pa = p16;
    Pa = r4 * Pa + p12;
    Pa = r4 * Pa + p8;
    Pa = r4 * Pa + p4;
    Pa = r4 * Pa + p0;
    
    Pb = p17;
    Pb = r4 * Pb + p13;
    Pb = r4 * Pb + p9;
    Pb = r4 * Pb + p5;
    Pb = r4 * Pb + p1;
    
    Pc = p18;
    Pc = r4 * Pc + p14;
    Pc = r4 * Pc + p10;
    Pc = r4 * Pc + p6;
    Pc = r4 * Pc + p2;
    
    Pd = p19;
    Pd = r4 * Pd + p15;
    Pd = r4 * Pd + p11;
    Pd = r4 * Pd + p7;
    Pd = r4 * Pd + p3;
    
    P = Pa + x2 * Pb + r2 * (Pc + x2 * Pd);
    
    P = x + x3 * P;

    return P;
}

/*-----------------------------------------------------------------
!!
!! The ArcSin[x] function is computed in terms of the ArcSin[y]
!! function (with |y|<0.5) according to the following reduction 
!! formula:
!!
!!
!! ( 0.5 < x )     ==>  Asin[x] =  pi/2  - 2 Asin[sqrt[(1-x)/2]]
!! ( -0.5 < x <.5) ==>  Asin[x] = Asin[x]
!! ( x < -0.5 )    ==>  Asin[x] = - pi/2 + 2 Asin[sqrt[(1+x)/2]] )
!! 
!!---------------------------------------------------------------*/

math_inline double asin( double x ) { 
    register double xx, y, absx, res;
    register double half, two;
    register double pihalflead, pihalftrail;

    _math_setdouble(  0x3ff9220000000000 , pihalflead  );
    _math_setdouble(  0xbed2aeef4b9ee59e , pihalftrail );

    half = 0.5;
    two  = 1.0 + 1.0;
    
    xx  = x;
    absx = xx;
    where ( xx < 0.0 ) {
      absx = - absx;
    }
    y = absx ;
    where ( absx > half ) {
      y = half - half * y;
      y = sqrt(y);
    }
    
    res = asin_k(y);
    
    where ( absx > half ) {
      res = pihalftrail - two * res;
      res = pihalflead  + res;
    }
    where ( xx < 0.0 ) {
      res = - res;
    }
 
    return res;
}

/*-----------------------------------------------------------------
!! The ArcCos[x] function is computed in terms of the ArcSin[y]
!! function (with |y|<0.5) according to the following reduction 
!! formula
!!
!! ( 0.5 < x )     ==>  Acos[x] = 2 Asin[sqrt[(1-x)/2]]
!! ( -0.5 < x <.5) ==>  Acos[x] = pi/2 - Asin[x]
!! ( x < -0.5 )    ==>  Acos[x] = 2 (pi/2 - Asin[sqrt[(1+x)/2]] )
!!---------------------------------------------------------------*/
math_inline double acos( double x ) {
    register double xx, y, absx, res;
    register double half, two;
    register double pihalflead, pihalftrail;

    _math_setdouble(  0x3ff9220000000000 , pihalflead  );
    _math_setdouble(  0xbed2aeef4b9ee59e , pihalftrail );

    half = 0.5;
    two  = 1.0 + 1.0;
    
    xx   =x;
    absx = xx;
    where ( xx < 0.0 ) {
      absx = - absx;
    }
    y = absx;
    where ( absx > half ) {
      y = half - half * absx;
      y = sqrt(y);
    }
    
    res = asin_k(y);
    
    where ( xx < 0.0 ) {
      res = - res;
    }
    
    where ( absx > half ) {
      res = two * res;
    }
    where ( absx <= half ) {
      res = pihalftrail - res;
      res = pihalflead  + res;
    }
    where ( xx + half < 0.0) {
      res = two * pihalftrail + res;
      res = two * pihalflead  + res;
    }

    return res;
}


#undef _math_setdouble

#endif

---