math_suppl.h

Go to the documentation of this file.

/*==================================================================
!! Routines needed for the apeNEXT C library math.h
!!          which are not standard to the C library
!!
!! $Id: math_suppl.h,v 2.10 2005/11/11 10:22:09 morinl 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:
!!   ->  isqrt(x)    <- 1/sqrt(x)
!!   ->  isqrt2(x) 
!!   ->  sqrt2(x) 
!!   ->  sin_1oct(x) <- sin(x), -pi/4<x<pi/4  
!!   ->  cos_1oct(x) <- cos(x), -pi/4<x<pi/4  
!!   ->  atan_k(x)   <- atan(x), -1<x<1 
!!   ->  asin_k(x)   <- asin(x), -1/2<x<1/2 
!! in addition:
!!   functions to which the following macros expand to (EL 05/09/2003)
!!   -> _signbit(x)
!!   -> _isfinite(x)
!!   -> _isinf(x)
!!   -> _isnan(x)
!! added (EL 17/10/2003)
!!   -> _isgreater(x,y)
!!   -> _isgreaterequal(x,y)
!!   -> _isless(x,y)
!!   -> _islessequal(x,y)
!!   -> _islessgreater(x,y)
!!==================================================================*/

#ifndef _MATH_H_              
#include <math.h>
#else
#ifndef __MATH_SUPPL_H__
#define __MATH_SUPPL_H__

#include <nlibc.h>

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

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

/*-----------------------------------------------------------------
!!
!! 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
!!
!!---------------------------------------------------------------*/

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double isqrt2( double x );
#else
#if defined(_uses_isqrt2_math_suppl_h) || !defined(__cflow_processed)
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( 0x3fe0000000000176, p0 );
    _math_setdouble( 0x3fd8000000000640, p1 );
    _math_setdouble( 0x3fd3fffffd11f1d0, p2 );
    _math_setdouble( 0x3fd17ffffbfd0ec7, p3 );
    _math_setdouble( 0x3fcf81774c2db32f, p4 );
    _math_setdouble( 0x3fcce1962db980e2, 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;
}
#endif
#endif // Has Main

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double sqrt2( double x ) ;
#else
#if defined(_uses_sqrt2_math_suppl_h) || !defined(__cflow_processed)
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( 0x3fe0000000000176, p0 );
    _math_setdouble( 0x3fd8000000000640, p1 );
    _math_setdouble( 0x3fd3fffffd11f1d0, p2 );
    _math_setdouble( 0x3fd17ffffbfd0ec7, p3 );
    _math_setdouble( 0x3fcf81774c2db32f, p4 );
    _math_setdouble( 0x3fcce1962db980e2, 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;
}
#endif
#endif // Has Main

/*-----------------------------------------------------------------
!!
!! Kernel function for the computation of Sin[x]
!!
!! Sin[x] = x ( 1 + x^2 Q[x^2] )    ;   -pi/4 < x < pi/4
!!
!!---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double sin_1oct( double xx );
#else
#if defined(_uses_sin_1oct_math_suppl_h) || !defined(__cflow_processed)
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;
    
}
#endif
#endif // Has Main

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

}
#endif
#endif // Has Main

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

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

/* -----------------------------------------------------
!!  signbit(x)  -> see Notes in math.h
!!  
!!  alternativ: ie = ileadone.L(xu.xint) & 63
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline int _signbit(double x);
#else
#if defined(_uses__signbit_math_suppl_h) || !defined(__cflow_processed)
math_inline int _signbit(double x)
{

    register unsigned int _mask;
    register int _sh,ie;
    register union {
      double xr;
      unsigned int xint;
    } xu;

    _mask=0x8000000000000000ULL;
    _sh  =63;

    xu.xr = x;
    ie = (_mask&xu.xint)>>_sh;

    return ie;
}
#endif
#endif // Has Main
/* -----------------------------------------------------
!!  isfinite(x)  -> see Notes in math.h
!!  version 1  clk cnt 110  run cnt 90
!!  version 2  clk cnt 115  run cnt 93
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline int _isfinite(double x);
#else
#if defined(_uses__isfinite_math_suppl_h) || !defined(__cflow_processed)
math_inline int _isfinite(double x) {

#ifdef __cflow_processed
# warning "_isfinite(): NaN and Infinity are not supported by ApeNEXT. Generates systematically an exception."
#endif

    register unsigned int _emask;
    register int _sh,ie,ij,_im;
    register union {
      double xr;
      unsigned int xint;
    } xu;

    _emask=0x7FF0000000000000ULL;
//    _sh  =52;
//    _im  =0x7ff;

    xu.xr = x;
/* version 1 */
    ie = _emask != (_emask&xu.xint); 

/*    ie = 0x7ff != ((_emask&xu.xint)>>_sh); */
/*  bit shift not necessary
    ii = (_emask&xu.xint)>>_sh;
    ie = 0;
    where (ii != _im) {
          ie = 1;
    }
*/
/* version 2 */
/*
    ii = (_emask&xu.xint);
    ie = 0;
    where (ii != _emask) {
          ie = 1;
    }
*/

    return ie;
}
#endif
#endif // Has Main
/* -----------------------------------------------------
!!  isinf(x)  -> see Notes in math.h
!! ---------------------------------------------------------------*/

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline int _isinf(double x);
#else
#if defined(_uses__isinf_math_suppl_h) || !defined(__cflow_processed)
math_inline int _isinf(double x) {

#ifdef __cflow_processed
# warning "_isinf(): NaN and Infinity are not supported by ApeNEXT. Generates systematically an exception."
#endif

    register unsigned int _emask,_mmask;
    register int ie,ii,ij;
    register union {
      double xr;
      unsigned int xint;
    } xu;

    _emask=0x7FF0000000000000ULL;
    _mmask=0x000fffffffffffffULL;

    xu.xr = x;
    ii = (_emask&xu.xint);
    ij = (_mmask&xu.xint);
    ie = 0;
    where ( (ii == _emask) && (ij == 0) ) {
          ie = 1;
    }

    return ie;
}
#endif
#endif // Has Main
/* -----------------------------------------------------
!!  isnan(x)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline int _isnan(double x) ;
#else
#if defined(_uses__isnan_math_suppl_h) || !defined(__cflow_processed)
math_inline int _isnan(double x) {

#ifdef __cflow_processed
# warning "_isnan(): NaN and Infinity are not supported by ApeNEXT. Generates systematically an exception."
#endif

    register unsigned int _emask,_mmask;
    register int ie,ii,ij;
    register union {
      double xr;
      unsigned int xint;
    } xu;

    _emask=0x7FF0000000000000ULL;
    _mmask=0x000fffffffffffffULL;

    xu.xr = x;
    ii = (_emask&xu.xint);
    ij = (_mmask&xu.xint);
    ie = 0;
    where ( (ii == _emask) && (ij != 0) ) {
          ie = 1;
    }

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

/* -----------------------------------------------------
!!  fpclassify(x)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline int _fpclassify(double x);
#else
#if defined(_uses__fpclassify_math_suppl_h) || !defined(__cflow_processed)
math_inline int _fpclassify(double x)
{
#ifdef __cflow_processed
# warning "_fpclassify(): NaN and Infinity are not supported by ApeNEXT. Generates systematically an exception."
#endif
    register unsigned int _emask,_mmask;
    register int ie,ii,ij;
    register union {
      double xr;
      unsigned int xint;
    } xu;


    _emask=0x7FF0000000000000ULL;
    _mmask=0x000fffffffffffffULL;

    xu.xr = x;
    ii = (_emask&xu.xint);
    ij = (_mmask&xu.xint);

/* isfinite = isnormal: => 1 !! */
    ie = _emask != ii;

/* isinf: => FP_INFINITE */
    where ( (ii == _emask) && (ij == 0) ) {
          ie = FP_INFINITE;
    }

/* isnan: => FP_NAN */
    where ( (ii == _emask) && (ij != 0) ) {
          ie = FP_NAN;
    }

/* iszero: => FP_ZERO */
    where ( ii == 0 ) {
          ie = FP_ZERO;
    }

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

/* -----------------------------------------------------
!!  isgreater(x,y)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline int _isgreater(double x, double y);
#else
#if defined(_uses__isgreater_math_suppl_h) || !defined(__cflow_processed)
math_inline int _isgreater(double x, double y) {

   register unsigned int _smask,_nsmask;
   register int sx,sy,ax,ay;
   register int ie;
   register double xr,yr;

   _smask  = 0x8000000000000000ull;
   _nsmask = 0x7fffffffffffffffull;

   xr = x;
   yr = y;

// sign of x
   asm("iand.L %0 %1 %2" : "=r" (sx) : "r" (_smask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (sx) : "r" (sx) );
// sign of y
   asm("iand.L %0 %1 %2" : "=r" (sy) : "r" (_smask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (sy) : "r" (sy) );

// abs value of x
   asm("iand.L %0 %1 %2" : "=r" (ax) : "r" (_nsmask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (ax) : "r" (ax) );
// abs value of y
   asm("iand.L %0 %1 %2" : "=r" (ay) : "r" (_nsmask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (ay) : "r" (ay) );
/* wenn ich den kram eh in asm programmieren muss, dann
   vielleicht auch bis zum bitteren ende, i.e. mit
   stack-manipulationsroutinen fuer die where cond.,
   der assembler der where cond. sieht naemlich nicht
   berauschend aus                                    */

   ie = 0;
   where (   ( (ax>ay) && (sx==0) ) \
           ||( (ax<ay) && (sy!=0) ) ) {
      ie = 1;
   }

   return ie;
}
#endif
#endif // Has Main
/* -----------------------------------------------------
!!  isgreaterequal(x,y)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double _isgreaterequal(double x, double y);
#else
#if defined(_uses__isgreaterequal_math_suppl_h) || !defined(__cflow_processed)
math_inline double _isgreaterequal(double x, double y) {

   register unsigned int _smask,_nsmask;
   register int sx,sy,ax,ay;
   register int ie;
   register double xr,yr;

   _smask  = 0x8000000000000000ull;
   _nsmask = 0x7fffffffffffffffull;

   xr = x;
   yr = y;

// sign of x
   asm("iand.L %0 %1 %2" : "=r" (sx) : "r" (_smask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (sx) : "r" (sx) );
// sign of y
   asm("iand.L %0 %1 %2" : "=r" (sy) : "r" (_smask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (sy) : "r" (sy) );

// abs value of x
   asm("iand.L %0 %1 %2" : "=r" (ax) : "r" (_nsmask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (ax) : "r" (ax) );
// abs value of y
   asm("iand.L %0 %1 %2" : "=r" (ay) : "r" (_nsmask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (ay) : "r" (ay) );
/* wenn ich den kram eh in asm programmieren muss, dann
   vielleicht auch bis zum bitteren ende, i.e. mit
   stack-manipulationsroutinen fuer die where cond.,
   der assembler der where cond. sieht naemlich nicht
   berauschend aus                                    */

   ie = 0;
   where (   ( (ax>=ay) && (sx==0) ) \
           ||( (ax<=ay) && (sy!=0) ) ) {
      ie = 1;
   }

   return ie;
}
#endif
#endif // Has Main
/* -----------------------------------------------------
!!  isless(x,y)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double _isless(double x, double y);
#else
#if defined(_uses__isless_math_suppl_h) || !defined(__cflow_processed)
math_inline double _isless(double x, double y) {

   register unsigned int _smask,_nsmask;
   register int sx,sy,ax,ay;
   register int ie;
   register double xr,yr;

   _smask  = 0x8000000000000000ull;
   _nsmask = 0x7fffffffffffffffull;

   xr = x;
   yr = y;

// sign of x
   asm("iand.L %0 %1 %2" : "=r" (sx) : "r" (_smask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (sx) : "r" (sx) );
// sign of y
   asm("iand.L %0 %1 %2" : "=r" (sy) : "r" (_smask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (sy) : "r" (sy) );

// abs value of x
   asm("iand.L %0 %1 %2" : "=r" (ax) : "r" (_nsmask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (ax) : "r" (ax) );
// abs value of y
   asm("iand.L %0 %1 %2" : "=r" (ay) : "r" (_nsmask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (ay) : "r" (ay) );
/* wenn ich den kram eh in asm programmieren muss, dann
   vielleicht auch bis zum bitteren ende, i.e. mit
   stack-manipulationsroutinen fuer die where cond.,
   der assembler der where cond. sieht naemlich nicht
   berauschend aus                                    */

   ie = 0;
   where (   ( (ax<ay) && (sy==0) ) \
           ||( (ax>ay) && (sx!=0) ) ) {
      ie = 1;
   }

   return ie;
}
#endif
#endif // Has Main
/* -----------------------------------------------------
!!  islessequal(x,y)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double _islessequal(double x, double y);
#else
#if defined(_uses__islessequal_math_suppl_h) || !defined(__cflow_processed)
math_inline double _islessequal(double x, double y) {

   register unsigned int _smask,_nsmask;
   register int sx,sy,ax,ay;
   register int ie;
   register double xr,yr;

   _smask  = 0x8000000000000000ull;
   _nsmask = 0x7fffffffffffffffull;

   xr = x;
   yr = y;

// sign of x
   asm("iand.L %0 %1 %2" : "=r" (sx) : "r" (_smask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (sx) : "r" (sx) );
// sign of y
   asm("iand.L %0 %1 %2" : "=r" (sy) : "r" (_smask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (sy) : "r" (sy) );

// abs value of x
   asm("iand.L %0 %1 %2" : "=r" (ax) : "r" (_nsmask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (ax) : "r" (ax) );
// abs value of y
   asm("iand.L %0 %1 %2" : "=r" (ay) : "r" (_nsmask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (ay) : "r" (ay) );
/* wenn ich den kram eh in asm programmieren muss, dann
   vielleicht auch bis zum bitteren ende, i.e. mit
   stack-manipulationsroutinen fuer die where cond.,
   der assembler der where cond. sieht naemlich nicht
   berauschend aus                                    */

   ie = 0;
   where (   ( (ax<=ay) && (sy==0) ) \
           ||( (ax>=ay) && (sx!=0) ) ) {
      ie = 1;
   }

   return ie;
}
#endif
#endif // Has Main
/* -----------------------------------------------------
!!  islessgreater(x,y)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double _islessgreater(double x, double y);
#else
#if defined(_uses__islessgreater_math_suppl_h) || !defined(__cflow_processed)
math_inline double _islessgreater(double x, double y) {

   register unsigned int _smask,_nsmask;
   register int sx,sy,ax,ay;
   register int ie;
   register double xr,yr;

   _smask  = 0x8000000000000000ull;
   _nsmask = 0x7fffffffffffffffull;

   xr = x;
   yr = y;

// sign of x
   asm("iand.L %0 %1 %2" : "=r" (sx) : "r" (_smask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (sx) : "r" (sx) );
// sign of y
   asm("iand.L %0 %1 %2" : "=r" (sy) : "r" (_smask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (sy) : "r" (sy) );

// abs value of x
   asm("iand.L %0 %1 %2" : "=r" (ax) : "r" (_nsmask), "r" (xr) );
   asm("lutcross.H %0 %1" : "=r" (ax) : "r" (ax) );
// abs value of y
   asm("iand.L %0 %1 %2" : "=r" (ay) : "r" (_nsmask), "r" (yr) );
   asm("lutcross.H %0 %1" : "=r" (ay) : "r" (ay) );
/* wenn ich den kram eh in asm programmieren muss, dann
   vielleicht auch bis zum bitteren ende, i.e. mit
   stack-manipulationsroutinen fuer die where cond.,
   der assembler der where cond. sieht naemlich nicht
   berauschend aus                                    */

   ie = 1;
   where ( (ax==ay) && (sx==sy) ) {
      ie = 0;
   }

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

/* -----------------------------------------------------
!!  isunordered(x,y)  -> see Notes in math.h
!! ---------------------------------------------------------------*/
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double _isunordered(double x, double y);
#else
#if defined(_uses__isunordered_math_suppl_h) || !defined(__cflow_processed)
math_inline double _isunordered(double x, double y) {

    register unsigned int _emask,_mmask;
    register int ie,ii,ij,ie1;
    register union {
      double xr;
      unsigned int xint;
    } xu;

    _emask=0x7FF0000000000000ULL;
    _mmask=0x000fffffffffffffULL;

    xu.xr = x;
    ii = (_emask&xu.xint);
    ij = (_mmask&xu.xint);
    ie1 = 0;
    where ( (ii == _emask) && (ij != 0) ) {
          ie1 = 1;
    }
    xu.xr = y;
    ii = (_emask&xu.xint);
    ij = (_mmask&xu.xint);
    ie = 0;
    where ( (ii == _emask) && (ij != 0) ) {
          ie = ie1;
    }

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

#endif //  __MATH_SUPPL_H__
#endif //  ifdef _MATH_H_

---