quadsim.c

Go to the documentation of this file.

/*

  Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.

  This program is free software; you can redistribute it and/or modify it
  under the terms of version 2 of the GNU General Public License as
  published by the Free Software Foundation.

  This program is distributed in the hope that it would be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  

  Further, this software is distributed without any warranty that it is
  free of the rightful claim of any third person regarding infringement 
  or the like.  Any license provided herein, whether implied or 
  otherwise, applies only to this software file.  Patent licenses, if 
  any, provided herein do not apply to combinations of this program with 
  other software, or any other product whatsoever.  

  You should have received a copy of the GNU General Public License along
  with this program; if not, write the Free Software Foundation, Inc., 59
  Temple Place - Suite 330, Boston MA 02111-1307, USA.

  Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
  Mountain View, CA 94043, or:

  http://www.sgi.com

  For further information regarding this notice, see:

  http://oss.sgi.com/projects/GenInfo/NoticeExplan

*/


/* =======================================================================
 * =======================================================================
 *
 *
 * =======================================================================
 * =======================================================================
 */



#include "quad.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <fp_class.h>
#include "defs.h"
#include "quadsim.h"

/* For fp_class */

#define DMANTWIDTH      52
#define DEXPWIDTH       11
#define DSIGNMASK       0x7fffffffffffffffll
#define DEXPMASK        0x800fffffffffffffll
#define DQNANBITMASK    0xfff7ffffffffffffll

#define MANTWIDTH       23
#define EXPWIDTH        8
#define SIGNMASK        0x7fffffff
#define EXPMASK         0x807fffff
#define QNANBITMASK     0xffbfffff

typedef union
{
        struct
        {
                UINT32 hi;
                UINT32 lo;
        } word;

        double  d;
} du;

static const du         m_twop31 =
{0xc1e00000,    0x00000000};

static const du         twop31m1 =
{0x41dfffff,    0xffc00000};

static const du         twop32m1 =
{0x41efffff,    0xffe00000};

static const du         twop52 =
{0x43300000,    0x00000000};

static const    du      twop62 =
{0x43d00000,    0x00000000};

static const    du      m_twop63 =
{0xc3e00000,    0x00000000};

static const    du      twop63 =
{0x43e00000,    0x00000000};

static const    du      twop64 =
{0x43f00000,    0x00000000};

static const    du      twopm916 =
{0x06b00000,    0x00000000};

/* Note: the name 'infinity' may conflict with <math.h> */
static const    du      myinfinity =
{0x7ff00000,    0x00000000};

extern  INT     c_q_le(QUAD x, QUAD y, INT *p_err );
extern  INT     c_q_ge(QUAD x, QUAD y, INT *p_err );
extern  INT32   c_ji_qint(QUAD x, INT *p_err );
extern  INT32   c_fp_class_q(QUAD x);
extern  UINT32  c_ji_quint(QUAD x, INT *p_err );
extern  INT64   c_ki_qint(QUAD x, INT *p_err );
extern  UINT64  c_ki_quint(QUAD x, INT *p_err );
extern  float   c_sngl_q(QUAD x, INT *p_err );
extern  double  c_dble_q(QUAD x, INT *p_err );
extern  QUAD    c_q_flotj(INT32 n, INT *p_err );
extern  QUAD    c_q_flotju(UINT32 n, INT *p_err );
extern  QUAD    c_q_flotk(INT64 n, INT *p_err );
extern  QUAD    c_q_flotku(UINT64 n, INT *p_err );
extern  QUAD    c_q_ext( float x, INT *p_err );
extern  QUAD    c_q_extd(double x, INT *p_err );
extern  QUAD    c_q_trunc(QUAD x, INT *p_err );
extern  double  trunc(double x);

#if defined(_GCC_NO_PRAGMAWEAK) || defined(__CYGWIN__)

#define c_q_trunc __c_q_trunc
#define c_q_ge __c_q_ge
#define c_ki_qint __c_ki_qint
#define c_q_le __c_q_le

#endif

/* quad to INT32 */

#pragma weak c_ji_qint = __c_ji_qint
#define c_ji_qint __c_ji_qint

INT32
c_ji_qint(QUAD x, INT *p_err )
{
QUAD    xi;
INT32   result;

        *p_err = 0;

        if ( x.hi != x.hi )
        {
                *p_err = 1;

                return ( (INT32)x.hi );
        }

        if ( fabs(x.hi) == myinfinity.d )
        {
                *p_err = 1;

                return ( (INT32)x.hi );
        }

        xi = c_q_trunc(x, p_err);

        if ( xi.hi > twop31m1.d )
                *p_err = 1;

        if ( xi.hi < m_twop31.d )
                *p_err = 1;

        result = (INT32)xi.hi;

        return ( result );
}

/* quad to UINT32 */

#pragma weak c_ji_quint = __c_ji_quint
#define c_ji_quint __c_ji_quint

UINT32
c_ji_quint(QUAD x, INT *p_err )
{
QUAD    xi;
INT64   n;
UINT32  result;

        *p_err = 0;

        if ( x.hi != x.hi )
        {
                *p_err = 1;

                return ( (UINT32)x.hi );
        }

        if ( fabs(x.hi) == myinfinity.d )
        {
                *p_err = 1;

                return ( (UINT32)x.hi );
        }

        xi = c_q_trunc(x, p_err);

        if ( xi.hi > twop32m1.d )
        {
                *p_err = 1;

                return ( (UINT32)x.hi );
        }

        if ( xi.hi < 0.0 )
        {
                *p_err = 1;

                return ( (UINT32)x.hi );
        }

        n = c_ki_qint(xi, p_err);

        result = n & 0xffffffff;

        return ( result );
}

/* quad to INT64 */

#pragma weak c_ki_qint = __c_ki_qint
#define c_ki_qint __c_ki_qint

INT64
c_ki_qint(QUAD x, INT *p_err )
{
QUAD    xi;
INT64   t;

        *p_err = 0;

        if ( x.hi != x.hi )
        {
                *p_err = 1;

                return ( (INT64)x.hi );
        }

        if ( fabs(x.hi) == myinfinity.d )
        {
                *p_err = 1;

                return ( (INT64)x.hi );
        }

        xi = c_q_trunc(x, p_err);

        if ( (xi.hi > twop63.d) ||
             ((xi.hi == twop63.d) && (xi.lo >= 0.0))
           )
        {
                /* we get overflow here, so just let the hardware do it */

                *p_err = 1;

                return ( (INT64)xi.hi );
        }

        if ( (xi.hi < m_twop63.d) ||
             ((xi.hi == m_twop63.d) && (xi.lo < 0.0))
           )
        {
                /* we get overflow here, so just let the hardware do it */

                if ( xi.hi == m_twop63.d )
                {
                        xi.hi *= 2.0;
                }

                *p_err = 1;

                return ( (INT64)xi.hi );
        }

        if ( fabs(xi.hi) > twop62.d )
        {
                /* add things up very carefully to avoid overflow */

                t = 0.5*xi.hi;
                return ( ((INT64)xi.lo) + t + t );
        }

        /* just do it */

        return ( (INT64)xi.hi + (INT64)xi.lo );
}

/* quad to UINT64 */

#pragma weak c_ki_quint = __c_ki_quint
#define c_ki_quint __c_ki_quint

UINT64
c_ki_quint(QUAD x, INT *p_err )
{
QUAD    xi;
double  z;
UINT64  t;

        *p_err = 0;

        if ( x.hi != x.hi )
        {
                *p_err = 1;

                return ( (UINT64)x.hi );
        }

        if ( fabs(x.hi) == myinfinity.d )
        {
                *p_err = 1;

                return ( (UINT64)x.hi );
        }

        xi = c_q_trunc(x, p_err);

        if ( (xi.hi > twop64.d) ||
             ((xi.hi == twop64.d) && (xi.lo > -1.0))
           )
        {
                *p_err = 1;

                return ( (UINT64)x.hi );
        }

        if ( xi.hi < 0.0 )
        {
                *p_err = 1;

                return ( (UINT64)x.hi );
        }

        if ( (xi.hi > twop63.d) ||
             (xi.hi == twop63.d) && (xi.lo >= 0.0)
           )
        {
                /* subtract 2**64       */

                z = xi.hi - twop64.d;

                /* normalize result */

                xi.hi = z + xi.lo;
                xi.lo = z - xi.hi + xi.lo;
        }

        if ( fabs(xi.hi) > twop62.d )
        {
                /* add things up very carefully to avoid overflow */

                t = 0.5*xi.hi;
                return ( ((UINT64)xi.lo) + t + t );
        }

        /* just do it */

        return ( (UINT64)xi.hi + (UINT64)xi.lo );
}

/* quad to float */

#pragma weak c_sngl_q = __c_sngl_q
#define c_sngl_q __c_sngl_q

float
c_sngl_q(QUAD x, INT *p_err )
{
        *p_err = 0;

        return ( (float)x.hi );
}

/* quad to double */

#pragma weak c_dble_q = __c_dble_q
#define c_dble_q __c_dble_q

double
c_dble_q(QUAD x, INT *p_err )
{
        *p_err = 0;

        return ( x.hi );
}

/* INT32 to quad */

#pragma weak c_q_flotj = __c_q_flotj
#define c_q_flotj __c_q_flotj

QUAD
c_q_flotj(INT32 n, INT *p_err )
{
QUAD    result;

        *p_err = 0;

        result.hi = n;
        result.lo = 0.0;
        return ( result );
}

/* UINT32 to quad */

#pragma weak c_q_flotju = __c_q_flotju
#define c_q_flotju __c_q_flotju

QUAD
c_q_flotju(UINT32 n, INT *p_err )
{
QUAD    result;

        *p_err = 0;

        result.hi = n;
        result.lo = 0.0;
        return ( result );
}

/* INT64 to quad */

#pragma weak c_q_flotk = __c_q_flotk
#define c_q_flotk __c_q_flotk

QUAD
c_q_flotk(INT64 n, INT *p_err )
{
INT64 m;
double  w, ww;
QUAD  result;

        *p_err = 0;

        m = (n >> 11);
        m <<= 11;

        w = m;

        ww = n - m;

        /* normalize result */

        result.hi = w + ww;
        result.lo = w - result.hi + ww;

        return ( result );
}

/* UINT64 to QUAD */

#pragma weak c_q_flotku = __c_q_flotku
#define c_q_flotku __c_q_flotku

QUAD
c_q_flotku(UINT64 n, INT *p_err )
{
UINT64 m;
double  w, ww;
QUAD  result;

        *p_err = 0;

        m = (n >> 11);
        m <<= 11;

        w = m;

        ww = n - m;

        /* normalize result */

        result.hi = w + ww;
        result.lo = w - result.hi + ww;

        return ( result );
}

/* float to quad */

#pragma weak c_q_ext = __c_q_ext
#define c_q_ext __c_q_ext

QUAD
c_q_ext( float x, INT *p_err )
{
QUAD    result;

        *p_err = 0;

        result.hi = x;
        result.lo = 0.0;
        return ( result );
}

/* double to quad */

#pragma weak c_q_extd = __c_q_extd
#define c_q_extd __c_q_extd

QUAD
c_q_extd(double x, INT *p_err )
{
QUAD    result;

        *p_err = 0;

        result.hi = x;
        result.lo = 0.0;
        return ( result );
}

/* truncates a long double, i.e. computes the integral part of
   a long double
*/

#pragma weak c_q_trunc = __c_q_trunc
#define c_q_trunc __c_q_trunc

QUAD
c_q_trunc(QUAD x, INT *p_err )
{
double  uhi, ulo;
QUAD    result;

        uhi = x.hi;
        ulo = x.lo;

        *p_err = 0;

        if ( uhi != uhi )
        {
                result.hi = uhi;
                result.lo = ulo;

                return ( result );
        }

        if ( uhi >= 0.0 )
        {
                if ( uhi < twop52.d )
                {
                        /* binary point occurs in uhi; truncate uhi to an integer
                        */

                        result.hi = trunc(uhi);

                        result.lo = 0.0;

                        if ( result.hi < uhi )
                                return ( result );

                        /* must adjust result by one if ulo < 0.0 */

                        if ( ulo < 0.0 )
                        {
                                result.hi -= 1.0;

                                return ( result );
                        }

                        return ( result );
                }
                else if ( fabs(ulo) < twop52.d )
                {
                        /* binary point occurs in ulo; truncate ulo to an integer
                        */

                        result.hi = uhi;

                        result.lo = trunc(ulo);

                        if ( result.lo > ulo )
                        {
                                result.lo -= 1.0;
                        }

                        return ( result );
                }

                /* arg is an integer */

                result.hi = uhi;
                result.lo = ulo;

                return ( result );
        }
        else
        {
                if ( fabs(uhi) < twop52.d )
                {
                        /* binary point occurs in uhi; truncate uhi to an integer
                        */

                        result.hi = trunc(uhi);

                        result.lo = 0.0;

                        if ( result.hi > uhi )
                                return ( result );

                        /* must adjust result by one if ulo > 0.0 */

                        if ( ulo > 0.0 )
                        {
                                result.hi += 1.0;

                                return ( result );
                        }

                        return ( result );
                }
                else if ( fabs(ulo) < twop52.d )
                {
                        /* binary point occurs in ulo; truncate ulo to an integer
                        */

                        result.hi = uhi;

                        result.lo = trunc(ulo);

                        if ( result.lo < ulo )
                        {
                                result.lo += 1.0;
                        }

                        return ( result );
                }

                /* arg is an integer */

                result.hi = uhi;
                result.lo = ulo;

                return ( result );
        }
}

/* ====================================================================
 *
 * FunctionName: c_fp_class_q
 *
 * Description: c_fp_class_q(u) returns the floating point class 
 *   to which u belongs, and should work equivalently to the
 *   functions declared in <fp_class.h>.
 *
 * ====================================================================
 */


#pragma weak c_fp_class_q = __c_fp_class_q
#define c_fp_class_q __c_fp_class_q

INT
fp_class_d( double x )
{
  UINT64 ll, exp, mantissa;
  INT32 sign;

  ll = *(UINT64*)&x;
  exp = (ll >> DMANTWIDTH);
  sign = (exp >> DEXPWIDTH);
  exp &= 0x7ff;
  mantissa = (ll & (DSIGNMASK & DEXPMASK));
  if ( exp == 0x7ff ) {
    /* result is an infinity, or a NaN */
    if ( mantissa == 0 )
      return ( (sign == 0) ? FP_POS_INF : FP_NEG_INF );
    else if ( mantissa & ~DQNANBITMASK )
      return ( FP_QNAN );
    else
      return ( FP_SNAN );
  }

  if ( exp == 0 ) {
    if ( mantissa == 0 )
      return ( (sign == 0) ? FP_POS_ZERO : FP_NEG_ZERO );
    else
      return ( (sign == 0) ? FP_POS_DENORM : FP_NEG_DENORM );
  }
  else
    return ( (sign == 0) ? FP_POS_NORM : FP_NEG_NORM );
}

INT32
c_fp_class_q(QUAD x)
{       
QUAD    y;
INT     err;
INT32   class;

        /* Notice that the definition of denormal for quad precision is
         * chosen so that normal quads have a full 107 bits precision
         * in their mantissas.
         */

        class = fp_class_d(x.hi);

        if ( (class != FP_POS_NORM) && (class != FP_NEG_NORM) )
                return ( class );

        class = fp_class_d(x.lo);

        if ( (class == FP_POS_DENORM) && (class == FP_NEG_DENORM) )
                return ( class );

        y.hi = twopm916.d;
        y.lo = 0.0;

        if ( c_q_ge(x, y, &err) )
                return ( FP_POS_NORM );

        y.hi = -y.hi;

        if ( c_q_le(x, y, &err) )
                return ( FP_POS_NORM );

        if ( x.hi > 0.0 )
                return ( FP_POS_DENORM );

        return ( FP_NEG_DENORM );
}