math.h

Go to the documentation of this file.

#ifndef _MATH_H_              
#define _MATH_H_

#include <nlibc.h>

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

#define _math_setdouble(U,D) { \
  register union { double d; unsigned i; } u; \
  asm("\tatr0h %0 <"#U">" : "=r" (u.i)); D=u.d; \
}
/*  old version 
#define _math_setdouble(U,D) { \
  register union { double d; unsigned i; } u; \
  u.i=U; D=u.d; \
}
*/

#include <stdint.h>
#include <limits.h>
#include <float.h>


/* ++++++++++++        macros       +++++++++++++++++++*/

#define INFINITY  0x1.0p32767

#define HUGE_VAL        INFINITY


#define HUGE_VALF  HUGE_VAL
#define HUGE_VALL  HUGE_VAL


#undef NAN

#define FP_INFINITE  2
#define FP_NAN       3
#define FP_NORMAL    1
#define FP_SUBNORMAL 1
#define FP_ZERO      4

#define FP_FAST_FMA 

#define FP_FAST_FMAF FP_FAST_FMA
#define FP_FAST_FMAL FP_FAST_FMA

#define FP_ILOGB0     -INT_MAX
#define FP_ILOGBNAN    INT_MAX

#define MATH_ERRNO  1
#define MATH_ERREXCEPT 2
#define math_errhandling 

#include "math_suppl.h"


/* -----------------------------------------------------------------------*/
#pragma STDC FP_CONTRACT 

/* --------------------------------------------------------------------*/
#define fpclassify(x) _fpclassify(x)


/*---------------------------------------------------------------------------*/
#define isfinite(x) _isfinite(x)

/* ----------------------------------------------------------------------------------------------*/
/* #define isinf(x) (int) ( !(0x7ff != ((0x7FF0000000000000ULL & (*((unsigned int *)(&x))))>>52) ) )
   der geht wenn auf integer zugewiesen, aber nicht in where */
#define isinf(x) _isinf(x)

/* ----------------------------------------------------------------------------------------------*/
#define isnan(x) _isnan(x)

/* ----------------------------------------------------------------*/
#define isnormal(x) _isfinite(x)

/* ----------------------------------------------------------------------------------------------*/
#define signbit(x) _signbit(x)

/* -----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double sqrt( double a );
#else
#if defined(_uses_sqrt_math_h) || !defined(__cflow_processed) 
math_inline double sqrt( double a )
{
    register double ris, tmp_a2;
    register double three_half, one_half;
    register double aR;

    _math_setdouble(0x3FE0000000000000,one_half);   // 0.5
    _math_setdouble(0x3FF8000000000000,three_half); // 1.5

    aR = a;
    ris = 0.0;
/*  where ( aR>0.0 ) {   // avoids exception if x<0  */
    where ( aR ) {       // original 
        asm("\tlutmove %0 $ZERO" : "=r" (ris));
        asm("\tlutisqrt.L %0 %1" : "=r" (ris) : "r" (aR));
    }
    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;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double log(double x);
#else
#if defined(_uses_log_math_h) || !defined(__cflow_processed)
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;

    asm("\n\t!!-- begin log()");
    _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("\tlutmove %0 $ZERO\n\tlutlogm.L %0 %1" : "=r" (x0) : "r" (xx));
    asm("\tlutmove %0 $ZERO\n\tlutloge.L %0 %1" : "=r" (e1) : "r" (xx));

    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;

    asm("\n\t!!-- end log()");
    return P;
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double exp(double x);
#else
#if defined(_uses_exp_math_h) || !defined(__cflow_processed)
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("\tlutmove %0 $ZERO" : "=r" (PP));
    asm("\tlutexp.L %0 %1" : "=r" (PP) : "r" (z));

    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;
}    
#endif
#endif // Has Main

/* ----------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double expm1(double x);
#else
#if defined(_uses_expm1_math_h) || !defined(__cflow_processed)
math_inline double expm1(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("\tlutmove %0 $ZERO" : "=r" (PP));
    asm("\tlutexp.L %0 %1" : "=r" (PP) : "r" (z));

    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;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double acos(double x);
#else
#if defined(_uses_acos_math_h) || !defined(__cflow_processed)
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;
} 
#endif
#endif // Has Main

/* -----------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double asin(double x);
#else
#if defined(_uses_asin_math_h) || !defined(__cflow_processed)
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;
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double atan(double x);
#else
#if defined(_uses_atan_math_h) || !defined(__cflow_processed)
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;    
}
#endif
#endif // Has Main

/* -----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double atan2(double y, double x);
#else
#if defined(_uses_atan2_math_h) || !defined(__cflow_processed)
math_inline double atan2(double y, double x)
{
    register double xx, absx, z, res, pi;
    register double pihalflead, pihalftrail;

    _math_setdouble(  0x3ff9220000000000 , pihalflead  );
    _math_setdouble(  0xbed2aeef4b9ee59e , pihalftrail );
    pi = 3.14159265358979323846;  // replace by pilead,pitrail ? /

    where (x == 0.0) {
       res = 0.0;

       where ( y > 0.0) { 
             res = 0.5 * pi;
       }
       where ( y < 0.0) { 
             res = -0.5 * pi;
       }

    } else {

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

    z = absx;
    where ( absx > 1.0 ) {
      z = 1.0 / absx;
    }

    res = atan_k( z );

    where ( absx > 1.0 ) {
      res = pihalftrail - res;
      res = pihalflead  + res;
    }

    where ( x < 0.0 ) {
       res = res - pi;
    }
    where (xx <= 0.0 ) {  // "<=" takes care of x<0, y=0
       res = -res;
    }

    }

    return res;    
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double cos(double x);
#else
#if defined(_uses_cos_math_h) || !defined(__cflow_processed)
math_inline double cos(double x)
{
    register double y, siny, cosy;
    register double res, 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;

    // Compute 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;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double sin(double x);
#else
#if defined(_uses_sin_math_h) || !defined(__cflow_processed)
math_inline double sin(double x)
{
    register double y, siny, cosy;
    register double res, 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;

    // -----------------------------------------
    // Compute 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;
} 
#endif
#endif // Has Main

/* ---------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double tan(double x);
#else
#if defined(_uses_tan_math_h) || !defined(__cflow_processed)
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;
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double acosh(double x);
#else
#if defined(_uses_acosh_math_h) || !defined(__cflow_processed)
math_inline double acosh(double x)
{
    register double xr, res;

    xr = x;
    res = log(xr+sqrt(xr*xr-1.0));

    return res;
}    
#endif
#endif // Has Main

/* -----------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double asinh(double x);
#else
#if defined(_uses_asinh_math_h) || !defined(__cflow_processed)
math_inline double asinh(double x)
{
    register double xr, res;

    xr = x;
    res = log(xr+sqrt(xr*xr+1.0));

    return res; 
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double atanh(double x);
#else
#if defined(_uses_atanh_math_h) || !defined(__cflow_processed)
math_inline double atanh(double x)
{
  register double xr, res;

  xr = x;
  res = 0.5*log( (1.0+xr)/(1.0-xr) );

  return res;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double cosh(double x);
#else
#if defined(_uses_cosh_math_h) || !defined(__cflow_processed)
math_inline double cosh(double x)
{
    register double xr, xm, expx, expxm, half, res;

/* version 1 */
/*
    half = 0.5;
    xr =   x;
    xm = - xr;

    expx  = expm1( xr );
    expxm = expm1( xm );

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

/* version 2 */
    half = 0.5;
    xr = x;
    expx = exp(xr);
    res = 1.0/expx;
    res = half * (expx+res);

    return res;
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double sinh(double x);
#else
#if defined(_uses_sinh_math_h) || !defined(__cflow_processed)
math_inline double sinh(double x)
{
    register double xr, xm, expx, expxm, half, res;

/* version 1 */
/* 
    half = 0.5;
    xr =   x;
    xm = - xr;

    expx  = expm1(xr);
    expxm = expm1(xm);

    res = half * (expx - expxm);
*/

/* version 1 */
    half = 0.5;
    xr =   x;
    expx  = exp(xr);
    res = 1.0/expx;
    res = half * (expx - res);

    return res;
}    
#endif
#endif // Has Main

/* -----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double tanh(double xxx);
#else
#if defined(_uses_tanh_math_h) || !defined(__cflow_processed)
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  = expm1(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;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double exp2(double x);
#else
#if defined(_uses_exp2_math_h) || !defined(__cflow_processed)
math_inline double exp2(double x)
{
    register double loge2,y,res;

    loge2 = 0.69314718055994530942;

    y = x*loge2;
    res = exp(y);

    return res;
}    
#endif
#endif // Has Main

/* ----------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double frexp(double value, int *eptr);
#else
#if defined(_uses_frexp_math_h) || !defined(__cflow_processed)
math_inline double frexp(double value, int *eptr)
{
 register int ptr;
 register uint64_t ix;
 register double res;
 register union {
   double zd;
   int zint;
 } z;

        z.zd = value;
                                                       // write_int(hx);
        ix = 0x7fffffffffffffffull & z.zint;           // strip sign
                                                       // write_int(ix);
/*      *eptr = 0;

        where (ix >= 0x7ff0000000000000 || (ix==0)) {
           return x;                                   // 0,inf,nan
        }
*/
        ptr = (ix >> 52)-1022;                         // exponent+1
                                                       // write_int(hx);
        z.zint = ( (z.zint & 0x800fffffffffffff) | 0x3fe0000000000000);
                                                       // sign plus mantissa
                                                       // normalize and divide by 2
                                                       // write_int(hx);
        *eptr = ptr;
        res   = z.zd;
        
  return res;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline int ilogb(double x);
#else
#if defined(_uses_ilogb_math_h) || !defined(__cflow_processed)
math_inline int ilogb(double x)
{
     register double xr;
     register uint64_t _mask;
     register int _sh,ie;

     _mask=0x7FF0000000000000ULL;
     _sh  =52;

     xr = x;
     asm("\tiand.L %0 %1 %2" : "=r" (xr) : "r" (_mask) , "r" (xr) );
     asm("\tilsh.L %0 %1 %2" : "=r" (ie) : "r" (xr) , "r" (_sh) );
     asm("\tlutcross.H %0 %1" : "=r" (ie) : "r" (ie) );
     ie = ie - 0x3ff;

     return ie;
}    
#endif
#endif // Has Main

/* ----------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double ldexp(double x, int n);
#else
#if defined(_uses_ldexp_math_h) || !defined(__cflow_processed)
math_inline double ldexp(double x, int n)
{
    register double xr = 0;
    register int nr,_sh;

    _sh = -52;

    nr = n+1023;
    asm("\tilsh.L %0 %1 %2" : "=r" (xr) : "r" (nr) , "r" (_sh) );

    xr = x*xr;

    return xr;
}    
#endif
#endif // Has Main

/* ----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double log10(double x);
#else
#if defined(_uses_log10_math_h) || !defined(__cflow_processed)
math_inline double log10(double x)
{
    register double log10e, res;

    log10e = 0.43429448190325182765;

    res = log(x)*log10e;

    return res;
}    
#endif
#endif // Has Main

/* -----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double log1p(double xxx);
#else
#if defined(_uses_log1p_math_h) || !defined(__cflow_processed)
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("\tlutmove %0 $ZERO" : "=r" (x0));
    asm("\tlutmove %0 $ZERO" : "=r" (e1));
    asm("\tlutlogm.L %0 %1" : "=r" (x0) : "r" (xx));
    asm("\tlutloge.L %0 %1" : "=r" (e1) : "r" (xx));

    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;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double log2(double x);
#else
#if defined(_uses_log2_math_h) || !defined(__cflow_processed)
math_inline double log2(double x)
{
    register double log2e, res, xr;

    _math_setdouble(  0x3ff71547652b82fe, log2e    );  

    xr  = x;
    res = log(xr)*log2e;

    return res;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if (FLT_RADIX==2)
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double logb(double x);
#else
#if defined(_uses_logb_math_h) || !defined(__cflow_processed)
math_inline double logb(double x)
{
     register double xr;
     register unsigned int _mask;
     register int _sh,ie;

     _mask=0x7FF0000000000000ULL;
     _sh  =52;

     xr = x;
     asm("\tiand.L %0 %1 %2" : "=r" (xr) : "r" (_mask) , "r" (xr) );
     asm("\tilsh.L %0 %1 %2" : "=r" (ie) : "r" (xr) , "r" (_sh) );
     xr = ie - 0x3ff;

     return xr;
}    
#endif
#endif // Has Main
#else
#error "FLT_RADIX = 2 requested"
#endif

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double modf(double value, double *dptr);
#else
#if defined(_uses_modf_math_h) || !defined(__cflow_processed)
math_inline double modf(double value, double *dptr)
{
  register double y,ptr;
  register uint64_t _emask,_mmask,_i;
  register int _exp,_sh;
  register union {
    double zd;
    int zint;
  } x,z;

  _emask = 0x7ff0000000000000ULL;
  _mmask = 0x000fffffffffffffULL;
  _sh    = 52;

  x.zd = value;
  _exp = ((x.zint&_emask)>>_sh)-1023;        // exponent of x
                                             // write_int(_exp);
  where (_exp < 52) {                        // chance of fractional part
       where (_exp < 0) {                    // |x| < 1
            ptr = 0.0;                       // correct for sign ? (+-0?)
            y = x.zd;
       } else {
            _i = _mmask>>_exp;
            where ( (x.zint&_i) == 0) {         // x is integer
                 ptr = x.zd;
                 y     = 0.0;                // correct sign ? (+-0?)
            } else {
                 z.zint = x.zint&(~_i);
                 ptr = z.zd;
                 y = x.zd-ptr;
            }
       }
  } else {                                   // no fractional part
       ptr = x.zd;
       y     = 0.0;
  }

  *dptr = ptr;
  return y;
}
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if (FLT_RADIX==2)
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double scalbn(double x, int n);
#else
#if defined(_uses_scalbn_math_h) || !defined(__cflow_processed)
math_inline double scalbn(double x, int n)
{
    register double xr;
    register int nr,_sh;

    _sh = -52;

    nr = n+1023;
    asm("\tilsh.L %0 %1 %2" : "=r" (xr) : "r" (nr) , "r" (_sh) );

    xr = x*xr;

    return xr;
}    
#endif
#endif // Has Main
#else
#error "FLT_RADIX = 2 requested"
#endif

/* ----------------------------------------------------------------*/
#if (FLT_RADIX==2)
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double scalbln(double x, long int n);
#else
#if defined(_uses_scalbln_math_h) || !defined(__cflow_processed)
math_inline double scalbln(double x, long int n)
{
    register double xr;
    register int nr,_sh;

    _sh = -52;

    nr = n+1023;
    asm("\tilsh.L %0 %1 %2" : "=r" (xr) : "r" (nr) , "r" (_sh) );

    xr = x*xr;

    return xr;
}    
#endif
#endif
#else
#error "FLT_RADIX = 2 requested"
#endif

/* ----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double cbrt(double x);
#else
#if defined(_uses_cbrt_math_h) || !defined(__cflow_processed)
math_inline double cbrt(double x)
{
  register double xr, xabs, res, one_third;
  _math_setdouble( 0x3fd5555555555555 , one_third);

  xr = x;
  xabs = xr;
  where (xr < 0.0) {
        xabs = -xr;
  }

  res = 0.0;
  where ( xr != 0.0) {
     res = exp(one_third*log(xabs));
  } 

  where (xr<0.0) {
        res = -res;
  } 
  
  return res;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double fabs(double x);
#else
#if defined(_uses_fabs_math_h) || !defined(__cflow_processed)
math_inline double fabs(double x)
{
  register double xr, res;
  register uint64_t _mask;

/* version 1 */
  _mask= 0x7fffffffffffffffULL;

  xr = x;
  asm("\tiand.L %0 %1 %2" : "=r" (xr) : "r" (_mask) , "r" (xr) );
  return xr;
 
/* version 2 */ 
/*
  xr = x;
  where( xr < 0.0 ) {
    xr = -xr;
  } 
  return xr;
*/

/* version 3 */ 
/*
  xr = x;
  asm("\tfnorm_ap %0 $F_ONE %1 $ZERO" : "=r" (res) : "r" (xr));
  return res;
*/

}    
#endif
#endif // Has Main

/* ----------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double hypot(double x, double y);
#else
#if defined(_uses_hypot_math_h) || !defined(__cflow_processed)
math_inline double hypot(double x, double y)
{
  register double res;

  res = sqrt(x*x+y*y);

  return res;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double erf(double x);
#else
#if defined(_uses_erf_math_h) || !defined(__cflow_processed) 
math_inline double erf(double x)
{
#ifdef __cflow_processed
  #error "Not yet implemented"
#endif
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double erfc(double x);
#else
#if defined(_uses_erfc_math_h) || !defined(__cflow_processed) 
math_inline double erfc(double x)
{
#ifdef __cflow_processed
  #error "Not yet implemented"
#endif
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double lgamma(double x);
#else
#if defined(_uses_lgamma_math_h) || !defined(__cflow_processed) 
math_inline double lgamma(double x)
{
#ifdef __cflow_processed
  #error "Not yet implemented"
#endif
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double tgamma(double x);
#else
#if defined(_uses_tgamma_math_h) || !defined(__cflow_processed) 
math_inline double tgamma(double x)
{
#ifdef __cflow_processed
  #error "Not yet implemented"
#endif
}    
#endif
#endif // Has Main

/* -----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double ceil(double x);
#else
#if defined(_uses_ceil_math_h) || !defined(__cflow_processed)
math_inline double ceil(double x)
{
    register double tmpx,res,nummag,delta;


    _math_setdouble(  0x4337ffffffffffff, nummag   );

    tmpx = x;
    res = ( (tmpx+nummag) - nummag); 
    where (res < x) {
          res = res + 1.0;
    }

/* old version
    _math_setdouble(  0x4330000000000000, nummag   );

    tmpx = x;
    where (x<0.0) {
          tmpx = -x;
    }
    res = ( ((tmpx+0.5)+nummag) - nummag)-1.0; 
    delta = tmpx - res;
    where (delta >= 1.0) {
          res = res + 1.0;
    }
    where (x<0.0) {
          res = -res;
    } 
    where (res < x) {
          res = res + 1.0;
    }
*/

    return res;
}    
#endif
#endif // Has Main

/* --------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double floor(double x);
#else
#if defined(_uses_floor_math_h) || !defined(__cflow_processed)
math_inline double floor(double x)
{
    register double tmpx,res,nummag;

//    _math_setdouble(  0x4337ffffffffffff, nummag   );
    _math_setdouble(  0x4330000000000000, nummag   );

    tmpx = x;
    res = (tmpx+nummag) - nummag; 
    where (res > x) {
          res = res - 1.0;
    }

/* old version
    _math_setdouble(  0x4330000000000000, nummag   );

    tmpx = x;
    where (x<0.0) {
          tmpx = -x;
    }
    res = ( ((tmpx+0.5)+nummag) - nummag)-1.0; 
    delta = tmpx - res;
    where (delta >= 1.0) {
          res = res + 1.0;
    }
    where (x<0.0) {
          res = -res;
    } 
    where (res > x) {
          res = res - 1.0;
    }
*/

    return res;
}    
#endif
#endif // Has Main

/* --------------------------------------------------------------------*/
#if (FLT_ROUNDS==1)
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double nearbyint(double x);
#else
#if defined(_uses_nearbyint_math_h) || !defined(__cflow_processed)
math_inline double nearbyint(double x)
{
     register double res;

     where (x>0.0) {
           res = floor(x);
           where ( (x-res) >= 0.5 ) {
           res = res + 1.0;
           }
     } else {
           res = ceil(x);
           where ( (res-x) >= 0.5 ) {
           res = res - 1.0;
           }
     }

     return res;
}
#endif
#endif // Has Main
#else
#error "FLT_ROUNDS = 1 expected"
#endif

/* -------------------------------------------------------------------*/
#if (FLT_ROUNDS==1)
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double rint(double x);
#else
#if defined(_uses_rint_math_h) || !defined(__cflow_processed)
math_inline double rint(double x)
{
     register double res;

     where (x>0.0) {
           res = floor(x);
           where ( (x-res) >= 0.5 ) {
           res = res + 1.0;
           }
     } else {
           res = ceil(x);
           where ( (res-x) >= 0.5 ) {
           res = res - 1.0;
           }
     }

     return res;
}    
#endif
#endif // Has Main
#else
#error "FLT_ROUNDS = 1 expected"
#endif

/* -------------------------------------------------------------------*/
#if (FLT_ROUNDS==1)
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline long int lrint(double x);
#else
#if defined(_uses_lrint_math_h) || !defined(__cflow_processed)
math_inline long int lrint(double x)
{
     register double res;
     register long int res_i;

     where (x>0.0) {
           res = floor(x);
           where ( (x-res) >= 0.5 ) {
           res = res + 1.0;
           }
     } else {
           res = ceil(x);
           where ( (res-x) >= 0.5 ) {
           res = res - 1.0;
           }
     }

     res_i = res;
     return res_i;
}    
#endif // Has Main
#endif
#else
#error "FLT_ROUNDS = 1 expected"
#endif

/* ----------------------------------------------------------------------*/
#define llrint(x) lrint(x)

/* ---------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double round(double x);
#else
#if defined(_uses_round_math_h) || !defined(__cflow_processed)
math_inline double round(double x)
{
     register double res,tmpx,nummag;

/*  not correct e.g. 16.5
    _math_setdouble(  0x4337fffffffffffe, nummag   );

    tmpx = x;           
    res = (tmpx+nummag) - nummag; 
*/

/*  old version */
     where (x>0.0) {
           res = floor(x);
           where ( (x-res) >= 0.5 ) {
           res = res + 1.0;
           }
     } else {
           res = ceil(x);
           where ( (res-x) >= 0.5 ) {
           res = res - 1.0;
           }
     }

     return res;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline long int lround(double x);
#else
#if defined(_uses_lround_math_h) || !defined(__cflow_processed)
math_inline long int lround(double x)
{
     register double res;
     register long int res_i;

     where (x>0.0) {
           res = floor(x);
           where ( (x-res) >= 0.5 ) {
           res = res + 1.0;
           }
     } else {
           res = ceil(x);
           where ( (res-x) >= 0.5 ) {
           res = res - 1.0;
           }
     }

     res_i = res;
     return res_i;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#define llround(x) lround(x)

/* --------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double trunc(double x);
#else
#if defined(_uses_trunc_math_h) || !defined(__cflow_processed)
math_inline double trunc(double x)
{
    register double tmpx,res,nummag,delta;

    //  -- rtc version --
    //  
    //      _math_setdouble(  0x4330000000000000, nummag   );
    //  
    //      tmpx = x;
    //      where (x<0.0) {
    //            tmpx = -x;
    //      }
    //      res = (((tmpx+0.5)+nummag) - nummag)-1.0; 
    //      delta = tmpx - res;
    //      where (delta >= 1.0) {
    //            res = res + 1.0;
    //      }
    //      where (x<0.0) {
    //            res = -res;
    //      } 


    // --  C version --
     register double xr,tmp;
     register uint64_t _emask;
     register int _sh,ie;

     _emask=0x7FF0000000000000ULL;
     _sh  =52;

     xr = x;
     // tmp=unshifted exponent
     asm("\tiand.L %0 %1 %2" : "=r" (tmp) : "r" (_emask) , "r" (xr) );
     // ie=shifted exponent
     asm("\tilsh.L %0 %1 %2" : "=r" (ie) : "r" (tmp) , "r" (_sh) );
     // integer exponent
     ie = ie - 0x3ff;

     // return 0 for -1 < x < 1
     if( ie < 0 ) {
        double rv;
        asm("\tiand %R0 %1 %2\n\tior %0 $ZERO %R0" : "=r" (rv) : "r" (x), "r" (0x8000000000000000ULL));
        return rv;
     }
     _sh = _sh - ie;
     asm("\tilsh.L %0 %1 %2" : "=r" (xr) : "r" (xr) , "r" (_sh) );
     _sh = -_sh;
     asm("\tilsh.L %0 %1 %2" : "=r" (res) : "r" (xr) , "r" (_sh) );

    return res;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double copysign(double x, double y);
#else
#if defined(_uses_copysign_math_h) || !defined(__cflow_processed)
math_inline double copysign(double x, double y)
{
  register uint64_t _smask;
  register double res;
  register union {
    double zd;
    int zint;
  } _x,_y,z;

  _smask = 0x8000000000000000ULL;

  _x.zd = x;
  _y.zd = y;
                                  // write_int(ix);
                                  // write_int(iy);
  _x.zint &= ~_smask;             // everything except sign of x
  _y.zint &= _smask;              // get sign of y
                                  // write_int(ix);
                                  // write_int(iy);
  _x.zint |= _y.zint;             // sign of y into sign of x
                                  // write_int(ix);

  z.zint = _x.zint;
  res = z.zd;

  return res;
}    
#endif
#endif // Has Main

/* ----------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double fmod(double x, double y);
#else
#if defined(_uses_fmod_math_h) || !defined(__cflow_processed)
math_inline double fmod(double x, double y)
{
      register double ratio,res;

      ratio = x/y;
      res = x - trunc(ratio) * y;

      return res;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double remainder(double x, double y);
#else
#if defined(_uses_remainder_math_h) || !defined(__cflow_processed)
math_inline double remainder(double x, double y)
{
   register double res, ratio, n, two, two_c;

/* old version: not always correct
   res = fabs(y)-fabs(fmod(x,y));
   res = copysign(res,x);
   res = -res;
*/
   two = 1.0 + 1.0;
   ratio = x/y;
   n = round(ratio);
   two_c = fmod(n,two);
//   write_double(n);
   where ( fabs(n-ratio) == 0.5 ) {
      n = n - two_c;
   }
//   write_double(n);
   res = x - n * y;

   return res;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double remquo(double x, double y, int *quo);
#else
#if defined(_uses_remquo_math_h) || !defined(__cflow_processed)
math_inline double remquo(double x, double y, int *quo)
{
  register double res,ratio,eight,dquo;
  register int iquo;
  eight = 8.0;

  res = remainder(x,y);

  ratio = trunc(x/y);
  dquo = fmod(ratio,eight);
  dquo = copysign(dquo,ratio);
  iquo = dquo;
  *quo = iquo;

  return res;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double nan(const char *tagp);
#else
#if defined(_uses_nan_math_h) || !defined(__cflow_processed)
math_inline double nan(const char *tagp)
{
#ifdef __cflow_processed
# warning "nan(): NaN and Infinity are not supported by ApeNEXT. Generates systematically an exception."
#endif
  return 0;
}    
#endif
#endif // Has Main

/* --------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double nextafter(double x, double y);
#else
#if defined(_uses_nextafter_math_h) || !defined(__cflow_processed)
math_inline double nextafter(double x, double y)
{
  register double xr,yr,res;
  register unsigned int ix,_lastbitmask;
  register union {
    double zd;
    int zint;
  } z;

  _lastbitmask = 0x0000000000000001ULL;

  xr = x;
  yr = y;
  z.zd = xr;
  ix   = z.zint;
                                    //  write_int(ix);

  where (xr != yr) {           
     where (x==0.0) {
        z.zint = _lastbitmask;
        res = copysign(z.zd,yr);
     } else {
        where (  ((xr>0.0)&&(yr>xr)) \
               ||((xr<0.0)&&(yr<xr)) ) {
             z.zint = ix + _lastbitmask;
        } 
        else {
             z.zint = ix - _lastbitmask;
        }
        res = z.zd;
     }
  } else {
        res = yr;
  }
// write_int(z.zint);     // one more digit in write_double useful

  return res;
}    
#endif
#endif // Has Main

/* -------------------------------------------------------------------*/
/* math_inline double nexttoward(double x, long double y) */
#define nexttoward(x,y) nextafter(x,y)

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double pow(double x, double y);
#else
#if defined(_uses_pow_math_h) || !defined(__cflow_processed)
math_inline double pow(double x, double y)
{
  register double tmp_y,tmp_x,res,sign,two;
  register int    tmp_i;
  register double tmp;
  register uint64_t _emask;
  register int _sh,ie;

  asm("\n\t!!-- begin pow()");
  res = 0.0;
//  two = 2.0;

/* old version */
/*
  where (x<0.0) {
        tmp_i = y;
        tmp_y = tmp_i;
        where (tmp_y==y) {
              tmp_x = -x;
              res = exp(y*log(tmp_x));
              sign = 1.0-2.0*fmod(y,two); 
              res = sign*res;
        } else {
        tmp_x = log(x);    // crash
        }
  } 
  where (x>0.0) {
        res =  exp( y * log(x) );
  }
*/

/* new version */
  sign = 1.0;
  tmp_x = x;
  tmp_y = y;
  where (tmp_x < 0.0) {
     tmp_x = -tmp_x;
/* version 1
     tmp_i = y;       //better to use trunc or round
     tmp_y = tmp_i;   // danger for large y
*/
/*   tmp_y = trunc(y); trunc introduced */
/*
       sign = sign - two*fmod(fabs(y),two);
       where (tmp_y != y) {
        sign = 0x7ff0100000000000ULL;
     } 
*/

/* version 2 */
     _emask=0x7FF0000000000000ULL;
     _sh  =52;

     // compute tmp=trunc(y/2)*2
     asm("\tiand.L %0 %1 %2" : "=r" (tmp) : "r" (_emask) , "r" (tmp_y) );
     asm("\tilsh.L %0 %1 %2" : "=r" (ie) : "r" (tmp) , "r" (_sh) );
     ie = ie - 0x400;
     _sh = _sh - ie;
     where( ie < 0 ){ 
        tmp = 0;
     } else {
        asm("\tilsh.L %0 %1 %2" : "=r" (tmp) : "r" (tmp_y) , "r" (_sh) );
        _sh = -_sh;
        asm("\tilsh.L %0 %1 %2" : "=r" (tmp) : "r" (tmp) , "r" (_sh) );
     }

     // check for evenness of tmp
     tmp = tmp_y - tmp;
     where (tmp != 0.0) {
        // y not even
        res = ((union { unsigned u; double d; }){.u=0x7ff0100000000000ULL}).d; // (NAN)  
        where ( (tmp == 1.0) || (tmp == -1.0) )
        {  // y odd
           sign = -1.0;
        }
        else
           tmp_x = 0.0;
     }
  }
/*   asm("rtm 0x1000 [ %0 %1 %2 %3 ]" : : "r" (sign), "r" (tmp), "r" (tmp_x), "r" (res) ); -- LM ????????????????? */
  where (tmp_x != 0.0) {
     res = sign*exp(y*log(tmp_x));
  }
  else
      where (y==0.0)
          res=1.0;
  
  asm("\t!!-- end pow()\n");
  return res;
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double fdim(double x, double y);
#else
#if defined(_uses_fdim_math_h) || !defined(__cflow_processed)
math_inline double fdim(double x, double y)
{
  register double res;

  where( x > y ) {
    res = x - y;
  } else {
    res = 0;
  }
  return res;
}    
#endif
#endif // Has Main

/* ------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double fmax(double x, double y);
#else
#if defined(_uses_fmax_math_h) || !defined(__cflow_processed)
math_inline double fmax(double x, double y)
{
  register double res;

  where( x > y ) {
    res = x;
  } else {
    res = y;
  }
  return res;
}    
#endif
#endif // Has Main

/* --------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double fmin(double x, double y);
#else
#if defined(_uses_fmin_math_h) || !defined(__cflow_processed)
math_inline double fmin(double x, double y)
{
  register double res;

  where( x < y ) {
    res = x;
  } else {
    res = y;
  }
  return res;
}    
#endif
#endif // Has Main

/* ---------------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double fma(double x, double y, double z);
#else
#if defined(_uses_fma_math_h) || !defined(__cflow_processed)
math_inline double fma(double x, double y, double z)
{
  register double res;

  res = x*y + z;

  return res;
}    
#endif
#endif // Has Main

/* --------------------------------------------------------------------*/
#define isgreater(x,y) _isgreater(x,y)

/* -----------------------------------------------------------------*/
#define isgreaterequal(x,y) _isgreaterequal(x,y)

/* -------------------------------------------------------------------*/
#define isless(x,y) _isless(x,y)

/* -------------------------------------------------------------------------*/
#define islessequal(x,y) _islessequal(x,y)

/* --------------------------------------------------------------------*/
#define islessgreater(x,y) _islessgreater(x,y)

/* ----------------------------------------------------------------*/
#define isunordered( x, y) _isunordered( x, y)

#undef _math_setdouble
#endif  /* _MATH_H_ */

---