---

gul_ieee.cpp

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/*
 *  implementation of missing ieee functions for windows platforms
 */

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */

#include "stdafx.h"

#if defined(__MINGW32__) || (defined(_MSC_VER))

#include "gul_types.h"

using gul::int32;
using gul::uint32;

#define GET_HIGH_WORD(w,x) w = ((int32 *)&(x))[1]
#define GET_LOW_WORD(w,x) w = ((int32 *)&(x))[0]
#define EXTRACT_WORDS(h,l,x) h = ((int32 *)&(x))[1]; l = ((int32 *)&(x))[0];
#define SET_HIGH_WORD(x,w) ((int32 *)&(x))[1] = w;

/*
 * copysign(double x, double y)
 * copysign(x,y) returns a value with the magnitude of x and
 * with the sign bit of y.
 */

GULAPI double copysign(double x, double y)
{
  uint32 hx,hy;
  GET_HIGH_WORD(hx,x);
  GET_HIGH_WORD(hy,y);
  SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
  return x;
}

/* ilogb(double x)
 * return the binary exponent of non-zero x
 * ilogb(0) = 0x80000001
 * ilogb(inf/NaN) = 0x7fffffff (no signal is raised)
 */

GULAPI int ilogb(double x)
{
  int32 hx,lx,ix;

  GET_HIGH_WORD(hx,x);
  hx &= 0x7fffffff;
  if(hx<0x00100000)
  {
    GET_LOW_WORD(lx,x);
    if((hx|lx)==0) 
      return 0x80000001;      /* ilogb(0) = 0x80000001 */
    else                      /* subnormal x */
    {
      if(hx==0)
      {
        for (ix = -1043; lx>0; lx<<=1) ix -=1;
      } 
      else 
      {
        for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
      }
      return ix;
    }
  }
  else if (hx<0x7ff00000) 
    return (hx>>20)-1023;

  return 0x7fffffff;
}

/* 
 * scalbn (double x, int n)
 * scalbn(x,n) returns x* 2**n  computed by  exponent  
 * manipulation rather than by actually performing an 
 * exponentiation or a multiplication.
 */

static const double
two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
hug   = 1.0e+300,
tiny   = 1.0e-300;

GULAPI double scalbn (double x, int n)
{
  int32 k,hx,lx;
  EXTRACT_WORDS(hx,lx,x);
  k = (hx&0x7ff00000)>>20;            /* extract exponent */
  if (k==0)                           /* 0 or subnormal x */
  {                             
    if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
    x *= two54; 
    GET_HIGH_WORD(hx,x);
    k = ((hx&0x7ff00000)>>20) - 54; 
    if (n< -50000) return tiny*x;   /*underflow*/
  }
  if (k==0x7ff) return x+x;             /* NaN or Inf */
  k = k+n; 
  if (k >  0x7fe) return hug*copysign(hug,x); /* overflow  */
  if (k > 0)                            /* normal result */
  { 
    SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); 
    return x;
  }
  if (k <= -54)
    if (n > 50000)      /* in case integer overflow in n+k */
      return hug*copysign(hug,x);       /*overflow*/
    else 
      return tiny*copysign(tiny,x);     /*underflow*/
  k += 54;                              /* subnormal result */
  SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
  return x*twom54;
}

#endif

---