---

guar_exact.h

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.
 */

#ifndef GUAR_EXACT_H
#define GUAR_EXACT_H

#include <math.h>
// #include <new.h>

namespace guar {

using gul::Assert;
using gul::InternalError;
using gul::ndebug;
using gul::Swap;

/* ---------------------------------------------------------
  Converts int-string to float
----------------------------------------------------------- */

template< class T >
void IntTo( const unsigned int na, const unsigned long *a, T& t )
{
  int i;

  if( na == 0 )
  {
    t = (T)0.0;
    return;
  }
  t = (T)a[na-1];
  
  for( i = na-2; i >= 0; i-- )
  {
    t = t*(T)(LL(0x100000000)) + (T)a[i];
  }  
}

/* ----------------------------------------------------------------
  Calculates the sum of two operands
------------------------------------------------------------------ */
void IntAdd( const unsigned int na, const unsigned long *a, 
             const unsigned int nb, const unsigned long *b, 
             unsigned int *nc, unsigned long *c );

/* ---------------------------------------------------------------
  Calculates c=a-b with length of c = max(length(a),length(b)).
  If a < b the result c is a negative number (two-complement, i.e.
  not(fabs(c))+1), and the fuction returns -1, else it returns 0.
------------------------------------------------------------------ */
int IntSub( const int unsigned na, const unsigned long *a, 
            const int unsigned nb, const unsigned long *b, 
            unsigned int *nc, unsigned long *c );
                 
/* ---------------------------------------------------------------
  Short Multiplication
---------------------------------------------------------------- */
unsigned long IntMultShort( 
         const unsigned int na, const unsigned long *a, const unsigned long b,
         unsigned long *c );
         
/* ----------------------------------------------------------------
  Calculates the product of two operands, slow but simple
  TODO: Multiplication with FFT
------------------------------------------------------------------ */
void IntMult( const unsigned int na, const unsigned long *a, 
              const unsigned int nb, const unsigned long *b,
              unsigned int *nc, unsigned long *c );

/* ---------------------------------------------------------
  Divide a int-string by a single int
----------------------------------------------------------- */
unsigned long IntDivShort( 
       const unsigned int na, const unsigned long *a, const unsigned long b, 
       unsigned int *nq, unsigned long *q );

/* ---------------------------------------------------------
  Divide a int-string by a int-string
----------------------------------------------------------- */
int IntDiv( const unsigned int na, const unsigned long *a, 
            const unsigned int nb, const unsigned long *b, 
            unsigned int *nq, unsigned long *q, 
            unsigned int *nr, unsigned long *r );

/* ---------------------------------------------------------
  Converts double to int-string
----------------------------------------------------------- */
int DoubleToInt( const double d, const unsigned int na, unsigned long *a,
                 unsigned int *retLen );

/* ---------------------------------------------------------
  Converts int-string to double
----------------------------------------------------------- */
double IntToDouble( const unsigned int na, const unsigned long *a );


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

struct Interval
{
  friend Interval operator+( const Interval &a, const Interval &b );
  friend Interval operator-( const Interval &a, const Interval &b );
  friend Interval operator*( const Interval &a, const Interval &b );
  friend Interval operator/( const Interval &a, const Interval &b );
  friend int Test( const Interval &a, const Interval &b );

  double         m_low;
  double         m_high;

  Interval() { m_low = 0.0; m_high = 0.0; }
  Interval(double val) { m_low = val; m_high = val; }
  Interval(double low, double high) { m_low = low; m_high = high; }
};

/* -------------------------------------------------------------
  wrapper class, because i don't want implicit double -> long
  conversion
---------------------------------------------------------------*/
class ULong
{
  unsigned long l;
public:
  explicit ULong( unsigned long ul ) : l(ul) { }
  operator unsigned long() const { return l; }
};

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

class rational
{
public:
  struct rational_rep
  {
    unsigned int   m_na;  /* number of words for a */
    unsigned long *m_a;   /* numerator   */

    unsigned int   m_nb;  /* number of words for b */
    unsigned long *m_b;   /* denominator */

    size_t         m_size_a;  /* allocation size in byte of a */ 
    size_t         m_size_b;  /* allocation size of byte of b */ 

    int            m_refcount;    /* Referencecounter */ 

    Interval       m_i;

    int            m_sign   : 2;  /* -1 = negative, 0 = positive or zero */
    unsigned int   m_bounds : 1;  /* lower and upper bound initialized */

    rational_rep()
    {
      m_refcount = 1;
      m_sign = 0;
      m_bounds = m_na = m_nb = 0;
      m_size_a = m_size_b = 0;
      m_a = m_b = NULL;
    }

    rational_rep( int sign, int na, unsigned long *a )
    {
      size_t bsize;

      m_refcount = 1;
      m_sign = 0;
      m_bounds = m_na = m_nb = 0;
      m_size_b = m_size_a = 0; 
      m_a = m_b = 0;

      if( na == 0 ) return;
      
      m_a = (unsigned long *)gust::PoolAlloc( na*sizeof(unsigned long), &bsize );
      if( m_a == NULL ) throw gul::PoolAllocError();
      for( int i=0; i<na; i++ ) m_a[i] = a[i];

      m_sign = sign;
      m_na = na;
      m_size_a = bsize;
    }
    
    rational_rep( int sign, unsigned long a )
    {
      size_t bsize;

      m_refcount = 1;
      m_sign = 0;
      m_bounds = m_na = m_nb = 0;
      m_size_a = m_size_b = 0;
      m_a = m_b = NULL;
      if( a == 0 ) return; 
     
      m_a = (unsigned long *)gust::PoolAlloc( sizeof(unsigned long), &bsize );
      if( m_a == NULL ) throw gul::PoolAllocError();
      m_a[0] = a;

      m_sign = sign;
      m_na = 1;
      m_size_a = bsize;
    }
 
            // this only reserves memory, but leaves it uninitialized
    rational_rep( int sign, int na, int nb )
    {
      size_t bsize;

      m_refcount = 1;
      m_sign = 0;
      m_bounds = m_na = m_nb = 0;
      m_size_a = m_size_b = 0; 
      m_a = m_b = NULL;

      if( na != 0 )       
      {
        m_sign = sign;
        m_a = (unsigned long *)gust::PoolAlloc( na*sizeof(unsigned long),&bsize);
        if( m_a == NULL ) throw gul::PoolAllocError();
        m_size_a = bsize;
        m_na = na;
      }
      if( nb != 0 )       
      {     
        m_b = (unsigned long *)gust::PoolAlloc( nb*sizeof(unsigned long),&bsize);
        if( m_b == NULL ) throw gul::PoolAllocError();
        m_size_b = bsize;
        m_nb = nb;
      }
    }

    ~rational_rep()
    {
      if( (m_size_a != 0) && (m_a != NULL) ) gust::PoolFree( m_a, m_size_a );      
      if( (m_size_b != 0) && (m_b != NULL) ) gust::PoolFree( m_b, m_size_b );
    }

    template< class T >
    void dump( T& num, T& den )
    { 
      if( m_na == 0 )
      {  num = (T)0.0; den = (T)1.0; return; }
      else
        IntTo<T>( m_na, m_a, num );

      if( m_sign ) 
        num = -num;

      if( m_nb == 0 )
        den = (T)1.0;
      else
        IntTo<T>( m_nb, m_b, den );  
    }

    void calc_bounds( void )
    { 
      Interval inum,iden;
      double num,den;

      dump( num, den );
      if( num != 0.0 )
      {
        inum.m_low = num - DBL_EPSILON*fabs(num);
        inum.m_high = num + DBL_EPSILON*fabs(num);
      }
      else
      {
        inum.m_low = inum.m_high = 0.0;
      }
      
      if( den != 1.0 )
      {
        iden.m_low = den - DBL_EPSILON*fabs(den);
        iden.m_high = den + DBL_EPSILON*fabs(den);
      }
      else
      {
        iden.m_low = iden.m_high = 1.0;
      }

      m_i = inum/iden;

      m_bounds = 1;
    }
   
    void negative( void )    /* for internal use only, no checks !!!! */ 
    {      
      if( m_na != 0 )
        m_sign ^= -1;
    }
    
    void reciprocal( void )  /* for internal use only, no checks !!!! */
    {
      size_t bsize;
      unsigned long *p;
      int i;
      
      if( m_na == 0 ) throw gul::Error();
      
      if( (m_na == 1) && (m_a[0] == 1) )
      {
        gust::PoolFree( m_a, m_size_a );
        m_na = 0;
        m_size_a = 0;
        m_a = NULL;
      }
      
      if( m_size_b == 0 )
      {
        m_b = (unsigned long *)gust::PoolAlloc( sizeof(unsigned long), &bsize );
        if( m_b == NULL ) throw gul::PoolAllocError();
        m_nb = 1;
        m_b[0] = 1;
        m_size_b = bsize;
      }

      i = m_na; m_na = m_nb; m_nb = i;
      bsize = m_size_a; m_size_a = m_size_b; m_size_b = bsize;
      p = m_a; m_a = m_b; m_b = p;
    }
  };

public:
  rational_rep   *m;
  
  GULAPI rational( void );
  GULAPI explicit rational( int na, unsigned long *a, int sign = 0 );
  GULAPI explicit rational( ULong a, int sign = 0 );
  GULAPI explicit rational( int sign, int size_a, int size_b );

  rational( const rational& other )
  {
    other.m->m_refcount++;
    m = other.m;
  }

  GULAPI rational& operator=( const rational& other );
  GULAPI ~rational();

  GULAPI rational get_copy() const;

  rational numerator() const
  { 
    if( m->m_na )
      return rational(m->m_na,m->m_a,m->m_sign);

    return rational();
  }
  rational denominator() const
  { 
    if( m->m_nb )
      return rational(m->m_nb,m->m_b,m->m_sign);

    return rational(ULong(1));
  }

  // GULAPI int to_int() const;
  GULAPI void calc_bounds( void ) const;
  GULAPI const Interval& get_bounds( void ) const;

  template< class T > GULAPI void dump( T& num, T& den ) const;
  template< class T > void dump( T& t ) const
  {
    T n,d;

    m->dump( n, d );
    
    t = n/d;
  } 
  double dump() const
  {
    double t;
    dump(t);
    return t;
  } 

  bool is_zero( void ) const { return(m->m_na == 0); }

  int  test( void ) const
  {
    if( m->m_sign ) return(-1);
    else if( m->m_na == 0 ) return(0);
    return(1);
  }

  GULAPI void div_mod( rational *q, rational *r ) const;

  GULAPI friend rational operator*( const rational& a, const rational& b );
  GULAPI friend rational operator+( const rational& A, const rational& B );
  GULAPI friend rational operator/( const rational& A, const rational& B );
  GULAPI friend rational operator-( const rational& A, const rational& B );
  // friend int gul::compare( const rational& A, const rational &B );
};

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

/* ---------------------------------------------------------------
  functions for intervall arithmetic
--------------------------------------------------------------- */

inline Interval operator+( const Interval &a, const Interval &b )
{
  Interval c;
  
  c.m_low = a.m_low + b.m_low;
  c.m_low = c.m_low - DBL_EPSILON*fabs(c.m_low);

  c.m_high = a.m_high + b.m_high;
  c.m_high = c.m_high + DBL_EPSILON*fabs(c.m_high);
    
  return( c );
}

inline Interval operator-( const Interval &a, const Interval &b )
{
  Interval c(-b.m_high,-b.m_low);
  
  return(a+c);
}

inline Interval operator/( const Interval &a, const Interval &b )
{
  if( (b.m_low < 0.0) && (b.m_high > 0.0) )
    throw gul::IntervalDivZeroError();

  Interval c(1.0/b.m_high,1.0/b.m_low);
  
  return(a*c);
}

// according to the Ansi C standard 'sqrt' has a rounding error less then 1 ULP
inline Interval Sqrt( const Interval &a )
{
  Interval c;
  
  c.m_low = sqrt(a.m_low);
  c.m_low = c.m_low - DBL_EPSILON*fabs(c.m_low);

  c.m_high = sqrt(a.m_high);
  c.m_high = c.m_high + DBL_EPSILON*fabs(c.m_high);
    
  return( c );
}

inline int Test( const Interval &a )
{
  if( a.m_high < 0.0 ) return(-1);
  if( a.m_low > 0.0 ) return(+1);

  return(0);   /* dunno :) */
}

inline int Test( const Interval &a, const Interval &b )
{
  Interval c = a-b;
  return Test(c);
}

inline double GetRatio( const Interval &a )
{
  double h,l;

  h = a.m_high;
  l = a.m_low;

  if( l < 0.0 )
  {
    if( h >= 0.0 ) return gul::rtr<double>::giant(); // relative error = inf

    h = -l;
    l = -h;
  }
   
  if( l == 0.0 )
  {
    if( h == 0.0 ) return 1.0;
    return gul::rtr<double>::giant();
  }

  return h/l;
}

/* ---------------------------------------------------------------
  functions for fractions arithmetic
--------------------------------------------------------------- */

/*
inline bool operator<( const rational& A, int b )
{
  rational B,C;
  unsigned long buf;
  
  if( b < 0 )
  {
    B.m->m_sign = -1;
    buf = -b;
    B.m->m_a = &buf;
    B.m->m_na = 1; 
  }
  else if( b > 0 )
  {
    buf = b;
    B.m->m_a = &buf;
    B.m->m_na = 1; 
  }
  
  C = A - B;

  return( C.m->m_sign != 0 );
}

inline bool operator<=( const rational& A, int b )
{
  rational B,C;
  unsigned long buf;
  
  if( b < 0 )
  {
    B.m->m_sign = -1;
    buf = -b;
    B.m->m_a = &buf;
    B.m->m_na = 1; 
  }
  else if( b > 0 )
  {
    buf = b;
    B.m->m_a = &buf;
    B.m->m_na = 1; 
  }
  
  C = A - B;

  if( C.m->m_sign || C.m->m_na == 0)
    return(true);

  return( false );
}

inline bool operator>=( const rational& A, int b )
{
  return( !(A < b) );
}
*/

/*
class point2<rational>
{
public:
  rational m_x;
  rational m_y;

  point2<rational>() { }
  ~point2<rational>() { }
  explicit point2<rational>( const point2<rational>& o ) : m_x(o.m_x), m_y(o.m_y) { }
  explicit point2<rational>( ULong x, ULong y ) : m_x(x), m_y(y) { }
  explicit point2<rational>( const rational& x, const rational& y ) : 
                                                        m_x(x), m_y(y) { }
};

class PrecPoint3D
{
public:
  rational m_x;
  rational m_y;
  rational m_z;

  PrecPoint3D() { }
  ~PrecPoint3D() { }
  PrecPoint3D( const PrecPoint3D& o ) 
                    : m_x(o.m_x), m_y(o.m_y), m_z(o.m_z) { }
  PrecPoint3D( ULong x, ULong y, ULong z ) 
                    : m_x(x), m_y(y), m_z(z) { }
  PrecPoint3D( 
     const rational& x, const rational& y, const rational& z ) 
                    : m_x(x), m_y(y), m_z(z) { }

  PrecPoint3D& operator=( const PrecPoint3D& o ) 
                  { m_x = o.m_x; m_y = o.m_y; m_z = o.m_z; return *this; }
};
*/

}

namespace gul {

template<> struct rtr< guar::rational >
{
  // these constants are needed in some templates, which i also want to 
  // instantiate for the rational number class, and  i do not want to write
  //  an operator which converts floats or ints to rationals (this is too
  // dangerous, because of the many (unwanted) implicit conversions, 
  // which c++ does)
  static guar::rational zero() { return guar::rational(); }
  static guar::rational one() { return guar::rational(guar::ULong(1)); }
};

template<>
inline int test( const guar::rational& a )
{
  return a.test();
}

#ifdef GUL_VECTOR_H
template<> inline point<guar::rational> cross_product( 
      const point2<guar::rational>& a, const point2<guar::rational>& b )
{
  point<guar::rational> c;
  c.z = a.x*b.y-b.x*a.y;
  return c;
}
#endif

template<>
inline int compare( const guar::rational& A, const guar::rational &B )
{
  guar::Interval i;
  
  if( A.m == B.m )
    return 0;

  if( !A.m->m_bounds ) A.calc_bounds();
  if( !B.m->m_bounds ) B.calc_bounds();  
  i = A.m->m_i - B.m->m_i;
  if( i.m_high < 0.0 ) return(-1); 
  else if( i.m_low > 0.0 ) return(1); 

  guar::rational c = A - B;
  if( c.m->m_sign ) return(-1);
  if( c.m->m_na == 0 ) return(0);

  return(1);
}

inline std::ostream& operator<<( std::ostream& s, const guar::rational& r )
{
  return s << r.dump();
} 

}

#endif

---