---

gul_types.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 GUL_TYPES_H
#define GUL_TYPES_H


#ifdef __GNUC__

  /*---------------------------------------------------------------------------
  BTW,  "gcc -dM -E - < /dev/null" shows symbols automatically defined by gcc
  ----------------------------------------------------------------------------*/

  #ifndef alloca
    #define alloca __builtin_alloca
  #endif

#endif

/***********************************************************************
  Configuration Options
***********************************************************************/
// define this only when doing memory leak debugging, very slow !
// #define POOLDBG
// define this only when debugging
// #define RANGECHECK

// at least sparc ultras require an alignment of 8 bytes, or
// else writing doubles causes a bus error. in any case,
// the alignment must be >= sizeof(void*)
#define ALIGNMENT 8
#define POOL_ALIGNMENT (ALIGNMENT > sizeof(void*) ? ALIGNMENT : sizeof(void*))

#ifdef KAI_BSD_HACK
  #include <sys/cdefs.h>
  #include <machine/ansi.h>
  #undef _BSD_WCHAR_T_
#endif

#include <cstddef>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <float.h>
#include <math.h>

#ifdef RANGECHECK
  #include <iostream>
#endif

#ifdef _MSC_VER
  #include <malloc.h>
  #ifdef USEDLL
    #ifdef CREATE_GUL_DLL
      #define GULAPI __declspec(dllexport)
    #else
      #define GULAPI __declspec(dllimport)
    #endif
  #else
    #define GULAPI
  #endif
#else
  #define GULAPI
#endif

#ifdef _MSC_VER
  GULAPI void *test_gul_dll_malloc( size_t s );
  GULAPI void test_gul_dll_free( void *s );

  // a normal function, not a template, onto which can be set a breakpoint
  GULAPI void gul_halt( void );
#endif

#if (defined(__MINGW32__)) || (defined(_MSC_VER))
  // still have to define this somewhere :((
  GULAPI double scalbn(double,int);
  GULAPI int ilogb(double);
  GULAPI long random(void);
  GULAPI void srandom(long);
#endif

#include "gul_error.h"
#include "gust_pool.h"
#include "gust_new.h"

// #if !defined(_MSC_VER) || defined(NEED_NEW_H)
//  #include <new.h>
// #endif

namespace gul
{
#ifndef NDEBUG
  const bool ndebug=false;
#else
  const bool ndebug=true;
#endif
 
#ifdef _MSC_VER
  typedef unsigned _int64 uint64;
  typedef _int64 int64;
  #define isinf(f) (!_finite((f)))
  // 64-Bit Integer constants
  #define LL(i) i##i64
#else
#ifdef __BORLANDC__
  typedef unsigned __int64 uint64;
  typedef __int64 int64;
  #define isinf(f) (!_finite((f)))
  // 64-Bit Integer constants
  #define LL(i) i##i64
#else
  typedef unsigned long long uint64;
  typedef long long int64;
  #ifdef __MINGW32__
//    #include <float.h>
//    #define isinf(f) (!_finite((f)))
  #endif
  // 64-Bit Integer constants
  #define LL(i) i##LL
#endif
#endif

typedef unsigned long uint32;
typedef long int32;

template< class T >
inline int compare( const T& a, const T& b )
{
  if( a < b ) 
    return -1;
  else if( a == b )
    return 0;
  return 1;
} 
template< class T >
inline int test( const T& a )
{
  if( a < 0 ) return -1;
  else if( a == 0 ) return 0;
  return 1;
} 

template< class T > struct rtr  /* rational number traits */
{ 
};     

template<> struct rtr< float >
{
  static float epsilon() { return FLT_EPSILON; }
  static float epsilon_inv() { return 1.0f/FLT_EPSILON; }
  static float zero_tol() { return 1.e-4f; }
  static size_t mantissa_length() { return 4; } // in multiples of 4
  static float maximum() { return FLT_MAX; }
  static float minimum() { return FLT_MIN; }
  static float pi() { return 3.14159265358979323846f; }
  static float pi_180() { return 3.14159265358979323846f/180.0f; }
  static float pi_180_inv() { return 180.0f/3.14159265358979323846f; }
  static float root2_2() { return 1.41421356237309514547f; }
  static float root2_3() { return 1.7320508076887719318f; }
  static float golden_r() { return 0.61803399f; }
  static float golden_c() { return 1.0f-golden_r(); }
  // i never understood for what this 'tiny' is needed in "Numerical Recipes",
  // but anyway here it is
  static float tiny() { return 1e-20f; }
  static float giant() { return 1.0f/tiny(); }
  // 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 float zero() { return 0.0f; }
  static float one() { return 1.0f; }

  static float  **m_BinCoeff;
  static int      m_BinCoeff_Pmax;
  static void     InitBinCoeff( const int Pmax );
  static void     ExitBinCoeff();
  static float BinCoeff( const int p, const int k )
  {  
    if( (k > p) || (k < 0) || (p < 0) ) return 0.0;
    if( p > m_BinCoeff_Pmax ) InitBinCoeff( p ); 
    return m_BinCoeff[p][k];
  }  
 #if (!defined(__MINGW32__)) && (!defined(__BORLANDC__)) && (!defined(sun)) && \
     (!defined(_MSC_VER))
  static float floor( const float a ) { return ::floorf(a); }
  static float ceil( const float a ) { return ::ceilf(a); }
  static float fabs( const float a ) { return ::fabsf(a); }
  static float sqrt( const float a ) { return ::sqrtf(a); }
  static float sin( const float a ) { return ::sinf(a); }
  static float cos( const float a ) { return ::cosf(a); }
  static float acos( const float a ) { return ::acosf(a); }
  // IEEE functions
  static int ilogb( const float a ) { return ::ilogbf(a); }
  static float scalbn( const float x, int n ) { return ::scalbnf(x,n); } 
 #else
  // windows runtime has no float functions !!!
  static float floor( const float a ) { return (float)::floor((double)a); }
  static float ceil( const float a ) { return (float)::ceil((double)a); }
  static float fabs( const float a ) { return (float)::fabs((double)a); }
  static float sqrt( const float a ) { return (float)::sqrt((double)a); }
  static float sin( const float a ) { return (float)::sin((double)a); }
  static float cos( const float a ) { return (float)::cos((double)a); }
  static float acos( const float a ) { return (float)::acos((double)a); }
  // IEEE functions
  static int ilogb( const float a ) { return ::ilogb((double)a); }
  static float scalbn(const float x, int n) 
                                    { return (float)::scalbn((double)x,n); } 
 #endif
  static float rad( const float d ) { return d * pi_180(); }
  static float deg( const float r ) { return r * pi_180_inv(); }
  static int cmp( const float a, const float b ) { return compare<float>(a,b);}

  static int id() { return(1); }
};

template<> struct rtr< double >
{
  static double epsilon() { return DBL_EPSILON; }
  static double epsilon_inv() { return 1.0/DBL_EPSILON; }
  static double zero_tol() { return 1.0e-9; }
  static size_t mantissa_length() { return 8; } // in multiples of 4
  static double maximum() { return DBL_MAX; }
  static double minimum() { return DBL_MIN; }
  static double pi() { return 3.14159265358979323846; }
  static double pi_180() { return 3.14159265358979323846/180.0; }
  static double pi_180_inv() { return 180.0/3.14159265358979323846; }
  static double root2_2() { return 1.41421356237309514547; }
  static double root2_3() { return 1.7320508076887719318; }
  static double golden_r() { return 0.61803399; }
  static double golden_c() { return 1.0-golden_r(); }
  // i never understood for what this 'tiny' is needed in "Numerical Recipes",
  // but anyway here it is
  static double tiny() { return 1e-40; }
  static double giant() { return 1.0f/tiny(); }
  // 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 float zero() { return 0.0; }
  static float one() { return 1.0; }

  static double  **m_BinCoeff;
  static int       m_BinCoeff_Pmax;
  static void      InitBinCoeff( const int Pmax );
  static void      ExitBinCoeff();
  static double BinCoeff( const int p, const int k )
  {  
    if( (k > p) || (k < 0) || (p < 0) ) return 0.0f;
    if( p > m_BinCoeff_Pmax ) InitBinCoeff( p ); 
    return m_BinCoeff[p][k];
  }  
  static double floor( const double a ) { return ::floor(a); }
  static double ceil( const double a ) { return ::ceil(a); }
  static double fabs( const double a ) { return ::fabs(a); }
  static double sqrt( const double a ) { return ::sqrt(a); }
  static double sin( const double a ) { return ::sin(a); }
  static double cos( const double a ) { return ::cos(a); }
  static double acos( const double a ) { return ::acos(a); }
  // IEEE functions
  static int ilogb( const double a ) { return ::ilogb(a); }
  static double scalbn( const double x, int n ) { return ::scalbn(x,n); } 

  static double rad( const double d ) { return d * pi_180(); }
  static double deg( const double r ) { return r * pi_180_inv(); }
  static int cmp( const double a, const double b ){return compare<double>(a,b);}

  static int id() { return(2); }
};

template<> struct rtr< int >
{
  static int id() { return(3); }
};

template<> struct rtr< unsigned int >
{
  static int id() { return(4); }
};


/*
// this tolerance structure is not longer needed, since from now on only
// the relative error is measured, when checking equality

template< class T > struct rtl  // tolerances for comparisons - struct
{ 
  T zero_tol;
  T zero_tol_v2;
  T zero_tol_v3;
  T zero_tol_v4;
  T point_coincidence_tol;
  T knot_coincidence_tol;
  T knot_uniqueness_tol;
    
  rtl() 
  {
    // this tolerances only make sense when the values are normalized
    zero_tol = rtr<T>::epsilon() * (T)16;
    zero_tol_v2 = zero_tol / rtr<T>::root2_2();
    zero_tol_v3 = zero_tol / rtr<T>::root2_3();
    zero_tol_v4 = zero_tol / (T)2;
    point_coincidence_tol = rtr<T>::epsilon() * (T)16;
    knot_coincidence_tol = rtr<T>::epsilon() * (T)16;
    knot_uniqueness_tol = rtr<T>::epsilon() * (T)16; 
  }

  rtl( T zt, T pct, T kct, T kut )
  {
    zero_tol = zt;
    zero_tol_v2 = zero_tol / rtr<T>::root2_2();
    zero_tol_v3 = zero_tol / rtr<T>::root2_3();
    zero_tol_v4 = zero_tol / (T)2;
    point_coincidence_tol = pct;
    knot_coincidence_tol = kct;
    knot_uniqueness_tol = kut;
  }

  rtl( const rtl& a )
  {
    zero_tol = a.zero_tol;
    zero_tol_v2 = a.zero_tol_v2;
    zero_tol_v3 = a.zero_tol_v3;
    zero_tol_v4 = a.zero_tol_v4;
    point_coincidence_tol = a.point_coincidence_tol;
    knot_coincidence_tol = a.knot_coincidence_tol;
    knot_uniqueness_tol = a.knot_uniqueness_tol;
  }
};     
*/

template< class T >
inline T Max( const T& a, const T& b ) { return a > b ? a : b; } 
template< class T >
inline T Min( const T& a, const T& b ) { return a < b ? a : b; } 

template< class T >
inline T Max( const T& a, const T& b, const T& c ) 
{ if((a > b) && (a > c)) return a; 
  else if(b > c) return b;
  return c; 
} 
template< class T >
inline T Min( const T& a, const T& b, const T& c ) 
{ if((a < b) && (a < c)) return a; 
  else if(b < c) return b;
  return c; 
} 
 
template< class T >
inline void Swap( T& a, T& b ) { T tmp=a; a=b; b=tmp; }
template< class T >
inline T Square( const T& a ) { return a*a; }
                                            // give 'a' the sign of 'b':
template< class T > inline T Sign( const T& a, const T& b ) 
{ return (b < 0.0 ? -rtr<T>::fabs(a) : rtr<T>::fabs(a)); }



/*-----------------------------------------------------------
  Union for some dirty tricks
-----------------------------------------------------------*/

union ptr_int_union
{
  void *p;
  int   i;
  void  set_i( int ii ) { p = 0; i = ii; } 
};

/* -------------------------------------------------------------------------
  Pointer with reference counter semantic
-------------------------------------------------------------------------- */

template < class T > 
class Ptr
{ 
  template< class M >
  struct Ptr_rep
  {
    T           *m_ptr;
    unsigned int m_refs;
    unsigned int m_type : 2;     /* storage class */
    unsigned int m_size : 30;
    unsigned int m_nelems;

    Ptr_rep( T *p, unsigned int i, unsigned int t, size_t s, unsigned int n) :
                     m_ptr(p), m_refs(i), m_type(t), m_nelems(n) 
    { 
      #ifdef RANGECHECK
        if( s > 0x3fffffff )
        {
          std::cout << "Ptr: array is too big\n";
          std::cout.flush();
          gul_halt();
        }
      #endif

      m_size = (unsigned int)s;
    }

    inline ~Ptr_rep()
    { 
      T *p = m_ptr;

      if( m_type != gust::dumb )
      {
        for( unsigned int j = 0; j < m_nelems; j++ )
        {
          p->~T();
          p++;
        }
      }
      gust::release((void *)m_ptr, m_type, m_size);
    }

    // this only decreases the counter for constructed elements
    void shrink( unsigned int n )
    {
      if( m_type != gust::dumb )
        for( int i = n; i < m_nelems; i++ ) m_ptr[i].~T();
      m_nelems = n;
    }

    // only for objects without constructors/destructors !!!!
    void realloc(unsigned int n) 
    { 
      if( m_type == heap )
      {
        m_ptr = (T *)::realloc( m_ptr, sizeof(T)*n );
        m_size = sizeof(T)*n;
        m_nelems = n;
      }
      else if( m_size < sizeof(T)*n )
        throw InternalError();  /* something went wrong */
    }
  };

  Ptr_rep<T>  *m_rep;

public:
  explicit Ptr( T *p = 0, unsigned int t = 0, size_t s = 0, 
                unsigned int n = 0 )
  { 
    size_t dummy;

    if( p != 0 )    
    {
      m_rep = (Ptr_rep<T> *)gust::PoolAlloc( sizeof(Ptr_rep<T>), &dummy );
      if( m_rep == 0 ) throw gul::AllocError();

      ::new(m_rep) Ptr_rep<T>(p,1,t,s,n);
    }
    else
        m_rep = 0;
  }

  Ptr( const Ptr& other )
  {
    if( other.m_rep ) other.m_rep->m_refs++;
    m_rep = other.m_rep;
  }

  Ptr& operator=( const Ptr& other )
  {
    if( other.m_rep ) other.m_rep->m_refs++;

    if( (m_rep != 0) && (--m_rep->m_refs == 0) ) 
    {
      m_rep->~Ptr_rep<T>();
      gust::PoolFree( m_rep, sizeof(Ptr_rep<T>) );      
    }

    m_rep = other.m_rep;

    return( *this );    
  } 

  ~Ptr() 
  {
    if( (m_rep != 0) && (--m_rep->m_refs == 0) ) 
    {
      m_rep->~Ptr_rep<T>();
      gust::PoolFree( m_rep, sizeof(Ptr_rep<T>) );      
    }
  }

  // gets the number of constructed elements in the area to which it points
  int nElems() const  
  {
    if( m_rep != 0 ) 
      return m_rep->m_nelems; 
    else 
      return 0;
  }
  const Ptr_rep<T> *rep( void ) const { return m_rep; }

  void reset( T *p, unsigned int r, unsigned int t, size_t s, 
              unsigned int n )
  {
    size_t dummy;
    
    if( m_rep == 0 )
    {
      m_rep = (Ptr_rep<T> *)gust::PoolAlloc( sizeof(Ptr_rep<T>), &dummy );
      if( m_rep == 0 ) throw gul::AllocError();
    }
    else
    {
      if( --(m_rep->m_refs) == 0 )
      {
        m_rep->~Ptr_rep<T>();
      }
      else
      {
        m_rep = (Ptr_rep<T> *)gust::PoolAlloc( sizeof(Ptr_rep<T>), &dummy );
        if( m_rep == 0 ) throw gul::AllocError();
      }
    }
    ::new(m_rep) Ptr_rep<T>(p,r,t,s,n);   
  }

  // Caution when using this ! The source p must not be changed in any way
  // as long the pointer uses it
  void use_pointer( T *p, unsigned int n )
  {
    size_t dummy;
    
    if( m_rep == 0 )
    {
      m_rep = (Ptr_rep<T> *)gust::PoolAlloc( sizeof(Ptr_rep<T>), &dummy );
      if( m_rep == 0 ) throw gul::AllocError();
    }
    else
    {
      if( --(m_rep->m_refs) == 0 )
      {
        m_rep->~Ptr_rep<T>();
      }
      else
      {
        m_rep = (Ptr_rep<T> *)gust::PoolAlloc( sizeof(Ptr_rep<T>), &dummy );
        if( m_rep == 0 ) throw gul::AllocError();
      }
    }
    ::new(m_rep) Ptr_rep<T>(p,1,gust::dumb,0,n);   
  }

  void reserve_global( unsigned int n )
  {
    size_t s;
    T *p = (T *)gust::reserve_global( n*sizeof(T), &s );
    
    reset( p, 1, gust::heap, s, n );
    for( size_t i = 0; i < n; i++ ) ::new(p++) T();  /* call constructors */
  }

  void reserve_pool( unsigned int n )
  {
    size_t s;
    T *p = (T *)gust::reserve_pool( n*sizeof(T), &s );
    
    reset( p, 1, gust::pool, s, n );
    for( size_t i = 0; i < n; i++ ) ::new(p++) T();  /* call constructors */
  }

  void reserve_place( void *buf, unsigned int n )
  {
    T *p = (T *)buf;
    reset( p, 1, gust::place, n*sizeof(T), n );      /* call constructors */
    for( size_t i = 0; i < n; i++ ) ::new(p++) T();
  }


  void reset()       //! delete all elements
  {
    if( m_rep == 0 ) return;

    if( --(m_rep->m_refs) == 0 )
    {
      m_rep->~Ptr_rep<T>();
      gust::PoolFree( m_rep, sizeof(Ptr_rep<T>) ); 
    }
    m_rep = 0;
  }
  // this only decreases the counter for constructed elements
  void shrink(unsigned int n)
  {
    if( m_rep == 0 ) throw gul::Error();
    m_rep->shrink(n);
  }
  T& operator*() const 
  { 
    #ifdef RANGECHECK
      if( nElems() <= 0 )
      {
        cout << "Ptr: pointer not initialized\n";
        cout.flush();
        gul_halt();
      }
    #endif
    return *(m_rep->m_ptr);
  }
  T& operator[](int i) const 
  {
    #ifdef RANGECHECK
      if( (i >= nElems()) || (i < 0) )
      {
        std::cout << "Ptr: accessing elements behind array bounds\n";
        std::cout.flush();
        gul_halt();
      }
      if( m_rep && m_rep->m_size )
      {
        if( (size_t)i*sizeof(T) >= m_rep->m_size )
        {
          std::cout << "Ptr: size and element count don't match\n";
          std::cout.flush();
          gul_halt();
        }
      }
    #endif
    return (m_rep->m_ptr)[i];
  }

  //  T* operator->() const { return m_rep->m_ptr; }
  T* get() const { return m_rep->m_ptr; }
  //  operator T*() const { return m_rep->m_ptr; }
 
  T* release()
  {
    if( (m_rep == 0) || (m_rep->m_refs != 1) )
      throw Error();
    T *p = m_rep->m_ptr;
    gust::PoolFree( m_rep, sizeof(Ptr_rep<T>) ); 
    m_rep = 0;
    return(p);
  }

  void realloc(unsigned int n) const  // low level, not for classes !
  { if(m_rep != 0) m_rep->realloc(n); else throw InternalError();  }
};

enum DomainDirection
{
  u_direction,
  v_direction
};

typedef int sub_range[2];


template< class T > class point;
template< class T > class point2;
template< class T > class point1;
template< class T > class hpoint;
template< class T > class hpoint2;
template< class T > class hpoint1;


template< class T >
class vec4
{
public:
  T m[4];
  T& operator[] (size_t i) { return m[i]; }

  /*
  template< class T2 > inline vec4<T>& operator=( const point<T2>& a );
  template< class T2 > inline vec4<T>& operator=( const point2<T2>& a );
  template< class T2 > inline vec4<T>& operator=( const point1<T2>& a );

  template< class T2 > inline vec4<T>& operator=( const hpoint<T2>& a );
  template< class T2 > inline vec4<T>& operator=( const hpoint2<T2>& a );
  template< class T2 > inline vec4<T>& operator=( const hpoint1<T2>& a );
  */
};

template< class T >
struct mat4x4
{
  vec4<T> m[4];    
  vec4<T>& operator[] (size_t i) { return m[i]; }
  /*---------------------------------------------------------------*//**
    creates 4x4 unit matrix                                           */
  /*------------------------------------------------------------------*/
  static mat4x4<T> identity()
  {
    mat4x4 a = { {{{(T)1,(T)0,(T)0,(T)0}},
                  {{(T)0,(T)1,(T)0,(T)0}},
                  {{(T)0,(T)0,(T)1,(T)0}},
                  {{(T)0,(T)0,(T)0,(T)1}}} };
    return a;
  }
  /*-----------------------------------------------------------------*//**
    create 4x4 rotation matrix about x axis.
    the angle theta (in radians) describes the clockwise rotation, if you
    look along the axis towards the origin                              */
  /*--------------------------------------------------------------------*/
  static mat4x4<T> rotate_x( T theta )
  {
    mat4x4 a = identity();
    T ct = rtr<T>::cos(theta);
    T st = rtr<T>::sin(theta);
    a[1][1] = a[2][2] = ct;
    a[1][2] = -st;
    a[2][1] = st;
    return a;
  }
  /*-----------------------------------------------------------------*//**
    create 4x4 rotation matrix about y axis.
    the angle theta (in radians) describes the clockwise rotation, if you
    look along the axis towards the origin                              */
  /*--------------------------------------------------------------------*/
  static mat4x4<T> rotate_y( T theta )
  {
    mat4x4 a = identity();
    T ct = rtr<T>::cos(theta);
    T st = rtr<T>::sin(theta);
    a[0][0] = a[2][2] = ct;
    a[0][2] = st;
    a[2][0] = -st;
    return a;
  }
  /*-----------------------------------------------------------------*//**
    create 4x4 rotation matrix about z axis.
    the angle theta (in radians) describes the clockwise rotation, if you
    look along the axis towards the origin                              */
  /*--------------------------------------------------------------------*/
  static mat4x4<T> rotate_z( T theta )
  {
    mat4x4 a = identity();
    T ct = rtr<T>::cos(theta);
    T st = rtr<T>::sin(theta);
    a[0][0] = a[1][1] = ct;
    a[0][1] = -st;
    a[1][0] = st;
    return a;
  }
};
  
template< class T >
struct bounding_box
{
  T x1;
  T x2;
  T y1;
  T y2;
  T z1;
  T z2;
};

/*------------------------------------------------------------------
   point types
-------------------------------------------------------------------*/

template< class T >
class point1
{
public:
  T x;
  static int dim() { return 1; }
  static bool hom() { return false; }
  T weight() const { return (T)1; }
  point1<T>() { }
  template< class T2 >
  point1<T>( const T2& ax )
  {
    x = (T)ax;
  }  
  inline point1<T>& operator+= ( const point1<T>& a )
  {
    x += a.x;
    return *this;
  }
  inline point1<T>& operator-= ( const point1<T>& a )
  {
    x -= a.x;
    return *this;
  }
  inline point1<T>& operator*= ( T& a )
  {
    x *= a;
    return *this;
  }
  const T& operator[] ( size_t i ) const { return (&x)[i]; }

  operator vec4<T>() const
  {
    vec4<T> v;
    v[0] = x; v[1] = 0; v[2] = 0; v[3] = 1;
    return v;
  }
};
template< class T >
class point2
{
public:
  T x;
  T y;
  static int dim() { return 2; }
  static bool hom() { return false; }
  T weight() const { return rtr<T>::one(); }
  point2<T>() { }
  template< class T2 >
  point2<T>( const T2& ax, const T2& ay )
  {
    x = (T)ax;
    y = (T)ay;
  }
  template< class T2 >
  point2<T>& operator=( const point2<T2>& a )
  { 
    x = (T)a.x;
    y = (T)a.y; 
    return *this;
  }
  inline point2<T>& operator+= ( const point2<T>& a )
  {
    x += a.x; y += a.y;
    return *this;
  }
  inline point2<T>& operator-= ( const point2<T>& a )
  {
    x -= a.x; y -= a.y;
    return *this;
  }
  inline point2<T>& operator*= ( T& a )
  {
    x *= a; y *= a;
    return *this;
  }
  const T& operator[] ( size_t i ) const { return (&x)[i]; }

  operator vec4<T>() const
  {
    vec4<T> v;
    v[0] = x; v[1] = y; v[2] = 0; v[3] = 1;
    return v;
  }
};
template< class T >
class point
{
public:
  T x;
  T y;
  T z;
  static int dim() { return 3; }
  static bool hom() { return false; }
  T weight() const { return rtr<T>::one(); }
  point<T>() { }
  template< class T2 >
  point<T>( const T2& ax, const T2& ay, const T2& az )
  {
    x = (T)ax;
    y = (T)ay; 
    z = (T)az; 
  }
  template< class T2 >
  point<T>& operator=( const point2<T2>& a )
  { 
    x = (T)a.x;
    y = (T)a.y; 
    z = (T)0; 
    return *this;
  }
  template< class T2 >
  point<T>& operator=( const point<T2>& a )
  { 
    x = (T)a.x;
    y = (T)a.y; 
    z = (T)a.z; 
    return *this;
  }
  inline point<T>& operator+= ( const point<T>& a )
  {
    x += a.x; y += a.y; z += a.z;
    return *this;
  }
  inline point<T>& operator-= ( const point<T>& a )
  {
    x -= a.x; y -= a.y; z -= a.z;
    return *this;
  }
  inline point<T>& operator*= ( T& a )
  {
    x *= a; y *= a; z *= a;
    return *this;
  }
  const T& operator[] ( size_t i ) const { return (&x)[i]; }

  operator vec4<T>() const
  {
    vec4<T> v;
    v[0] = x; v[1] = y; v[2] = z; v[3] = 1;
    return v;
  }
};

template< class T >
class hpoint1
{
public:
  T x;
  T w;
  static int dim() { return 2; }
  static bool hom() { return true; }
  T weight() const { return w; }
  hpoint1<T>() { }
  template< class T2 >
  hpoint1<T>( const T2& ax, const T2& aw )
  {
    x = (T)ax;
    w = (T)aw; 
  }
  inline hpoint1<T>& operator+= ( const hpoint1<T>& a )
  {
    x += a.x; w += a.w;
    return *this;
  }
  inline hpoint1<T>& operator-= ( const hpoint1<T>& a )
  {
    x -= a.x; w -= a.w;
    return *this;
  }
  inline hpoint1<T>& operator*= ( T& a )
  {
    x *= a; w *= a;
    return *this;
  }
  const T& operator[] ( size_t i ) const { return (&x)[i]; }

  operator vec4<T>() const
  {
    vec4<T> v;
    v[0] = x; v[1] = 0; v[2] = 0; v[3] = w;
    return v;
  }
};
template< class T >
class hpoint2
{
public:
  T x;
  T y;
  T w;
  static int dim() { return 3; }
  static bool hom() { return true; }
  T weight() const { return w; }
  hpoint2<T>() { }
  template< class T2 >
  hpoint2<T>( const T2& ax, const T2& ay, const T2& aw )
  {
    x = (T)ax;
    y = (T)ay; 
    w = (T)aw; 
  }
  inline hpoint2<T>& operator+= ( const hpoint2<T>& a )
  {
    x += a.x; y += a.y; w += a.w;
    return *this;
  }
  inline hpoint2<T>& operator-= ( const hpoint2<T>& a )
  {
    x -= a.x; y -= a.y; w -= a.w;
    return *this;
  }
  inline hpoint2<T>& operator*= ( T& a )
  {
    x *= a; y *= a; w *= a;
    return *this;
  }
  const T& operator[] ( size_t i ) const { return (&x)[i]; }

  operator vec4<T>() const
  {
    vec4<T> v;
    v[0] = x; v[1] = y; v[2] = 0; v[3] = w;
    return v;
  }
};
template< class T >
class hpoint
{
public:
  T x;
  T y;
  T z;
  T w;
  static int dim() { return 4; }
  static bool hom() { return true; }
  T weight() const { return w; }
  hpoint<T>() { }
  template< class T2 >
  hpoint<T>( const T2& ax, const T2& ay, const T2& az, const T2& aw )
  {
    x = (T)ax;
    y = (T)ay; 
    z = (T)az; 
    w = (T)aw; 
  }
  template< class T2 >
  hpoint<T>& operator=( const point<T2>& a )
  { 
    x = (T)a.x;
    y = (T)a.y; 
    z = (T)a.z; 
    w = (T)1; 
    return *this; 
  }
  template< class T2 >
  hpoint<T>& operator=( const point2<T2>& a )
  { 
    x = (T)a.x;
    y = (T)a.y; 
    z = (T)0; 
    w = (T)1; 
    return *this; 
  }
  template< class T2 >
  hpoint<T>& operator=( const hpoint2<T2>& a )
  { 
    x = (T)a.x;
    y = (T)a.y; 
    z = (T)0; 
    w = (T)a.w; 
    return *this;
  }
  inline hpoint<T>& operator+= ( const hpoint<T>& a )
  {
    x += a.x; y += a.y; z += a.z; w += a.w;
    return *this;
  }
  inline hpoint<T>& operator-= ( const hpoint<T>& a )
  {
    x -= a.x; y -= a.y; z -= a.z; w -= a.w;
    return *this;
  }
  inline hpoint<T>& operator*= ( T& a )
  {
    x *= a; y *= a; z *= a; w *= a;
    return *this;
  }
  const T& operator[] ( size_t i ) const { return (&x)[i]; }

  operator vec4<T>() const
  {
    vec4<T> v;
    v[0] = x; v[1] = y; v[2] = z; v[3] = w;
    return v;
  }
};

/*
template< class T > template< class T2 >
inline vec4<T>& vec4<T>::operator=( const point<T2>& a )
{
  m[0] = (T)a.x; m[1] = (T)a.y; m[2] = (T)a.z; m[3] = (T)1;
  return *this;
}
template< class T > template< class T2 >
inline vec4<T>& vec4<T>::operator=( const point2<T2>& a )
{
  m[0] = (T)a.x; m[1] = (T)a.y; m[2] = (T)0; m[3] = (T)1;
  return *this;
}
template< class T > template< class T2 >
inline vec4<T>& vec4<T>::operator=( const point1<T2>& a )
{
  m[0] = (T)a.x; m[1] = (T)0; m[2] = (T)0; m[3] = (T)1;
  return *this;
}

template< class T > template< class T2 >
inline vec4<T>& vec4<T>::operator=( const hpoint<T2>& a )
{
  m[0] = (T)a.x; m[1] = (T)a.y; m[2] = (T)a.z; m[3] = (T)a.w;
  return *this;
}
template< class T > template< class T2 >
inline vec4<T>& vec4<T>::operator=( const hpoint2<T2>& a )
{
  m[0] = (T)a.x; m[1] = (T)a.y; m[2] = (T)0; m[3] = (T)a.w;
  return *this;
}
template< class T > template< class T2 >
inline vec4<T>& vec4<T>::operator=( const hpoint1<T2>& a )
{
  m[0] = (T)a.x; m[1] = (T)0; m[2] = (T)0; m[3] = (T)a.w;
  return *this;
}
*/

template< class T >
struct line
{
  point<T> P1;
  point<T> P2;
  const point<T>& operator[] ( size_t i ) const { return (&P1)[i]; }
};
template< class T >
struct line2
{
  point2<T> P1;
  point2<T> P2;
  const point2<T>& operator[] ( size_t i ) const { return (&P1)[i]; }
};

template< class T >
struct triangle
{
  point<T> P1;
  point<T> P2;
  point<T> P3;
};
template< class T >
struct triangle2
{
  point2<T> P1;
  point2<T> P2;
  point2<T> P3;
};
struct itriangle
{
  int v[3]; 
};

template< class EP >
struct polyline
{
  int      n;
  Ptr<EP>  P;
};


template< class T >
struct knot_vec
{
  int         m;
  Ptr< T >    U;
};
template< class HP >
struct point_vec
{
  int              n;
  Ptr< HP >        Pw;
};
template< class T, class HP >
struct curve
{
  int             p;
  knot_vec<T>     knt;
  point_vec<HP>   cpt;
};

template< class HP >
struct point_net
{  
  int               nu;
  int               nv;
  Ptr< Ptr< HP > >  Pw;
};
template< class T, class HP > 
struct surface
{
  int            pu;
  int            pv;
  knot_vec<T>    knu;
  knot_vec<T>    knv;
  point_net<HP>  net;
};  
                                        /* simple list templates */
template< class T > class ListNode
{
public:
  T el;
  ListNode *next;
  ListNode *prev;

  ListNode() { }
  ListNode( const T& a ) : el(a) { }
    
  void *operator new( size_t s )
  {
    size_t dummy;
    void *p = gust::PoolAlloc( s, &dummy );
    if( p == NULL ) throw PoolAllocError();
    return(p);
  }
  void operator delete( void *p, size_t s )
  {
    gust::PoolFree( p, s );
  }
};

template< class T > class List   // simple List
{
public:
  T   *head;
  int  nElems;
  
  List() { head = NULL; nElems = 0; }
  ~List() { DeleteElems(); }

  T *First(void) const { return head; }
  T *Last(void) const
  {
    T *n0 = 0;
    T *n = head;
    while(n) { n0 = n; n = n->next; }
    return n0;
  }
  bool IsEmpty(void) const { return(nElems == 0); }

  void Append( T *node )
  {
    node->next = head; 
    node->prev = NULL;
    if( head != NULL ) head->prev = node; 
    head = node;
    nElems++;
  }

  void Remove( T *node )
  {
    if( node->prev != NULL ) node->prev->next = node->next;
    if( node->next != NULL ) node->next->prev = node->prev;
    if( head == node ) head = node->next;
    nElems--;
  }
  void DeleteElems( void )
  {
    T *n = head, *nnext;
    while( n != NULL ) { nnext = n->next; delete n; n = nnext; }
    head = NULL; nElems = 0;
  }
  void ReleaseElems( void )
  {
    head = NULL; nElems = 0;
  }
  
  List<T>& operator+=( List<T>& other )
  {
    T *olast;

    if( &other == this ) return *this;

    if( head ) 
    {
      olast = other.Last();
      if( olast ) 
      {
        olast->next = head;
        head->prev = olast;
        head = other.First();
        nElems += other.nElems;
        other.ReleaseElems();
      }
    }
    else
    {
      head = other.First();
      nElems = other.nElems;
      other.ReleaseElems();
    }
    return *this;
  }
  List<T>& operator=( List<T>& other )
  {
    if( this == &other ) return *this;

    DeleteElems();
    return ((*this) += other);
  }
  List( List<T>& other )
  {
    head = NULL; nElems = 0;
    (*this) += other;
  }

// private:                      // don't allow simple copying
//   List( const List<T>& other );
};

template< class T >
struct ListNodeInfo
{
  List<T> *m_L;
  T       *m_n;

  ListNodeInfo( List<T> *L, T *n ) { m_L = L; m_n = n; }  
};

template< class T > class List2   // simple List, bidirectional
{
public:
  T   *tail;
  T   *head;
  int  nElems;
  
  List2() { head = NULL; tail = NULL; nElems = 0; }
  ~List2() { DeleteElems(); }

  T *First(void) const { return head; }
  T *Last(void) const  { return tail; }
  
  bool IsEmpty(void) const { return(nElems == 0); }
  
  void AppendLeft( T *node )
  {
    node->next = head; 
    node->prev = NULL;
    if( head != NULL ) { head->prev = node; } else { tail = node; } 
    head = node;
    nElems++;
  }
  void AppendRight( T *node )
  {
    node->prev = tail; 
    node->next = NULL;
    if( tail != NULL ) { tail->next = node; } else { head = node; } 
    tail = node;
    nElems++;
  }
  void Remove( T *node )
  {
    if( node->prev != NULL ) node->prev->next = node->next;
    if( node->next != NULL ) node->next->prev = node->prev;
    if( head == node ) head = node->next;
    if( tail == node ) tail = node->prev;
    nElems--;
  }
  void DeleteElems( void )
  {
    T *n = head, *nnext;
    while( n != NULL ) { nnext = n->next; delete n; n = nnext; }
    head = tail = NULL; nElems = 0;
  }
  void ReleaseElems( void )
  {
    head = tail = NULL; nElems = 0;
  }

  List2<T>& operator+=( List2<T>& other )
  {
    T *olast;

    if( &other == this ) return *this;

    if( head ) 
    {
      olast = other.Last();
      if( olast ) 
      {
        olast->next = head;
        head->prev = olast;
        head = other.First();  // tail remains the same
        nElems += other.nElems;
        other.ReleaseElems();
      }
    }
    else
    {
      head = other.First();
      tail = other.Last();
      nElems = other.nElems;
      other.ReleaseElems();
    }
    return *this;
  }
  List2<T>& operator=( List2<T>& other )
  {
    if( this == &other ) return *this;

    DeleteElems();
    return ((*this) += other);
  }
  // caution, this emptys the input list
  List2( List2<T>& other )
  {
    head = NULL; tail = NULL; nElems = 0;
    (*this) += other;
  }

// private:                      // don't allow simple copying
//   List2( const List2<T>& other );
};

                              /* simple stack template */
template< class T >     
struct Stack
{
  T     *m_sp;
  Ptr<T> m_buf;

  Stack()  { m_sp = 0; } 
  ~Stack() { }

  void init(unsigned int maxelems) 
  { 
    m_buf.reserve_pool(maxelems); 
    m_sp = &m_buf[0]; 
  } 
  void reset() { m_sp = &m_buf[0]; }
  unsigned int release( T **elems ) 
  { 
    unsigned int i = (unsigned int)(m_sp - &m_buf[0]);
    *elems = m_buf.release(); 
    m_sp = 0;
    return i;
  }
  void push( const T& el ) { *m_sp = el; m_sp++; }
  T pop() { m_sp--; return *m_sp; }
};

}

#include "guma_sorting.h"

namespace gul {

using guma::BBTNode;
using guma::PtrHeapSort;
using guma::InitBBT;
using guma::InsertElementIntoBBT;
using guma::DeleteElementFromBBT;
using guma::SearchElementInBBT;
using guma::SearchSuccessorInBBT;
using guma::SearchPredecessorInBBT;

/* ----------------------------------------------------------------------
  Balanced Binary Tree (allocates/deallocates memory for keys)
------------------------------------------------------------------------ */
template< class T >
struct Map
{
  struct Node
  {
    BBTNode *m_node;
    Node() : m_node(0) { }
    Node(BBTNode *node) : m_node(node) { }
    BBTNode *bbtnode() { return m_node; }
    T *key() { return (T *)m_node->key; }
    bool IsEmpty() { return m_node == 0; }
    Node left() { return Node(m_node->left); }
    Node right() { return Node(m_node->right); }
  };

  BBTNode *root;

  Map()
  {
    size_t dummy;

    root = (BBTNode *)gust::PoolAlloc( sizeof(BBTNode), &dummy );
    if( root == NULL ) throw PoolAllocError();
    (*((int *)&root->parent)) = 1;
    InitBBT(root);   // must leave parent field untouched !!!!
  }
  ~Map()
  {
    // misusing the unused parent field of the special root node as reference
    // counter (to have at least pointer copy semantics)

    if( --(*((int *)&root->parent)) == 0 )
    {
      DeleteElems();
      gust::PoolFree( root, sizeof(BBTNode) );   
    }
  }
  Map( const Map& other )
  {
    (*((int *)&other->root->parent))++; 
    root = other->root;
  }
  Map& operator=( const Map& other )
  {
    (*((int *)&other->root->parent))++;
    if( --(*((int *)&root->parent)) == 0 )
    {
      DeleteElems();
      gust::PoolFree( root, sizeof(BBTNode) );   
    }
    root = other->root;
  }

  Node Head( void ) { return(Node(root->right)); }
  int Height( void ) { return( *((int *)(&root->left)) ); }
  bool IsEmpty( void ) { return root->right == 0; }
  bool IsRemoved( Node node ) 
  { return node.bbtnode()->parent == 0; }

  Node NewNode()
  {
    BBTNode *node;
    size_t dummy;
    
    node = (BBTNode *)gust::PoolAlloc( sizeof(BBTNode), &dummy );
    if( node == NULL ) throw PoolAllocError();
    node->key = gust::PoolAlloc( sizeof(T), &dummy );
    if( node->key == NULL ) throw PoolAllocError();
    ::new(node->key) T();

    return(Node(node));
  }
  template< class T1 >
  Node NewNode( T1& key )
  {
    BBTNode *node;
    size_t dummy;
    
    node = (BBTNode *)gust::PoolAlloc( sizeof(BBTNode), &dummy );
    if( node == NULL ) throw PoolAllocError();
    node->key = gust::PoolAlloc( sizeof(T), &dummy );
    if( node->key == NULL ) throw PoolAllocError();
    ::new(node->key) T(key);

    return(Node(node));
  }
  //
  // same as above, but for objects which can only be passed by value
  // (e.g. something like "int(4)")
  //
  template< class T1 >
  Node NewNodeV( T1 key )
  {
    BBTNode *node;
    size_t dummy;
    
    node = (BBTNode *)gust::PoolAlloc( sizeof(BBTNode), &dummy );
    if( node == NULL ) throw PoolAllocError();
    node->key = gust::PoolAlloc( sizeof(T), &dummy );
    if( node->key == NULL ) throw PoolAllocError();
    ::new(node->key) T(key);

    return(Node(node));
  }
  T *ReleaseKey( Node node )
  {
    T *t = 0;
    if( !node.IsEmpty() )
    {
      t = node.key();
      node.bbtnode()->key = 0;
    }
    return t;
  }
  void FreeNode( Node node )
  {
    if( !node.IsEmpty() )
    {
      if( node.key() != 0 )
      {
        node.key()->~T();
        gust::PoolFree( node.key(), sizeof(T) );
      }
      gust::PoolFree( node.bbtnode(), sizeof(BBTNode) );   
    }
  }
  void do_destroy( Node node )
  {
    if( node.IsEmpty() ) return;
    
    if( !node.left().IsEmpty() ) do_destroy( node.left() );
    if( !node.right().IsEmpty() ) do_destroy( node.right() );

    FreeNode( node );
  }
  void DeleteElems( void )
  {
    do_destroy( Head() );
    InitBBT(root);
  }
  void InsertNode( Node element, Node where, int side )
  {
    InsertElementIntoBBT( element.bbtnode(), where.bbtnode(), side, root );
  }
  void RemoveNode( Node element )
  {
    BBTNode *el = element.bbtnode();
    DeleteElementFromBBT( el, root );
    el->parent = 0; // to indicate that node isn't in the tree anymore
  }
  template< class K >
  int SearchNode( K *key, int (*compfunc)(K *, T *), Node *element )
  {
    BBTNode *el;
    int side = SearchElementInBBT(root, (void *)key, 
                                  (int (*)(void *, void *))compfunc, &el );
    *element = Node(el);
    return side;
  }
  Node Successor( Node element )
  {
    return Node( SearchSuccessorInBBT(element.bbtnode(), root) );
  }
  Node Predecessor( Node element )
  {
    return Node( SearchPredecessorInBBT(element.bbtnode(), root) );
  }
  Node First( void )  // deliver first element in symmetric order
  { 
    BBTNode *node = NULL, *next = root->right;

    while( next != NULL )
    {
      node = next;
      next = node->left;
    }
    
    return( Node(node) );
  }
  Node Last( void )  // deliver last element in symmetric order
  { 
    BBTNode *node = NULL, *next = root->right;

    while( next != NULL )
    {
      node = next;
      next = node->right;
    } 
    return( Node(node) );
  }
  void do_convert_to_array( int *nEl, Node node, Ptr<T *> *A )
  {
    if( !node.IsEmpty() )
    {
      do_convert_to_array( nEl, node.left(), A );
      (*A)[*nEl] = node.key();
      (*nEl)++;
      do_convert_to_array( nEl, node.right(), A );
    }
  }
  int ConvertToArray( Ptr<T *> *retKey ) 
  {
    int nEl = 0, maxEl = 1 << Height();
    
    (*retKey).reserve_pool( maxEl );
    
    do_convert_to_array( &nEl, Head(), retKey );

    return nEl;
  }
};

/* ----------------------------------------------------------------------
  Balanced Binary Tree (doesn't allocates/deallocates memory for keys)
------------------------------------------------------------------------ */
template< class T >
struct RefMap
{
  struct Node
  {
    BBTNode *m_node;
    Node() : m_node(0) { }
    Node(BBTNode *node) : m_node(node) { }
    BBTNode *bbtnode() { return m_node; }
    T *key() { return (T *)m_node->key; }
    bool IsEmpty() { return m_node == 0; }
    Node left() { return Node(m_node->left); }
    Node right() { return Node(m_node->right); }
  };

  BBTNode *root;

  RefMap()
  {
    size_t dummy;
    
    root = (BBTNode *)gust::PoolAlloc( sizeof(BBTNode), &dummy );
    if( root == NULL ) throw PoolAllocError();
    (*((int *)&root->parent)) = 1;
    InitBBT(root);   // must leave parent field untouched !!!!
  }
  ~RefMap()
  {
    // misusing the unused parent field of the special root node as reference
    // counter (to have at least pointer copy semantics)

    if( --(*((int *)&root->parent)) == 0 )
    {
      DeleteElems();
      gust::PoolFree( root, sizeof(BBTNode) );   
    }
  }
  RefMap( const RefMap& other )
  {
    (*((int *)&other->root->parent))++; 
    root = other->root;
  }
  RefMap& operator=( const RefMap& other )
  {
    (*((int *)&other->root->parent))++;
    if( --(*((int *)&root->parent)) == 0 )
    {
      DeleteElems();
      gust::PoolFree( root, sizeof(BBTNode) );   
    }
    root = other->root;
  }
  Node Head( void ) { return(Node(root->right)); }
  int Height( void ) { return( *((int *)(&root->left)) ); }
  bool IsEmpty( void ) { return root->right == 0; }
  bool IsRemoved( Node node ) 
  { return node.bbtnode()->parent == 0; }

  Node NewNode( T *key )
  {
    BBTNode *node;
    size_t dummy;
    
    node = (BBTNode *)gust::PoolAlloc( sizeof(BBTNode), &dummy );
    if( node == NULL ) throw PoolAllocError();
    node->key = key;

    return(Node(node));
  }
  void FreeNode( Node node )
  {
    if( !node.IsEmpty() )
      gust::PoolFree( node.bbtnode(), sizeof(BBTNode) );   
  }
  void do_destroy( Node node )
  {
    if( node.IsEmpty() ) return;
    
    if( !node.left().IsEmpty() ) do_destroy( node.left() );
    if( !node.right().IsEmpty() ) do_destroy( node.right() );

    gust::PoolFree( node.bbtnode(), sizeof(BBTNode) );   
  }
  void DeleteElems( void )
  {
    do_destroy( Head() );
    InitBBT(root);
  }
  void InsertNode( Node element, Node where, int side )
  {
    InsertElementIntoBBT( element.bbtnode(), where.bbtnode(), side, root );
  }
  void RemoveNode( Node element )
  {
    BBTNode *el = element.bbtnode();
    DeleteElementFromBBT( el, root );
    el->parent = 0; // to indicate that node isn't in the tree anymore
  }
  template< class K >
  int SearchNode( K *key, int (*compfunc)(K *, T *), Node *element )
  {
    BBTNode *el;
    int side = SearchElementInBBT(root, (void *)key, 
                                  (int (*)(void *, void *))compfunc, &el );
    *element = Node(el);
    return side;
  }
  Node Successor( Node element )
  {
    return Node( SearchSuccessorInBBT(element.bbtnode(), root) );
  }
  Node Predecessor( Node element )
  {
    return Node( SearchPredecessorInBBT(element.bbtnode(), root) );
  }
  Node First( void )  // deliver first element in symmetric order
  { 
    BBTNode *node = NULL, *next = root->right;

    while( next != NULL )
    {
      node = next;
      next = node->left;
    } 
    return( Node(node) );
  }
  Node Last( void )  // deliver last element in symmetric order
  { 
    BBTNode *node = NULL, *next = root->right;

    while( next != NULL )
    {
      node = next;
      next = node->right;
    } 
    return( Node(node) );
  }
  template< class U >
  void do_convert_to_array( int *nEl, Node node, Ptr<U> *A )
  {
    if( !node.IsEmpty() )
    {
      do_convert_to_array( nEl, node.left(), A );
      (*A)[*nEl] = (U)node.key();
      (*nEl)++;
      do_convert_to_array( nEl, node.right(), A );
    }
  }
  template< class U >
  int ConvertToArray( Ptr<U> *retKey ) 
  {
    int nEl = 0, maxEl = 1 << Height();
    
    (*retKey).reserve_pool( maxEl );
    
    do_convert_to_array<U>( &nEl, Head(), retKey );

    return nEl;
  }
};

// base class for objects which are reserved in the pool 

class pool_object
{
public:
  void *operator new( size_t s )
  {
    size_t dummy;
    void *p = gust::PreferPoolAlloc( s, &dummy );
    if( p == NULL ) throw PoolAllocError();
    return(p);
  }
  
  void operator delete( void *p, size_t s )
  {
    gust::PreferPoolFree( p, s );
  }
};

// these structs are used for storing records in a kdtree

// structure which can be used for records with any number of dimensions
template< class T >
class kdarray : public pool_object
{
public:
  T *m_base;
  kdarray() { }
  kdarray( T *base ) { m_base = base; }
  T& operator[]( size_t i) const { return m_base[i]; } 
};

// structure which can be used to attach additional data to records
template< class T, class EP >
struct kdpoint : public pool_object
{
  EP m_P;
  void *m_id;  // for casting the base class to derived classes,
               // or for tagging additional data to a record
  kdpoint( const EP& P, void *id ) { m_P = P; m_id = id; }
  const T& operator[] ( size_t i ) const { return m_P[i]; }
};


template< class T >
class poolalloc
{
public:
  typedef T         value_type;
  typedef size_t    size_type;
  typedef ptrdiff_t difference_type;
  typedef T*        pointer;
  typedef const T*  const_pointer;
  typedef T&        reference;
  typedef const T&  const_reference;
  
  pointer address( reference r ) const { return &r; }
  const_pointer address( const_reference r ) const { return &r; } 

  poolalloc() throw() { }
  template< class U > poolalloc( const poolalloc<U>& ) throw() { }
  ~poolalloc() throw() { }
  
  inline pointer allocate( size_type n, const void *h = 0 ); 
  inline void deallocate( pointer p, size_type n );
  
  void construct( pointer p, const T& val ) { new(p) T(val); }
  void destroy( pointer p ) { p->~T(); }
  size_type max_size() const throw() { return ((size_t)-1); }
  
  template< class U >
  struct rebind { typedef poolalloc<U> other; };
};

template< class T >
inline bool operator==( const poolalloc<T>&, const poolalloc<T>& ) throw() 
{ return true; }
template< class T >
inline bool operator!=( const poolalloc<T>&, const poolalloc<T>& ) throw() 
{ return false; }

template<> class poolalloc<void>
{
public:
  typedef void*       pointer;
  typedef const void* const_pointer;
  typedef void        value_type;

  template< class U >
  struct rebind { typedef poolalloc<U> other; }; 
};

template< class T >
inline T *poolalloc<T>::allocate( size_type n, poolalloc<void>::const_pointer h )
{
  size_t dummy;
  T *p = (T *)gust::PreferPoolAlloc( n*sizeof(T), &dummy );
  if( p == 0 ) throw gul::AllocError();
  return p;
}
template< class T >
inline void poolalloc<T>::deallocate( pointer p, size_type n )
{
  gust::PreferPoolFree( p, n*sizeof(T) );
}

/*-------------------------------------------------------------------------
  this is a hack to get the 'basic_string' class of Gnu C 2.95.2 working with
  PoolAlloc (in this implementation 'basic_string' doesn't uses the
  std::allocator !!!! ) 
-------------------------------------------------------------------------*/

class pool_bytealloc
{
public:
  static void *allocate( size_t n )
  {
    size_t dummy;
    void *p = gust::PreferPoolAlloc( n, &dummy );
    if( p == 0 ) throw gul::AllocError();
    return p;
  }
  static void deallocate( void *p, size_t n )
  {
    gust::PreferPoolFree( p, n );
  }
};



}

#endif

---