/*******************************************************************************
+
+  LEDA 3.5
+
+  floatingpoint.h
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/

// This is the file floatingpoint.h written by
// Christoph Burnikel (contact burnikel@mpi-sb.mpg.de)
// It contains declarations of non-standard basic functions for
// hardware-specific floating point arithmetic with doubles.
// The corresponding implementations are in floatingpoint.c

// We distinguish between two types of hardware:
// BIG_ENDIAN    machines (double starts with high word)
// LITTLE_ENDIAN machines (double starts with low  word)
// Please find out what type of machine you are using
// For LITTLE ENDIAN machines use flag 
//         -DLITTLE_ENDIAN_MACHINE 
// in compilation

// You can try the test program floatingpoint_test.c to 
// see whether this works fine.
// Note that we test by default in floatingpoint.c whether the
// current platform agrees with the used flags -DLITTLE_ENDIAN_MACHINE
// or -DWORD_LENGTH_64.
// For example, if on an old LITTLE_ENDIAN machine like 
// linux with a x86 processor the flag -DLITTLE_ENDIAN_MACHINE is omitted,
// this is automatically detected and an error message is 
// displayed.

#ifndef LEDA_FLOATINGPOINT_H
#define LEDA_FLOATINGPOINT_H

#include <LEDA/system.h>
#include <math.h>

// -------------------------------------------------------------------
// the special values +-0, +-00, NaN as double constants
// -------------------------------------------------------------------

extern const double double_min;
extern const double NaN_double;
extern const double pInf_double;
extern const double nInf_double;
extern const double pZero_double;
extern const double nZero_double;


// -------------------------------------------------------------------
// some commonly used functions for double arithmetic
// -------------------------------------------------------------------

// Note: Many compilers like sunpro do not know bool (!) and hence 
// we leave all functions return "long" instead of "bool"

inline long sign_bit_set(const double& x);
// tells whether the sign bit of x is set

inline void clear_sign_bit(double& x);
// clears the sign bit of x to zero (making x nonnegative)

inline long is_not_underflowed(const double& x);
// tells whether the biased exponent of x is nonzero

inline long is_finite(const double& x);
// tells whether the biased exponent of x is smaller than 2047

inline long is_infinite(const double& x);
// tells whether the biased exponent of x is equal to 2047

inline long is_single_precision(const double& x);
// tells whether the last 27 bit of x's significant are zero

inline long is_nan(const double& x);
// tells whether the last 27 bit of x's significant are zero

inline double float_abs(const double& x);
// returns the absolute value of x

double power_two(long l);
// returns 2^l as a double but works only if -1022 <= l <= 1023

int float_sign(const double& x);
// returns the sign of x

double compose_parts(int sign_1, unsigned long exp_11,
          unsigned long most_sig_20, unsigned long least_sig_32);
// returns a double composed of the following parts:
// 1 sign bit                       : sign_1
// 11 bits of the exponent          : exponent_11
// higher 20 bit of the significant : most_sig_20
// lower  32 bit of the significant : least_sig_32

void read_parts(const double& a,
               int& sign_1, long& exp_11,
               unsigned long& most_sig_20, unsigned long& least_sig_32);
// disassembles a into the four parts described in the function
// compose_parts (hence this is the inverse of compose_parts)

long get_biased_exponent(const double& x);
// returns the exponent of x in unsigned long format, biased by +1023

long get_unbiased_exponent(const double& x);
// read the unbiased exponent of x in unsigned long format

long binary_equal(const double& x, const double& y);
// says whether x and y are equal as binaries


// ------------------------------------------------------------------------
// Now the definitions of inline functions follow
// ------------------------------------------------------------------------

#ifndef WORD_LENGTH_64
static const unsigned long signbit_mask          = 0x80000000;
static const unsigned long clear_signbit_mask    = 0x7fffffff;
static const unsigned long exponent_mask         = 0x7ff00000;
static const unsigned long single_precision_mask = 0x07ffffff;
static const unsigned long high_mantissa_mask    = 0x000fffff;
#else
static const unsigned long signbit_mask          = 0x8000000000000000;
static const unsigned long clear_signbit_mask    = 0x7fffffffffffffff;
static const unsigned long exponent_mask         = 0x7ff0000000000000;
static const unsigned long single_precision_mask = 0x0000000007ffffff;
static const unsigned long mantissa_mask         = 0x000fffffffffffff;
#endif

#ifndef LITTLE_ENDIAN_MACHINE

//-------------------------------------------------------------------------
// code for big-endian architecture (default)
// ------------------------------------------------------------------------

inline long sign_bit_set(const double& x)
{
   return *((unsigned long*) &x) & signbit_mask;
}

inline void clear_sign_bit(double& x)
{
   *((unsigned long*) &x) &= clear_signbit_mask;
}

inline double float_abs(const double& y)
{
   double x = y;
   *((unsigned long*) &x) &= clear_signbit_mask;
   return x;
}

inline long is_not_underflowed(const double& x)
{
   return *((unsigned long*) &x) & exponent_mask;
}

inline long is_finite(const double& x)
{
   unsigned long exponent_part = *((unsigned long*) &x) & exponent_mask;
   return exponent_part ^ exponent_mask;
}

inline long is_infinite(const double& x)
{
   unsigned long exponent_part = *((unsigned long*) &x) & exponent_mask;
   return (exponent_part == exponent_mask);
}

inline long is_nan(const double& x)
{
   unsigned long* ptr = (unsigned long*) &x ;
   unsigned long exponent_part = *ptr & exponent_mask;
#ifndef WORD_LENGTH_64
   return 
    (exponent_part == exponent_mask && (*(ptr+1) || *ptr & high_mantissa_mask));
#else
   return
    (exponent_part == exponent_mask && (*ptr & mantissa_mask));
#endif
}

inline long is_single_precision(const double& x)
{
   unsigned long* low_part = (unsigned long*) &x;
#ifndef WORD_LENGTH_64
   return !(*++low_part & single_precision_mask);
#else
   return !(*low_part   & single_precision_mask);
#endif
}

#else

// ------------------------------------------------------------------------
// code for little-endian architecture (only on some 32 bit machines)
// ------------------------------------------------------------------------

inline long sign_bit_set(const double& x)
{
   unsigned long* p = (unsigned long*) &x;
   return *(++p) & signbit_mask;
}

inline void clear_sign_bit(double& x)
{
   unsigned long* p = (unsigned long*) &x;
   *(++p) &= clear_signbit_mask;
}

inline double float_abs(const double& y)
{
  double x = y;
  unsigned long* p = (unsigned long*) &x;
   *(++p) &= clear_signbit_mask;
  return x;
}

inline long is_not_underflowed(const double& x)
{
   unsigned long* p = (unsigned long*) &x;
   return *(++p) & exponent_mask;
}

inline long is_finite(const double& x)
{
   unsigned long* p = (unsigned long*) &x;
   unsigned long exponent_part = *(++p) & exponent_mask;
   return exponent_part ^ exponent_mask;
}

inline long is_infinite(const double& x)
{
   unsigned long* p = (unsigned long*) &x;
   unsigned long exponent_part = *(++p) & exponent_mask;
   return (exponent_part == exponent_mask);
}

inline long is_nan(const double& x)
{
   unsigned long* ptr = (unsigned long*) &x ;
   unsigned long exponent_part = *(ptr+1) & exponent_mask;
   return ( (exponent_part == exponent_mask) 
         && (*ptr || *(ptr+1) & high_mantissa_mask));
}

inline long is_single_precision(const double& x)
{
   unsigned long* low_part = (unsigned long*) &x;
   return !(*low_part & 0x07ffffff);
}

#endif

#endif
