Standard C



August 01, 1991
URL:http://www.drdobbs.com/standard-c/184402404

August 1991/Standard C/Listing 1

Listing 1

/* _Dtest function - IEEE 754 version */
#include "xmath.h"

short _Dtest (double *px)
       {/*      categorize *px */
       unsigned short *ps = (unsigned short *)px;
       short xchar = (ps[_D0] & _DMASK) >> _DOFF;

       if (xchar == _DMAX)     /* NaN or INF */
              return (ps[_D0] & _DFRAC || ps[_D1]
                     || ps[_D2] || ps[_D3] ? NAN : INF);
       else if (0 < xchar || ps[_D0] & _DFRAC
              || ps[_D1] || ps[_D2] || ps[_D3])
              return (FINITE);        /* finite */
       else
              return (0);     /* zero */
       }

/* End of File */
August 1991/Standard C/Listing 10

Listing 10

/* _Dnorm function - IEEE 754 version */
#include "xmath.h"

short _Dnorm(unsigned short *ps)
       { /* normalize double fraction */
       short xchar;
       unsigned short sign = ps[_D0] & _DSIGN;

       xchar = 0;
       if ((ps[_D0] &= _DFRAC) != 0 || ps[_D1]
              || ps[_D2] || ps[_D3])
              {       /* nonzero, scale */
              for (; ps[_D0] == 0; xchar -= 16)
                     {       /* shift left by 16 */
                     ps[_D0] = ps[_D1], ps[_D1] = ps[_D2];
                     ps[_D2] = ps[_D3], ps[_D3] = 0;
                     }
              for (; ps[_D0] < 1<< D0FF; -xchar)
                     {       /* shift left by 1 */
                     ps[_D0] = ps[_D0] << 1 | ps[_D1] >> 15;
                     ps[_D1] = ps[_D1] << 1 | ps[_D2] >> 15;
                     ps[_D2] = ps[_D2] << 1 | ps[_D3] >> 15;
                     ps[_D3] <<= 1;
                     }
              for (; 1<<_D0FF+1 <= ps[_D0]; ++xchar)
                     {       /* shift right by 1 */
                     ps[_D3] = ps[_D3] >> 1 | ps[_D2] << 15;
                     ps[_D2] = ps[_D2] >> 1 | ps[_D1] << 15;
                     ps[_D1] = ps[_D1] >> 1 | ps[_D0] << 15;
                     ps[_D0] >>= 1;
                     }
              ps[_D0] &= _DFRAC;
              }
       ps[_D0] |= sign;
       return (xchar);
       }

/* End of File */
August 1991/Standard C/Listing 11

Listing 11

/* fmod function */
#include "xmath.h"

double (fmod)(double x, double y)
       {       /* compute fmod(x, y) */
       const short errx = _Dtest(&x);
       const short erry = _Dtest(&y);

       if (errx == NAN || erry == NAN || errx == INF || erry == 0)
              {       /* fmod undefined */
              errno = EDOM;
              return (errx == NAN ? x : erry == NAN ? y : _Nan._D);
              }
       else if (errx == 0 || erry == INF)



              return (x); /* fmod(0,nonzero) or fmod(finite,INF) */
       else
              {       /* fmod(finite,finite) */
              double t;
              short n, neg, ychar;

              if (y < 0.0)
                     y = -y;
              if (x < 0.0)
                     x = -x, neg = 1;
              else
                     neg = 0;
              for (t = y, _Dunscale(&ychar, &t), n = 0; ; )
                     {       /* subtract |y| until |x| y| */
                     short xchar;

                     t = x;
                     if (n < 0 || _Dunscale(&xchar, &t) == 0
                            || (n = xchar - ychar) 0)
                            return (neg ? -x : x);
                     for (; 0 <= n; -n)
                            {       /* try to subtract |y|*2^n */
                            t = y, _Dscale(&t, n);
                            if (t <= x)
                                   {
                                   x -= t;
                                   break;
                                   }
                            }
                     }
              }

August 1991/Standard C/Listing 2

Listing 2

/* ceil function */
#include "xmath.h"

double (ceil)(double x)
       {       /* compute ceil(x) */
       return (_Dint(&x, 0) < 0 && 0.0 < x ? x + 1.0 : x);
       }

/* End of File */
August 1991/Standard C/Listing 3

Listing 3

/* floor function */
#include "xmath.h"

double (floor)(double x)
       {       /* compute floor(x) */
       return (_Dint(&x, 0) < 0 && x < 0.0 ? x - 1.0 : x);
       }

/* End of File */
August 1991/Standard C/Listing 4

Listing 4

/* _Dint function - IEEE 754 version */
#include "xmath.h"

short _Dint(double *px, short xexp)
       {       /* test and drop (scaled) fraction bits */
       unsigned short *ps = (unsigned short *)px;
       unsigned short frac = ps[_D0] & _DFRAC
              || ps[_D1] || ps[_D2] || ps[_D3];
       short xchar = (ps[_D0] & _DMASK) >> _DOFF;

       if (xchar == 0 && !frac)
              return (0);     /* zero */
       else if (xchar != _DMAX)
              ;       /* finite */
       else if (!frac)
              return (INF);
       else
              {       /* NaN */
              errno = EDOM;
              return (NAN);
              }
       xchar = (_DBIAS+48+_DOFF+1) - xchar - xexp;
       if (xchar < 0)
              return (0);     /* no frac bits to drop */
       else if ((48+_DOFF)  xchar)
              {       /* all frac bits */
              ps[_D0] = 0, ps[_D1] : 0;
              ps[_D2] = 0, ps[_D3] = 0;
              return (FINITE);
              }
       else
              {       /* strip out frac bits */
              static const unsigned short mask[] = {
                     0x0000, 0x0001, 0x0003, 0x0007,
                     0x000f, 0x001f, 0x003f, 0x007f,
                     0x00ff, 0x01ff, 0x03ff, 0x07ff,
                     0x0fff, 0x1fff, 0x3fff, 0x7fff};
              static const size_t sub[] = {_D3, _D2, _D1, _D0};

              frac = mask[xchar & 0xf];
              xchar >>= 4;
              frac &= ps[sub[xchar]];
              ps[sub[xchar]] ^= frac;
              switch (xchar)
                     {      /* cascade through! */
              case 3:
                     frac |= ps[_D1], ps[_D1] = 0;
              case 2:
                     frac |= ps[_D2], ps[_D2] = 0;
              case 1:
                     frac |= ps[_D3], ps[_D3] : 0;
                     }
              return (frac ? FINITE : 0);
              }
       }

/* End of File */
August 1991/Standard C/Listing 5

Listing 5

/* modf function */
#include "xmath.h"

double (modf)(double x, double *pint)
       {       /* compute modf(x, &intpart) */
       *pint = x;
       switch (_Dint(pint, 0))
              {       /* test for special codes */
       case NAN:
              return (x);
       case INF:
       case 0:
              return (0.0);
       default:        /* finite */
              return (x - *pint);
              }
       }

/* End of File */
August 1991/Standard C/Listing 6

Listing 6

/* frexp function */
#include "xmath.h"

double (frexp)(double x, int *pexp)
       {       /* compute frexp(x, &i) */
       short binexp;

       switch (_Dunscale(&binexp, &x))
              {       /* test for special codes */
       case NAN:
       case INF:
              errno = EDOM;
              *pexp = 0;
              return (x);
       case 0:
              *pexp = 0;
              return (0.0);
       default:        /* finite */
              *pexp = binexp;
              return (x);
              }
       }

/* End of File */
August 1991/Standard C/Listing 7

Listing 7

/* ldexp function */
#include "xmath.h"

double (ldexp)(double x, int xexp)
       {       /* compute ldexp(x, xexp) */
       switch (_Dtest(&x))
              {       /* test for special codes */
       case NAN:
              errno : EDOM;
              break;
       case INF:
              errno = ERANGE;
              break;
       case 0:
              break;
       default:         /* finite */
              if (0 <= _Dscale(&x, xexp))
                     errno = ERANGE;
              }
       return (x);
       }

/* End of File */
August 1991/Standard C/Listing 8

Listing 8

/* _Dunscale function - IEEE 754 version */
#include "xmath.h"

short _Dunscale(short *pex, double *px)
       {       /* separate *px to 1/2 <= |frac| < 1 and 2^*pex */
       unsigned short *ps = (unsigned short *)px;
       short xchar = (ps[_D0] & _DMASK) >> _D0FF;

       if (xchar == _DMAX)
              {       /* NaN or INF */
              *pex = 0;
              return (ps[_D0] & _DFRAC || ps[_D1]
                     || ps[_D2] || ps[_D3] ? NAN : INF);
              }
       else if (0 < xchar || (xchar = _Dnorm(ps)) != 0)
              {       /* finite, reduce to [1/2, 1) */
              ps[_D0] = ps[_D0] & _DMASK | _DBIAS << _D0FF;
              *pex = xchar - _DBIAS;
              return (FINITE);
              }
       else
              {       /* zero */
              *pex = 0;
              return (0);
              }
       }

/* End of File */

Terms of Service | Privacy Statement | Copyright © 2024 UBM Tech, All rights reserved.