#include <stdlib.h>
#include <string.h>

#include "rpplugin.h"
#include "rpdbgerr.h"

#include <rtieeef.h>


#if (!defined(DOXYGEN))
static const char rcsid[] __RWUNUSED__ = 
    "@@@@(#)$Id: rtieeef.c,v 1.16 2001/07/18 14:45:57 johns Exp $";
#endif /* (!defined(DOXYGEN)) */


/**
 * \ingroup rtieeef
 * \page rtieeefoverview RtIEEEf Toolkit Overview
 *
 * Functions for IEEE 754 standard single precsision floating point arithmetic.
 * See also
 * - OpenBSD code
 *   -# ftp://ftp.openbsd.org/pub/OpenBSD/src/lib/libm/
 * - IEEE Standard 754 Floating Point Numbers
 *   -# http://www.research.microsoft.com/~hollasch/cgindex/coding/ieeefloat.html
 * - IEEE standard 754 for binary floating-point arithmetic
 *   -# http://cch.loria.fr/documentation/IEEE754/
 * - 754-1985 IEEE Standard for Binary Floating-Point Arithmetic 1985 
 *   -# http://shop.ieee.org/store/product.asp?prodno=SS10116
 */

/*
 * from: @(#)fdlibm.h 5.1 93/09/24
 */

/* The original fdlibm code used statements like:
        n0 = ((*(int*)&one)>>29)^1;             * index of high word *
        ix0 = *(n0+(int*)&x);                   * high word of x *
        ix1 = *((1-n0)+(int*)&x);               * low word of x *
   to dig two 32 bit words out of the 64 bit IEEE floating point
   value.  That is non-ANSI, and, moreover, the gcc instruction
   scheduler gets it wrong.  We instead use the following macros.
   Unlike the original code, we determine the endianness at compile
   time, not at run time; I don't see much benefit to selecting
   endianness at run time.  */

/* A union which permits us to convert between a double and two 32 bit
   ints.  */

/*
 * The arm32 port is little endian except for the FP word order which is
 * big endian.
 */

typedef struct rtieee_parts_bigendian rtieee_parts_bigendian;
struct rtieee_parts_bigendian
{
    unsigned int        msw;
    unsigned int        lsw;
};

typedef struct rtieee_parts_littleendian rtieee_parts_littleendian;
struct rtieee_parts_littleendian
{
    unsigned int        lsw;
    unsigned int        msw;
};

#if (defined(rwLITTLEENDIAN))
typedef struct rtieee_parts_littleendian rtieee_parts;
#else /* (defined(rwLITTLEENDIAN)) */
typedef struct rtieee_parts_bigendian rtieee_parts;
#endif /* (defined(rwLITTLEENDIAN)) */

typedef union rtieee_double_shape_type rtieee_double_shape_type;
union rtieee_double_shape_type
{
    double              value;
    rtieee_parts        parts;
};

/* Get two 32 bit ints from a double.  */

#define EXTRACT_WORDS(ix0,ix1,d)                \
do {                                            \
  rtieee_double_shape_type ew_u;                \
  ew_u.value = (d);                             \
  (ix0) = ew_u.parts.msw;                       \
  (ix1) = ew_u.parts.lsw;                       \
} while (0)

/* Get the more significant 32 bit int from a double.  */

#define GET_HIGH_WORD(i,d)                      \
do {                                            \
  rtieee_double_shape_type gh_u;                \
  gh_u.value = (d);                             \
  (i) = gh_u.parts.msw;                         \
} while (0)

/* Get the less significant 32 bit int from a double.  */

#define GET_LOW_WORD(i,d)                       \
do {                                            \
  rtieee_double_shape_type gl_u;                \
  gl_u.value = (d);                             \
  (i) = gl_u.parts.lsw;                         \
} while (0)

/* Set a double from two 32 bit ints.  */

#define INSERT_WORDS(d,ix0,ix1)                 \
do {                                            \
  rtieee_double_shape_type iw_u;                \
  iw_u.parts.msw = (ix0);                       \
  iw_u.parts.lsw = (ix1);                       \
  (d) = iw_u.value;                             \
} while (0)

/* Set the more significant 32 bits of a double from an int.  */

#define SET_HIGH_WORD(d,v)                      \
do {                                            \
  rtieee_double_shape_type sh_u;                \
  sh_u.value = (d);                             \
  sh_u.parts.msw = (v);                         \
  (d) = sh_u.value;                             \
} while (0)

/* Set the less significant 32 bits of a double from an int.  */

#define SET_LOW_WORD(d,v)                       \
do {                                            \
  rtieee_double_shape_type sl_u;                \
  sl_u.value = (d);                             \
  sl_u.parts.lsw = (v);                         \
  (d) = sl_u.value;                             \
} while (0)

/* A union which permits us to convert between a float and a 32 bit
   int.  */

typedef union rtieee_float_shape_type rtieee_float_shape_type;
union rtieee_float_shape_type
{
    float               value;
    unsigned int        word;
};

/* Get a 32 bit int from a float.  */

#define GET_FLOAT_WORD(i,d)                     \
do {                                            \
  rtieee_float_shape_type gf_u;                 \
  gf_u.value = (d);                             \
  (i) = gf_u.word;                              \
} while (0)

/* Set a float from a 32 bit int.  */

#define SET_FLOAT_WORD(d,i)                     \
do {                                            \
  rtieee_float_shape_type sf_u;                 \
  sf_u.word = (i);                              \
  (d) = sf_u.value;                             \
} while (0)

/************************** constants *********************************/

#define   C ( (float)5.4285717010e-01 )
#define   C1 ( (float)4.1666667908e-02 )
#define   C2 ( (float)-1.3888889225e-03 )
#define   C3 ( (float)2.4801587642e-05 )
#define   C4 ( (float)-2.7557314297e-07 )
#define   C5 ( (float)2.0875723372e-09 )
#define   C6 ( (float)-1.1359647598e-11 )
#define   D ( (float)-7.0530611277e-01 )
#define   E ( (float)1.4142856598e+00 )
#define   F ( (float)1.6071428061e+00 )
#define   G ( (float)3.5714286566e-01 )
#define   L1 ( (float)6.0000002384e-01 )
#define   L2 ( (float)4.2857143283e-01 )
#define   L3 ( (float)3.3333334327e-01 )
#define   L4 ( (float)2.7272811532e-01 )
#define   L5 ( (float)2.3066075146e-01 )
#define   L6 ( (float)2.0697501302e-01 )
#define   Lg1 ( (float)6.6666668653e-01 )
#define   Lg2 ( (float)4.0000000596e-01 )
#define   Lg3 ( (float)2.8571429849e-01 )
#define   Lg4 ( (float)2.2222198546e-01 )
#define   Lg5 ( (float)1.8183572590e-01 )
#define   Lg6 ( (float)1.5313838422e-01 )
#define   Lg7 ( (float)1.4798198640e-01 )
#define   Lp1 ( (float)6.6666668653e-01 )
#define   Lp2 ( (float)4.0000000596e-01 )
#define   Lp3 ( (float)2.8571429849e-01 )
#define   Lp4 ( (float)2.2222198546e-01 )
#define   Lp5 ( (float)1.8183572590e-01 )
#define   Lp6 ( (float)1.5313838422e-01 )
#define   Lp7 ( (float)1.4798198640e-01 )
#define   P1 ( (float)1.6666667163e-01 )
#define   P2 ( (float)-2.7777778450e-03 )
#define   P3 ( (float)6.6137559770e-05 )
#define   P4 ( (float)-1.6533901999e-06 )
#define   P5 ( (float)4.1381369442e-08 )
#define   Q1 ( (float)-3.3333335072e-02 )
#define   Q2 ( (float)1.5873016091e-03 )
#define   Q3 ( (float)-7.9365076090e-05 )
#define   Q4 ( (float)4.0082177293e-06 )
#define   Q5 ( (float)-2.0109921195e-07 )
#define   R02 ( (float)1.5625000000e-02 )
#define   R03 ( (float)-1.8997929874e-04 )
#define   R04 ( (float)1.8295404516e-06 )
#define   R05 ( (float)-4.6183270541e-09 )
#define   S01 ( (float)1.5619102865e-02 )
#define   S02 ( (float)1.1692678527e-04 )
#define   S03 ( (float)5.1354652442e-07 )
#define   S04 ( (float)1.1661400734e-09 )
#define   S1 ( (float)-1.6666667163e-01 )
#define   S2 ( (float)8.3333337680e-03 )
#define   S3 ( (float)-1.9841270114e-04 )
#define   S4 ( (float)2.7557314297e-06 )
#define   S5 ( (float)-2.5050759689e-08 )
#define   S6 ( (float)1.5896910177e-10 )
#define   cp ( (float)9.6179670095e-01 )
#define   cp_h ( (float)9.6179199219e-01 )
#define   cp_l ( (float)4.7017383622e-06 )
#define   efx ( (float)1.2837916613e-01 )
#define   efx8 ( (float)1.0270333290e+00 )
#define   erx ( (float)8.4506291151e-01 )
#define   half ( (float)5.0000000000e-01 )
#define   invln2 ( (float)1.4426950216e+00 )
#define   invpio2 ( (float)6.3661980629e-01 )
#define   invsqrtpi ( (float)5.6418961287e-01 )
#define   ivln10 ( (float)4.3429449201e-01 )
#define   ivln2 ( (float)1.4426950216e+00 )
#define   ivln2_h ( (float)1.4426879883e+00 )
#define   ivln2_l ( (float)7.0526075433e-06 )
#define   lg2 ( (float)6.9314718246e-01 )
#define   lg2_h ( (float)6.93145752e-01 )
#define   lg2_l ( (float)1.42860654e-06 )
#define   ln2 ( (float)6.9314718246e-01 )
#define   ln2_hi ( (float)6.9313812256e-01 )
#define   ln2_lo ( (float)9.0580006145e-06 )
#define   log10_2hi ( (float)3.0102920532e-01 )
#define   log10_2lo ( (float)7.9034151668e-07 )
#define   o_threshold ( (float)8.8721679688e+01 )
#define   one ( (float)1.0 )
#define   ovt ( (float)4.2995665694e-08 )
#define   pS0 ( (float)1.6666667163e-01 )
#define   pS1 ( (float)-3.2556581497e-01 )
#define   pS2 ( (float)2.0121252537e-01 )
#define   pS3 ( (float)-4.0055535734e-02 )
#define   pS4 ( (float)7.9153501429e-04 )
#define   pS5 ( (float)3.4793309169e-05 )
#define   pa0 ( (float)-2.3621185683e-03 )
#define   pa1 ( (float)4.1485610604e-01 )
#define   pa2 ( (float)-3.7220788002e-01 )
#define   pa3 ( (float)3.1834661961e-01 )
#define   pa4 ( (float)-1.1089469492e-01 )
#define   pa5 ( (float)3.5478305072e-02 )
#define   pa6 ( (float)-2.1663755178e-03 )
#define   pi ( (float)3.1415925026e+00 )
#define   pi_lo ( (float)1.5099578832e-07 )
#define   pi_o_2 ( (float)1.5707963705e+00 )
#define   pi_o_4 ( (float)7.8539818525e-01 )
#define   pio2_1 ( (float)1.5707855225e+00 )
#define   pio2_1t ( (float)1.0804334124e-05 )
#define   pio2_2 ( (float)1.0804273188e-05 )
#define   pio2_2t ( (float)6.0770999344e-11 )
#define   pio2_3 ( (float)6.0770943833e-11 )
#define   pio2_3t ( (float)6.1232342629e-17 )
#define   pio2_hi ( (float)1.5707962513e+00 )
#define   pio2_lo ( (float)7.5497894159e-08 )
#define   pio4 ( (float)7.8539812565e-01 )
#define   pio4_hi ( (float)7.8539818525e-01 )
#define   pio4lo ( (float)3.7748947079e-08 )
#define   pp0 ( (float)1.2837916613e-01 )
#define   pp1 ( (float)-3.2504209876e-01 )
#define   pp2 ( (float)-2.8481749818e-02 )
#define   pp3 ( (float)-5.7702702470e-03 )
#define   pp4 ( (float)-2.3763017452e-05 )
#define   qS1 ( (float)-2.4033949375e+00 )
#define   qS2 ( (float)2.0209457874e+00 )
#define   qS3 ( (float)-6.8828397989e-01 )
#define   qS4 ( (float)7.7038154006e-02 )
#define   qa1 ( (float)1.0642088205e-01 )
#define   qa2 ( (float)5.4039794207e-01 )
#define   qa3 ( (float)7.1828655899e-02 )
#define   qa4 ( (float)1.2617121637e-01 )
#define   qa5 ( (float)1.3637083583e-02 )
#define   qa6 ( (float)1.1984500103e-02 )
#define   qq1 ( (float)3.9791721106e-01 )
#define   qq2 ( (float)6.5022252500e-02 )
#define   qq3 ( (float)5.0813062117e-03 )
#define   qq4 ( (float)1.3249473704e-04 )
#define   qq5 ( (float)-3.9602282413e-06 )
#define   r00 ( (float)-6.2500000000e-02 )
#define   r01 ( (float)1.4070566976e-03 )
#define   r02 ( (float)-1.5995563444e-05 )
#define   r03 ( (float)4.9672799207e-08 )
#define   ra0 ( (float)-9.8649440333e-03 )
#define   ra1 ( (float)-6.9385856390e-01 )
#define   ra2 ( (float)-1.0558626175e+01 )
#define   ra3 ( (float)-6.2375331879e+01 )
#define   ra4 ( (float)-1.6239666748e+02 )
#define   ra5 ( (float)-1.8460508728e+02 )
#define   ra6 ( (float)-8.1287437439e+01 )
#define   ra7 ( (float)-9.8143291473e+00 )
#define   rb0 ( (float)-9.8649431020e-03 )
#define   rb1 ( (float)-7.9928326607e-01 )
#define   rb2 ( (float)-1.7757955551e+01 )
#define   rb3 ( (float)-1.6063638306e+02 )
#define   rb4 ( (float)-6.3756646729e+02 )
#define   rb5 ( (float)-1.0250950928e+03 )
#define   rb6 ( (float)-4.8351919556e+02 )
#define   s01 ( (float)1.9153760746e-02 )
#define   s02 ( (float)1.8594678841e-04 )
#define   s03 ( (float)1.1771846857e-06 )
#define   s04 ( (float)5.0463624390e-09 )
#define   s05 ( (float)1.2354227016e-11 )
#define   sa1 ( (float)1.9651271820e+01 )
#define   sa2 ( (float)1.3765776062e+02 )
#define   sa3 ( (float)4.3456588745e+02 )
#define   sa4 ( (float)6.4538726807e+02 )
#define   sa5 ( (float)4.2900814819e+02 )
#define   sa6 ( (float)1.0863500214e+02 )
#define   sa7 ( (float)6.5702495575e+00 )
#define   sa8 ( (float)-6.0424413532e-02 )
#define   sb1 ( (float)3.0338060379e+01 )
#define   sb2 ( (float)3.2579251099e+02 )
#define   sb3 ( (float)1.5367296143e+03 )
#define   sb4 ( (float)3.1998581543e+03 )
#define   sb5 ( (float)2.5530502930e+03 )
#define   sb6 ( (float)4.7452853394e+02 )
#define   sb7 ( (float)-2.2440952301e+01 )
#define   shuge ( (float)1.0e37 )
#define   tiny ( (float)1.0e-30 )
#define   tpi ( (float)6.3661974669e-01 )
#define   two ( (float)2.0000000000e+00 )
#define   two24 ( (float)16777216.0 )
#define   two25 ( (float)3.3554432000e+07 )
#define   two8 ( (float)2.5600000000e+02 )
#define   twom100 ( (float)7.8886090522e-31 )
#define   twom25 ( (float)2.9802322388e-08 )
#define   twon8 ( (float)3.9062500000e-03 )
#define   u00 ( (float)-7.3804296553e-02 )
#define   u01 ( (float)1.7666645348e-01 )
#define   u02 ( (float)-1.3818567619e-02 )
#define   u03 ( (float)3.4745343146e-04 )
#define   u04 ( (float)-3.8140706238e-06 )
#define   u05 ( (float)1.9559013964e-08 )
#define   u06 ( (float)-3.9820518410e-11 )
#define   u_threshold ( (float)-1.0397208405e+02 )
#define   v01 ( (float)1.2730483897e-02 )
#define   v02 ( (float)7.6006865129e-05 )
#define   v03 ( (float)2.5915085189e-07 )
#define   v04 ( (float)4.4111031494e-10 )
static const float  rhuge = ((float) 1.0000000000e+30);
static const float  zero = ((float) 0.0);

#define  B1 ( (unsigned int)709958130 )
#define  B2 ( (unsigned int)642849266 )

static const float  PIo2[] =    /* array */
  { (float) 1.5703125000e+00, (float) 4.5776367188e-04,
    (float) 2.5987625122e-05, (float) 7.5437128544e-08,
    (float) 6.0026650317e-11, (float) 7.3896444519e-13,
    (float) 5.3845816694e-15, (float) 5.6378512969e-18,
    (float) 8.3009228831e-20, (float) 3.2756352257e-22,
    (float) 6.3331015649e-25
};
static const float  PS2[5] =    /* array */
  { (float) 2.2220300674e+01, (float) 1.3620678711e+02,
    (float) 2.7047027588e+02, (float) 1.5387539673e+02,
    (float) 1.4657617569e+01
};
static const float  PS3[5] =    /* array */
  { (float) 3.5856033325e+01, (float) 3.6151397705e+02,
    (float) 1.1936077881e+03, (float) 1.1279968262e+03,
    (float) 1.7358093262e+02
};
static const float  PS5[5] =    /* array */
  { (float) 6.0753936768e+01, (float) 1.0512523193e+03,
    (float) 5.9789707031e+03, (float) 9.6254453125e+03,
    (float) 2.4060581055e+03
};
static const float  QS2[6] =    /* array */
  { (float) 3.0365585327e+01, (float) 2.6934811401e+02,
    (float) 8.4478375244e+02, (float) 8.8293585205e+02,
    (float) 2.1266638184e+02, (float) -5.3109550476e+00
};
static const float  QS3[6] =    /* array */
  { (float) 4.8758872986e+01, (float) 7.0968920898e+02,
    (float) 3.7041481934e+03, (float) 6.4604252930e+03,
    (float) 2.5163337402e+03, (float) -1.4924745178e+02
};
static const float  TWO23[2] =  /* array */
  { (float) 8.3886080000e+06, (float) -8.3886080000e+06 };
static const float  T[] =       /* array */
  { (float) 3.3333334327e-01, (float) 1.3333334029e-01,
    (float) 5.3968254477e-02, (float) 2.1869488060e-02,
    (float) 8.8632395491e-03, (float) 3.5920790397e-03,
    (float) 1.4562094584e-03, (float) 5.8804126456e-04,
    (float) 2.4646313977e-04, (float) 7.8179444245e-05,
    (float) 7.1407252108e-05, (float) -1.8558637748e-05,
    (float) 2.5907305826e-05
};
static const float  U0[5] =     /* array */
  { (float) -1.9605709612e-01, (float) 5.0443872809e-02,
    (float) -1.9125689287e-03, (float) 2.3525259166e-05,
    (float) -9.1909917899e-08
};
static const float  V0[5] =     /* array */
  { (float) 1.9916731864e-02, (float) 2.0255257550e-04,
    (float) 1.3560879779e-06, (float) 6.2274145840e-09,
    (float) 1.6655924903e-11
};
static const float  Zero[] =    /* array */
  { (float) 0.0, (float) -0.0 };
static const float  aT[] =      /* array */
  { (float) 3.3333334327e-01, (float) -2.0000000298e-01,
    (float) 1.4285714924e-01, (float) -1.1111110449e-01,
    (float) 9.0908870101e-02, (float) -7.6918758452e-02,
    (float) 6.6610731184e-02, (float) -5.8335702866e-02,
    (float) 4.9768779427e-02, (float) -3.6531571299e-02,
    (float) 1.6285819933e-02
};
static const float  atanhi[] =  /* array */
  { (float) 4.6364760399e-01, (float) 7.8539812565e-01,
    (float) 9.8279368877e-01, (float) 1.5707962513e+00
};
static const float  atanlo[] =  /* array */
  { (float) 5.0121582440e-09, (float) 3.7748947079e-08,
    (float) 3.4473217170e-08, (float) 7.5497894159e-08
};
static const float  bp[] =      /* array */
  { (float) 1.0, (float) 1.5 };
static const float  dp_h[] =    /* array */
  { (float) 0.0, (float) 5.84960938e-01 };
static const float  dp_l[] =    /* array */
  { (float) 0.0, (float) 1.56322085e-06 };
static const float  halF[2] =   /* array */
  { (float) 0.5, (float) -0.5 };
static const float  ln2HI[2] =  /* array */
  { (float) 6.9313812256e-01, (float) -6.9313812256e-01 };
static const float  ln2LO[2] =  /* array */
  { (float) 9.0580006145e-06, (float) -9.0580006145e-06 };
static const float  pR2[6] =    /* array */
  { (float) -8.8753431271e-08, (float) -7.0303097367e-02,
    (float) -1.4507384300e+00, (float) -7.6356959343e+00,
    (float) -1.1193166733e+01, (float) -3.2336456776e+00
};
static const float  pR3[6] =    /* array */
  { (float) -2.5470459075e-09, (float) -7.0311963558e-02,
    (float) -2.4090321064e+00, (float) -2.1965976715e+01,
    (float) -5.8079170227e+01, (float) -3.1447946548e+01
};
static const float  pR5[6] =    /* array */
  { (float) -1.1412546255e-11, (float) -7.0312492549e-02,
    (float) -4.1596107483e+00, (float) -6.7674766541e+01,
    (float) -3.3123129272e+02, (float) -3.4643338013e+02
};
static const float  pR8[6] =    /* array */
  { (float) 0.0000000000e+00, (float) -7.0312500000e-02,
    (float) -8.0816707611e+00, (float) -2.5706311035e+02,
    (float) -2.4852163086e+03, (float) -5.2530439453e+03
};
static const float  pS8[5] =    /* array */
  { (float) 1.1653436279e+02, (float) 3.8337448730e+03,
    (float) 4.0597855469e+04, (float) 1.1675296875e+05,
    (float) 4.7627726562e+04
};
static const float  pr2[6] =    /* array */
  { (float) 1.0771083225e-07, (float) 1.1717621982e-01,
    (float) 2.3685150146e+00, (float) 1.2242610931e+01,
    (float) 1.7693971634e+01, (float) 5.0735230446e+00
};
static const float  pr3[6] =    /* array */
  { (float) 3.0250391081e-09, (float) 1.1718686670e-01,
    (float) 3.9329774380e+00, (float) 3.5119403839e+01,
    (float) 9.1055007935e+01, (float) 4.8559066772e+01
};
static const float  pr5[6] =    /* array */
  { (float) 1.3199052094e-11, (float) 1.1718749255e-01,
    (float) 6.8027510643e+00, (float) 1.0830818176e+02,
    (float) 5.1763616943e+02, (float) 5.2871520996e+02
};
static const float  pr8[6] =    /* array */
  { (float) 0.0000000000e+00, (float) 1.1718750000e-01,
    (float) 1.3239480972e+01, (float) 4.1205184937e+02,
    (float) 3.8747453613e+03, (float) 7.9144794922e+03
};
static const float  ps2[5] =    /* array */
  { (float) 2.1436485291e+01, (float) 1.2529022980e+02,
    (float) 2.3227647400e+02, (float) 1.1767937469e+02,
    (float) 8.3646392822e+00
};
static const float  ps3[5] =    /* array */
  { (float) 3.4791309357e+01, (float) 3.3676245117e+02,
    (float) 1.0468714600e+03, (float) 8.9081134033e+02,
    (float) 1.0378793335e+02
};
static const float  ps5[5] =    /* array */
  { (float) 5.9280597687e+01, (float) 9.9140142822e+02,
    (float) 5.3532670898e+03, (float) 7.8446904297e+03,
    (float) 1.5040468750e+03
};
static const float  ps8[5] =    /* array */
  { (float) 1.1420736694e+02, (float) 3.6509309082e+03,
    (float) 3.6956207031e+04, (float) 9.7602796875e+04,
    (float) 3.0804271484e+04
};
static const float  qR2[6] =    /* array */
  { (float) 1.5044444979e-07, (float) 7.3223426938e-02,
    (float) 1.9981917143e+00, (float) 1.4495602608e+01,
    (float) 3.1666231155e+01, (float) 1.6252708435e+01
};
static const float  qR3[6] =    /* array */
  { (float) 4.3774099900e-09, (float) 7.3241114616e-02,
    (float) 3.3442313671e+00, (float) 4.2621845245e+01,
    (float) 1.7080809021e+02, (float) 1.6673394775e+02
};
static const float  qR5[6] =    /* array */
  { (float) 1.8408595828e-11, (float) 7.3242180049e-02,
    (float) 5.8356351852e+00, (float) 1.3511157227e+02,
    (float) 1.0272437744e+03, (float) 1.9899779053e+03
};
static const float  qR8[6] =    /* array */
  { (float) 0.0000000000e+00, (float) 7.3242187500e-02,
    (float) 1.1768206596e+01, (float) 5.5767340088e+02,
    (float) 8.8591972656e+03, (float) 3.7014625000e+04
};
static const float  qS5[6] =    /* array */
  { (float) 8.2776611328e+01, (float) 2.0778142090e+03,
    (float) 1.8847289062e+04, (float) 5.6751113281e+04,
    (float) 3.5976753906e+04, (float) -5.3543427734e+03
};
static const float  qS8[6] =    /* array */
  { (float) 1.6377603149e+02, (float) 8.0983447266e+03,
    (float) 1.4253829688e+05, (float) 8.0330925000e+05,
    (float) 8.4050156250e+05, (float) -3.4389928125e+05
};
static const float  qr2[6] =    /* array */
  { (float) -1.7838172539e-07, (float) -1.0251704603e-01,
    (float) -2.7522056103e+00, (float) -1.9663616180e+01,
    (float) -4.2325313568e+01, (float) -2.1371921539e+01
};
static const float  qr3[6] =    /* array */
  { (float) -5.0783124372e-09, (float) -1.0253783315e-01,
    (float) -4.6101160049e+00, (float) -5.7847221375e+01,
    (float) -2.2824453735e+02, (float) -2.1921012878e+02
};
static const float  qr5[6] =    /* array */
  { (float) -2.0897993405e-11, (float) -1.0253904760e-01,
    (float) -8.0564479828e+00, (float) -1.8366960144e+02,
    (float) -1.3731937256e+03, (float) -2.6124443359e+03
};
static const float  qr8[6] =    /* array */
  { (float) 0.0000000000e+00, (float) -1.0253906250e-01,
    (float) -1.6271753311e+01, (float) -7.5960174561e+02,
    (float) -1.1849806641e+04, (float) -4.8438511719e+04
};
static const float  qs2[6] =    /* array */
  { (float) 2.9533363342e+01, (float) 2.5298155212e+02,
    (float) 7.5750280762e+02, (float) 7.3939318848e+02,
    (float) 1.5594900513e+02, (float) -4.9594988823e+00
};
static const float  qs3[6] =    /* array */
  { (float) 4.7665153503e+01, (float) 6.7386511230e+02,
    (float) 3.3801528320e+03, (float) 5.5477290039e+03,
    (float) 1.9031191406e+03, (float) -1.3520118713e+02
};
static const float  qs5[6] =    /* array */
  { (float) 8.1276550293e+01, (float) 1.9917987061e+03,
    (float) 1.7468484375e+04, (float) 4.9851425781e+04,
    (float) 2.7948074219e+04, (float) -4.7191835938e+03
};
static const float  qs8[6] =    /* array */
  { (float) 1.6139537048e+02, (float) 7.8253862305e+03,
    (float) 1.3387534375e+05, (float) 7.1965775000e+05,
    (float) 6.6660125000e+05, (float) -2.9449025000e+05
};
static const int    init_jk[] = /* array */
  { 4, 7, 9 };
static const int    npio2_hw[] = /* array */
  { 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300,
    0x4116cb00, 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300,
    0x418a3a00, 0x4196cb00, 0x41a35c00, 0x41afed00, 0x41bc7e00,
    0x41c90f00, 0x41d5a000, 0x41e23100, 0x41eec200, 0x41fb5300,
    0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00, 0x421d1400,
    0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
    0x4242c700, 0x42490f00
};
static const int    two_over_pi[] = /* array */
  { 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC, 0x27, 0x57,
    0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62, 0x95, 0x99, 0x3C,
    0x43, 0x90, 0x41, 0xFE, 0x51, 0x63, 0xAB, 0xDE, 0xBB, 0xC5,
    0x61, 0xB7, 0x24, 0x6E, 0x3A, 0x42, 0x4D, 0xD2, 0xE0, 0x06,
    0x49, 0x2E, 0xEA, 0x09, 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB,
    0x1C, 0xB1, 0x29, 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E,
    0xBB, 0x44, 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E,
    0x41, 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
    0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8, 0x97,
    0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11, 0x8B, 0x5A,
    0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF, 0x27, 0xCB, 0x09,
    0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E, 0x5F, 0xEA, 0x2D, 0x75,
    0x27, 0xBA, 0xC7, 0xEB, 0xE5, 0xF1, 0x7B, 0x3D, 0x07, 0x39,
    0xF7, 0x8A, 0x52, 0x92, 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F,
    0x8D, 0x5D, 0x08, 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B,
    0xAB, 0xF0, 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9,
    0xE3, 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
    0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80, 0x4D,
    0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA, 0x73, 0xA8,
    0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B
};

/***********************************************************************/

int                 _rtIEEEferrno;

/***********************************************************************/
static float
pzerof(float x)
{
    const float        *p = (const float *) NULL;
    const float        *q = (const float *) NULL;
    float               z, r, s;
    int                 ix;

    RWFUNCTION(RWSTRING("pzerof"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;
    if (ix >= 0x41000000)
    {
        p = pR8;
        q = pS8;
    }
    else if (ix >= 0x40f71c58)
    {
        p = pR5;
        q = PS5;
    }
    else if (ix >= 0x4036db68)
    {
        p = pR3;
        q = PS3;
    }
    else if (ix >= 0x40000000)
    {
        p = pR2;
        q = PS2;
    }
    z = one / (x * x);
    r = p[0] + z * (p[1] +
                    z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
    s = one + z * (q[0] +
                   z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
    RWRETURN(one + r / s);
}

static float
qzerof(float x)
{
    const float        *p = (const float *) NULL;
    const float        *q = (const float *) NULL;
    float               s, r, z;
    int                 ix;

    RWFUNCTION(RWSTRING("qzerof"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;
    if (ix >= 0x41000000)
    {
        p = qR8;
        q = qS8;
    }
    else if (ix >= 0x40f71c58)
    {
        p = qR5;
        q = qS5;
    }
    else if (ix >= 0x4036db68)
    {
        p = qR3;
        q = QS3;
    }
    else if (ix >= 0x40000000)
    {
        p = qR2;
        q = QS2;
    }
    z = one / (x * x);
    r = p[0] + z * (p[1] +
                    z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
    s = one + z * (q[0] +
                   z * (q[1] +
                        z * (q[2] +
                             z * (q[3] + z * (q[4] + z * q[5])))));
    RWRETURN((-(float) .125 + r / s) / x);
}

static float
ponef(float x)
{
    const float        *p = (const float *) NULL;
    const float        *q = (const float *) NULL;
    float               z, r, s;
    int                 ix;

    RWFUNCTION(RWSTRING("ponef"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;
    if (ix >= 0x41000000)
    {
        p = pr8;
        q = ps8;
    }
    else if (ix >= 0x40f71c58)
    {
        p = pr5;
        q = ps5;
    }
    else if (ix >= 0x4036db68)
    {
        p = pr3;
        q = ps3;
    }
    else if (ix >= 0x40000000)
    {
        p = pr2;
        q = ps2;
    }
    z = one / (x * x);
    r = p[0] + z * (p[1] +
                    z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
    s = one + z * (q[0] +
                   z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
    RWRETURN(one + r / s);
}

static float
qonef(float x)
{
    const float        *p = (const float *) NULL;
    const float        *q = (const float *) NULL;
    float               s, r, z;
    int                 ix;

    RWFUNCTION(RWSTRING("qonef"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;
    if (ix >= 0x40200000)
    {
        p = qr8;
        q = qs8;
    }
    else if (ix >= 0x40f71c58)
    {
        p = qr5;
        q = qs5;
    }
    else if (ix >= 0x4036db68)
    {
        p = qr3;
        q = qs3;
    }
    else if (ix >= 0x40000000)
    {
        p = qr2;
        q = qs2;
    }
    z = one / (x * x);
    r = p[0] + z * (p[1] +
                    z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
    s = one + z * (q[0] +
                   z * (q[1] +
                        z * (q[2] +
                             z * (q[3] + z * (q[4] + z * q[5])))));
    RWRETURN(((float) .375 + r / s) / x);
}

/***********************************************************************/

/* 
 * Standard conformance (non-IEEE) on exception cases.
 * Mapping:
 *      1 -- acos(|x|>1)
 *      2 -- asin(|x|>1)
 *      3 -- atan2(+-0,+-0)
 *      4 -- hypot overflow
 *      5 -- cosh overflow
 *      6 -- exp overflow
 *      7 -- exp underflow
 *      8 -- y0(0)
 *      9 -- y0(-ve)
 *      10-- y1(0)
 *      11-- y1(-ve)
 *      12-- yn(0)
 *      13-- yn(-ve)
 *      14-- lgamma(finite) overflow
 *      15-- lgamma(-integer)
 *      16-- log(0)
 *      17-- log(x<0)
 *      18-- log10(0)
 *      19-- log10(x<0)
 *      20-- pow(0.0,0.0)
 *      21-- pow(x,y) overflow
 *      22-- pow(x,y) underflow
 *      23-- pow(0,negative) 
 *      24-- pow(neg,non-integral)
 *      25-- sinh(finite) overflow
 *      26-- sqrt(negative)
 *      27-- fmod(x,0)
 *      28-- remainder(x,0)
 *      29-- acosh(x<1)
 *      30-- atanh(|x|>1)
 *      31-- atanh(|x|=1)
 *      32-- scalb overflow
 *      33-- scalb underflow
 *      34-- j0(|x|>X_TLOSS)
 *      35-- y0(x>X_TLOSS)
 *      36-- j1(|x|>X_TLOSS)
 *      37-- y1(x>X_TLOSS)
 *      38-- jn(|x|>X_TLOSS, n)
 *      39-- yn(x>X_TLOSS, n)
 *      40-- gamma(finite) overflow
 *      41-- gamma(-integer)
 *      42-- pow(NaN,0.0)
 */

double
_rtIEEEstandard(double x, double __RWUNUSED__ y, int type)
{
    RWFUNCTION(RWSTRING("_rtIEEEstandard"));

    if (type)
    {
        switch (type)
        {
            case 1:           /* acos(|x|>1) */
                break;
            case 2:           /* asin(|x|>1) */
                break;
            case 3:           /* atan2(+-0,+-0) */
                break;
            case 4:           /* hypot overflow */
                break;
            case 5:           /* cosh overflow */
                break;
            case 6:           /* exp overflow */
                break;
            case 7:           /* exp underflow */
                break;
            case 8:           /* y0(0) */
                break;
            case 9:           /* y0(-ve) */
                break;
            case 10:          /* y1(0) */
                break;
            case 11:          /* y1(-ve) */
                break;
            case 12:          /* yn(0) */
                break;
            case 13:          /* yn(-ve) */
                break;
            case 14:          /* lgamma(finite) overflow */
                break;
            case 15:          /* lgamma(-integer) */
                break;
            case 16:          /* log(0) */
                break;
            case 17:          /* log(x<0) */
                break;
            case 18:          /* log10(0) */
                break;
            case 19:          /* log10(x<0) */
                break;
            case 20:          /* pow(0.0,0.0) */
                break;
            case 21:          /* pow(x,y) overflow */
                break;
            case 22:          /* pow(x,y) underflow */
                break;
            case 23:          /* pow(0,negative)  */
                break;
            case 24:          /* pow(neg,non-integral) */
                break;
            case 25:          /* sinh(finite) overflow */
                break;
            case 26:          /* sqrt(negative) */
                break;
            case 27:          /* fmod(x,0) */
                break;
            case 28:          /* remainder(x,0) */
                break;
            case 29:          /* acosh(x<1) */
                break;
            case 30:          /* atanh(|x|>1) */
                break;
            case 31:          /* atanh(|x|=1) */
                break;
            case 32:          /* scalb overflow */
                break;
            case 33:          /* scalb underflow */
                break;
            case 34:          /* j0(|x|>X_TLOSS) */
                break;
            case 35:          /* y0(x>X_TLOSS) */
                break;
            case 36:          /* j1(|x|>X_TLOSS) */
                break;
            case 37:          /* y1(x>X_TLOSS) */
                break;
            case 38:          /* jn(|x|>X_TLOSS, n) */
                break;
            case 39:          /* yn(x>X_TLOSS, n) */
                break;
            case 40:          /* gamma(finite) overflow */
                break;
            case 41:          /* gamma(-integer) */
                break;
            case 42:          /* pow(NaN,0.0) */
                break;

        }
    }

    RWRETURN(x);
}

float
_rtIEEEacosf(float x)
{
    float               z, p, q, r, w, s, c, df;
    int                 hx, ix;

    RWFUNCTION(RWSTRING("_rtIEEEacosf"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix == 0x3f800000)
    {                          /* |x|==1 */
        if (hx > 0)
            RWRETURN(0.0);     /* acos(1) = 0  */
        else
            RWRETURN(pi + (float) 2.0 * pio2_lo); /* acos(-1)= pi */
    }
    else if (ix > 0x3f800000)
    {                          /* |x| >= 1 */
        RWRETURN((x - x) / (x - x)); /* acos(|x|>1) is NaN */
    }
    if (ix < 0x3f000000)
    {                          /* |x| < 0.5 */
        if (ix <= 0x23000000)
            RWRETURN(pio2_hi + pio2_lo); /*if|x|<2**-57 */
        z = x * x;
        p = z * (pS0 +
                 z * (pS1 +
                      z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
        q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
        r = p / q;
        RWRETURN(pio2_hi - (x - (pio2_lo - x * r)));
    }
    else if (hx < 0)
    {                          /* x < -0.5 */
        z = (one + x) * (float) 0.5;
        p = z * (pS0 +
                 z * (pS1 +
                      z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
        q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
        s = RtIEEEsqrtf(z);
        r = p / q;
        w = r * s - pio2_lo;
        RWRETURN(pi - (float) 2.0 * (s + w));
    }
    else
    {                          /* x > 0.5 */
        int                 idf;

        z = (one - x) * (float) 0.5;
        s = RtIEEEsqrtf(z);
        df = s;
        GET_FLOAT_WORD(idf, df);
        SET_FLOAT_WORD(df, idf & 0xfffff000);
        c = (z - df * df) / (s + df);
        p = z * (pS0 +
                 z * (pS1 +
                      z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
        q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
        r = p / q;
        w = r * s + c;
        RWRETURN((float) 2.0 * (df + w));
    }
}

float
_rtIEEEacoshf(float x)
{
    float               t;
    int                 hx;

    RWFUNCTION(RWSTRING("_rtIEEEacoshf"));

    GET_FLOAT_WORD(hx, x);
    if (hx < 0x3f800000)
    {                          /* x < 1 */
        RWRETURN((x - x) / (x - x));
    }
    else if (hx >= 0x4d800000)
    {                          /* x > 2**28 */
        if (hx >= 0x7f800000)
        {                      /* x is inf of NaN */
            RWRETURN(x + x);
        }
        else
            RWRETURN(RtIEEElogf(x) + ln2); /* acosh(rhuge)=log(2x) */
    }
    else if (hx == 0x3f800000)
    {
        RWRETURN(0.0);         /* acosh(1) = 0 */
    }
    else if (hx > 0x40000000)
    {                          /* 2**28 > x > 2 */
        t = x * x;
        RWRETURN(RtIEEElogf((float) 2.0 * x -
                            one / (x + RtIEEEsqrtf(t - one))));
    }
    else
    {                          /* 1<x<2 */
        t = x - one;
        RWRETURN(RtIEEElog1pf
                 (t + RtIEEEsqrtf((float) 2.0 * t + t * t)));
    }
}

float
_rtIEEEasinf(float x)
{
    float               t = 0, w, p, q, c, r, s;
    int                 hx, ix;

    RWFUNCTION(RWSTRING("_rtIEEEasinf"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix == 0x3f800000)
    {
        /* asin(1)=+-pi/2 with inexact */
        RWRETURN(x * pio2_hi + x * pio2_lo);
    }
    else if (ix > 0x3f800000)
    {                          /* |x|>= 1 */
        RWRETURN((x - x) / (x - x)); /* asin(|x|>1) is NaN */
    }
    else if (ix < 0x3f000000)
    {                          /* |x|<0.5 */
        if (ix < 0x32000000)
        {                      /* if |x| < 2**-27 */
            if (rhuge + x > one)
                RWRETURN(x);   /* return x with inexact if x!=0 */
        }
        else
            t = x * x;
        p = t * (pS0 +
                 t * (pS1 +
                      t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
        q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
        w = p / q;
        RWRETURN(x + x * w);
    }
    /* 1> |x|>= 0.5 */
    w = one - RtIEEEfabsf(x);
    t = w * (float) 0.5;
    p = t * (pS0 +
             t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
    q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
    s = RtIEEEsqrtf(t);
    if (ix >= 0x3F79999A)
    {                          /* if |x| > 0.975 */
        w = p / q;
        t = pio2_hi - ((float) 2.0 * (s + s * w) - pio2_lo);
    }
    else
    {
        int                 iw;

        w = s;
        GET_FLOAT_WORD(iw, w);
        SET_FLOAT_WORD(w, iw & 0xfffff000);
        c = (t - w * w) / (s + w);
        r = p / q;
        p = (float) 2.0    *s * r - (pio2_lo - (float) 2.0 * c);
        q = pio4_hi - (float) 2.0 *w;

        t = pio4_hi - (p - q);
    }
    if (hx > 0)
        RWRETURN(t);
    else
        RWRETURN(-t);
}

float
_rtIEEEatanhf(float x)
{
    float               t;
    int                 hx, ix;

    RWFUNCTION(RWSTRING("_rtIEEEatanhf"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix > 0x3f800000)       /* |x|>1 */
        RWRETURN((x - x) / (x - x));
    if (ix == 0x3f800000)
        RWRETURN(x / zero);
    if (ix < 0x31800000 && (rhuge + x) > zero)
        RWRETURN(x);           /* x<2**-28 */
    SET_FLOAT_WORD(x, ix);
    if (ix < 0x3f000000)
    {                          /* x < 0.5 */
        t = x + x;
        t = (float) 0.5    *RtIEEElog1pf(t + t * x / (one - x));
    }
    else
        t = (float) 0.5    *RtIEEElog1pf((x + x) / (one - x));

    if (hx >= 0)
        RWRETURN(t);
    else
        RWRETURN(-t);
}

float
_rtIEEEcosf(float x, float y)
{
    float               a, hz, z, r, qx;
    int                 ix;

    RWFUNCTION(RWSTRING("_rtIEEEcosf"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;          /* ix = |x|'s high word */
    if (ix < 0x32000000)
    {                          /* if x < 2**27 */
        if (((int) x) == 0)
            RWRETURN(one);     /* generate inexact */
    }
    z = x * x;
    r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
    if (ix < 0x3e99999a)       /* if |x| < 0.3 */
        RWRETURN(one - ((float) 0.5 * z - (z * r - x * y)));
    else
    {
        if (ix > 0x3f480000)
        {                      /* x > 0.78125 */
            qx = (float) 0.28125;
        }
        else
        {
            SET_FLOAT_WORD(qx, ix - 0x01000000); /* x/4 */
        }
        hz = (float) 0.5   *z - qx;

        a = one - qx;
        RWRETURN(a - (hz - (z * r - x * y)));
    }
}

int
_rtIEEErem_pio2f(float *x, float *y, int e0, int nx, int prec,
                 const int *ipio2)
{
    int                 jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k,
        m, q0, ih;
    float               z, fw, f[20], fq[20], q[20];

    RWFUNCTION(RWSTRING("_rtIEEErem_pio2f"));

    /* initialize jk */
    jk = init_jk[prec];
    jp = jk;

    /* determine jx,jv,q0, note that 3>q0 */
    jx = nx - 1;
    jv = (e0 - 3) / 8;
    if (jv < 0)
        jv = 0;
    q0 = e0 - 8 * (jv + 1);

    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
    j = jv - jx;
    m = jx + jk;
    for (i = 0; i <= m; i++, j++)
        f[i] = (j < 0) ? zero : (float) ipio2[j];

    /* compute q[0],q[1],...q[jk] */
    for (i = 0; i <= jk; i++)
    {
        for (j = 0, fw = 0.0; j <= jx; j++)
            fw += x[j] * f[jx + i - j];
        q[i] = fw;
    }

    jz = jk;
  recompute:
    /* distill q[] into iq[] reversingly */
    for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
    {
        fw = (float) ((int) (twon8 * z));
        iq[i] = (int) (z - two8 * fw);
        z = q[j - 1] + fw;
    }

    /* compute n */
    z = RtIEEEscalbnf(z, q0);  /* actual value of z */
    z -= (float) 8.0   *RtIEEEfloorf(z * (float) 0.125); /* trim off integer >= 8 */

    n = (int) z;
    z -= (float) n;
    ih = 0;
    if (q0 > 0)
    {                          /* need iq[jz-1] to determine n */
        i = (iq[jz - 1] >> (8 - q0));
        n += i;
        iq[jz - 1] -= i << (8 - q0);
        ih = iq[jz - 1] >> (7 - q0);
    }
    else if (q0 == 0)
        ih = iq[jz - 1] >> 8;
    else if (z >= (float) 0.5)
        ih = 2;

    if (ih > 0)
    {                          /* q > 0.5 */
        n += 1;
        carry = 0;
        for (i = 0; i < jz; i++)
        {                      /* compute 1-q */
            j = iq[i];
            if (carry == 0)
            {
                if (j != 0)
                {
                    carry = 1;
                    iq[i] = 0x100 - j;
                }
            }
            else
                iq[i] = 0xff - j;
        }
        if (q0 > 0)
        {                      /* rare case: chance is 1 in 12 */
            switch (q0)
            {
                case 1:
                    iq[jz - 1] &= 0x7f;
                    break;
                case 2:
                    iq[jz - 1] &= 0x3f;
                    break;
            }
        }
        if (ih == 2)
        {
            z = one - z;
            if (carry != 0)
                z -= RtIEEEscalbnf(one, q0);
        }
    }

    /* check if recomputation is needed */
    if (z == zero)
    {
        j = 0;
        for (i = jz - 1; i >= jk; i--)
            j |= iq[i];
        if (j == 0)
        {                      /* need recomputation */
            for (k = 1; iq[jk - k] == 0; k++) ; /* k = no. of terms needed */

            for (i = jz + 1; i <= jz + k; i++)
            {                  /* add q[jz+1] to q[jz+k] */
                f[jx + i] = (float) ipio2[jv + i];
                for (j = 0, fw = 0.0; j <= jx; j++)
                    fw += x[j] * f[jx + i - j];
                q[i] = fw;
            }
            jz += k;
            goto recompute;
        }
    }

    /* chop off zero terms */
    if (z == (float) 0.0)
    {
        jz -= 1;
        q0 -= 8;
        while (iq[jz] == 0)
        {
            jz--;
            q0 -= 8;
        }
    }
    else
    {                          /* break z into 8-bit if necessary */
        z = RtIEEEscalbnf(z, -q0);
        if (z >= two8)
        {
            fw = (float) ((int) (twon8 * z));
            iq[jz] = (int) (z - two8 * fw);
            jz += 1;
            q0 += 8;
            iq[jz] = (int) fw;
        }
        else
            iq[jz] = (int) z;
    }

    /* convert integer "bit" chunk to floating-point value */
    fw = RtIEEEscalbnf(one, q0);
    for (i = jz; i >= 0; i--)
    {
        q[i] = fw * (float) iq[i];
        fw *= twon8;
    }

    /* compute PIo2[0,...,jp]*q[jz,...,0] */
    for (i = jz; i >= 0; i--)
    {
        for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
            fw += PIo2[k] * q[i + k];
        fq[jz - i] = fw;
    }

    /* compress fq[] into y[] */
    switch (prec)
    {
        case 0:
            fw = 0.0;
            for (i = jz; i >= 0; i--)
                fw += fq[i];
            y[0] = (ih == 0) ? fw : -fw;
            break;
        case 1:
        case 2:
            fw = 0.0;
            for (i = jz; i >= 0; i--)
                fw += fq[i];
            y[0] = (ih == 0) ? fw : -fw;
            fw = fq[0] - fw;
            for (i = 1; i <= jz; i++)
                fw += fq[i];
            y[1] = (ih == 0) ? fw : -fw;
            break;
        case 3:               /* painful */
            for (i = jz; i > 0; i--)
            {
                fw = fq[i - 1] + fq[i];
                fq[i] += fq[i - 1] - fw;
                fq[i - 1] = fw;
            }
            for (i = jz; i > 1; i--)
            {
                fw = fq[i - 1] + fq[i];
                fq[i] += fq[i - 1] - fw;
                fq[i - 1] = fw;
            }
            for (fw = 0.0, i = jz; i >= 2; i--)
                fw += fq[i];
            if (ih == 0)
            {
                y[0] = fq[0];
                y[1] = fq[1];
                y[2] = fw;
            }
            else
            {
                y[0] = -fq[0];
                y[1] = -fq[1];
                y[2] = -fw;
            }
    }
    RWRETURN(n & 7);
}

float
_rtIEEEsinf(float x, float y, int iy)
{
    float               z, r, v;
    int                 ix;

    RWFUNCTION(RWSTRING("_rtIEEEsinf"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;          /* high word of x */
    if (ix < 0x32000000)       /* |x| < 2**-27 */
    {
        if ((int) x == 0)
            RWRETURN(x);
    }                          /* generate inexact */
    z = x * x;
    v = z * x;
    r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
    if (iy == 0)
        RWRETURN(x + v * (S1 + z * r));
    else
        RWRETURN(x - ((z * (half * y - v * r) - y) - v * S1));
}

float
_rtIEEEtanf(float x, float y, int iy)
{
    float               z, r, v, w, s;
    int                 ix, hx;

    RWFUNCTION(RWSTRING("_rtIEEEtanf"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;      /* high word of |x| */
    if (ix < 0x31800000)       /* x < 2**-28 */
    {
        if ((int) x == 0)
        {                      /* generate inexact */
            if ((ix | (iy + 1)) == 0)
                RWRETURN(one / RtIEEEfabsf(x));
            else
                RWRETURN((iy == 1) ? x : -one / x);
        }
    }
    if (ix >= 0x3f2ca140)
    {                          /* |x|>=0.6744 */
        if (hx < 0)
        {
            x = -x;
            y = -y;
        }
        z = pio4 - x;
        w = pio4lo - y;
        x = z + w;
        y = 0.0;
    }
    z = x * x;
    w = z * z;
    /* Break x^5*(T[1]+x^2*T[2]+...) into
     *   x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
     *   x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
     */
    r = T[1] + w * (T[3] +
                    w * (T[5] + w * (T[7] + w * (T[9] + w * T[11]))));
    v = z * (T[2] +
             w * (T[4] +
                  w * (T[6] + w * (T[8] + w * (T[10] + w * T[12])))));
    s = z * x;
    r = y + z * (s * (r + v) + y);
    r += T[0] * s;
    w = x + r;
    if (ix >= 0x3f2ca140)
    {
        v = (float) iy;
        RWRETURN((float) (1 - ((hx >> 30) & 2)) *
                 (v - (float) 2.0 * (x - (w * w / (w + v) - r))));
    }
    if (iy == 1)
        RWRETURN(w);
    else
    {                          /* if allow error up to 2 ulp, 
                                * simply return -1.0/(x+r) here */
        /*  compute -1.0/(x+r) accurately */
        float               a, t;
        int                 i;

        z = w;
        GET_FLOAT_WORD(i, z);
        SET_FLOAT_WORD(z, i & 0xfffff000);
        v = r - (z - x);       /* z+v = r+x */
        t = a = -(float) 1.0 / w; /* a = -1.0/w */
        GET_FLOAT_WORD(i, t);
        SET_FLOAT_WORD(t, i & 0xfffff000);
        s = (float) 1.0 + t * z;
        RWRETURN(t + a * (s + t * v));
    }
}

float
_rtIEEEatan2f(float y, float x)
{
    float               z;
    int                 k, m, hx, hy, ix, iy;

    RWFUNCTION(RWSTRING("_rtIEEEatan2f"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    GET_FLOAT_WORD(hy, y);
    iy = hy & 0x7fffffff;
    if ((ix > 0x7f800000) || (iy > 0x7f800000)) /* x or y is NaN */
        RWRETURN(x + y);
    if (hx == 0x3f800000)
        RWRETURN((float) RtIEEEatanf(y)); /* x=1.0 */
    m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */

    /* when y = 0 */
    if (iy == 0)
    {
        switch (m)
        {
            case 0:
            case 1:
                RWRETURN(y);   /* atan(+-0,+anything)=+-0 */
            case 2:
                RWRETURN(pi + tiny); /* atan(+0,-anything) = pi */
            case 3:
                RWRETURN(-pi - tiny); /* atan(-0,-anything) =-pi */
        }
    }
    /* when x = 0 */
    if (ix == 0)
        RWRETURN((hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny);

    /* when x is INF */
    if (ix == 0x7f800000)
    {
        if (iy == 0x7f800000)
        {
            switch (m)
            {
                case 0:
                    RWRETURN(pi_o_4 + tiny); /* atan(+INF,+INF) */
                case 1:
                    RWRETURN(-pi_o_4 - tiny); /* atan(-INF,+INF) */
                case 2:
                    RWRETURN((float) 3.0 * pi_o_4 + tiny); /*atan(+INF,-INF) */

                case 3:
                    RWRETURN((float) -3.0 * pi_o_4 - tiny); /*atan(-INF,-INF) */
            }
        }
        else
        {
            switch (m)
            {
                case 0:
                    RWRETURN(zero); /* atan(+...,+INF) */
                case 1:
                    RWRETURN(-zero); /* atan(-...,+INF) */
                case 2:
                    RWRETURN(pi + tiny); /* atan(+...,-INF) */
                case 3:
                    RWRETURN(-pi - tiny); /* atan(-...,-INF) */
            }
        }
    }
    /* when y is INF */
    if (iy == 0x7f800000)
        RWRETURN((hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny);

    /* compute y/x */
    k = (iy - ix) >> 23;
    if (k > 60)
        z = pi_o_2 + (float) 0.5 *pi_lo; /* |y/x| >  2**60 */

    else if (hx < 0 && k < -60)
        z = 0.0;               /* |y|/x < -2**60 */
    else
        z = (float) RtIEEEatanf(RtIEEEfabsf(y / x)); /* safe to do y/x */
    switch (m)
    {
        case 0:
            RWRETURN(z);       /* atan(+,+) */
        case 1:
            {
                unsigned int        zh;

                GET_FLOAT_WORD(zh, z);
                SET_FLOAT_WORD(z, zh ^ 0x80000000);
            }
            RWRETURN(z);       /* atan(-,+) */
        case 2:
            RWRETURN(pi - (z - pi_lo)); /* atan(+,-) */
        default:              /* case 3 */
            RWRETURN((z - pi_lo) - pi); /* atan(-,-) */
    }
}

/************************************************************************/

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEcoshf
 * computes the hyperbolic cosine
 * \param x input
 * \return the hyperbolic cosine
 */
float
RtIEEEcoshf(float x)
{
    float               t, w;
    int                 ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEcoshf"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;

    /* x is INF or NaN */
    if (ix >= 0x7f800000)
        RWRETURN(x * x);

    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
    if (ix < 0x3eb17218)
    {
        t = RtIEEEexpm1f(RtIEEEfabsf(x));
        w = one + t;
        if (ix < 0x24000000)
            RWRETURN(w);       /* cosh(tiny) = 1 */
        RWRETURN(one + (t * t) / (w + w));
    }

    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2 */
    if (ix < 0x41b00000)
    {
        t = RtIEEEexpf(RtIEEEfabsf(x));
        RWRETURN(half * t + half / t);
    }

    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
    if (ix < 0x42b17180)
        RWRETURN(half * RtIEEEexpf(RtIEEEfabsf(x)));

    /* |x| in [log(maxdouble), overflowthresold] */
    if (ix <= 0x42b2d4fc)
    {
        w = RtIEEEexpf(half * RtIEEEfabsf(x));
        t = half * w;
        RWRETURN(t * w);
    }

    /* |x| > overflowthresold, cosh(x) overflow */
    RWRETURN(rhuge * rhuge);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEfmodf
 * returns the floating-point remainder of x / y.
 * \param x numerator
 * \param y denominator
 * \return the remainder.
 */
float
RtIEEEfmodf(float x, float y)
{
    int                 n, hx, hy, hz, ix, iy, sx, i;

    RWAPIFUNCTION(RWSTRING("RtIEEEfmodf"));

    GET_FLOAT_WORD(hx, x);
    GET_FLOAT_WORD(hy, y);
    sx = hx & 0x80000000;      /* sign of x */
    hx ^= sx;                  /* |x| */
    hy &= 0x7fffffff;          /* |y| */

    /* purge off exception values */
    if (hy == 0 || (hx >= 0x7f800000) || /* y=0,or x not finite */
        (hy > 0x7f800000))     /* or y is NaN */
        RWRETURN((x * y) / (x * y));
    if (hx < hy)
        RWRETURN(x);           /* |x|<|y| return x */
    if (hx == hy)
        RWRETURN(Zero[(unsigned int) sx >> 31]); /* |x|=|y| return x*0 */

    /* determine ix = ilogb(x) */
    if (hx < 0x00800000)
    {                          /* subnormal x */
        for (ix = -126, i = (hx << 8); i > 0; i <<= 1)
            ix -= 1;
    }
    else
        ix = (hx >> 23) - 127;

    /* determine iy = ilogb(y) */
    if (hy < 0x00800000)
    {                          /* subnormal y */
        for (iy = -126, i = (hy << 8); i >= 0; i <<= 1)
            iy -= 1;
    }
    else
        iy = (hy >> 23) - 127;

    /* set up {hx,lx}, {hy,ly} and align y to x */
    if (ix >= -126)
        hx = 0x00800000 | (0x007fffff & hx);
    else
    {                          /* subnormal x, shift x to normal */
        n = -126 - ix;
        hx = hx << n;
    }
    if (iy >= -126)
        hy = 0x00800000 | (0x007fffff & hy);
    else
    {                          /* subnormal y, shift y to normal */
        n = -126 - iy;
        hy = hy << n;
    }

    /* fix point fmod */
    n = ix - iy;
    while (n--)
    {
        hz = hx - hy;
        if (hz < 0)
        {
            hx = hx + hx;
        }
        else
        {
            if (hz == 0)       /* return sign(x)*0 */
                RWRETURN(Zero[(unsigned int) sx >> 31]);
            hx = hz + hz;
        }
    }
    hz = hx - hy;
    if (hz >= 0)
    {
        hx = hz;
    }

    /* convert back to floating value and restore the sign */
    if (hx == 0)               /* return sign(x)*0 */
        RWRETURN(Zero[(unsigned int) sx >> 31]);
    while (hx < 0x00800000)
    {                          /* normalize x */
        hx = hx + hx;
        iy -= 1;
    }
    if (iy >= -126)
    {                          /* normalize output */
        hx = ((hx - 0x00800000) | ((iy + 127) << 23));
        SET_FLOAT_WORD(x, hx | sx);
    }
    else
    {                          /* subnormal output */
        n = -126 - iy;
        hx >>= n;
        SET_FLOAT_WORD(x, hx | sx);
        x *= one;              /* create necessary signal */
    }
    RWRETURN(x);               /* exact output */
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEhypotf
 * computes the hypoteneuse sqrt(x*x+y*y)
 * in such a way that underflow will not happen, and overflow
 * occurs only if the final result deserves it.
 * \param x 
 * \param y
 * \return the hypoteneuse.
 */
float
RtIEEEhypotf(float x, float y)
{
    float               a = x, b = y, t1, t2, y1, y2, w;
    int                 j, k, ha, hb;

    RWAPIFUNCTION(RWSTRING("RtIEEEhypotf"));

    GET_FLOAT_WORD(ha, x);
    ha &= 0x7fffffff;
    GET_FLOAT_WORD(hb, y);
    hb &= 0x7fffffff;
    if (hb > ha)
    {
        a = y;
        b = x;
        j = ha;
        ha = hb;
        hb = j;
    }
    else
    {
        a = x;
        b = y;
    }
    SET_FLOAT_WORD(a, ha);     /* a <- |a| */
    SET_FLOAT_WORD(b, hb);     /* b <- |b| */
    if ((ha - hb) > 0xf000000)
    {
        RWRETURN(a + b);
    }                          /* x/y > 2**30 */
    k = 0;
    if (ha > 0x58800000)
    {                          /* a>2**50 */
        if (ha >= 0x7f800000)
        {                      /* Inf or NaN */
            w = a + b;         /* for sNaN */
            if (ha == 0x7f800000)
                w = a;
            if (hb == 0x7f800000)
                w = b;
            RWRETURN(w);
        }
        /* scale a and b by 2**-60 */
        ha -= 0x5d800000;
        hb -= 0x5d800000;
        k += 60;
        SET_FLOAT_WORD(a, ha);
        SET_FLOAT_WORD(b, hb);
    }
    if (hb < 0x26800000)
    {                          /* b < 2**-50 */
        if (hb <= 0x007fffff)
        {                      /* subnormal b or 0 */
            if (hb == 0)
                RWRETURN(a);
            SET_FLOAT_WORD(t1, 0x3f000000); /* t1=2^126 */
            b *= t1;
            a *= t1;
            k -= 126;
        }
        else
        {                      /* scale a and b by 2^60 */
            ha += 0x5d800000;  /* a *= 2^60 */
            hb += 0x5d800000;  /* b *= 2^60 */
            k -= 60;
            SET_FLOAT_WORD(a, ha);
            SET_FLOAT_WORD(b, hb);
        }
    }
    /* medium size a and b */
    w = a - b;
    if (w > b)
    {
        SET_FLOAT_WORD(t1, ha & 0xfffff000);
        t2 = a - t1;
        w = RtIEEEsqrtf(t1 * t1 - (b * (-b) - t2 * (a + t1)));
    }
    else
    {
        a = a + a;
        SET_FLOAT_WORD(y1, hb & 0xfffff000);
        y2 = b - y1;
        SET_FLOAT_WORD(t1, ha + 0x00800000);
        t2 = a - t1;
        w = RtIEEEsqrtf(t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)));
    }
    if (k != 0)
    {
        SET_FLOAT_WORD(t1, 0x3f800000 + (k << 23));
        RWRETURN(t1 * w);
    }
    else
        RWRETURN(w);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEj0f
 * computes the Bessel function of the first kind of order 0
 * \param x input
 * \return the Bessel function value
 */
float
RtIEEEj0f(float x)
{
    float               z, s, c, ss, cc, r, u, v;
    int                 hx;
    unsigned int        ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEj0f"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix >= 0x7f800000)
        RWRETURN(one / (x * x));
    x = RtIEEEfabsf(x);
    if (ix >= 0x40000000)
    {                          /* |x| >= 2.0 */
        s = RtIEEEsinf(x);
        c = RtIEEEcosf(x);
        ss = s - c;
        cc = s + c;
        if (ix < 0x7f000000)
        {                      /* make sure x+x not overflow */
            z = -RtIEEEcosf(x + x);
            if ((s * c) < zero)
                cc = z / ss;
            else
                ss = z / cc;
        }
        /*
         * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
         * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
         */
        if (ix > 0x80000000)
            z = (invsqrtpi * cc) / RtIEEEsqrtf(x);
        else
        {
            u = pzerof(x);
            v = qzerof(x);
            z = invsqrtpi * (u * cc - v * ss) / RtIEEEsqrtf(x);
        }
        RWRETURN(z);
    }
    if (ix < 0x39000000)
    {                          /* |x| < 2**-13 */
        if (rhuge + x > one)
        {                      /* raise inexact if x != 0 */
            if (ix < 0x32000000)
                RWRETURN(one); /* |x|<2**-27 */
            else
                RWRETURN(one - (float) 0.25 * x * x);
        }
    }
    z = x * x;
    r = z * (R02 + z * (R03 + z * (R04 + z * R05)));
    s = one + z * (S01 + z * (S02 + z * (S03 + z * S04)));
    if (ix < 0x3F800000)
    {                          /* |x| < 1.00 */
        RWRETURN(one + z * ((float) -0.25 + (r / s)));
    }
    else
    {
        u = (float) 0.5    *x;

        RWRETURN(((one + u) * (one - u) + z * (r / s)));
    }
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEy0f
 * computes the linearly independent
 * Bessel function of the second kind of order
 * 0 for the
 * positive integer value
 * expressed as a double.
 * \param x input
 * \return the Bessel function value
 */
float
RtIEEEy0f(float x)
{
    float               z, s, c, ss, cc, u, v;
    int                 hx;
    unsigned int        ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEy0f"));

    GET_FLOAT_WORD(hx, x);
    ix = 0x7fffffff & hx;
    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
    if (ix >= 0x7f800000)
        RWRETURN(one / (x + x * x));
    if (ix == 0)
        RWRETURN(-one / zero);
    if (hx < 0)
        RWRETURN(zero / zero);
    if (ix >= 0x40000000)
    {                          /* |x| >= 2.0 */
        /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
         * where x0 = x-pi/4
         *      Better formula:
         *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
         *                      =  1/sqrt(2) * (sin(x) + cos(x))
         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
         *                      =  1/sqrt(2) * (sin(x) - cos(x))
         * To avoid cancellation, use
         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
         * to compute the worse one.
         */
        s = RtIEEEsinf(x);
        c = RtIEEEcosf(x);
        ss = s - c;
        cc = s + c;
        /*
         * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
         * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
         */
        if (ix < 0x7f000000)
        {                      /* make sure x+x not overflow */
            z = -RtIEEEcosf(x + x);
            if ((s * c) < zero)
                cc = z / ss;
            else
                ss = z / cc;
        }
        if (ix > 0x80000000)
            z = (invsqrtpi * ss) / RtIEEEsqrtf(x);
        else
        {
            u = pzerof(x);
            v = qzerof(x);
            z = invsqrtpi * (u * ss + v * cc) / RtIEEEsqrtf(x);
        }
        RWRETURN(z);
    }
    if (ix <= 0x32000000)
    {                          /* x < 2**-27 */
        RWRETURN((u00 + tpi * RtIEEElogf(x)));
    }
    z = x * x;
    u = u00 + z * (u01 +
                   z * (u02 +
                        z * (u03 + z * (u04 + z * (u05 + z * u06)))));
    v = one + z * (v01 + z * (v02 + z * (v03 + z * v04)));
    RWRETURN((u / v + tpi * (RtIEEEj0f(x) * RtIEEElogf(x))));
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEj1f
 * computes the
 * Bessel function of the first kind of order
 * 1 for the input value.
 * \param x input
 * \return the Bessel function value
 */
float
RtIEEEj1f(float x)
{
    float               z, s, c, ss, cc, r, u, v, y;
    int                 hx;
    unsigned int        ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEj1f"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix >= 0x7f800000)
        RWRETURN(one / x);
    y = RtIEEEfabsf(x);
    if (ix >= 0x40000000)
    {                          /* |x| >= 2.0 */
        s = RtIEEEsinf(y);
        c = RtIEEEcosf(y);
        ss = -s - c;
        cc = s - c;
        if (ix < 0x7f000000)
        {                      /* make sure y+y not overflow */
            z = RtIEEEcosf(y + y);
            if ((s * c) > zero)
                cc = z / ss;
            else
                ss = z / cc;
        }
        /*
         * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
         * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
         */
        if (ix > 0x80000000)
            z = (invsqrtpi * cc) / RtIEEEsqrtf(y);
        else
        {
            u = ponef(y);
            v = qonef(y);
            z = invsqrtpi * (u * cc - v * ss) / RtIEEEsqrtf(y);
        }
        if (hx < 0)
            RWRETURN(-z);
        else
            RWRETURN(z);
    }
    if (ix < 0x32000000)
    {                          /* |x|<2**-27 */
        if (rhuge + x > one)
            RWRETURN((float) 0.5 * x); /* inexact if x!=0 necessary */
    }
    z = x * x;
    r = z * (r00 + z * (r01 + z * (r02 + z * r03)));
    s = one + z * (s01 + z * (s02 + z * (s03 + z * (s04 + z * s05))));
    r *= x;
    RWRETURN((x * (float) 0.5 + r / s));
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEy1f
 * computes the linearly independent
 * Bessel function of the second kind of order
 * 1 for the positive integer value (expressed as a float)
 * \param x input
 * \return the Bessel function value

 */
float
RtIEEEy1f(float x)
{
    float               z, s, c, ss, cc, u, v;
    int                 hx, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEy1f"));

    GET_FLOAT_WORD(hx, x);
    ix = 0x7fffffff & hx;
    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
    if (ix >= 0x7f800000)
        RWRETURN(one / (x + x * x));
    if (ix == 0)
        RWRETURN(-one / zero);
    if (hx < 0)
        RWRETURN(zero / zero);
    if (ix >= 0x40000000)
    {                          /* |x| >= 2.0 */
        s = RtIEEEsinf(x);
        c = RtIEEEcosf(x);
        ss = -s - c;
        cc = s - c;
        if (ix < 0x7f000000)
        {                      /* make sure x+x not overflow */
            z = RtIEEEcosf(x + x);
            if ((s * c) > zero)
                cc = z / ss;
            else
                ss = z / cc;
        }
        /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
         * where x0 = x-3pi/4
         *      Better formula:
         *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
         *                      =  1/sqrt(2) * (sin(x) - cos(x))
         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
         *                      = -1/sqrt(2) * (cos(x) + sin(x))
         * To avoid cancellation, use
         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
         * to compute the worse one.
         */
        if (ix > 0x48000000)
            z = (invsqrtpi * ss) / RtIEEEsqrtf(x);
        else
        {
            u = ponef(x);
            v = qonef(x);
            z = invsqrtpi * (u * ss + v * cc) / RtIEEEsqrtf(x);
        }
        RWRETURN(z);
    }
    if (ix <= 0x24800000)
    {                          /* x < 2**-54 */
        RWRETURN((-tpi / x));
    }
    z = x * x;
    u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * U0[4])));
    v = one + z * (V0[0] +
                   z * (V0[1] + z * (V0[2] + z * (V0[3] + z * V0[4]))));
    RWRETURN((x * (u / v) +
              tpi * (RtIEEEj1f(x) * RtIEEElogf(x) - one / x)));
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEynf
 * computes the
 * Bessel function of the second kind for the integer order
 * n for the positive integer value x expressed as a double.
 * \param n order
 * \param x positive integer value
 * \return the Bessel function value
 */
float
RtIEEEynf(int n, float x)
{
    int                 i, hx, ix, ib;
    int                 sign;
    float               a, b, temp;

    RWAPIFUNCTION(RWSTRING("RtIEEEynf"));

    GET_FLOAT_WORD(hx, x);
    ix = 0x7fffffff & hx;
    /* if Y(n,NaN) is NaN */
    if (ix > 0x7f800000)
        RWRETURN(x + x);
    if (ix == 0)
        RWRETURN(-one / zero);
    if (hx < 0)
        RWRETURN(zero / zero);
    sign = 1;
    if (n < 0)
    {
        n = -n;
        sign = 1 - ((n & 1) << 1);
    }
    if (n == 0)
        RWRETURN((RtIEEEy0f(x)));
    if (n == 1)
        RWRETURN((sign * RtIEEEy1f(x)));
    if (ix == 0x7f800000)
        RWRETURN(zero);

    a = RtIEEEy0f(x);
    b = RtIEEEy1f(x);
    /* quit if b is -inf */
    GET_FLOAT_WORD(ib, b);
    for (i = 1; i < n && ib != (int) 0xff800000; i++)
    {
        temp = b;
        b = ((float) (i + i) / x) * b - a;
        GET_FLOAT_WORD(ib, b);
        a = temp;
    }
    if (sign > 0)
        RWRETURN(b);
    else
        RWRETURN(-b);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEElog10f
 * computes the value of the logarithm of argument
 * x to base 10.
 * \param x input
 * \return the logarithm to base 10
 */
float
RtIEEElog10f(float x)
{
    float               y, z;
    int                 i, k, hx;

    RWAPIFUNCTION(RWSTRING("RtIEEElog10f"));

    GET_FLOAT_WORD(hx, x);

    k = 0;
    if (hx < 0x00800000)
    {                          /* x < 2**-126  */
        if ((hx & 0x7fffffff) == 0)
            RWRETURN(-two25 / zero); /* log(+-0)=-inf */
        if (hx < 0)
            RWRETURN((x - x) / zero); /* log(-#) = NaN */
        k -= 25;
        x *= two25;            /* subnormal number, scale up x */
        GET_FLOAT_WORD(hx, x);
    }
    if (hx >= 0x7f800000)
        RWRETURN(x + x);
    k += (hx >> 23) - 127;
    i = ((unsigned int) k & 0x80000000) >> 31;
    hx = (hx & 0x007fffff) | ((0x7f - i) << 23);
    y = (float) (k + i);
    SET_FLOAT_WORD(x, hx);
    z = y * log10_2lo + ivln10 * RtIEEElogf(x);
    RWRETURN(z + y * log10_2hi);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEElogf
 * computes the value of the natural logarithm of argument x .
 * \param x input
 * \return the natual logarithm.
 */
float
RtIEEElogf(float x)
{
    float               hfsq, f, s, z, R, w, t1, t2, dk;
    int                 k, ix, i, j;

    RWAPIFUNCTION(RWSTRING("RtIEEElogf"));

    GET_FLOAT_WORD(ix, x);

    k = 0;
    if (ix < 0x00800000)
    {                          /* x < 2**-126  */
        if ((ix & 0x7fffffff) == 0)
            RWRETURN(-two25 / zero); /* log(+-0)=-inf */
        if (ix < 0)
            RWRETURN((x - x) / zero); /* log(-#) = NaN */
        k -= 25;
        x *= two25;            /* subnormal number, scale up x */
        GET_FLOAT_WORD(ix, x);
    }
    if (ix >= 0x7f800000)
        RWRETURN(x + x);
    k += (ix >> 23) - 127;
    ix &= 0x007fffff;
    i = (ix + (0x95f64 << 3)) & 0x800000;
    SET_FLOAT_WORD(x, ix | (i ^ 0x3f800000)); /* normalize x or x/2 */
    k += (i >> 23);
    f = x - (float) 1.0;
    if ((0x007fffff & (15 + ix)) < 16)
    {                          /* |f| < 2**-20 */
        if (f == zero)
        {
            if (k == 0)
            {
                RWRETURN(zero);
            }
            else
            {
                dk = (float) k;
                RWRETURN(dk * ln2_hi + dk * ln2_lo);
            }
        }

        R = f * f * ((float) 0.5 - (float) 0.33333333333333333 * f);
        if (k == 0)
            RWRETURN(f - R);
        else
        {
            dk = (float) k;
            RWRETURN(dk * ln2_hi - ((R - dk * ln2_lo) - f));
        }
    }
    s = f / ((float) 2.0 + f);
    dk = (float) k;
    z = s * s;
    i = ix - (0x6147a << 3);
    w = z * z;
    j = (0x6b851 << 3) - ix;
    t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
    t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
    i |= j;
    R = t2 + t1;
    if (i > 0)
    {
        hfsq = (float) 0.5 *f * f;

        if (k == 0)
            RWRETURN(f - (hfsq - s * (hfsq + R)));
        else
            RWRETURN(dk * ln2_hi -
                     ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f));
    }
    else
    {
        if (k == 0)
            RWRETURN(f - s * (f - R));
        else
            RWRETURN(dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f));
    }
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEpowf
 * computes the value
 * of x
 * to the exponent y .
 * \param x input
 * \param y exponent
 * \return x^y
 */
float
RtIEEEpowf(float x, float y)
{
    float               z, ax, z_h, z_l, p_h, p_l;
    float               y1, t1, t2, r, s, t, u, v, w;
    int                 i, j, k, yisint, n;
    int                 hx, hy, ix, iy, is;

    RWAPIFUNCTION(RWSTRING("RtIEEEpowf"));

    GET_FLOAT_WORD(hx, x);
    GET_FLOAT_WORD(hy, y);
    ix = hx & 0x7fffffff;
    iy = hy & 0x7fffffff;

    /* y==zero: x**0 = 1 */
    if (iy == 0)
        RWRETURN(one);

    /* +-NaN return x+y */
    if (ix > 0x7f800000 || iy > 0x7f800000)
        RWRETURN(x + y);

    /* determine if y is an odd int when x < 0
     * yisint = 0 ... y is not an integer
     * yisint = 1 ... y is an odd int
     * yisint = 2 ... y is an even int
     */
    yisint = 0;
    if (hx < 0)
    {
        if (iy >= 0x4b800000)
            yisint = 2;        /* even integer y */
        else if (iy >= 0x3f800000)
        {
            k = (iy >> 23) - 0x7f; /* exponent */
            j = iy >> (23 - k);
            if ((j << (23 - k)) == iy)
                yisint = 2 - (j & 1);
        }
    }

    /* special value of y */
    if (iy == 0x7f800000)
    {                          /* y is +-inf */
        if (ix == 0x3f800000)
            RWRETURN(y - y);   /* inf**+-1 is NaN */
        else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
            RWRETURN((hy >= 0) ? y : zero);
        else                   /* (|x|<1)**-,+inf = inf,0 */
            RWRETURN((hy < 0) ? -y : zero);
    }
    if (iy == 0x3f800000)
    {                          /* y is  +-1 */
        if (hy < 0)
            RWRETURN(one / x);
        else
            RWRETURN(x);
    }
    if (hy == 0x40000000)
        RWRETURN(x * x);       /* y is  2 */
    if (hy == 0x3f000000)
    {                          /* y is  0.5 */
        if (hx >= 0)           /* x >= +0 */
            RWRETURN(RtIEEEsqrtf(x));
    }

    ax = RtIEEEfabsf(x);
    /* special value of x */
    if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000)
    {
        z = ax;                /*x is +-0,+-inf,+-1 */
        if (hy < 0)
            z = one / z;       /* z = (1/|x|) */
        if (hx < 0)
        {
            if (((ix - 0x3f800000) | yisint) == 0)
            {
                z = (z - z) / (z - z); /* (-1)**non-int is NaN */
            }
            else if (yisint == 1)
                z = -z;        /* (x<0)**odd = -(|x|**odd) */
        }
        RWRETURN(z);
    }

    /* (x<0)**(non-int) is NaN */
    if (((((unsigned int) hx >> 31) - 1) | yisint) == 0)
        RWRETURN((x - x) / (x - x));

    /* |y| is rhuge */
    if (iy > 0x4d000000)
    {                          /* if |y| > 2**27 */
        /* over/underflow if x is not close to one */
        if (ix < 0x3f7ffff8)
            RWRETURN((hy < 0) ? rhuge * rhuge : tiny * tiny);
        if (ix > 0x3f800007)
            RWRETURN((hy > 0) ? rhuge * rhuge : tiny * tiny);
        /* now |1-x| is tiny <= 2**-20, suffice to compute 
         * log(x) by x-x^2/2+x^3/3-x^4/4 */
        t = x - 1;             /* t has 20 trailing zeros */
        w = (t * t) * ((float) 0.5 -
                       t * ((float) 0.333333333333 - t * (float) 0.25));
        u = ivln2_h * t;       /* ivln2_h has 16 sig. bits */
        v = t * ivln2_l - w * ivln2;
        t1 = u + v;
        GET_FLOAT_WORD(is, t1);
        SET_FLOAT_WORD(t1, is & 0xfffff000);
        t2 = v - (t1 - u);
    }
    else
    {
        float               s2, s_h, s_l, t_h, t_l;

        n = 0;
        /* take care subnormal number */
        if (ix < 0x00800000)
        {
            ax *= two24;
            n -= 24;
            GET_FLOAT_WORD(ix, ax);
        }
        n += ((ix) >> 23) - 0x7f;
        j = ix & 0x007fffff;
        /* determine interval */
        ix = j | 0x3f800000;   /* normalize ix */
        if (j <= 0x1cc471)
            k = 0;             /* |x|<sqrt(3/2) */
        else if (j < 0x5db3d7)
            k = 1;             /* |x|<sqrt(3)   */
        else
        {
            k = 0;
            n += 1;
            ix -= 0x00800000;
        }
        SET_FLOAT_WORD(ax, ix);

        /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
        u = ax - bp[k];        /* bp[0]=1.0, bp[1]=1.5 */
        v = one / (ax + bp[k]);
        s = u * v;
        s_h = s;
        GET_FLOAT_WORD(is, s_h);
        SET_FLOAT_WORD(s_h, is & 0xfffff000);
        /* t_h=ax+bp[k] High */
        SET_FLOAT_WORD(t_h,
                       ((ix >> 1) | 0x20000000) + 0x0040000 +
                       (k << 21));
        t_l = ax - (t_h - bp[k]);
        s_l = v * ((u - s_h * t_h) - s_h * t_l);
        /* compute log(ax) */
        s2 = s * s;
        r = s2 * s2 * (L1 +
                       s2 * (L2 +
                             s2 * (L3 +
                                   s2 * (L4 + s2 * (L5 + s2 * L6)))));
        r += s_l * (s_h + s);
        s2 = s_h * s_h;
        t_h = (float) 3.0 + s2 + r;
        GET_FLOAT_WORD(is, t_h);
        SET_FLOAT_WORD(t_h, is & 0xfffff000);
        t_l = r - ((t_h - (float) 3.0) - s2);
        /* u+v = s*(1+...) */
        u = s_h * t_h;
        v = s_l * t_h + t_l * s;
        /* 2/(3log2)*(s+...) */
        p_h = u + v;
        GET_FLOAT_WORD(is, p_h);
        SET_FLOAT_WORD(p_h, is & 0xfffff000);
        p_l = v - (p_h - u);
        z_h = cp_h * p_h;      /* cp_h+cp_l = 2/(3*log2) */
        z_l = cp_l * p_h + p_l * cp + dp_l[k];
        /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
        t = (float) n;
        t1 = (((z_h + z_l) + dp_h[k]) + t);
        GET_FLOAT_WORD(is, t1);
        SET_FLOAT_WORD(t1, is & 0xfffff000);
        t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
    }

    s = one;                   /* s (sign of result -ve**odd) = -1 else = 1 */
    if (((((unsigned int) hx >> 31) - 1) | (yisint - 1)) == 0)
        s = -one;              /* (-ve)**(odd int) */

    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
    GET_FLOAT_WORD(is, y);
    SET_FLOAT_WORD(y1, is & 0xfffff000);
    p_l = (y - y1) * t1 + y * t2;
    p_h = y1 * t1;
    z = p_l + p_h;
    GET_FLOAT_WORD(j, z);
    if (j > 0x43000000)        /* if z > 128 */
    {
        RWRETURN(s * rhuge * rhuge); /* overflow */
    }
    else if (j == 0x43000000)
    {                          /* if z == 128 */
        if (p_l + ovt > z - p_h)
        {
            RWRETURN(s * rhuge * rhuge); /* overflow */
        }

    }
    else if ((j & 0x7fffffff) > 0x43160000) /* z <= -150 */
    {
        RWRETURN(s * tiny * tiny); /* underflow */
    }
    else if (j == (int) 0xc3160000)
    {                          /* z == -150 */
        if (p_l <= z - p_h)
            RWRETURN(s * tiny * tiny); /* underflow */
    }
    /*
     * compute 2**(p_h+p_l)
     */
    i = j & 0x7fffffff;
    k = (i >> 23) - 0x7f;
    n = 0;
    if (i > 0x3f000000)
    {                          /* if |z| > 0.5, set n = [z+0.5] */
        n = j + (0x00800000 >> (k + 1));
        k = ((n & 0x7fffffff) >> 23) - 0x7f; /* new k for n */
        SET_FLOAT_WORD(t, n & ~(0x007fffff >> k));
        n = ((n & 0x007fffff) | 0x00800000) >> (23 - k);
        if (j < 0)
            n = -n;
        p_h -= t;
    }
    t = p_l + p_h;
    GET_FLOAT_WORD(is, t);
    SET_FLOAT_WORD(t, is & 0xfffff000);
    u = t * lg2_h;
    v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
    z = u + v;
    w = v - (z - u);
    t = z * z;
    t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
    r = (z * t1) / (t1 - two) - (w + z * w);
    z = one - (r - z);
    GET_FLOAT_WORD(j, z);
    j += (n << 23);
    if ((j >> 23) <= 0)
        z = RtIEEEscalbnf(z, n); /* subnormal output */
    else
        SET_FLOAT_WORD(z, j);
    RWRETURN(s * z);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEremainderf
 * is recommended by the IEEE 754 standard.
 * This returns the fractional remainder part of
 * the quotienty formed buy the input values
 * \param numerator
 * \param denominator
 * \return remainder
 */
float
RtIEEEremainderf(float x, float p)
{
    int                 hx, hp;
    unsigned int        sx;
    float               p_half;

    RWAPIFUNCTION(RWSTRING("RtIEEEremainderf"));

    GET_FLOAT_WORD(hx, x);
    GET_FLOAT_WORD(hp, p);
    sx = hx & 0x80000000;
    hp &= 0x7fffffff;
    hx &= 0x7fffffff;

    /* purge off exception values */
    if (hp == 0)
        RWRETURN((x * p) / (x * p)); /* p = 0 */
    if ((hx >= 0x7f800000) ||  /* x not finite */
        ((hp > 0x7f800000)))   /* p is NaN */
        RWRETURN((x * p) / (x * p));

    if (hp <= 0x7effffff)
        x = RtIEEEfmodf(x, p + p); /* now x < 2p */
    if ((hx - hp) == 0)
        RWRETURN(zero * x);
    x = RtIEEEfabsf(x);
    p = RtIEEEfabsf(p);
    if (hp < 0x01000000)
    {
        if (x + x > p)
        {
            x -= p;
            if (x + x >= p)
                x -= p;
        }
    }
    else
    {
        p_half = (float) 0.5 *p;

        if (x > p_half)
        {
            x -= p;
            if (x >= p_half)
                x -= p;
        }
    }
    GET_FLOAT_WORD(hx, x);
    SET_FLOAT_WORD(x, hx ^ sx);
    RWRETURN(x);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEscalbf
 * increments the exponent of the first by the integer part of second argment
 * \param x input
 * \param fn exponent delta
 * \return value modified float
 */
float
RtIEEEscalbf(float x, float fn)
{
    RWAPIFUNCTION(RWSTRING("RtIEEEscalbf"));

    if (RtIEEEisnanf(x) || RtIEEEisnanf(fn))
        RWRETURN(x * fn);
    if (!RtIEEEfinitef(fn))
    {
        if (fn > (float) 0.0)
            RWRETURN(x * fn);
        else
            RWRETURN(x / (-fn));
    }
    if (RtIEEErintf(fn) != fn)
        RWRETURN((fn - fn) / (fn - fn));
    if (fn > (float) 65000.0)
        RWRETURN(RtIEEEscalbnf(x, 65000));
    if (-fn > (float) 65000.0)
        RWRETURN(RtIEEEscalbnf(x, -65000));
    RWRETURN(RtIEEEscalbnf(x, (int) fn));
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEsinhf
 * computes the hyperbolic sine of x .
 * \param x input
 * \return hyperbolic sine
 */
float
RtIEEEsinhf(float x)
{
    float               t, w, h;
    int                 ix, jx;

    RWAPIFUNCTION(RWSTRING("RtIEEEsinhf"));

    GET_FLOAT_WORD(jx, x);
    ix = jx & 0x7fffffff;

    /* x is INF or NaN */
    if (ix >= 0x7f800000)
        RWRETURN(x + x);

    h = 0.5;
    if (jx < 0)
        h = -h;
    /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
    if (ix < 0x41b00000)
    {                          /* |x|<22 */
        if (ix < 0x31800000)   /* |x|<2**-28 */
            if (shuge + x > one)
                RWRETURN(x);   /* sinh(tiny) = tiny with inexact */
        t = RtIEEEexpm1f(RtIEEEfabsf(x));
        if (ix < 0x3f800000)
            RWRETURN(h * ((float) 2.0 * t - t * t / (t + one)));
        RWRETURN(h * (t + t / (t + one)));
    }

    /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
    if (ix < 0x42b17180)
        RWRETURN(h * RtIEEEexpf(RtIEEEfabsf(x)));

    /* |x| in [log(maxdouble), overflowthresold] */
    if (ix <= 0x42b2d4fc)
    {
        w = RtIEEEexpf((float) 0.5 * RtIEEEfabsf(x));
        t = h * w;
        RWRETURN(t * w);
    }

    /* |x| > overflowthresold, sinh(x) overflow */
    RWRETURN(x * shuge);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEsqrtf
 * computes
 * the non-negative square root of x.
 * \param x input
 * \return the non-negative square root
 */
float
RtIEEEsqrtf(float x)
{
    float               z;
    int                 sign = (int) 0x80000000;
    int                 ix, s, q, m, t, i;
    unsigned int        r;

    RWAPIFUNCTION(RWSTRING("RtIEEEsqrtf"));

    GET_FLOAT_WORD(ix, x);

    /* take care of Inf and NaN */
    if ((ix & 0x7f800000) == 0x7f800000)
    {
        RWRETURN(x * x + x);   /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
                                * * sqrt(-inf)=sNaN */
    }
    /* take care of zero */
    if (ix <= 0)
    {
        if ((ix & (~sign)) == 0)
            RWRETURN(x);       /* sqrt(+-0) = +-0 */
        else if (ix < 0)
            RWRETURN((x - x) / (x - x)); /* sqrt(-ve) = sNaN */
    }
    /* normalize x */
    m = (ix >> 23);
    if (m == 0)
    {                          /* subnormal x */
        for (i = 0; (ix & 0x00800000) == 0; i++)
            ix <<= 1;
        m -= i - 1;
    }
    m -= 127;                  /* unbias exponent */
    ix = (ix & 0x007fffff) | 0x00800000;
    if (m & 1)                 /* odd m, double x to make it even */
        ix += ix;
    m >>= 1;                   /* m = [m/2] */

    /* generate sqrt(x) bit by bit */
    ix += ix;
    q = s = 0;                 /* q = sqrt(x) */
    r = 0x01000000;            /* r = moving bit from right to left */

    while (r != 0)
    {
        t = s + r;
        if (t <= ix)
        {
            s = t + r;
            ix -= t;
            q += r;
        }
        ix += ix;
        r >>= 1;
    }

    /* use floating add to find out rounding direction */
    if (ix != 0)
    {
        z = one - tiny;        /* trigger inexact flag */
        if (z >= one)
        {
            z = one + tiny;
            if (z > one)
                q += 2;
            else
                q += (q & 1);
        }
    }
    ix = (q >> 1) + 0x3f000000;
    ix += (m << 23);
    SET_FLOAT_WORD(z, ix);
    RWRETURN(z);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEasinhf
 * computes the inverse hyperbolic sine
 * of the real
 * argument
 * \param x input
 * \return the inverse hyperbolic sine
 */
float
RtIEEEasinhf(float x)
{
    float               t, w;
    int                 hx, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEasinhf"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix >= 0x7f800000)
        RWRETURN(x + x);       /* x is inf or NaN */
    if (ix < 0x31800000)
    {                          /* |x|<2**-28 */
        if (rhuge + x > one)
            RWRETURN(x);       /* return x inexact except 0 */
    }
    if (ix > 0x4d800000)
    {                          /* |x| > 2**28 */
        w = RtIEEElogf(RtIEEEfabsf(x)) + ln2;
    }
    else if (ix > 0x40000000)
    {                          /* 2**28 > |x| > 2.0 */
        t = RtIEEEfabsf(x);
        w = RtIEEElogf((float) 2.0 * t +
                       one / (RtIEEEsqrtf(x * x + one) + t));
    }
    else
    {                          /* 2.0 > |x| > 2**-28 */
        t = x * x;
        w = RtIEEElog1pf(RtIEEEfabsf(x) +
                         t / (one + RtIEEEsqrtf(one + t)));
    }
    if (hx > 0)
        RWRETURN(w);
    else
        RWRETURN(-w);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEatanf
 * compute the principal value of the arc tangent of x
 * \param x input
 * \return the principal value of the arc tangent 
 */
float
RtIEEEatanf(float x)
{
    float               w, s1, s2, z;
    int                 ix, hx, id;

    RWAPIFUNCTION(RWSTRING("RtIEEEatanf"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix >= 0x50800000)
    {                          /* if |x| >= 2^34 */
        if (ix > 0x7f800000)
            RWRETURN(x + x);      /* NaN */
        if (hx > 0)
            RWRETURN(atanhi[3] + atanlo[3]);
        else
            RWRETURN(-atanhi[3] - atanlo[3]);
    }
    if (ix < 0x3ee00000)
    {                          /* |x| < 0.4375 */
        if (ix < 0x31000000)
        {                      /* |x| < 2^-29 */
            if (rhuge + x > one)
                RWRETURN(x);   /* raise inexact */
        }
        id = -1;
    }
    else
    {
        x = RtIEEEfabsf(x);
        if (ix < 0x3f980000)
        {                      /* |x| < 1.1875 */
            if (ix < 0x3f300000)
            {                  /* 7/16 <=|x|<11/16 */
                id = 0;
                x = ((float) 2.0 * x - one) / ((float) 2.0 + x);
            }
            else
            {                  /* 11/16<=|x|< 19/16 */
                id = 1;
                x = (x - one) / (x + one);
            }
        }
        else
        {
            if (ix < 0x401c0000)
            {                  /* |x| < 2.4375 */
                id = 2;
                x = (x - (float) 1.5) / (one + (float) 1.5 * x);
            }
            else
            {                  /* 2.4375 <= |x| < 2^66 */
                id = 3;
                x = -(float) 1.0 / x;
            }
        }
    }
    /* end of argument reduction */
    z = x * x;
    w = z * z;
    /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
    s1 = z * (aT[0] +
              w * (aT[2] +
                   w * (aT[4] +
                        w * (aT[6] + w * (aT[8] + w * aT[10])))));
    s2 = w * (aT[1] +
              w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
    if (id < 0)
        RWRETURN(x - x * (s1 + s2));
    else
    {
        z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x);
        RWRETURN((hx < 0) ? -z : z);
    }
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEcbrtf
 * computes
 * the cube root of x .
 * \param x input
 * \return the cube root
 */
float
RtIEEEcbrtf(float x)
{
    float               r, s, t;
    int                 hx;
    unsigned int        sign;
    unsigned int        high;

    RWAPIFUNCTION(RWSTRING("RtIEEEcbrtf"));

    GET_FLOAT_WORD(hx, x);
    sign = hx & 0x80000000;    /* sign= sign(x) */
    hx ^= sign;
    if (hx >= 0x7f800000)
        RWRETURN((x + x));     /* cbrt(NaN,INF) is itself */
    if (hx == 0)
        RWRETURN((x));         /* cbrt(0) is itself */

    SET_FLOAT_WORD(x, hx);     /* x <- |x| */
    /* rough cbrt to 5 bits */
    if (hx < 0x00800000)       /* subnormal number */
    {
        SET_FLOAT_WORD(t, 0x4b800000); /* set t= 2**24 */
        t *= x;
        GET_FLOAT_WORD(high, t);
        SET_FLOAT_WORD(t, high / 3 + B2);
    }
    else
        SET_FLOAT_WORD(t, hx / 3 + B1);

    /* new cbrt to 23 bits */
    r = t * t / x;
    s = C + r * t;
    t *= G + F / (s + E + D / s);

    /* retore the sign bit */
    GET_FLOAT_WORD(high, t);
    SET_FLOAT_WORD(t, high | sign);
    RWRETURN((t));
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEceilf
 * returns the smallest integral value
 * greater than or equal to x .
 * \param x input
 * \return ceiling
 */
float
RtIEEEceilf(float x)
{
    int                 i0, j0;
    unsigned int        i;

    RWAPIFUNCTION(RWSTRING("RtIEEEceilf"));

    GET_FLOAT_WORD(i0, x);
    j0 = ((i0 >> 23) & 0xff) - 0x7f;
    if (j0 < 23)
    {
        if (j0 < 0)
        {                      /* raise inexact if x != 0 */
            if (rhuge + x > (float) 0.0)
            {                  /* return 0*sign(x) if |x|<1 */
                if (i0 < 0)
                {
                    i0 = 0x80000000;
                }
                else if (i0 != 0)
                {
                    i0 = 0x3f800000;
                }
            }
        }
        else
        {
            i = (0x007fffff) >> j0;
            if ((i0 & i) == 0)
                RWRETURN(x);   /* x is integral */
            if (rhuge + x > (float) 0.0)
            {                  /* raise inexact flag */
                if (i0 > 0)
                    i0 += (0x00800000) >> j0;
                i0 &= (~i);
            }
        }
    }
    else
    {
        if (j0 == 0x80)
            RWRETURN(x + x);   /* inf or NaN */
        else
            RWRETURN(x);       /* x is integral */
    }
    SET_FLOAT_WORD(x, i0);
    RWRETURN(x);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEcopysignf
 * is recommended by the IEEE 754 standard, and
 * returns x
 * with its sign changed to that of y 
 * \param x input
 * \param y target sign
 * \return sign-modified value
 */
float
RtIEEEcopysignf(float x, float y)
{
    unsigned int        ix, iy;

    RWAPIFUNCTION(RWSTRING("RtIEEEcopysignf"));

    GET_FLOAT_WORD(ix, x);
    GET_FLOAT_WORD(iy, y);
    SET_FLOAT_WORD(x, (ix & 0x7fffffff) | (iy & 0x80000000));
    RWRETURN(x);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEcosf
 * computes the cosine of x (measured in radians).
 * \param x input in radians
 * \return the cosine 
 * \warning A large magnitude argument may yield a result with little or no
 * significance.
 */
float
RtIEEEcosf(float x)
{
    float               y[2], z = 0.0;
    int                 n, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEcosf"));

    GET_FLOAT_WORD(ix, x);

    /* |x| ~< pi/4 */
    ix &= 0x7fffffff;
    if (ix <= 0x3f490fd8)
        RWRETURN(_rtIEEEcosf(x, z));

    /* cos(Inf or NaN) is NaN */
    else if (ix >= 0x7f800000)
        RWRETURN(x - x);

    /* argument reduction needed */
    else
    {
        n = RtIEEErem_pio2f(x, y);
        switch (n & 3)
        {
            case 0:
                RWRETURN(_rtIEEEcosf(y[0], y[1]));
            case 1:
                RWRETURN(-_rtIEEEsinf(y[0], y[1], 1));
            case 2:
                RWRETURN(-_rtIEEEcosf(y[0], y[1]));
            default:
                RWRETURN(_rtIEEEsinf(y[0], y[1], 1));
        }
    }
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEerff
 * calculates the error function of x .
 * \param x input
 * \return the error function.
 */
float
RtIEEEerff(float x)
{
    int                 hx, ix, i;
    float               R, S, P, Q, s, y, z, r;

    RWAPIFUNCTION(RWSTRING("RtIEEEerff"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix >= 0x7f800000)
    {                          /* erf(nan)=nan */
        i = ((unsigned int) hx >> 31) << 1;
        RWRETURN((float) (1 - i) + one / x); /* erf(+-inf)=+-1 */
    }

    if (ix < 0x3f580000)
    {                          /* |x|<0.84375 */
        if (ix < 0x31800000)
        {                      /* |x|<2**-28 */
            if (ix < 0x04000000)
                /*avoid underflow */
                RWRETURN((float) 0.125 * ((float) 8.0 * x + efx8 * x));

            RWRETURN(x + efx * x);
        }
        z = x * x;
        r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
        s = one + z * (qq1 +
                       z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
        y = r / s;
        RWRETURN(x + x * y);
    }
    if (ix < 0x3fa00000)
    {                          /* 0.84375 <= |x| < 1.25 */
        s = RtIEEEfabsf(x) - one;
        P = pa0 + s * (pa1 +
                       s * (pa2 +
                            s * (pa3 +
                                 s * (pa4 + s * (pa5 + s * pa6)))));
        Q = one + s * (qa1 +
                       s * (qa2 +
                            s * (qa3 +
                                 s * (qa4 + s * (qa5 + s * qa6)))));
        if (hx >= 0)
            RWRETURN(erx + P / Q);
        else
            RWRETURN(-erx - P / Q);
    }
    if (ix >= 0x40c00000)
    {                          /* inf>|x|>=6 */
        if (hx >= 0)
            RWRETURN(one - tiny);
        else
            RWRETURN(tiny - one);
    }
    x = RtIEEEfabsf(x);
    s = one / (x * x);
    if (ix < 0x4036DB6E)
    {                          /* |x| < 1/0.35 */
        R = ra0 + s * (ra1 +
                       s * (ra2 +
                            s * (ra3 +
                                 s * (ra4 +
                                      s * (ra5 +
                                           s * (ra6 + s * ra7))))));
        S = one + s * (sa1 +
                       s * (sa2 +
                            s * (sa3 +
                                 s * (sa4 +
                                      s * (sa5 +
                                           s * (sa6 +
                                                s * (sa7 +
                                                     s * sa8)))))));
    }
    else
    {                          /* |x| >= 1/0.35 */
        R = rb0 + s * (rb1 +
                       s * (rb2 +
                            s * (rb3 +
                                 s * (rb4 + s * (rb5 + s * rb6)))));
        S = one + s * (sb1 +
                       s * (sb2 +
                            s * (sb3 +
                                 s * (sb4 +
                                      s * (sb5 +
                                           s * (sb6 + s * sb7))))));
    }
    GET_FLOAT_WORD(ix, x);
    SET_FLOAT_WORD(z, ix & 0xfffff000);
    r = RtIEEEexpf(-z * z -
                   (float) 0.5625) * RtIEEEexpf((z - x) * (z +
                                                           x) + R / S);
    if (hx >= 0)
        RWRETURN(one - r / x);
    else
        RWRETURN(r / x - one);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEerfcf
 *  calculates the complementary error function of x 
 * \param x input
 * \return the complementary error function.
 */
float
RtIEEEerfcf(float x)
{
    int                 hx, ix;
    float               R, S, P, Q, s, y, z, r;

    RWAPIFUNCTION(RWSTRING("RtIEEEerfcf"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix >= 0x7f800000)
    {                          /* erfc(nan)=nan */
        /* erfc(+-inf)=0,2 */
        RWRETURN((float) (((unsigned int) hx >> 31) << 1) + one / x);
    }

    if (ix < 0x3f580000)
    {                          /* |x|<0.84375 */
        if (ix < 0x23800000)   /* |x|<2**-56 */
            RWRETURN(one - x);
        z = x * x;
        r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
        s = one + z * (qq1 +
                       z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
        y = r / s;
        if (hx < 0x3e800000)
        {                      /* x<1/4 */
            RWRETURN(one - (x + x * y));
        }
        else
        {
            r = x * y;
            r += (x - half);
            RWRETURN(half - r);
        }
    }
    if (ix < 0x3fa00000)
    {                          /* 0.84375 <= |x| < 1.25 */
        s = RtIEEEfabsf(x) - one;
        P = pa0 + s * (pa1 +
                       s * (pa2 +
                            s * (pa3 +
                                 s * (pa4 + s * (pa5 + s * pa6)))));
        Q = one + s * (qa1 +
                       s * (qa2 +
                            s * (qa3 +
                                 s * (qa4 + s * (qa5 + s * qa6)))));
        if (hx >= 0)
        {
            z = one - erx;
            RWRETURN(z - P / Q);
        }
        else
        {
            z = erx + P / Q;
            RWRETURN(one + z);
        }
    }
    if (ix < 0x41e00000)
    {                          /* |x|<28 */
        x = RtIEEEfabsf(x);
        s = one / (x * x);
        if (ix < 0x4036DB6D)
        {                      /* |x| < 1/.35 ~ 2.857143 */
            R = ra0 + s * (ra1 +
                           s * (ra2 +
                                s * (ra3 +
                                     s * (ra4 +
                                          s * (ra5 +
                                               s * (ra6 + s * ra7))))));
            S = one + s * (sa1 +
                           s * (sa2 +
                                s * (sa3 +
                                     s * (sa4 +
                                          s * (sa5 +
                                               s * (sa6 +
                                                    s * (sa7 +
                                                         s * sa8)))))));
        }
        else
        {                      /* |x| >= 1/.35 ~ 2.857143 */
            if (hx < 0 && ix >= 0x40c00000)
                RWRETURN(two - tiny); /* x < -6 */
            R = rb0 + s * (rb1 +
                           s * (rb2 +
                                s * (rb3 +
                                     s * (rb4 + s * (rb5 + s * rb6)))));
            S = one + s * (sb1 +
                           s * (sb2 +
                                s * (sb3 +
                                     s * (sb4 +
                                          s * (sb5 +
                                               s * (sb6 + s * sb7))))));
        }
        GET_FLOAT_WORD(ix, x);
        SET_FLOAT_WORD(z, ix & 0xfffff000);
        r = RtIEEEexpf(-z * z - (float) 0.5625) *
            RtIEEEexpf((z - x) * (z + x) + R / S);
        if (hx > 0)
            RWRETURN(r / x);
        else
            RWRETURN(two - r / x);
    }
    else
    {
        if (hx > 0)
            RWRETURN(tiny * tiny);
        else
            RWRETURN(two - tiny);
    }
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEfabsf
 * computes the absolute value of a floating-point number x .
 * \param x input
 * \return the absolute value
 */
float
RtIEEEfabsf(float x)
{
    unsigned int        ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEfabsf"));

    GET_FLOAT_WORD(ix, x);
    SET_FLOAT_WORD(x, ix & 0x7fffffff);
    RWRETURN(x);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEfinitef
 * is recommended by the IEEE 754 standard, and 
 * returns the value 1 when the input value is finite, or 0 otherwise
 * \param x input
 * \return 1 when the input value is finite, or 0 otherwise
 */
int
RtIEEEfinitef(float x)
{
    int                 ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEfinitef"));

    GET_FLOAT_WORD(ix, x);
    RWRETURN((int) ((unsigned int) ((ix & 0x7fffffff) -
                                    0x7f800000) >> 31));
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEfloorf
 * returns the largest integral value less than or equal to x .
 * \param x input
 * \return floor
 */
float
RtIEEEfloorf(float x)
{
    int                 i0, j0;
    unsigned int        i;

    RWAPIFUNCTION(RWSTRING("RtIEEEfloorf"));

    GET_FLOAT_WORD(i0, x);
    j0 = ((i0 >> 23) & 0xff) - 0x7f;
    if (j0 < 23)
    {
        if (j0 < 0)
        {                      /* raise inexact if x != 0 */
            if (rhuge + x > (float) 0.0)
            {                  /* return 0*sign(x) if |x|<1 */
                if (i0 >= 0)
                {
                    i0 = 0;
                }
                else if ((i0 & 0x7fffffff) != 0)
                {
                    i0 = 0xbf800000;
                }
            }
        }
        else
        {
            i = (0x007fffff) >> j0;
            if ((i0 & i) == 0)
                RWRETURN(x);   /* x is integral */
            if (rhuge + x > (float) 0.0)
            {                  /* raise inexact flag */
                if (i0 < 0)
                    i0 += (0x00800000) >> j0;
                i0 &= (~i);
            }
        }
    }
    else
    {
        if (j0 == 0x80)
            RWRETURN(x + x);   /* inf or NaN */
        else
            RWRETURN(x);       /* x is integral */
    }
    SET_FLOAT_WORD(x, i0);
    RWRETURN(x);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEilogbf
 * is recommended by the IEEE 754 standard.
 * returns x's exponent in integer format.
 * \param x inpout value
 * \return exponent in integer format.
 */
int
RtIEEEilogbf(float x)
{
    int                 hx, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEilogbf"));

    GET_FLOAT_WORD(hx, x);
    hx &= 0x7fffffff;
    if (hx < 0x00800000)
    {
        if (hx == 0)
            RWRETURN(0x80000001); /* ilogb(0) = 0x80000001 */
        else                   /* subnormal x */
            for (ix = -126, hx <<= 8; hx > 0; hx <<= 1)
                ix -= 1;
        RWRETURN(ix);
    }
    else if (hx < 0x7f800000)
        RWRETURN((hx >> 23) - 127);
    else
        RWRETURN(0x7fffffff);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEElog1pf
 * computes
 * the value of log(1+x) accurately even for tiny argument x .
 * \param x input
 * \return the value of log(1+x)
 */
float
RtIEEElog1pf(float x)
{
    float               hfsq, f = 0, c = 0, s, z, R, u;
    int                 k, hx, hu = 0, ax;

    RWAPIFUNCTION(RWSTRING("RtIEEElog1pf"));

    GET_FLOAT_WORD(hx, x);
    ax = hx & 0x7fffffff;

    k = 1;
    if (hx < 0x3ed413d7)
    {                          /* x < 0.41422  */
        if (ax >= 0x3f800000)
        {                      /* x <= -1.0 */
            if (x == (float) -1.0)
                RWRETURN(-two25 / zero); /* log1p(-1)=+inf */
            else
                RWRETURN((x - x) / (x - x)); /* log1p(x<-1)=NaN */
        }
        if (ax < 0x31000000)
        {                      /* |x| < 2**-29 */
            if (two25 + x > zero /* raise inexact */
                && ax < 0x24800000) /* |x| < 2**-54 */
                RWRETURN(x);
            else
                RWRETURN(x - x * x * (float) 0.5);
        }
        if (hx > 0 || hx <= ((int) 0xbe95f61f))
        {
            k = 0;
            f = x;
            hu = 1;
        }                      /* -0.2929<x<0.41422 */
    }
    if (hx >= 0x7f800000)
        RWRETURN(x + x);
    if (k != 0)
    {
        if (hx < 0x5a000000)
        {
            u = (float) 1.0 + x;
            GET_FLOAT_WORD(hu, u);
            k = (hu >> 23) - 127;
            /* correction term */
            c = (k > 0) ? (float) 1.0 - (u - x) : x - (u - (float) 1.0);
            c /= u;
        }
        else
        {
            u = x;
            GET_FLOAT_WORD(hu, u);
            k = (hu >> 23) - 127;
            c = 0;
        }
        hu &= 0x007fffff;
        if (hu < 0x3504f7)
        {
            SET_FLOAT_WORD(u, hu | 0x3f800000); /* normalize u */
        }
        else
        {
            k += 1;
            SET_FLOAT_WORD(u, hu | 0x3f000000); /* normalize u/2 */
            hu = (0x00800000 - hu) >> 2;
        }
        f = u - (float) 1.0;
    }
    hfsq = (float) 0.5 *f * f;

    if (hu == 0)
    {                          /* |f| < 2**-20 */
        if (f == zero)
        {
            if (k == 0)
            {
                RWRETURN(zero);
            }
            else
            {
                c += k * ln2_lo;
                RWRETURN(k * ln2_hi + c);
            }
        }

        R = hfsq * ((float) 1.0 - (float) 0.66666666666666666 * f);
        if (k == 0)
            RWRETURN(f - R);
        else
            RWRETURN(k * ln2_hi - ((R - (k * ln2_lo + c)) - f));
    }
    s = f / ((float) 2.0 + f);
    z = s * s;
    R = z * (Lp1 +
             z * (Lp2 +
                  z * (Lp3 +
                       z * (Lp4 + z * (Lp5 + z * (Lp6 + z * Lp7))))));
    if (k == 0)
        RWRETURN(f - (hfsq - s * (hfsq + R)));
    else
        RWRETURN(k * ln2_hi -
                 ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f));
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEElogbf
 * returns x's exponent as a real
 * \param x input
 * \return exponent as a real
 */
float
RtIEEElogbf(float x)
{
    int                 ix;

    RWAPIFUNCTION(RWSTRING("RtIEEElogbf"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;          /* high |x| */
    if (ix == 0)
        RWRETURN((float) -1.0 / RtIEEEfabsf(x));
    if (ix >= 0x7f800000)
        RWRETURN(x * x);
    if ((ix >>= 23) == 0)      /* IEEE 754 logb */
        RWRETURN(-126.0);
    else
        RWRETURN((float) (ix - 127));
}

/**
 * \ingroup rtieeef
 * \ref RtIEEErintf
 *  returns the integral value (represented as a real)
 * nearest to x according to the prevailing rounding mode.
 * \param x input 
 * \return rounded value
 */
float
RtIEEErintf(float x)
{
    int                 i0, j0, sx;
    unsigned int        i, i1;
    float               w, t;

    RWAPIFUNCTION(RWSTRING("RtIEEErintf"));

    GET_FLOAT_WORD(i0, x);
    sx = (i0 >> 31) & 1;
    j0 = ((i0 >> 23) & 0xff) - 0x7f;
    if (j0 < 23)
    {
        if (j0 < 0)
        {
            if ((i0 & 0x7fffffff) == 0)
                RWRETURN(x);
            i1 = (i0 & 0x07fffff);
            i0 &= 0xfff00000;
            i0 |= ((i1 | -((signed int)i1)) >> 9) & 0x400000;
            SET_FLOAT_WORD(x, i0);
            w = TWO23[sx] + x;
            t = w - TWO23[sx];
            GET_FLOAT_WORD(i0, t);
            SET_FLOAT_WORD(t, (i0 & 0x7fffffff) | (sx << 31));
            RWRETURN(t);
        }
        else
        {
            i = (0x007fffff) >> j0;
            if ((i0 & i) == 0)
                RWRETURN(x);   /* x is integral */
            i >>= 1;
            if ((i0 & i) != 0)
                i0 = (i0 & (~i)) | ((0x100000) >> j0);
        }
    }
    else
    {
        if (j0 == 0x80)
            RWRETURN(x + x);   /* inf or NaN */
        else
            RWRETURN(x);       /* x is integral */
    }
    SET_FLOAT_WORD(x, i0);
    w = TWO23[sx] + x;
    RWRETURN(w - TWO23[sx]);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEscalbnf
 * increments the exponent of the first by the second argment
 * \param x input
 * \param n exponent delta
 * \return value modified float
 */
float
RtIEEEscalbnf(float x, int n)
{
    int                 k, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEscalbnf"));

    GET_FLOAT_WORD(ix, x);
    k = (ix & 0x7f800000) >> 23; /* extract exponent */
    if (k == 0)
    {                          /* 0 or subnormal x */
        if ((ix & 0x7fffffff) == 0)
            RWRETURN(x);       /* +-0 */
        x *= two25;
        GET_FLOAT_WORD(ix, x);
        k = ((ix & 0x7f800000) >> 23) - 25;
        if (n < -50000)
            RWRETURN(tiny * x); /*underflow */
    }
    if (k == 0xff)
        RWRETURN(x + x);       /* NaN or Inf */
    k = k + n;
    if (k > 0xfe)
        RWRETURN(rhuge * RtIEEEcopysignf(rhuge, x)); /* overflow  */
    if (k > 0)                 /* normal result */
    {
        SET_FLOAT_WORD(x, (ix & 0x807fffff) | (k << 23));
        RWRETURN(x);
    }
    if (k <= -25)
    {
        if (n > 50000)         /* in case integer overflow in n+k */
        {
            RWRETURN(rhuge * RtIEEEcopysignf(rhuge, x)); /*overflow */
        }
        else
        {
            RWRETURN(tiny * RtIEEEcopysignf(tiny, x)); /*underflow */
        }
    }

    k += 25;                   /* subnormal result */
    SET_FLOAT_WORD(x, (ix & 0x807fffff) | (k << 23));
    RWRETURN(x * twom25);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEsignificandf
 * returns the significand or mantissa of the input value
 * \param x input
 * \return the significand or mantissa
 */
float
RtIEEEsignificandf(float x)
{
    RWAPIFUNCTION(RWSTRING("RtIEEEsignificandf"));

    RWRETURN(RtIEEEscalbf(x, (float) -RtIEEEilogbf(x)));
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEsinf
 * computes the sine of x (measured in radians).
 * \param x input in radians
 * \return the sine 
 * \warning A large magnitude argument may yield a result with little or no
 * significance.
 */
float
RtIEEEsinf(float x)
{
    float               y[2], z = 0.0;
    int                 n, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEsinf"));
    GET_FLOAT_WORD(ix, x);

    /* |x| ~< pi/4 */
    ix &= 0x7fffffff;
    if (ix <= 0x3f490fd8)
        RWRETURN(_rtIEEEsinf(x, z, 0));

    /* sin(Inf or NaN) is NaN */
    else if (ix >= 0x7f800000)
        RWRETURN(x - x);

    /* argument reduction needed */
    else
    {
        n = RtIEEErem_pio2f(x, y);
        switch (n & 3)
        {
            case 0:
                RWRETURN(_rtIEEEsinf(y[0], y[1], 1));
            case 1:
                RWRETURN(_rtIEEEcosf(y[0], y[1]));
            case 2:
                RWRETURN(-_rtIEEEsinf(y[0], y[1], 1));
            default:
                RWRETURN(-_rtIEEEcosf(y[0], y[1]));
        }
    }
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEtanf
 * computes the tangent of x (measured in radians).
 * \param x input in radians
 * \return the tangent 
 * \warning A large magnitude argument may yield a result with little or no
 * significance.
 */
float
RtIEEEtanf(float x)
{
    float               y[2], z = 0.0;
    int                 n, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEtanf"));

    GET_FLOAT_WORD(ix, x);

    /* |x| ~< pi/4 */
    ix &= 0x7fffffff;
    if (ix <= 0x3f490fda)
        RWRETURN(_rtIEEEtanf(x, z, 1));

    /* tan(Inf or NaN) is NaN */
    else if (ix >= 0x7f800000)
        RWRETURN(x - x);       /* NaN */

    /* argument reduction needed */
    else
    {
        n = RtIEEErem_pio2f(x, y);
        RWRETURN(_rtIEEEtanf(y[0], y[1], 1 - ((n & 1) << 1))); /*   1 -- n even
                                                                * * * -1 -- n odd */
    }
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEtanhf
 * computes thehyperbolic tangent of x (measured in radians).
 * \param x input in radians
 * \return thehyperbolic tangent 
 * \warning A large magnitude argument may yield a result with little or no
 * significance.
 */
float
RtIEEEtanhf(float x)
{
    float               t, z;
    int                 jx, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEtanhf"));

    GET_FLOAT_WORD(jx, x);
    ix = jx & 0x7fffffff;

    /* x is INF or NaN */
    if (ix >= 0x7f800000)
    {
        if (jx >= 0)
            RWRETURN(one / x + one); /* tanh(+-inf)=+-1 */
        else
            RWRETURN(one / x - one); /* tanh(NaN) = NaN */
    }

    /* |x| < 22 */
    if (ix < 0x41b00000)
    {                          /* |x|<22 */
        if (ix < 0x24000000)   /* |x|<2**-55 */
            RWRETURN(x * (one + x)); /* tanh(small) = small */
        if (ix >= 0x3f800000)
        {                      /* |x|>=1  */
            t = RtIEEEexpm1f(two * RtIEEEfabsf(x));
            z = one - two / (t + two);
        }
        else
        {
            t = RtIEEEexpm1f(-two * RtIEEEfabsf(x));
            z = -t / (t + two);
        }
        /* |x| > 22, return +-1 */
    }
    else
    {
        z = one - tiny;        /* raised inexact flag */
    }
    RWRETURN((jx >= 0) ? z : -z);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEacosf
 * compute the principal value of the arc cosine of x
 * \param x input
 * \return the principal value of the arc cosine 
 */
float
RtIEEEacosf(float x)            /* wrapper acosf */
{
    float               z;

    RWAPIFUNCTION(RWSTRING("RtIEEEacosf"));

    z = _rtIEEEacosf(x);
    if ( /* _LIB_VERSION == _IEEE_ || */ RtIEEEisnanf(x))
        RWRETURN(z);
    if (RtIEEEfabsf(x) > (float) 1.0)
    {
        /* RtIEEEacosf(|x|>1) */
        RWRETURN((float) _rtIEEEstandard((double) x, (double) x, 101));
    }
    else
        RWRETURN(z);

}

/**
 * \ingroup rtieeef
 * \ref RtIEEEacoshf
 * compute the inverse hyperbolic cosine
 * of the real
 * argument x.
 * \param x input
 * \return the inverse hyperbolic cosine
 */
float
RtIEEEacoshf(float x)           /* wrapper acoshf */
{
    float               z;

    RWAPIFUNCTION(RWSTRING("RtIEEEacoshf"));

    z = _rtIEEEacoshf(x);
    if ( /* _LIB_VERSION == _IEEE_ || */ RtIEEEisnanf(x))
        RWRETURN(z);
    if (x < (float) 1.0)
    {
        /* acosh(x<1) */
        RWRETURN((float) _rtIEEEstandard((double) x, (double) x, 129));
    }
    else
        RWRETURN(z);

}

/**
 * \ingroup rtieeef
 * \ref RtIEEEasinf
 * compute the principal value of the arc sine of x
 * \param x input
 * \return the principal value of the arc sine 
 */
float
RtIEEEasinf(float x)            /* wrapper asinf */
{
    float               z;

    RWAPIFUNCTION(RWSTRING("RtIEEEasinf"));

    z = _rtIEEEasinf(x);
    if ( /* _LIB_VERSION == _IEEE_ || */ RtIEEEisnanf(x))
        RWRETURN(z);
    if (RtIEEEfabsf(x) > (float) 1.0)
    {
        /* RtIEEEasinf(|x|>1) */
        RWRETURN((float) _rtIEEEstandard((double) x, (double) x, 102));
    }
    else
        RWRETURN(z);

}

/**
 * \ingroup rtieeef
 * \ref RtIEEEatan2f
 * computes the principal value of the arc tangent of y/x ,
 * using the signs of both arguments to determine the quadrant of
 * the return value.
 * \param y exp(R)*sin(theta)
 * \param x exp(R)*cos(theta)
 * \return theta
 */
float
RtIEEEatan2f(float y, float x)  /* wrapper atan2f */
{
    float               z;

    RWAPIFUNCTION(RWSTRING("RtIEEEatan2f"));

    z = _rtIEEEatan2f(y, x);
    if ( /* _LIB_VERSION == _IEEE_ || */ RtIEEEisnanf(x)
        || RtIEEEisnanf(y))
        RWRETURN(z);
    if (x == (float) 0.0 && y == (float) 0.0)
    {
        /* atan2f(+-0,+-0) */
        RWRETURN((float) _rtIEEEstandard((double) y, (double) x, 103));
    }
    else
        RWRETURN(z);

}

/**
 * \ingroup rtieeef
 * \ref RtIEEEatanhf
 * computes the inverse hyperbolic tangent
 * of the real
 * argument x.
 * \param x input
 * \return the inverse hyperbolic tangent
 */
float
RtIEEEatanhf(float x)           /* wrapper atanhf */
{
    float               z, y;

    RWAPIFUNCTION(RWSTRING("RtIEEEatanhf"));

    z = _rtIEEEatanhf(x);
    if ( /* _LIB_VERSION == _IEEE_ || */ RtIEEEisnanf(x))
        RWRETURN(z);
    y = RtIEEEfabsf(x);
    if (y >= (float) 1.0)
    {
        if (y > (float) 1.0)
            /* atanhf(|x|>1) */
            RWRETURN((float)
                     _rtIEEEstandard((double) x, (double) x, 130));
        else
            /* atanhf(|x|==1) */
            RWRETURN((float)
                     _rtIEEEstandard((double) x, (double) x, 131));
    }
    else
        RWRETURN(z);

}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEexpf
 * computes the exponential value of the given argument x .
 * \param x input
 * \return the exponential value 
 */
float
RtIEEEexpf(float x)             /* default IEEE double exp */
{
    float               y, hi = 0, lo = 0, c, t;
    int                 k = 0, xsb;
    unsigned int        hx;

    RWAPIFUNCTION(RWSTRING("RtIEEEexpf"));

    GET_FLOAT_WORD(hx, x);
    xsb = (hx >> 31) & 1;      /* sign bit of x */
    hx &= 0x7fffffff;          /* high word of |x| */

    /* filter out non-finite argument */
    if (hx >= 0x42b17218)
    {                          /* if |x|>=88.721... */
        if (hx > 0x7f800000)
            RWRETURN(x + x);   /* NaN */
        if (hx == 0x7f800000)
            RWRETURN((float) ((xsb == 0) ? x : 0.0)); /* exp(+-inf)={inf,0} */
        if (x > o_threshold)
            RWRETURN(rhuge * rhuge); /* overflow */
        if (x < u_threshold)
            RWRETURN(twom100 * twom100); /* underflow */
    }

    /* argument reduction */
    if (hx > 0x3eb17218)
    {                          /* if  |x| > 0.5 ln2 */
        if (hx < 0x3F851592)
        {                      /* and |x| < 1.5 ln2 */
            hi = x - ln2HI[xsb];
            lo = ln2LO[xsb];
            k = 1 - xsb - xsb;
        }
        else
        {
            k = (int) (invln2 * x + halF[xsb]);
            t = (float) k;
            hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */
            lo = t * ln2LO[0];
        }
        x = hi - lo;
    }
    else if (hx < 0x31800000)
    {                          /* when |x|<2**-28 */
        if (rhuge + x > one)
            RWRETURN(one + x); /* trigger inexact */
    }
    else
        k = 0;

    /* x is now in primary range */
    t = x * x;
    c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
    if (k == 0)
        RWRETURN(one - ((x * c) / (c - (float) 2.0) - x));
    else
        y = one - ((lo - (x * c) / ((float) 2.0 - c)) - hi);
    if (k >= -125)
    {
        unsigned int        hy;

        GET_FLOAT_WORD(hy, y);
        SET_FLOAT_WORD(y, hy + (k << 23)); /* add k to y's exponent */
        RWRETURN(y);
    }
    else
    {
        unsigned int        hy;

        GET_FLOAT_WORD(hy, y);
        SET_FLOAT_WORD(y, hy + ((k + 100) << 23)); /* add k to y's exponent */
        RWRETURN(y * twom100);
    }
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEEjnf
 * computes the
 * Bessel function of the first kind of the integer order n
 * for the real value x .
 * \param n order
 * \param x input
 * \return Bessel function value.
 */
float
RtIEEEjnf(int n, float x)
{
    int                 i, hx, ix, sgn;
    float               a, b, temp, di;
    float               z, w;

    RWAPIFUNCTION(RWSTRING("RtIEEEjnf"));

    /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
     * Thus, J(-n,x) = J(n,-x)
     */
    GET_FLOAT_WORD(hx, x);
    ix = 0x7fffffff & hx;
    /* if J(n,NaN) is NaN */
    if (ix > 0x7f800000)
        RWRETURN(x + x);
    if (n < 0)
    {
        n = -n;
        x = -x;
        hx ^= 0x80000000;
    }
    if (n == 0)
        RWRETURN((RtIEEEj0f(x)));
    if (n == 1)
        RWRETURN((RtIEEEj1f(x)));
    sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */
    x = RtIEEEfabsf(x);
    if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */
        b = zero;
    else if ((float) n <= x)
    {
        /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
        a = RtIEEEj0f(x);
        b = RtIEEEj1f(x);
        for (i = 1; i < n; i++)
        {
            temp = b;
            b = b * ((float) (i + i) / x) - a; /* avoid underflow */
            a = temp;
        }
    }
    else
    {
        if (ix < 0x30800000)
        {                      /* x < 2**-29 */
            /* x is tiny, return the first Taylor expansion of J(n,x) 
             * J(n,x) = 1/n!*(x/2)^n  - ...
             */
            if (n > 33)        /* underflow */
                b = zero;
            else
            {
                temp = x * (float) 0.5;
                b = temp;
                for (a = one, i = 2; i <= n; i++)
                {
                    a *= (float) i; /* a = n! */
                    b *= temp; /* b = (x/2)^n */
                }
                b = b / a;
            }
        }
        else
        {
            /* use backward recurrence */
            /*    x      x^2      x^2       
             *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
             *   2n  - 2(n+1) - 2(n+2)
             *
             *    1      1        1       
             *  (for large x)   =  ----  ------   ------   .....
             *   2n   2(n+1)   2(n+2)
             *   -- - ------ - ------ - 
             *    x     x         x
             *
             * Let w = 2n/x and h=2/x, then the above quotient
             * is equal to the continued fraction:
             *      1
             * = -----------------------
             *         1
             *    w - -----------------
             *     1
             *          w+h - ---------
             *         w+2h - ...
             *
             * To determine how many terms needed, let
             * Q(0) = w, Q(1) = w(w+h) - 1,
             * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
             * When Q(k) > 1e4 good for single 
             * When Q(k) > 1e9 good for double 
             * When Q(k) > 1e17 good for quadruple 
             */
            /* determine k */
            float               t, v;
            float               q0, q1, h, tmp;
            int                 k, m;

            w = (n + n) / (float) x;
            h = (float) 2.0 / (float) x;
            q0 = w;
            z = w + h;
            q1 = w * z - (float) 1.0;
            k = 1;
            while (q1 < (float) 1.0e9)
            {
                k += 1;
                z += h;
                tmp = z * q1 - q0;
                q0 = q1;
                q1 = tmp;
            }
            m = n + n;
            for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
                t = one / (i / x - t);
            a = t;
            b = one;
            /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
             *  Hence, if n*(log(2n/x)) > ...
             *  single 8.8722839355e+01
             *  double 7.09782712893383973096e+02
             *  long double 1.1356523406294143949491931077970765006170e+04
             *  then recurrent value may overflow and the result is 
             *  likely underflow to zero
             */
            tmp = (float) n;
            v = two / x;
            tmp = tmp * RtIEEElogf(RtIEEEfabsf(v * tmp));
            if (tmp < (float) 8.8721679688e+01)
            {
                for (i = n - 1, di = (float) (i + i); i > 0; i--)
                {
                    temp = b;
                    b *= di;
                    b = b / x - a;
                    a = temp;
                    di -= two;
                }
            }
            else
            {
                for (i = n - 1, di = (float) (i + i); i > 0; i--)
                {
                    temp = b;
                    b *= di;
                    b = b / x - a;
                    a = temp;
                    di -= two;
                    /* scale b to avoid spurious overflow */
                    if (b > (float) 1e10)
                    {
                        a /= b;
                        t /= b;
                        b = one;
                    }
                }
            }
            b = (t * RtIEEEj0f(x) / b);
        }
    }
    if (sgn == 1)
        RWRETURN(-b);
    else
        RWRETURN(b);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEexpm1f
 *  computes the value of exp(x) - 1 accurately even for tiny argument x .
 * \param x input
 * \return exp(x) - 1 
 */
float
RtIEEEexpm1f(float x)
{
    float               y, hi, lo, c = 0, t, e, hxs, hfx, r1;
    int                 k, xsb;
    unsigned int        hx;

    RWAPIFUNCTION(RWSTRING("RtIEEEexpm1f"));

    GET_FLOAT_WORD(hx, x);
    xsb = hx & 0x80000000;     /* sign bit of x */
    if (xsb == 0)
        y = x;
    else
        y = -x;                /* y = |x| */
    hx &= 0x7fffffff;          /* high word of |x| */

    /* filter out rhuge and non-finite argument */
    if (hx >= 0x4195b844)
    {
        /* if |x|>=27*ln2 */
        if (hx >= 0x42b17218)
        {
            /* if |x|>=88.721... */
            if (hx > 0x7f800000)
                RWRETURN(x + x); /* NaN */
            if (hx == 0x7f800000)
                RWRETURN((float) ((xsb == 0) ? x : -1.0)); /* exp(+-inf)=
                                                            * * * {
                                                            * * * inf,-1
                                                            * * * }
                                                            */
            if (x > o_threshold)
                RWRETURN(rhuge * rhuge); /* overflow */

        }

        if (xsb != 0)
        {
            /* x < -27*ln2, return -1.0 with inexact */
            if (x + tiny < (float) 0.0) /* raise inexact */
                RWRETURN(tiny - one); /* return -1 */

        }

    }

    /* argument reduction */
    if (hx > 0x3eb17218)
    {
        /* if  |x| > 0.5 ln2 */
        if (hx < 0x3F851592)
        {
            /* and |x| < 1.5 ln2 */
            if (xsb == 0)
            {
                hi = x - ln2_hi;
                lo = ln2_lo;
                k = 1;
            }
            else
            {
                hi = x + ln2_hi;
                lo = -ln2_lo;
                k = -1;
            }

        }
        else
        {
            k = (int) (invln2 * x +
                       ((xsb == 0) ? (float) 0.5 : (float) -0.5));
            t = (float) k;
            hi = x - t * ln2_hi; /* t*ln2_hi is exact here */
            lo = t * ln2_lo;

        }

        x = hi - lo;
        c = (hi - x) - lo;

    }

    else if (hx < 0x33000000)
    {
        /* when |x|<2**-25, return x */
        t = rhuge + x;         /* return x with inexact flags when x!=0 */
        RWRETURN(x - (t - (rhuge + x)));
    }

    else
        k = 0;

    /* x is now in primary range */
    hfx = (float) 0.5  *x;

    hxs = x * hfx;
    r1 = one + hxs * (Q1 +
                      hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
    t = (float) 3.0 - r1 * hfx;
    e = hxs * ((r1 - t) / ((float) 6.0 - x * t));
    if (k == 0)
    {
        RWRETURN(x - (x * e - hxs));

    }                          /* c is 0 */
    else
    {

        e = (x * (e - c) - c);
        e -= hxs;
        if (k == -1)
            RWRETURN((float) 0.5 * (x - e) - (float) 0.5);

        if (k == 1)
        {

            if (x < (float) -0.25)
            {

                RWRETURN(-(float) 2.0 * (e - (x + (float) 0.5)));
            }

            else
            {

                RWRETURN(one + (float) 2.0 * (x - e));
            }

        }

        if (k <= -2 || k > 56)
        {
            /* suffice to return exp(x)-1 */
            int                 i;

            y = one - (e - x);
            GET_FLOAT_WORD(i, y);
            SET_FLOAT_WORD(y, i + (k << 23)); /* add k to y's exponent */
            RWRETURN(y - one);

        }

        t = one;
        if (k < 23)
        {

            int                 i;

            SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000 >> k)); /* t=1-2^-k */
            y = t - (e - x);
            GET_FLOAT_WORD(i, y);
            SET_FLOAT_WORD(y, i + (k << 23)); /* add k to y's exponent */

        }
        else
        {

            int                 i;

            SET_FLOAT_WORD(t, ((0x7f - k) << 23)); /* 2^-k */
            y = x - (e + t);
            y += one;
            GET_FLOAT_WORD(i, y);
            SET_FLOAT_WORD(y, i + (k << 23)); /* add k to y's exponent */

        }

    }

    RWRETURN(y);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEfrexpf 
 * returns the mantissa and the exponent of the input value
 * \param x input
 * \param eptr receives exponent
 * \return mantissa 
 */
float
RtIEEEfrexpf(float x, int *eptr)
{
    int                 hx, ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEfrexpf"));

    GET_FLOAT_WORD(hx, x);
    ix = 0x7fffffff & hx;
    *eptr = 0;
    if (ix >= 0x7f800000 || (ix == 0))
        RWRETURN(x);           /* 0,inf,nan */
    if (ix < 0x00800000)
    {                          /* subnormal */
        x *= two25;
        GET_FLOAT_WORD(hx, x);
        ix = hx & 0x7fffffff;
        *eptr = -25;
    }
    *eptr += (ix >> 23) - 126;
    hx = (hx & 0x807fffff) | 0x3f000000;
    *(int *) &x = hx;
    RWRETURN(x);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEisinf 
 * is recommended by the IEEE 754 standard, and 
 * returns the value 1 when the double input value is infinite, or 0 otherwise
 * \param x input
 * \return 1 when the input value is infinite, or 0 otherwise
 */
int
RtIEEEisinf(double x)
{
    int                 hx, lx;

    RWAPIFUNCTION(RWSTRING("RtIEEEisinf"));

    EXTRACT_WORDS(hx, lx, x);
    hx &= 0x7fffffff;
    hx ^= 0x7ff00000;
    hx |= lx;
    RWRETURN((hx == 0));
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEisinff 
 * is recommended by the IEEE 754 standard, and 
 * returns the value 1 when the float input value is infinite, or 0 otherwise
 * \param x input
 * \return 1 when the input value is infinite, or 0 otherwise
 */
int
RtIEEEisinff(float x)
{
    int                 ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEisinff"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;
    ix ^= 0x7f800000;
    RWRETURN((ix == 0));
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEisnanf
 * is recommended by the IEEE 754 standard, and 
 * returns the value 1 when the float input value is not a number -- e.g.
 * the resulty of dividing zero by itself --  or 0 otherwise
 * \param x input
 * \return 1 when the input value is not a number, or 0 otherwise
 */
int
RtIEEEisnanf(float x)
{
    int                 ix;

    RWAPIFUNCTION(RWSTRING("RtIEEEisnanf"));

    GET_FLOAT_WORD(ix, x);
    ix &= 0x7fffffff;
    ix = 0x7f800000 - ix;
    RWRETURN((int) (((unsigned int) (ix)) >> 31));
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEldexpf 
 * increments the exponent of the first by the second argment
 * \param x input
 * \param n exponent delta
 * \return value modified float
 */
float
RtIEEEldexpf(float value, int exp)
{
    RWAPIFUNCTION(RWSTRING("RtIEEEldexpf"));

    if (!RtIEEEfinitef(value) || value == (float) 0.0)
        RWRETURN(value);
    value = RtIEEEscalbnf(value, exp);
    if (!RtIEEEfinitef(value) || value == (float) 0.0)
        _rtIEEEferrno = RTIEEEF_ERANGE;
    RWRETURN(value);
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEmodff
 * returns the signed fractional portion of x. 
 * \param x input
 * \param intptr receives integer portion
 * \return the signed fractional portion
 */
float
RtIEEEmodff(float x, float *iptr)
{
    int                 i0, j0;
    unsigned int        i;

    RWAPIFUNCTION(RWSTRING("RtIEEEmodff"));

    GET_FLOAT_WORD(i0, x);
    j0 = ((i0 >> 23) & 0xff) - 0x7f; /* exponent of x */
    if (j0 < 23)
    {                          /* integer part in x */
        if (j0 < 0)
        {                      /* |x|<1 */
            SET_FLOAT_WORD(*iptr, i0 & 0x80000000); /* *iptr = +-0 */
            RWRETURN(x);
        }
        else
        {
            i = (0x007fffff) >> j0;
            if ((i0 & i) == 0)
            {                  /* x is integral */
                unsigned int        ix;

                *iptr = x;
                GET_FLOAT_WORD(ix, x);
                SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */
                RWRETURN(x);
            }
            else
            {
                SET_FLOAT_WORD(*iptr, i0 & (~i));
                RWRETURN(x - *iptr);
            }
        }
    }
    else
    {                          /* no fraction part */
        unsigned int        ix;

        *iptr = x * one;
        GET_FLOAT_WORD(ix, x);
        SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */
        RWRETURN(x);
    }
}

/**
 * \ingroup rtieeef
 * \ref RtIEEEnextafterf
 * returns the next machine representable number from
 * x in direction y .
 * \param x input
 * \param y direction
 * \return the next machine representable number
 */
float
RtIEEEnextafterf(float x, float y)
{
    int                 hx, hy, ix, iy;

    RWAPIFUNCTION(RWSTRING("RtIEEEnextafterf"));

    GET_FLOAT_WORD(hx, x);
    GET_FLOAT_WORD(hy, y);
    ix = hx & 0x7fffffff;      /* |x| */
    iy = hy & 0x7fffffff;      /* |y| */

    if ((ix > 0x7f800000) ||   /* x is nan */
        (iy > 0x7f800000))     /* y is nan */
        RWRETURN(x + y);
    if (x == y)
        RWRETURN(x);           /* x=y, return x */
    if (ix == 0)
    {                          /* x == 0 */
        SET_FLOAT_WORD(x, (hy & 0x80000000) | 1); /* return +-minsubnormal */
        y = x * x;
        if (y == x)
            RWRETURN(y);
        else
            RWRETURN(x);       /* raise underflow flag */
    }
    if (hx >= 0)
    {                          /* x > 0 */
        if (hx > hy)
        {                      /* x > y, x -= ulp */
            hx -= 1;
        }
        else
        {                      /* x < y, x += ulp */
            hx += 1;
        }
    }
    else
    {                          /* x < 0 */
        if (hy >= 0 || hx > hy)
        {                      /* x < y, x -= ulp */
            hx -= 1;
        }
        else
        {                      /* x > y, x += ulp */
            hx += 1;
        }
    }
    hy = hx & 0x7f800000;
    if (hy >= 0x7f800000)
        RWRETURN(x + x);       /* overflow  */
    if (hy < 0x00800000)
    {                          /* underflow */
        y = x * x;
        if (y != x)
        {                      /* raise underflow flag */
            SET_FLOAT_WORD(y, hx);
            RWRETURN(y);
        }
    }
    SET_FLOAT_WORD(x, hx);
    RWRETURN(x);
}

/**
 * \ingroup rtIEEEf
 * \ref RtIEEErem_pio2f 
 * returns the remainder of x rem pi/2 in y[0]+y[1] 
 * \param x input 
 * \param y pointer to interval
 * \return the remainder
 */
int
RtIEEErem_pio2f(float x, float *y)
{
    float               z, w, t, r, fn;
    float               tx[3];
    int                 e0, i, j, nx, n, ix, hx;

    RWAPIFUNCTION(RWSTRING("RtIEEErem_pio2f"));

    GET_FLOAT_WORD(hx, x);
    ix = hx & 0x7fffffff;
    if (ix <= 0x3f490fd8)      /* |x| ~<= pi/4 , no need for reduction */
    {
        y[0] = x;
        y[1] = 0;
        RWRETURN(0);
    }
    if (ix < 0x4016cbe4)
    {                          /* |x| < 3pi/4, special case with n=+-1 */
        if (hx > 0)
        {
            z = x - pio2_1;
            if ((ix & 0xfffffff0) != 0x3fc90fd0)
            {                  /* 24+24 bit pi OK */
                y[0] = z - pio2_1t;
                y[1] = (z - y[0]) - pio2_1t;
            }
            else
            {                  /* near pi/2, use 24+24+24 bit pi */
                z -= pio2_2;
                y[0] = z - pio2_2t;
                y[1] = (z - y[0]) - pio2_2t;
            }
            RWRETURN(1);
        }
        else
        {                      /* negative x */
            z = x + pio2_1;
            if ((ix & 0xfffffff0) != 0x3fc90fd0)
            {                  /* 24+24 bit pi OK */
                y[0] = z + pio2_1t;
                y[1] = (z - y[0]) + pio2_1t;
            }
            else
            {                  /* near pi/2, use 24+24+24 bit pi */
                z += pio2_2;
                y[0] = z + pio2_2t;
                y[1] = (z - y[0]) + pio2_2t;
            }
            RWRETURN(-1);
        }
    }
    if (ix <= 0x43490f80)
    {                          /* |x| ~<= 2^7*(pi/2), medium size */
        t = RtIEEEfabsf(x);
        n = (int) (t * invpio2 + half);
        fn = (float) n;
        r = t - fn * pio2_1;
        w = fn * pio2_1t;      /* 1st round good to 40 bit */
        if (n < 32
            && (ix & 0xffffff00) != (unsigned int) npio2_hw[n - 1])
        {
            y[0] = r - w;      /* quick check no cancellation */
        }
        else
        {
            unsigned int        high;

            j = ix >> 23;
            y[0] = r - w;
            GET_FLOAT_WORD(high, y[0]);
            i = j - ((high >> 23) & 0xff);
            if (i > 8)
            {                  /* 2nd iteration needed, good to 57 */
                t = r;
                w = fn * pio2_2;
                r = t - w;
                w = fn * pio2_2t - ((t - r) - w);
                y[0] = r - w;
                GET_FLOAT_WORD(high, y[0]);
                i = j - ((high >> 23) & 0xff);
                if (i > 25)
                {              /* 3rd iteration need, 74 bits acc */
                    t = r;     /* will cover all possible cases */
                    w = fn * pio2_3;
                    r = t - w;
                    w = fn * pio2_3t - ((t - r) - w);
                    y[0] = r - w;
                }
            }
        }
        y[1] = (r - y[0]) - w;
        if (hx < 0)
        {
            y[0] = -y[0];
            y[1] = -y[1];
            RWRETURN(-n);
        }
        else
            RWRETURN(n);
    }
    /* 
     * all other (large) arguments
     */
    if (ix >= 0x7f800000)
    {                          /* x is inf or NaN */
        y[0] = y[1] = x - x;
        RWRETURN(0);
    }
    /* set z = scalbn(|x|,ilogb(x)-7) */
    e0 = (ix >> 23) - 134;     /* e0 = ilogb(z)-7; */
    SET_FLOAT_WORD(z, ix - ((int) (e0 << 23)));
    for (i = 0; i < 2; i++)
    {
        tx[i] = (float) ((int) (z));
        z = (z - tx[i]) * two8;
    }
    tx[2] = z;
    nx = 3;
    while (tx[nx - 1] == zero)
        nx--;                  /* skip zero term */
    n = _rtIEEErem_pio2f(tx, y, e0, nx, 2, two_over_pi);
    if (hx < 0)
    {
        y[0] = -y[0];
        y[1] = -y[1];
        RWRETURN(-n);
    }
    RWRETURN(n);
}
