<commit_before>/*
 *    This file is part of HMMoC-BFloat-Algebra, a tiny C++ library supplying
 *    the Algebra<> template, for polymorphic representation of non-negative
 *    floating-point values in either native (eg double), log-transformed, or
 *    BFloat representations.  BFloat is the Buoyant Float type, allowing
 *    fixed-precision representation at an arbitrary precision range.  The
 *    library is based on code from: <NAME>. HMMoC—a compiler for hidden
 *    Markov models. Bioinformatics (2007) 23(18): 2485-2487. The code was
 *    derived by <NAME> and <NAME> from <NAME>, <NAME>. An
 *    expectation maximization algorithm for training hidden substitution
 *    models.  J Mol Biol. 2002 Apr 12;317(5):753-64.  You may use at will,
 *    subject to the license (Apache v2.0), but *please cite those two papers* in
 *    your documentation and publications associated with uses of this library.
 *    Thank you!
 *
 *
 *    Copyright (C) 2011 by <NAME>, <NAME> Cancer Research
 *    Center.
 *
 *    This file is based on algebras.h in HMMoC 1.3, a hidden Markov model
 *    compiler.  Copyright (C) 2007 by <NAME>, Oxford University.
 *
 *    Licensed under the Apache License, Version 2.0 (the "License");
 *    you may not use this file except in compliance with the License.
 *    You may obtain a copy of the License at
 *    
 *        http://www.apache.org/licenses/LICENSE-2.0
 *    
 *    Unless required by applicable law or agreed to in writing, software
 *    distributed under the License is distributed on an "AS IS" BASIS,
 *    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *    See the License for the specific language governing permissions and
 *    limitations under the License.
\*/
//
// Algebra.hpp - extended real types
//
// <NAME>, June 1, 2011
// <NAME>, 27/8/04
//

#ifndef HMMOC_BFLOAT_ALGEBRA_HPP
#define HMMOC_BFLOAT_ALGEBRA_HPP

// MS compatible compilers support #pragma once
#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif

#if defined(_MSC_VER) && (_MSC_VER >= 1310)
#  pragma warning (disable : 4675) // suppress ADL warning
#endif

#include <assert.h>

// Do under- and overflow checking in BFloat? (undef for no)
#ifndef NDEBUG
#define BFLOAT_CHECK_UOFLOW
// Do extra checks in operator- to diagnose negative values? (set to 0 for no)
#define ALGEBRA_CHECK_NEGATIVE
#endif // NDEBUG

#ifndef NBOOST_SERIALIZATION // #define NBOOST_SERIALIZATION if you don't want to compile in the boost serialization library
#include <boost/serialization/nvp.hpp>
#include <boost/serialization/utility.hpp>
#include <boost/serialization/version.hpp>
#endif // NBOOST_SERIALIZATION

// TODO: PUT BACK ?  See below (search for seqan)...
//#include <seqan/basic.h>

#include <iostream>
#include <cmath> // for isinf and isnan
#include <limits>

// You can control some inlining recommentations to the compiler using these:
#define ALGEBRA_INLINE_ARITHMETIC inline
#define ALGEBRA_INLINE_CAST inline
#define ALGEBRA_INLINE_CAST_COMPARE inline
#define ALGEBRA_INLINE_CAST_ARITHMETIC inline

// typedefs
typedef float BFMantissa;
//typedef double BFMantissa;
const BFMantissa cBFloatRange = 20282409603651670423947251286016.0;  // 2.03e+31; 2^104
const BFMantissa cBFloatRangeInv = 1.0/cBFloatRange;
// <NAME> 6/7/7: need to typecast to avoid compiler warnings about imprecise FP representations
const BFMantissa cBFloatRangeSqrt    = (BFMantissa)1.0e+18;          // Value between square root of the exponent, and the exponent
const BFMantissa cBFloatRangeInvSqrt = (BFMantissa)1.0e-18;          // Square of this should still be representable, with full mantissa!
const BFMantissa logcBFloatRange     = ::log(cBFloatRange);
const int cBFloatDigits              = 7;                 // Number of significant digits (7 for floats, 16 for doubles?)
const int cBFloatInfinity            = 1000000000;        // Tiniest number representable is cBFloatRangeInv ^ BFloatInfinity
// PTE 15/9/08
const BFMantissa cBFloatEpsilon      = ( BFMantissa )1.193e-07;
const int cBFloatConvTableSize       = 100;               // This includes many zero entries, it makes additions a bit faster
const int cBFloatDoubleConvTableSize = 50;                // Table size for bfloat -> double conversion; cBFloatRange^(-size/2) is double 0
const double cLogspaceEpsilon        = 36.7; // TODO: Test & verify

namespace hmmoc {

// PTE 14/9/08
// To help disambiguate ::log( double ) from hmmoc::log(Algebra)
template <typename other_type>
inline
other_type
log ( other_type const & to_be_logged )
{
  return ::log( to_be_logged );
} // log( other_type const & )

// Forward decls
class Logspace;
class DoubleRealspace;
class LongDoubleRealspace;
class FloatRealspace;
class BFloat;
static inline std::ostream& bfloat_print ( std::ostream& out, const BFloat& x );
static inline void assertBFloatValidity ( BFloat const & bf );
static inline void BFloatNormalise ( BFloat& a );

//
// BFloats: more buoyant floats.
//
// NOTE: For now BFloats can represent only positive values.
//
// struct{ float + int } is 8 bytes; nice size makes noticable speed difference
//
// For BFloats, infinity is represented as ( e == cBFloatInfinity, and we never let f be std::numeric_limits<BFMantissa>::infinity() ).  0 is represented as ( f == e == 0 ), and NaN as ( f != f ).
class BFloat {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  inline
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( f );
    ar & BOOST_SERIALIZATION_NVP( e );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  static BFMantissa* aConversionLookup;              // used by addition
  static double* aDoubleConversionLookup;            // used by Value()
  BFMantissa f;
  int e;
 public:
  inline
  BFloat ( BFMantissa const & iF, int const & iE, bool normalise = false ) :
    f( iF ), e( iE )
  {
    if( normalise ) {
      BFloatNormalise( *this );
    } else {
      assertBFloatValidity( *this );
    }
  };

  inline
  BFloat () : f( 0.0 ), e( 0 ) {};

  inline
  BFloat ( double const & prob ) { *this = prob; }

  inline
  BFloat ( BFloat const & copy_from ) : f( copy_from.f ), e( copy_from.e ) {
    assertBFloatValidity( *this );
  };

  // Defined later, after these others are defined.
  inline
  BFloat ( Logspace const & copy_from );
  // Defined later, after these others are defined.
  inline
  BFloat ( DoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline
  BFloat ( LongDoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline
  BFloat ( FloatRealspace const & copy_from );

  inline
  ~BFloat () {};

  inline
  BFloat &
  operator= ( BFloat const & copy_from )
  {
    f = copy_from.f;
    e = copy_from.e;

    return *this;
  } // operator=( BFloat const & )

  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( Logspace const & copy_from );
  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( DoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( LongDoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( FloatRealspace const & copy_from );

  inline BFloat &
  operator= ( double prob )
  {
    if( prob != prob ) {
      // NaN
      // For now, don't let it be set to NaN.
      assert( prob == prob ); // NOT NaN
      f = prob;
      e = 0;
    } else
    if( prob == std::numeric_limits<double>::infinity() ) {
      f = 1.0;
      e = cBFloatInfinity;
    } else
    // Simplistic double-to-BFloat conversion - can be slow if 'standard' numbers get very large/small
    if( prob < std::numeric_limits<double>::min() ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
      if( prob < 0.0 ) {
        std::cerr << "BFloat: Negative number: " << prob << std::endl;
        // TODO: ?
        assert( prob >= 0.0 );
      }
#endif
      f = 0.0;
      e = 0;
    } else {
      f = 0.0;
      e = 0;
      while( prob > cBFloatRangeSqrt ) {
        prob *= cBFloatRangeInv;
        e++;
      }
      while( prob < cBFloatRangeInvSqrt ) {
        prob *= cBFloatRange;
        e--;
      }
      f = prob; 
    }
    assertBFloatValidity( *this );
    return *this;
  } // operator=( double const & )

  inline double Value () const { 
    if( e == 0 ) {
      return f;
    }
    if( e == cBFloatInfinity ) {
      return std::numeric_limits<double>::infinity();
    }
    if( f != f ) {
      return std::numeric_limits<double>::quiet_NaN();
    }
    if (abs(e) < cBFloatDoubleConvTableSize/2) {
      return (double)f * aDoubleConversionLookup[ e + cBFloatDoubleConvTableSize/2 ];
    } else if (e < cBFloatDoubleConvTableSize/2) {
      return 0.0;
    } else {
      return (double)f * exp((double)e * logcBFloatRange);
    }
  } // double Value()

  // sets to zero
  inline
  void clear () {
    f=0.0; e=0;
  }

  inline
  bool isZero () const
  {
    return ( f == 0.0 );
  } // isZero() const

  inline
  bool isOne () const
  {
    return ( f == 1.0 );
  } // isOne() const

  inline
  bool isInfinity () const
  {
    return ( e == cBFloatInfinity );
  } // isInfinity() const

  inline
  bool isNaN () const
  {
    return ( f != f );
  } // isNaN() const

}; // End class BFloat


//
// dummy class to initialise BFloat lookup table
//
class _BFloatInitialize {
public:
  _BFloatInitialize();
};


//
// implementations of BFloat calculations
//

static inline void assertBFloatValidity ( BFloat const & bf )
{
  // Never negative
  assert( bf.f >= 0.0 );

  // TODO: REMOVE?  NaN should be ok to represent, but we might want to raise a red flag...
  assert( !bf.isNaN() );

  // we represent infinity as ( e == cBFloatInfinity ), not as ( f == std::numeric_limits<BFMantissa>::infinity() ).
  assert( bf.f != std::numeric_limits<BFMantissa>::infinity() );

  if( bf.isZero() ) {
    // 0 is represented as both e and f being 0
    assert( bf.e == 0 );
    //} else if( !bf.isInfinity() ) {
    // f must be in this range (AFTER NORMALIZATION)
    //assert( bf.f <= cBFloatRangeSqrt ); 
    //assert( bf.f >= cBFloatRangeInvSqrt ); 
  }
    
  return;
} // static assertBFloatValidity( BFloat const & )

// Normalization of BFloat result of a single operation
#ifdef BFLOAT_CHECK_UOFLOW
static inline void BFloatNormalise ( BFloat& a )
     //#define BFloatNormalise(a)
{
  if (a.f > cBFloatRangeSqrt) {
    a.f *= cBFloatRangeInv;
    a.e++;
  } else if (a.f < cBFloatRangeInvSqrt) {
    if( a.f == 0.0 ) {
      a.e = 0; 
    } else {
      a.f *= cBFloatRange;
      a.e--;
    }
  }
  if( a.e > cBFloatInfinity ) {
    std::cerr << "BFloat: Overflow" << std::endl;
    a.e = cBFloatInfinity;
  } else if( a.e < -cBFloatInfinity ) {
    std::cerr << "BFloat: Underflow" << std::endl;
    a.e = 0;
    a.f = 0.0;
  }
  // TODO: REMOVE
#ifdef ALGEBRA_CHECK_NEGATIVE
    if( a.f < 0 ) {
      std::cerr << "BFloat: Negative number: ";
      bfloat_print( std::cerr, a );
      std::cerr << std::endl;
      std::cerr << "\tin BFloatNormalise( BFMantissa & ) [BFLOAT_CHECK_UOFLOW is true]" << std::endl;
    }
#endif // ALGEBRA_CHECK_NEGATIVE
#ifndef NDEBUG
  assertBFloatValidity( a );
#endif // NDEBUG
}; // BFloatNormalise ( BFloat & ) // when BFLOAT_CHECK_UOFLOW
#else // if BFLOAT_CHECK_UOFLOW .. else ..
static inline void BFloatNormDown ( BFloat& a ) { 
  a.f *= cBFloatRangeInv;
  a.e++;
  assert( a.f <= cBFloatRangeSqrt );
}
static inline void BFloatNormUp ( BFloat& a ) { 
  if (a.f == 0.0) {
    a.e = 0;
  } else {
    a.f *= cBFloatRange;
    a.e--;
    assert( a.f >= cBFloatRangeInvSqrt );
  }
}
static inline void BFloatNormalise ( BFloat& a )
     //#define BFloatNormalise(a) 
{
  if (a.f > cBFloatRangeSqrt) {
    BFloatNormDown(a);
  } else if (a.f < cBFloatRangeInvSqrt) {
    BFloatNormUp(a);
  }
#ifdef ALGEBRA_CHECK_NEGATIVE
    if( a.f < 0 ) {
      std::cerr << "BFloat: Negative number: ";
      bfloat_print( std::cerr, a );
      std::cerr << std::endl;
      std::cerr << "\tin BFloatNormalise( BFMantissa & ) [BFLOAT_CHECK_UOFLOW is false]" << std::endl;
    }
#endif // ALGEBRA_CHECK_NEGATIVE
#ifndef NDEBUG
  assertBFloatValidity( a );
#endif // NDEBUG
};  // BFloatNormalise ( BFloat & ) // when NOT BFLOAT_CHECK_UOFLOW
#endif // End if BFLOAT_CHECK_UOFLOW .. else ..

static inline void DoubleNormalise ( double& f, int& e )
{
  // comparing to 0.0 here fails, because the comparison is done
  // using higher-precision doubles, but the subsequent while-loop
  // uses true doubles, resulting in an infinite loop. (G.L. 3/9/07)
  if( f < std::numeric_limits<double>::min() ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if( f < 0.0 ) {
      std::cerr << "BFloat: Negative number: " << f << std::endl;
      // TODO: ?
      assert( f >= 0.0 );
    }
#endif // ALGEBRA_CHECK_NEGATIVE
    f = 0.0; 
    e = 0;
  } else {
    while (f > (double)cBFloatRangeSqrt) {
      f *= (double)cBFloatRangeInv;
      e++;
    }
    while (f < (double)cBFloatRangeInvSqrt) {
      f *= (double)cBFloatRange;
      e--;
    }
  }
 }; // DoubleNormalise( double &, int & )

// Logarithm of a BFloat
static inline double bfloat_doublelog ( const BFloat& a ) { return a.e*logcBFloatRange+::log(a.f); }

// BFloat exp of a double -- could be a teensy tiny value.
static inline BFloat bfloat_doubleexp ( double iA ) 
{
  if( iA == std::numeric_limits<double>::infinity() ) {
    // exp( inf ) is inf
    return BFloat( 1.0, cBFloatInfinity );
  }
  if( iA == -std::numeric_limits<double>::infinity() ) {
    // exp( -inf ) is 0
    return 0.0;
  }
  int iE = (int)floor( iA / ::log(cBFloatRange) );
  iA -= iE * ::log(cBFloatRange);
  return BFloat( ::exp(iA), iE, true );
} // static bfloat_doubleexp( double )

// Returns a double value - or underflow/overflow if it does not fit.
static inline double bfloat2double ( const BFloat & bfloat) { return bfloat.Value(); }

/// TODO: PUT BACK
//template <typename T>
//inline static double
//toDouble ( T const & v )
//{
//  return seqan::convert<T, double>( v );
//} // toDouble( T const & )

inline static double
toDouble ( BFloat const & v )
{
  return bfloat2double( v );
} // toDouble( BFloat const & )

inline static long double
toLongDouble ( BFloat const & v )
{
  return bfloat2double( v );
} // toLongDouble( BFloat const & )

inline static double
toLogDouble ( BFloat const & v )
{
  return bfloat_doublelog( v );
} // toLogDouble( BFloat const & )

static inline BFloat double2bfloat ( double prob )
{
  return prob;
}

static inline BFloat bfloat_pr_product ( const BFloat& a, const BFloat& b )
{ 
  //assertBFloatValidity( a );
  //assertBFloatValidity( b );
  if( a.isZero() ) {
    return 0.0;
  }
  if( b.isZero() ) {
    return 0.0;
  }
  return BFloat( a.f*b.f, a.e+b.e, true );
} // static bfloat_pr_product( BFloat const &, BFloat const & )

static inline BFloat bfloat_pr_double_product ( const BFloat& a, double const & b )
{ 
  //assertBFloatValidity( a );
  if( a.isZero() ) {
    return 0.0;
  }
  if( b == 0.0 ) {
    return 0.0;
  }
  register double mantisse = a.f*b;
  int exponent = a.e;
  DoubleNormalise(mantisse, exponent);
  return BFloat( mantisse, exponent, false );
} // bfloat_pr_double_product ( BFloat const &, double const & )

static inline void bfloat_pr_product_accum ( BFloat& a, const BFloat& b )
{
  if( a.isZero() ) {
    // It's already 0.
    return;
  }
  if( b.isZero() ) {
    // Multiplying by 0..
    a.clear();
    return;
  }
  a.f *= b.f; a.e += b.e; 
  BFloatNormalise( a );
} // static bfloat_pr_product_accum( BFloat &, BFloat const & )

static inline void bfloat_pr_double_product_accum (  BFloat& a, double const & b )
{ 
  if( a.isZero() ) {
    // It's already 0.
    return;
  }
  if( b == 0.0 ) {
    // Multiplying by 0..
    a.clear();
    return;
  }
  register double mantisse = a.f*b;
  DoubleNormalise( mantisse, a.e );
  a.f = mantisse;
  assertBFloatValidity( a );
  return;
} // static bfloat_pr_double_product_accum( BFloat &, double const & )

// PTE 14/9/08
static inline BFloat bfloat_pr_power (const BFloat& a, const BFloat& b) 
{ 
  return BFloat( ::pow(a.f,b.f), (a.e*b.e), true );
}

static inline BFloat bfloat_pr_double_power ( const BFloat& a, double const & b )
{
  // TODO: test, make sure it's safe.  Probably it's not.
  register double mantisse = ::pow(a.f, b);
  int exponent = a.e;
  DoubleNormalise(mantisse, exponent);
  return BFloat(mantisse, exponent);
}

static inline BFloat bfloat_pr_quotient ( const BFloat& a, const BFloat& b )
{
  if( b.isZero() ) {
    return std::numeric_limits<double>::quiet_NaN();
  } else if( b.isInfinity() ) {
    return 0.0;
  }
  return BFloat( a.f/b.f, a.e-b.e, true );
} // bfloat_pr_product ( BFloat const &, BFloat const & )
  
static inline void bfloat_pr_quotient_accum ( BFloat& a, const BFloat& b ) 
{ 
  if( b.isZero() ) {
    a = std::numeric_limits<double>::quiet_NaN();
  } else if( b.isInfinity() ) {
    a = 0.0;
  } else {
    a.f /= b.f; 
    a.e -= b.e;
  }
  BFloatNormalise( a ); 
} // bfloat_pr_quotient_accum( BFloat &, BFloat const & )

static inline BFloat bfloat_pr_sum ( const BFloat& a, const BFloat& b )
{
  if( b.isZero() ) {
    assertBFloatValidity( a );
    return a;
  }
  if( a.isZero() ) {
    assertBFloatValidity( b );
    return b;
  }
  if( b.isInfinity() || a.isInfinity() ) {
    return std::numeric_limits<double>::infinity();
  }

  BFloat r;
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize) {
      assertBFloatValidity( a );
      return a;
    } else {
      r.f = a.f + b.f * BFloat::aConversionLookup[ a.e - b.e ];
      r.e = a.e;
      assertBFloatValidity( r );
      // TODO: REMOVE?
      //BFloatNormalise( r );
      return r;
    }
  } else {
    if (a.e <= b.e - cBFloatConvTableSize) {
      assertBFloatValidity( b );
      return b;
    } else {
      r.f = b.f + a.f * BFloat::aConversionLookup[ b.e - a.e ];
      r.e = b.e;
      // TODO: REMOVE?
      //BFloatNormalise( r );
      return r;
    }
  }
} // bfloat_pr_sum ( BFloat const &, BFloat const & )
 
static inline void bfloat_pr_sum_accum ( BFloat& a, const BFloat& b) 
{
  //assertBFloatValidity( a );
  //assertBFloatValidity( b );
  if( b.isZero() ) {
    // Adding 0 does nothing.
    return;
  }
  if( a.isInfinity() ) {
    // Adding to infinity does nothing.
    return;
  }
  if( a.isZero() ) {
    // Adding to 0 is an assignment.
    a.e = b.e;
    a.f = b.f;
    assertBFloatValidity( a );
    return;
  }
  if( b.isInfinity() ) {
    // Adding infinity makes a infinite.
    a.e = b.e;
    a.f = b.f;
    return;
  }
  if (a.e >= b.e) {
    if (a.e < b.e + cBFloatConvTableSize)
      a.f += b.f * BFloat::aConversionLookup[ a.e - b.e ];
  } else {
    if (a.e > b.e - cBFloatConvTableSize) {
      a.f = b.f + a.f * BFloat::aConversionLookup[ b.e - a.e ];
      a.e = b.e;
    } else {
      // Assign b to a.
      a.f = b.f;
      a.e = b.e;
    }
  }
  // TODO: REMOVE?
  //BFloatNormalise( a );
} // bfloat_pr_sum_accum( BFloat &, BFloat const & )

static inline void bfloat_print ( const BFloat& x ) 
{
  bfloat_print( std::cout, x );
}
static inline std::ostream& bfloat_print ( std::ostream& out, const BFloat& x ) 
{
  static const double log10 = ::log(10.0);
  static const double maxmantisse = 10.0 * (1.0 - 0.55 * exp(-cBFloatDigits * log10));
  //out.setf(ios::fixed,ios::floatfield);
  out.precision( cBFloatDigits );
  if( x.isInfinity() ) {
    out << 1.0 << "e+Inf";
  } else if( x.isNaN() ) {
    out << "NaN";
  } else if( x.isZero() ) {
    out << "0.0";
  } else if( x.isOne() ) {
    out << "1.0";
  } else if( x.e == 0 ) {
    out << x.f;
  } else {
    double iM = (::log(x.f) + logcBFloatRange*(double)x.e) / log10;
    //long iExp = long(floor(iM));
    // PTE 16/9/08
    double iExp = floor(iM);
    iM = exp((iM - iExp) * log10);
    if (iM > maxmantisse) {
      iExp += 1;
      iM = 1.0;
    }
    out << iM << ( iExp<0 ? "e" : "e+" ) << iExp;
  }
  //out.setf(ios::fixed,ios::floatfield);  // default
  // TODO: MAGIC #
  out.precision( 6 );           // default
  return out;
} // bfloat_print ( ostream &, BFloat const & )

// forward decl.
static inline bool bfloat_equal ( const BFloat& a, const BFloat& b );

// PTE 14/9/08
static inline BFloat bfloat_pr_diff (const BFloat& a, const BFloat& b) 
{
  if( b.isZero() ) {
    // subtracting zero does nothing.
    assertBFloatValidity( a );
    return a;
  } else if( a.isZero() ) {
    // a is 0.  Negative!
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff( zero, ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    return a; // 0
  } // End if a is 0
  if( ( a.e == b.e ) && ( a.f <= b.f ) ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if(
      //( (a.e == b.e) && ( ( a.f - b.f ) <= -cBFloatEpsilon ) )
      ( (a.e == b.e) && ( ( a.f - b.f ) <= -(100.0f*cBFloatEpsilon) ) )
    ) {
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff( ";
      bfloat_print( std::cerr, a );
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
      std::cerr << "\t\t(a.e - b.e) is zero" << std::endl;
      std::cerr << "\t\t(a.f - b.f) is " << ( a.f - b.f ) << std::endl;
      return 0.0;
    } // End if it's negative ..
#endif // ALGEBRA_CHECK_NEGATIVE
    return 0.0;
  }
  if( (a.e > b.e) || ( (a.e == b.e) && ( ( a.f - b.f ) > -cBFloatEpsilon ) ) ) {
    // a positive result...
    if( a.e >= b.e + cBFloatConvTableSize ) {
      assertBFloatValidity( a );
      return a;
    } else {
      return BFloat( a.f - b.f * BFloat::aConversionLookup[ a.e - b.e ], a.e, true );
    }
  } else {
    // It could still be positive
    if( a.e > ( b.e - cBFloatConvTableSize ) ) {
      BFMantissa tmp_f = ( a.f * BFloat::aConversionLookup[ b.e - a.e ] ) - b.f;
      if( tmp_f > 0 ) {
        return BFloat( tmp_f, b.e, true );
      } // else negative
    }

    if( bfloat_equal( a, b ) ) {
      return 0.0;
    }

    // a negative result..
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff( ";
      bfloat_print( std::cerr, a );
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
      std::cerr << "\t\t(a.e - b.e) is " << ( a.e - b.e ) << std::endl;
      std::cerr << "\t\t(a.f - b.f) is " << ( a.f - b.f ) << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    return 0.0;

    //if (a.e <= b.e - cBFloatConvTableSize)
    //  return b;
    //else
    //  return BFloat( b.f - a.f * BFloat::aConversionLookup[ b.e - a.e ], b.e );
  } // End if positive .. else negative ..
} // bfloat_pr_diff ( const BFloat &, const BFloat & )
 
// PTE 14/9/08
static inline void bfloat_pr_diff_accum ( BFloat& a, const BFloat& b) 
{
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( b.isZero() ) {
    // subtracting zero does nothing.
    return;
  }
  if( a.isZero() ) {
    // a is 0.  Negative!
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff_accum( zero, ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    return;
  } // End if a is 0
  if( ( a.e == b.e ) && ( a.f <= b.f ) ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if(
      ( ( a.f - b.f ) <= -cBFloatEpsilon )
    ) {
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff_accum( ";
      bfloat_print( std::cerr, a );
      std::cerr << " [ a.f = " << a.f << ", a.e = " << a.e << " ]";
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " [ b.f = " << b.f << ", b.e = " << b.e << " ]";
      std::cerr << " )" << std::endl;
      a.f = 0.0;
      a.e = 0;
      return;
    } // End if it's negative.
#endif // ALGEBRA_CHECK_NEGATIVE
    a.f = 0.0;
    a.e = 0;
    return;
  }
  if( (a.e > b.e) || ( (a.e == b.e) && ( ( a.f - b.f ) > -cBFloatEpsilon ) ) ) {
    // a positive result...
    if( a.e < b.e + cBFloatConvTableSize ) {
      a.f -= b.f * BFloat::aConversionLookup[ a.e - b.e ];
      BFloatNormalise( a );
    }
  } else {
    // It could still be positive
    if( a.e > b.e - cBFloatConvTableSize ) {
      if( a.f * BFloat::aConversionLookup[ b.e - a.e ] > b.f ) {
        a.f = a.f * BFloat::aConversionLookup[ b.e - a.e ] - b.f;
        a.e = b.e;
        BFloatNormalise( a );
        return;
      } // else negative..
    }

    // a negative result..
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff_accum( ";
      bfloat_print( std::cerr, a );
      std::cerr << " [ a.f = " << a.f << ", a.e = " << a.e << " ]";
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " [ b.f = " << b.f << ", b.e = " << b.e << " ]";
      std::cerr << " )" << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    a.f = 0.0;
    a.e = 0;
  } // End if positive .. else negative ..
} // bfloat_pr_diff_accum ( BFloat &, const BFloat & )

static inline bool bfloat_less ( const BFloat& a, const BFloat& b) 
{
  // ASSUMPTION: a and b are non-negative.
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( a.isZero() ) {
    return ( !b.isZero() );
  }
  if( b.isZero() ) {
    return false;
  }
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize)
      return false;
    else
      return a.f < b.f * BFloat::aConversionLookup[ a.e - b.e ];
  }
  if (a.e <= b.e - cBFloatConvTableSize)
    return true;
  else
    return a.f * BFloat::aConversionLookup[ b.e - a.e ] < b.f;
} // static bfloat_less ( BFloat const &, BFloat const & )
  
static inline bool bfloat_equal ( const BFloat& a, const BFloat& b )
{
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( a.isZero() ) {
    return ( b.isZero() );
  }
  if( b.isZero() ) {
    return false;
  }
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize)
      return false;
    else
      return a.f == b.f * BFloat::aConversionLookup[ a.e - b.e ];
  }
  if (a.e <= b.e - cBFloatConvTableSize)
    return false;
  else
    return a.f * BFloat::aConversionLookup[ b.e - a.e ] == b.f;
} // static bfloat_equal( BFloat const &, BFloat const & )

static inline bool bfloat_lessequal ( const BFloat& a, const BFloat& b )
{
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( a.isZero() ) {
    return true;
  }
  if( b.isZero() ) {
    return false;
  }
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize)
      return false;
    else
      return a.f <= b.f * BFloat::aConversionLookup[ a.e - b.e ];
  }
  if (a.e <= b.e - cBFloatConvTableSize)
    return true;
  else
    return a.f * BFloat::aConversionLookup[ b.e - a.e ] <= b.f;
} // static bfloat_lessequal( BFloat const &, BFloat const & )


//
// Wrapper to allow BFloats to be used by Algebra template
//
struct BFloatMethods
{
  typedef BFloat Value;
  static inline double to_double (BFloat const & iX) { return bfloat2double(iX); }
  static inline BFloat from_double (double iP) { return double2bfloat(iP); }
  static inline BFloat pmul ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_product(iX,iY); }
  static inline BFloat pmuldouble ( BFloat const & iX, double iY) { return bfloat_pr_double_product(iX,iY); }
  // PTE 14/9/08
  static inline BFloat ppow ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_power(iX,iY); }
  // PTE 14/9/08
  static inline BFloat ppowdouble ( BFloat const & iX, double iY) { return bfloat_pr_double_power(iX,iY); }
  static inline BFloat pdiv ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_quotient(iX,iY); }
  static inline BFloat psum ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_sum(iX,iY); }
  static inline BFloat pdiff ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_diff(iX,iY); }
  static inline BFloat doubleexp ( double iX ) { return bfloat_doubleexp(iX); }
  static inline double doublelog ( BFloat const & iX ) { return bfloat_doublelog(iX); }
  static inline void pmulacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_product_accum(iX,iY); }
  static inline void pmulaccdouble ( BFloat& iX, double iY ) { bfloat_pr_double_product_accum(iX,iY); }
  static inline void pdivacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_quotient_accum(iX,iY); }
  static inline void psumacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_sum_accum(iX,iY); }
  static inline void pdiffacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_diff_accum(iX,iY); }
  static inline bool less ( BFloat const & iX, BFloat const & iY ) { return bfloat_less(iX,iY); }
  static inline bool equal ( BFloat const & iX, BFloat const & iY ) { return bfloat_equal(iX,iY); }
  static inline bool lessequal ( BFloat const & iX, BFloat const & iY ) { return bfloat_lessequal(iX,iY); }
  static inline std::ostream& print ( std::ostream& iOut, BFloat const & iX ) { return bfloat_print( iOut, iX ); }
};


//
// Wrapper to use Algebra with a double or a float.
//
// PTE 15/9/08
template <typename RealspaceType>
struct RealspaceMethods
{
  typedef RealspaceType Value;
  static inline double to_double ( Value const & iX ) { return (double)iX.x; }
  static inline Value from_double ( double const & iP ) { return Value(iP); }
  static inline Value pmul ( Value const & iX, Value const & iY ) { return iX.x*iY.x; }
  static inline Value pmuldouble ( Value const & iX, double const & iY ) { return iX.x*iY; }
  // PTE 14/9/08
  static inline Value ppow ( Value const & iX, Value const & iY ) { return ::pow(iX.x,iY.x); }
  // PTE 14/9/08
  static inline Value ppowdouble ( Value const & iX, double const & iY ) { return ::pow(iX.x,iY); }
  static inline Value pdiv ( Value const & iX, Value const & iY ) { return iX.x/iY.x; }
  static inline Value psum ( Value const & iX, Value const & iY ) { return iX.x+iY.x; }
  static inline Value pdiff ( Value const & iX, Value const & iY ) { return iX.x-iY.x; }
  static inline Value doubleexp ( double const & iX ) { return Value(exp(iX)); }
  static inline double doublelog ( Value const & iX ) { return (double)log(iX.x); }
  static inline void pmulacc ( Value& iX, Value const & iY ) { iX.x*=iY.x; }
  static inline void pmulaccdouble ( Value& iX, double const & iY ) { iX.x*=static_cast<double>(Value(iY)); }
  static inline void pdivacc ( Value& iX, Value const & iY ) { iX.x /= iY.x; }
  static inline void psumacc ( Value& iX, Value const & iY ) { iX.x += iY.x; }
  static inline void pdiffacc ( Value& iX, Value const & iY ) { iX.x -= iY.x; }
  static inline bool less ( Value const & iX, Value const & iY ) { return ( !equal( iX, iY ) && (iX.x<iY.x) ); } //{ return iX.x<iY.x; }
  // TODO: Implement a numeric_limits<DoubleRealspace> that inherits from numeric_limits<double>, and likewise for FloatRealspace, etc.
  //static inline bool equal ( Value const & iX, Value const & iY ) { return ( ( iX.x == iY.x ) || ( ( iX.x > iY.x ) ? ( ( iX.x-iY.x ) <= std::numeric_limits<RealspaceType>::epsilon() ) : ( ( iY.x-iX.x ) <= std::numeric_limits<RealspaceType>::epsilon() ) ) ); } //{ return iX.x==iY.x; }
  //static inline bool equal ( Value const & iX, Value const & iY ) { return ( ( iX.x == iY.x ) || ( ( iX.x > iY.x ) ? ( ( iX.x-iY.x ) <= std::numeric_limits<float>::epsilon() ) : ( ( iY.x-iX.x ) <= std::numeric_limits<float>::epsilon() ) ) ); } //{ return iX.x==iY.x; }
  static inline bool equal ( Value const & iX, Value const & iY ) { return ( ( iX.x == iY.x ) || ( ( iX.x > iY.x ) ? ( ( iX.x-iY.x ) <= std::numeric_limits<typename Value::XType>::epsilon() ) : ( ( iY.x-iX.x ) <= std::numeric_limits<typename Value::XType>::epsilon() ) ) ); } //{ return iX.x==iY.x; }
  static inline bool lessequal ( Value const & iX, Value const & iY ) { return iX.x<=iY.x; }
  static inline std::ostream& print ( std::ostream& iOut, Value const & iX ) { if( iX.x < 0 ) { iOut << '-'; return bfloat_print( iOut, double2bfloat( -iX.x ) ); } return bfloat_print( iOut, double2bfloat(iX.x) ); }
}; // End class RealspaceMethods

// PTE 15/9/08
class DoubleRealspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  typedef double XType;
  double x;

 public:
  DoubleRealspace () {}
  DoubleRealspace ( double x ) : x (x) {}
  template <typename other_type>
  DoubleRealspace ( other_type const & from_val )
  {
    x = toDouble( from_val );
  }
  //operator double& () { return x; }
  operator double () const { return x; }
  void clear () {
    x = 0.0;
  }
  DoubleRealspace &
  operator= ( double const & d )
  {
    x = d;
    return *this;
  } // operator=( double )
  template <typename other_type>
  DoubleRealspace &
  operator= ( other_type const & other_val )
  {
    x = toDouble( other_val );

    return *this;
  } // operator= ( other_type const & )
  
  inline double Value () const {
    return x;
  } // double Value()

  bool isZero () const
  {
    return ( x == 0.0 );
  } // isZero() const

  bool isOne () const
  {
    return ( x == 1.0 );
  } // isOne() const

  bool isInfinity () const
  {
    return std::isinf( x );
  } // isInfinity() const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN() const

}; // End class DoubleRealspace

// Added by <NAME> 15/9/08
inline static double
toDouble ( DoubleRealspace const & v )
{
  return v.x;
} // toDouble( DoubleRealspace const & )

inline static long double
toLongDouble ( DoubleRealspace const & v )
{
  return v.x;
} // toLongDouble( DoubleRealspace const & )

inline static double
toLogDouble ( DoubleRealspace const & v )
{
  return log( v.x );
} // toLogDouble( DoubleReaspace const & )

///// mark
// PTE 15/9/08
class LongDoubleRealspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  typedef long double XType;
  long double x;

 public:
  LongDoubleRealspace () {}
  LongDoubleRealspace ( double x ) : x (x) {}
  LongDoubleRealspace ( long double x ) : x (x) {}
  LongDoubleRealspace ( BFloat const & from_val )
  {
    x = toLongDouble( from_val );
  }
  template <typename other_type>
  LongDoubleRealspace ( other_type const & from_val )
  {
    x = toLongDouble( from_val );
  }
  //operator double& () { return x; }
  operator double () const { return static_cast<double>( x ); }
  operator long double () const { return x; }
  void clear () {
    x = 0.0;
  }
  LongDoubleRealspace &
  operator= ( long double const & d )
  {
    x = d;
    return *this;
  } // operator=( long double )
  template <typename other_type>
  LongDoubleRealspace &
  operator= ( other_type const & other_val )
  {
    x = toLongDouble( other_val );

    return *this;
  } // operator= ( other_type const & )
  
  inline long double Value () const {
    return x;
  } // long double Value()

  bool isZero () const
  {
    return ( x == static_cast<long double>( 0.0 ) );
  } // isZero() const

  bool isOne () const
  {
    return ( x == static_cast<long double>( 1.0 ) );
  } // isOne() const

  bool isInfinity () const
  {
    return std::isinf( x );
  } // isInfinity() const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN() const

}; // End class LongDoubleRealspace

// Added by <NAME> 15/9/08
inline static double
toDouble ( LongDoubleRealspace const & v )
{
  return static_cast<double>( v.x );
} // toDouble( LongDoubleRealspace const & )

inline static double
toLongDouble ( LongDoubleRealspace const & v )
{
  return v.x;
} // toDouble( LongDoubleRealspace const & )

inline static double
toLogDouble ( LongDoubleRealspace const & v )
{
  return log( v.x );
} // toLogDouble( DoubleReaspace const & )
///// endmark

///////////////////////
class FloatRealspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  typedef float XType;
  float x;

 public:
  FloatRealspace () {}
  FloatRealspace ( float x ) : x (x) {}
  FloatRealspace ( double x ) : x (x) {}
  template <typename other_type>
  FloatRealspace ( other_type const & from_val )
  {
    x = toDouble( from_val );
  }
  //operator double& () { return x; }
  operator double () const { return x; }
  void clear () {
    x = 0.0;
  }
  FloatRealspace &
  operator= ( float const & d )
  {
    x = d;
    return *this;
  } // operator=( float )
  FloatRealspace &
  operator= ( double const & d )
  {
    x = d;
    return *this;
  } // operator=( double )
  template <typename other_type>
  FloatRealspace &
  operator= ( other_type const & other_val )
  {
    x = toDouble( other_val );

    return *this;
  } // operator= ( other_type const & )
  
  inline double Value () const {
    return x;
  } // double Value()

  bool isZero () const
  {
    return ( x == 0.0f );
  } // isZero() const

  bool isOne () const
  {
    return ( x == 1.0f );
  } // isOne() const

  bool isInfinity () const
  {
    return std::isinf( x );
  } // isInfinity() const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN() const

}; // End class FloatRealspace

// Added by <NAME> 15/9/08
inline static double
toDouble ( FloatRealspace const & v )
{
  return v.x;
} // toDouble( FloatRealspace const & )

inline static long double
toLongDouble ( FloatRealspace const & v )
{
  return v.x;
} // toLongDouble( FloatRealspace const & )

inline static double
toLogDouble ( FloatRealspace const & v )
{
  return log( v.x );
} // toLogDouble( FloatReaspace const & )

///////////////////////
//
// Simple log-space numbers - don't use, except possibly for Viterbi
//
template <typename LogspaceType>
struct LogspaceMethods
{
  typedef LogspaceType Value;
  static inline double to_double ( Value const & iX ) {
    if( iX.isInfinity() ) {
      return std::numeric_limits<double>::infinity();
    } else if( iX.isNaN() ) {
      return std::numeric_limits<double>::quiet_NaN();
    } else if( iX.isZero() ) {
      // log( 0 ) = -inf
      return 0.0;
    } // else {
    return ::exp(iX.x);
    //}
  } // to_double( Value const & )

  static inline Value from_double ( double const & iP ) {
    if( iP == std::numeric_limits<double>::infinity() ) {
      return Value( iP );
    } else if( iP != iP /*NaN*/ ) {
      // ?
      assert( iP == iP );
      return Value( iP );
    } else if( iP == 0.0 ) {
      // log( 0 ) = -inf
      return Value( -std::numeric_limits<double>::infinity() );
    } else {
      return Value( ::log( iP ) );
    }
  } // from_double( double const & )

  static inline Value pmul ( Value const & iX, Value const & iY ) { return iX.x+iY.x; }
  static inline Value pmuldouble ( Value const & iX, double const & iY) { return iX.x+::log(iY); }
  // PTE 14/9/08
  static inline Value ppow ( Value const & iX, Value const & iY) { return iX.x*iY.x; }
  // PTE 14/9/08
  static inline Value ppowdouble ( Value const & iX, double const & iY) { return iX.x*::log(iY); }
  static inline Value pdiv ( Value const & iX, Value const & iY) { return iX.x-iY.x; }
  static inline Value psum ( Value const & iX, Value const & iY) { return logspace_add(iX.x,iY.x); }
  static inline Value pdiff ( Value const & iX, Value const & iY) { return logspace_diff(iX.x,iY.x); }
  static inline Value doubleexp ( double const & iX ) { return iX; }
  static inline double doublelog ( Value const & iX ) { return iX.x; }
  static inline void pmulacc ( Value& iX, Value const & iY) { iX.x+=iY.x; }
  static inline void pmulaccdouble ( Value& iX, double const & iY) { iX.x+=::log(iY); }
  static inline void pdivacc ( Value& iX, Value const & iY) { iX.x -= iY.x; }
  static inline void psumacc ( Value& iX, Value const & iY) { iX.x = logspace_add(iX.x,iY.x); }
  static inline void pdiffacc ( Value& iX, Value const & iY) { iX.x = logspace_diff(iX.x,iY.x); }
  static inline bool less ( Value const & iX, Value const & iY) { return iX.x<iY.x; }
  static inline bool equal ( Value const & iX, Value const & iY) { return iX.x==iY.x; }
  static inline bool lessequal ( Value const & iX, Value const & iY) { return iX.x<=iY.x; }
  static inline std::ostream& print ( std::ostream& iOut, Value const & iX ) {
    if( ( iX.x != iX.x ) /* NaN */ ) {
      iOut << iX.x;
      return iOut;
    } else if( iX.x == std::numeric_limits<double>::infinity() ) {
      iOut << iX.x;
      return iOut;
    } else if( iX.x == -std::numeric_limits<double>::infinity() ) {
      // -inf is log( 0 )
      iOut << "0.0";
      return iOut;
    }
    //iOut << iX.x;
    //return iOut;
    return bfloat_print( iOut, bfloat_doubleexp(iX.x) );
  }
};

class Logspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

  // PTE 14/9/08
 public:
  double x;

 public:
  Logspace () {}
  // TODO: ERE I AM.  Can't do this.  Need to expect arg double is already logged.
  Logspace ( double const & other_val ) {
    //operator=( other_val );
    x = other_val;
  }
  Logspace ( Logspace const & other_logspace )
  {
    x = other_logspace.x;
  }
  Logspace ( LongDoubleRealspace const & other_val )
  {
    operator=( other_val );
  } // <init>( LongDoubleRealspace const & )
  Logspace ( FloatRealspace const & other_val )
  {
    operator=( other_val );
  } // <init>( FloatRealspace const & )
  Logspace ( BFloat const & other_val )
  {
    operator=( other_val );
  } // <init>( BFloat const & )
  //operator double&(){ return x; }
  // PTE 15/9/08
  //operator double () const { return x; } //return ::exp(x); }
  void clear() {
    // Set to log( 0 ).
    x = -std::numeric_limits<double>::infinity();
  }

  Logspace &
  operator= ( Logspace const & other_logspace )
  {
    x = other_logspace.x;

    return *this;
  } // operator= ( Logspace const & )
  Logspace &
  operator= ( LongDoubleRealspace const & other_val )
  {
    return operator=( toDouble( other_val ) );
  } // operator= ( LongDoubleRealspace const & )
  Logspace &
  operator= ( FloatRealspace const & other_val )
  {
    return operator=( toDouble( other_val ) );
  } // operator= ( FloatRealspace const & )
  Logspace &
  operator= ( BFloat const & other_val )
  {
    if( other_val.isZero() ) {
      x = -std::numeric_limits<double>::infinity();
    } else if( other_val.isInfinity() ) {
      x = std::numeric_limits<double>::infinity();
    } else {
      x = toLogDouble( other_val );
    }
    return *this;
  } // operator= ( BFloat const & )
  Logspace &
  operator= ( double const & other_val )
  {
    if( other_val == std::numeric_limits<double>::infinity() ) {
      x = other_val;
    } else if( other_val != other_val /* NaN */ ) {
      // ?
      assert( other_val == other_val );
      x = other_val;
    } else if( other_val == 0.0 ) {
      // log( 0 ) = -inf
      x = -std::numeric_limits<double>::infinity();
    } else {
      x = ::log( other_val );
    }

    return *this;
  } // operator= ( double const & )

  bool isZero () const
  {
    return ( x == -std::numeric_limits<double>::infinity() );
  } // isZero () const

  bool isOne () const
  {
    return ( x == 0 );
  } // isOne () const

  bool isInfinity () const
  {
    return ( x == std::numeric_limits<double>::infinity() );
  } // isInfinity () const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN () const

  //friend double logspace_add( double iX, double iY );
  //friend double logspace_diff( double iX, double iY );
  //friend bool LogspaceMethods<Logspace>::less( Logspace iX, Logspace iY);
}; // End class Logspace

// OLD.
//inline Logspace logspace_addsmall( Logspace iX, Logspace iY ) {
//  if (iX.x - iY.x > 36.7) return iX;
//  return Logspace( iX.x + ::log(1.0+exp(iY.x-iX.x)) );
//}

// Paul changed from "iX>iY" to "iX.x>iY.x" to work around a g++ "i686-apple-darwin8-g++-4.0.1 (GCC) 4.0.1 (Apple Computer, Inc. build 5367)" complaint of an "ambiguous overload for operator>".
// PTE 15/9/08
inline double
logspace_add ( double lhs, double const & rhs )
{
  static const double log2 = ::log( 2.0 );
  if( lhs == rhs ) {
    lhs += log2;
  } else if( lhs > rhs ) {
    lhs += ::log1p(::exp(rhs-lhs));
  } else if( lhs < rhs ) {
    lhs = rhs + ::log1p(::exp(lhs-rhs));
  }
  return lhs;
} // logspace_add( double, double const & )

inline double
logspace_diff ( double lhs, double const & rhs )
{
  if( lhs <= rhs ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if(
      ( ( lhs - rhs ) < -cLogspaceEpsilon )
    ) {
      std::cerr << "Logspace: Negative number" << std::endl;
      std::cerr << "\tin logspace_diff( " << lhs << ", " << rhs << " )" << std::endl;
    }
#endif // ALGEBRA_CHECK_NEGATIVE
    return -std::numeric_limits<double>::infinity();
  }
  // assert( lhs > rhs )
  return ( lhs + ::log1p(-exp(rhs-lhs)) );
} // logspace_diff( double, double const & )

// Added by <NAME> 15/9/08
inline static double
toDouble ( Logspace const & v )
{
  return ::exp( v.x );
} // toDouble( Logspace const & )

inline static long double
toLongDouble ( Logspace const & v )
{
  return ::exp( static_cast<long double>( v.x ) );
} // toLongDouble( Logspace const & )

inline static double
toLogDouble ( Logspace const & v )
{
  return v.x;
} // toLogDouble( Logspace const & )

inline
BFloat::
  BFloat ( Logspace const & copy_from ) {
    *this = bfloat_doubleexp( copy_from.x );
  };
inline
BFloat::
  BFloat ( LongDoubleRealspace const & copy_from ) {
    *this = copy_from.x;
  };
inline
BFloat::
  BFloat ( DoubleRealspace const & copy_from ) {
    *this = copy_from.x;
  };
inline
BFloat::
  BFloat ( FloatRealspace const & copy_from ) {
    *this = copy_from.x;
  };
inline
BFloat::BFloat &
BFloat::
  operator= ( Logspace const & copy_from ) {
    *this = bfloat_doubleexp( copy_from.x );
    return *this;
  };
inline
BFloat::BFloat &
BFloat::
  operator= ( DoubleRealspace const & copy_from ) {
    *this = copy_from.x;
    return *this;
  };
inline
BFloat::BFloat &
BFloat::
  operator= ( LongDoubleRealspace const & copy_from ) {
  *this = static_cast<double>( copy_from.x );
    return *this;
  };
inline
BFloat::BFloat &
BFloat::
  operator= ( FloatRealspace const & copy_from ) {
    *this = copy_from.x;
    return *this;
  };


//
// Algebra class - Wrapper for overloading all arithmetic operators to use a different algebra.
//
// <NAME>, 19/3/03
// Based on logprob.h by by <NAME>.
//

template <class AlgebraMethods>
class Algebra {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive><commit_msg>Slight syntax modifications for compilation under C++ v. 4.5  (current Mac version is 4.2)<commit_after>/*
 *    This file is part of HMMoC-BFloat-Algebra, a tiny C++ library supplying
 *    the Algebra<> template, for polymorphic representation of non-negative
 *    floating-point values in either native (eg double), log-transformed, or
 *    BFloat representations.  BFloat is the Buoyant Float type, allowing
 *    fixed-precision representation at an arbitrary precision range.  The
 *    library is based on code from: <NAME>. HMMoC—a compiler for hidden
 *    Markov models. Bioinformatics (2007) 23(18): 2485-2487. The code was
 *    derived by <NAME> and <NAME> from <NAME>, <NAME>. An
 *    expectation maximization algorithm for training hidden substitution
 *    models.  J Mol Biol. 2002 Apr 12;317(5):753-64.  You may use at will,
 *    subject to the license (Apache v2.0), but *please cite those two papers* in
 *    your documentation and publications associated with uses of this library.
 *    Thank you!
 *
 *
 *    Copyright (C) 2011 by <NAME>, <NAME> Cancer Research
 *    Center.
 *
 *    This file is based on algebras.h in HMMoC 1.3, a hidden Markov model
 *    compiler.  Copyright (C) 2007 by <NAME>, Oxford University.
 *
 *    Licensed under the Apache License, Version 2.0 (the "License");
 *    you may not use this file except in compliance with the License.
 *    You may obtain a copy of the License at
 *    
 *        http://www.apache.org/licenses/LICENSE-2.0
 *    
 *    Unless required by applicable law or agreed to in writing, software
 *    distributed under the License is distributed on an "AS IS" BASIS,
 *    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *    See the License for the specific language governing permissions and
 *    limitations under the License.
\*/
//
// Algebra.hpp - extended real types
//
// <NAME>, June 1, 2011
// <NAME>, 27/8/04
//

#ifndef HMMOC_BFLOAT_ALGEBRA_HPP
#define HMMOC_BFLOAT_ALGEBRA_HPP

// MS compatible compilers support #pragma once
#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif

#if defined(_MSC_VER) && (_MSC_VER >= 1310)
#  pragma warning (disable : 4675) // suppress ADL warning
#endif

#include <assert.h>

// Do under- and overflow checking in BFloat? (undef for no)
#ifndef NDEBUG
#define BFLOAT_CHECK_UOFLOW
// Do extra checks in operator- to diagnose negative values? (set to 0 for no)
#define ALGEBRA_CHECK_NEGATIVE
#endif // NDEBUG

#ifndef NBOOST_SERIALIZATION // #define NBOOST_SERIALIZATION if you don't want to compile in the boost serialization library
#include <boost/serialization/nvp.hpp>
#include <boost/serialization/utility.hpp>
#include <boost/serialization/version.hpp>
#endif // NBOOST_SERIALIZATION

// TODO: PUT BACK ?  See below (search for seqan)...
//#include <seqan/basic.h>

#include <iostream>
#include <cmath> // for isinf and isnan
#include <limits>

//TAH 11/01/2011  this defines uint32_t, used below
#include <stdint.h>


// You can control some inlining recommentations to the compiler using these:
#define ALGEBRA_INLINE_ARITHMETIC inline
#define ALGEBRA_INLINE_CAST inline
#define ALGEBRA_INLINE_CAST_COMPARE inline
#define ALGEBRA_INLINE_CAST_ARITHMETIC inline

// typedefs
typedef float BFMantissa;
//typedef double BFMantissa;
const BFMantissa cBFloatRange = 20282409603651670423947251286016.0;  // 2.03e+31; 2^104
const BFMantissa cBFloatRangeInv = 1.0/cBFloatRange;
// <NAME> 6/7/7: need to typecast to avoid compiler warnings about imprecise FP representations
const BFMantissa cBFloatRangeSqrt    = (BFMantissa)1.0e+18;          // Value between square root of the exponent, and the exponent
const BFMantissa cBFloatRangeInvSqrt = (BFMantissa)1.0e-18;          // Square of this should still be representable, with full mantissa!
const BFMantissa logcBFloatRange     = ::log(cBFloatRange);
const int cBFloatDigits              = 7;                 // Number of significant digits (7 for floats, 16 for doubles?)
const int cBFloatInfinity            = 1000000000;        // Tiniest number representable is cBFloatRangeInv ^ BFloatInfinity
// PTE 15/9/08
const BFMantissa cBFloatEpsilon      = ( BFMantissa )1.193e-07;
const int cBFloatConvTableSize       = 100;               // This includes many zero entries, it makes additions a bit faster
const int cBFloatDoubleConvTableSize = 50;                // Table size for bfloat -> double conversion; cBFloatRange^(-size/2) is double 0
const double cLogspaceEpsilon        = 36.7; // TODO: Test & verify

namespace hmmoc {

// PTE 14/9/08
// To help disambiguate ::log( double ) from hmmoc::log(Algebra)
template <typename other_type>
inline
other_type
log ( other_type const & to_be_logged )
{
  return ::log( to_be_logged );
} // log( other_type const & )

// Forward decls
class Logspace;
class DoubleRealspace;
class LongDoubleRealspace;
class FloatRealspace;
class BFloat;
static inline std::ostream& bfloat_print ( std::ostream& out, const BFloat& x );
static inline void assertBFloatValidity ( BFloat const & bf );
static inline void BFloatNormalise ( BFloat& a );

//
// BFloats: more buoyant floats.
//
// NOTE: For now BFloats can represent only positive values.
//
// struct{ float + int } is 8 bytes; nice size makes noticable speed difference
//
// For BFloats, infinity is represented as ( e == cBFloatInfinity, and we never let f be std::numeric_limits<BFMantissa>::infinity() ).  0 is represented as ( f == e == 0 ), and NaN as ( f != f ).
class BFloat {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  inline
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( f );
    ar & BOOST_SERIALIZATION_NVP( e );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  static BFMantissa* aConversionLookup;              // used by addition
  static double* aDoubleConversionLookup;            // used by Value()
  BFMantissa f;
  int e;
 public:
  inline
  BFloat ( BFMantissa const & iF, int const & iE, bool normalise = false ) :
    f( iF ), e( iE )
  {
    if( normalise ) {
      BFloatNormalise( *this );
    } else {
      assertBFloatValidity( *this );
    }
  };

  inline
  BFloat () : f( 0.0 ), e( 0 ) {};

  inline
  BFloat ( double const & prob ) { *this = prob; }

  inline
  BFloat ( BFloat const & copy_from ) : f( copy_from.f ), e( copy_from.e ) {
    assertBFloatValidity( *this );
  };

  // Defined later, after these others are defined.
  inline
  BFloat ( Logspace const & copy_from );
  // Defined later, after these others are defined.
  inline
  BFloat ( DoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline
  BFloat ( LongDoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline
  BFloat ( FloatRealspace const & copy_from );

  inline
  ~BFloat () {};

  inline
  BFloat &
  operator= ( BFloat const & copy_from )
  {
    f = copy_from.f;
    e = copy_from.e;

    return *this;
  } // operator=( BFloat const & )

  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( Logspace const & copy_from );
  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( DoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( LongDoubleRealspace const & copy_from );
  // Defined later, after these others are defined.
  inline BFloat &
  operator= ( FloatRealspace const & copy_from );

  inline BFloat &
  operator= ( double prob )
  {
    if( prob != prob ) {
      // NaN
      // For now, don't let it be set to NaN.
      assert( prob == prob ); // NOT NaN
      f = prob;
      e = 0;
    } else
    if( prob == std::numeric_limits<double>::infinity() ) {
      f = 1.0;
      e = cBFloatInfinity;
    } else
    // Simplistic double-to-BFloat conversion - can be slow if 'standard' numbers get very large/small
    if( prob < std::numeric_limits<double>::min() ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
      if( prob < 0.0 ) {
        std::cerr << "BFloat: Negative number: " << prob << std::endl;
        // TODO: ?
        assert( prob >= 0.0 );
      }
#endif
      f = 0.0;
      e = 0;
    } else {
      f = 0.0;
      e = 0;
      while( prob > cBFloatRangeSqrt ) {
        prob *= cBFloatRangeInv;
        e++;
      }
      while( prob < cBFloatRangeInvSqrt ) {
        prob *= cBFloatRange;
        e--;
      }
      f = prob; 
    }
    assertBFloatValidity( *this );
    return *this;
  } // operator=( double const & )

  inline double Value () const { 
    if( e == 0 ) {
      return f;
    }
    if( e == cBFloatInfinity ) {
      return std::numeric_limits<double>::infinity();
    }
    if( f != f ) {
      return std::numeric_limits<double>::quiet_NaN();
    }
    if (abs(e) < cBFloatDoubleConvTableSize/2) {
      return (double)f * aDoubleConversionLookup[ e + cBFloatDoubleConvTableSize/2 ];
    } else if (e < cBFloatDoubleConvTableSize/2) {
      return 0.0;
    } else {
      return (double)f * exp((double)e * logcBFloatRange);
    }
  } // double Value()

  // sets to zero
  inline
  void clear () {
    f=0.0; e=0;
  }

  inline
  bool isZero () const
  {
    return ( f == 0.0 );
  } // isZero() const

  inline
  bool isOne () const
  {
    return ( f == 1.0 );
  } // isOne() const

  inline
  bool isInfinity () const
  {
    return ( e == cBFloatInfinity );
  } // isInfinity() const

  inline
  bool isNaN () const
  {
    return ( f != f );
  } // isNaN() const

}; // End class BFloat


//
// dummy class to initialise BFloat lookup table
//
class _BFloatInitialize {
public:
  _BFloatInitialize();
};


//
// implementations of BFloat calculations
//

static inline void assertBFloatValidity ( BFloat const & bf )
{
  // Never negative
  assert( bf.f >= 0.0 );

  // TODO: REMOVE?  NaN should be ok to represent, but we might want to raise a red flag...
  assert( !bf.isNaN() );

  // we represent infinity as ( e == cBFloatInfinity ), not as ( f == std::numeric_limits<BFMantissa>::infinity() ).
  assert( bf.f != std::numeric_limits<BFMantissa>::infinity() );

  if( bf.isZero() ) {
    // 0 is represented as both e and f being 0
    assert( bf.e == 0 );
    //} else if( !bf.isInfinity() ) {
    // f must be in this range (AFTER NORMALIZATION)
    //assert( bf.f <= cBFloatRangeSqrt ); 
    //assert( bf.f >= cBFloatRangeInvSqrt ); 
  }
    
  return;
} // static assertBFloatValidity( BFloat const & )

// Normalization of BFloat result of a single operation
#ifdef BFLOAT_CHECK_UOFLOW
static inline void BFloatNormalise ( BFloat& a )
     //#define BFloatNormalise(a)
{
  if (a.f > cBFloatRangeSqrt) {
    a.f *= cBFloatRangeInv;
    a.e++;
  } else if (a.f < cBFloatRangeInvSqrt) {
    if( a.f == 0.0 ) {
      a.e = 0; 
    } else {
      a.f *= cBFloatRange;
      a.e--;
    }
  }
  if( a.e > cBFloatInfinity ) {
    std::cerr << "BFloat: Overflow" << std::endl;
    a.e = cBFloatInfinity;
  } else if( a.e < -cBFloatInfinity ) {
    std::cerr << "BFloat: Underflow" << std::endl;
    a.e = 0;
    a.f = 0.0;
  }
  // TODO: REMOVE
#ifdef ALGEBRA_CHECK_NEGATIVE
    if( a.f < 0 ) {
      std::cerr << "BFloat: Negative number: ";
      bfloat_print( std::cerr, a );
      std::cerr << std::endl;
      std::cerr << "\tin BFloatNormalise( BFMantissa & ) [BFLOAT_CHECK_UOFLOW is true]" << std::endl;
    }
#endif // ALGEBRA_CHECK_NEGATIVE
#ifndef NDEBUG
  assertBFloatValidity( a );
#endif // NDEBUG
}; // BFloatNormalise ( BFloat & ) // when BFLOAT_CHECK_UOFLOW
#else // if BFLOAT_CHECK_UOFLOW .. else ..
static inline void BFloatNormDown ( BFloat& a ) { 
  a.f *= cBFloatRangeInv;
  a.e++;
  assert( a.f <= cBFloatRangeSqrt );
}
static inline void BFloatNormUp ( BFloat& a ) { 
  if (a.f == 0.0) {
    a.e = 0;
  } else {
    a.f *= cBFloatRange;
    a.e--;
    assert( a.f >= cBFloatRangeInvSqrt );
  }
}
static inline void BFloatNormalise ( BFloat& a )
     //#define BFloatNormalise(a) 
{
  if (a.f > cBFloatRangeSqrt) {
    BFloatNormDown(a);
  } else if (a.f < cBFloatRangeInvSqrt) {
    BFloatNormUp(a);
  }
#ifdef ALGEBRA_CHECK_NEGATIVE
    if( a.f < 0 ) {
      std::cerr << "BFloat: Negative number: ";
      bfloat_print( std::cerr, a );
      std::cerr << std::endl;
      std::cerr << "\tin BFloatNormalise( BFMantissa & ) [BFLOAT_CHECK_UOFLOW is false]" << std::endl;
    }
#endif // ALGEBRA_CHECK_NEGATIVE
#ifndef NDEBUG
  assertBFloatValidity( a );
#endif // NDEBUG
};  // BFloatNormalise ( BFloat & ) // when NOT BFLOAT_CHECK_UOFLOW
#endif // End if BFLOAT_CHECK_UOFLOW .. else ..

static inline void DoubleNormalise ( double& f, int& e )
{
  // comparing to 0.0 here fails, because the comparison is done
  // using higher-precision doubles, but the subsequent while-loop
  // uses true doubles, resulting in an infinite loop. (G.L. 3/9/07)
  if( f < std::numeric_limits<double>::min() ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if( f < 0.0 ) {
      std::cerr << "BFloat: Negative number: " << f << std::endl;
      // TODO: ?
      assert( f >= 0.0 );
    }
#endif // ALGEBRA_CHECK_NEGATIVE
    f = 0.0; 
    e = 0;
  } else {
    while (f > (double)cBFloatRangeSqrt) {
      f *= (double)cBFloatRangeInv;
      e++;
    }
    while (f < (double)cBFloatRangeInvSqrt) {
      f *= (double)cBFloatRange;
      e--;
    }
  }
 }; // DoubleNormalise( double &, int & )

// Logarithm of a BFloat
static inline double bfloat_doublelog ( const BFloat& a ) { return a.e*logcBFloatRange+::log(a.f); }

// BFloat exp of a double -- could be a teensy tiny value.
static inline BFloat bfloat_doubleexp ( double iA ) 
{
  if( iA == std::numeric_limits<double>::infinity() ) {
    // exp( inf ) is inf
    return BFloat( 1.0, cBFloatInfinity );
  }
  if( iA == -std::numeric_limits<double>::infinity() ) {
    // exp( -inf ) is 0
    return 0.0;
  }
  int iE = (int)floor( iA / ::log(cBFloatRange) );
  iA -= iE * ::log(cBFloatRange);
  return BFloat( ::exp(iA), iE, true );
} // static bfloat_doubleexp( double )

// Returns a double value - or underflow/overflow if it does not fit.
static inline double bfloat2double ( const BFloat & bfloat) { return bfloat.Value(); }

/// TODO: PUT BACK
//template <typename T>
//inline static double
//toDouble ( T const & v )
//{
//  return seqan::convert<T, double>( v );
//} // toDouble( T const & )

inline static double
toDouble ( BFloat const & v )
{
  return bfloat2double( v );
} // toDouble( BFloat const & )

inline static long double
toLongDouble ( BFloat const & v )
{
  return bfloat2double( v );
} // toLongDouble( BFloat const & )

inline static double
toLogDouble ( BFloat const & v )
{
  return bfloat_doublelog( v );
} // toLogDouble( BFloat const & )

static inline BFloat double2bfloat ( double prob )
{
  return prob;
}

static inline BFloat bfloat_pr_product ( const BFloat& a, const BFloat& b )
{ 
  //assertBFloatValidity( a );
  //assertBFloatValidity( b );
  if( a.isZero() ) {
    return 0.0;
  }
  if( b.isZero() ) {
    return 0.0;
  }
  return BFloat( a.f*b.f, a.e+b.e, true );
} // static bfloat_pr_product( BFloat const &, BFloat const & )

static inline BFloat bfloat_pr_double_product ( const BFloat& a, double const & b )
{ 
  //assertBFloatValidity( a );
  if( a.isZero() ) {
    return 0.0;
  }
  if( b == 0.0 ) {
    return 0.0;
  }
  register double mantisse = a.f*b;
  int exponent = a.e;
  DoubleNormalise(mantisse, exponent);
  return BFloat( mantisse, exponent, false );
} // bfloat_pr_double_product ( BFloat const &, double const & )

static inline void bfloat_pr_product_accum ( BFloat& a, const BFloat& b )
{
  if( a.isZero() ) {
    // It's already 0.
    return;
  }
  if( b.isZero() ) {
    // Multiplying by 0..
    a.clear();
    return;
  }
  a.f *= b.f; a.e += b.e; 
  BFloatNormalise( a );
} // static bfloat_pr_product_accum( BFloat &, BFloat const & )

static inline void bfloat_pr_double_product_accum (  BFloat& a, double const & b )
{ 
  if( a.isZero() ) {
    // It's already 0.
    return;
  }
  if( b == 0.0 ) {
    // Multiplying by 0..
    a.clear();
    return;
  }
  register double mantisse = a.f*b;
  DoubleNormalise( mantisse, a.e );
  a.f = mantisse;
  assertBFloatValidity( a );
  return;
} // static bfloat_pr_double_product_accum( BFloat &, double const & )

// PTE 14/9/08
static inline BFloat bfloat_pr_power (const BFloat& a, const BFloat& b) 
{ 
  return BFloat( ::pow(a.f,b.f), (a.e*b.e), true );
}

static inline BFloat bfloat_pr_double_power ( const BFloat& a, double const & b )
{
  // TODO: test, make sure it's safe.  Probably it's not.
  register double mantisse = ::pow(a.f, b);
  int exponent = a.e;
  DoubleNormalise(mantisse, exponent);
  return BFloat(mantisse, exponent);
}

static inline BFloat bfloat_pr_quotient ( const BFloat& a, const BFloat& b )
{
  if( b.isZero() ) {
    return std::numeric_limits<double>::quiet_NaN();
  } else if( b.isInfinity() ) {
    return 0.0;
  }
  return BFloat( a.f/b.f, a.e-b.e, true );
} // bfloat_pr_product ( BFloat const &, BFloat const & )
  
static inline void bfloat_pr_quotient_accum ( BFloat& a, const BFloat& b ) 
{ 
  if( b.isZero() ) {
    a = std::numeric_limits<double>::quiet_NaN();
  } else if( b.isInfinity() ) {
    a = 0.0;
  } else {
    a.f /= b.f; 
    a.e -= b.e;
  }
  BFloatNormalise( a ); 
} // bfloat_pr_quotient_accum( BFloat &, BFloat const & )

static inline BFloat bfloat_pr_sum ( const BFloat& a, const BFloat& b )
{
  if( b.isZero() ) {
    assertBFloatValidity( a );
    return a;
  }
  if( a.isZero() ) {
    assertBFloatValidity( b );
    return b;
  }
  if( b.isInfinity() || a.isInfinity() ) {
    return std::numeric_limits<double>::infinity();
  }

  BFloat r;
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize) {
      assertBFloatValidity( a );
      return a;
    } else {
      r.f = a.f + b.f * BFloat::aConversionLookup[ a.e - b.e ];
      r.e = a.e;
      assertBFloatValidity( r );
      // TODO: REMOVE?
      //BFloatNormalise( r );
      return r;
    }
  } else {
    if (a.e <= b.e - cBFloatConvTableSize) {
      assertBFloatValidity( b );
      return b;
    } else {
      r.f = b.f + a.f * BFloat::aConversionLookup[ b.e - a.e ];
      r.e = b.e;
      // TODO: REMOVE?
      //BFloatNormalise( r );
      return r;
    }
  }
} // bfloat_pr_sum ( BFloat const &, BFloat const & )
 
static inline void bfloat_pr_sum_accum ( BFloat& a, const BFloat& b) 
{
  //assertBFloatValidity( a );
  //assertBFloatValidity( b );
  if( b.isZero() ) {
    // Adding 0 does nothing.
    return;
  }
  if( a.isInfinity() ) {
    // Adding to infinity does nothing.
    return;
  }
  if( a.isZero() ) {
    // Adding to 0 is an assignment.
    a.e = b.e;
    a.f = b.f;
    assertBFloatValidity( a );
    return;
  }
  if( b.isInfinity() ) {
    // Adding infinity makes a infinite.
    a.e = b.e;
    a.f = b.f;
    return;
  }
  if (a.e >= b.e) {
    if (a.e < b.e + cBFloatConvTableSize)
      a.f += b.f * BFloat::aConversionLookup[ a.e - b.e ];
  } else {
    if (a.e > b.e - cBFloatConvTableSize) {
      a.f = b.f + a.f * BFloat::aConversionLookup[ b.e - a.e ];
      a.e = b.e;
    } else {
      // Assign b to a.
      a.f = b.f;
      a.e = b.e;
    }
  }
  // TODO: REMOVE?
  //BFloatNormalise( a );
} // bfloat_pr_sum_accum( BFloat &, BFloat const & )

static inline void bfloat_print ( const BFloat& x ) 
{
  bfloat_print( std::cout, x );
}
static inline std::ostream& bfloat_print ( std::ostream& out, const BFloat& x ) 
{
  static const double log10 = ::log(10.0);
  static const double maxmantisse = 10.0 * (1.0 - 0.55 * exp(-cBFloatDigits * log10));
  //out.setf(ios::fixed,ios::floatfield);
  out.precision( cBFloatDigits );
  if( x.isInfinity() ) {
    out << 1.0 << "e+Inf";
  } else if( x.isNaN() ) {
    out << "NaN";
  } else if( x.isZero() ) {
    out << "0.0";
  } else if( x.isOne() ) {
    out << "1.0";
  } else if( x.e == 0 ) {
    out << x.f;
  } else {
    double iM = (::log(x.f) + logcBFloatRange*(double)x.e) / log10;
    //long iExp = long(floor(iM));
    // PTE 16/9/08
    double iExp = floor(iM);
    iM = exp((iM - iExp) * log10);
    if (iM > maxmantisse) {
      iExp += 1;
      iM = 1.0;
    }
    out << iM << ( iExp<0 ? "e" : "e+" ) << iExp;
  }
  //out.setf(ios::fixed,ios::floatfield);  // default
  // TODO: MAGIC #
  out.precision( 6 );           // default
  return out;
} // bfloat_print ( ostream &, BFloat const & )

// forward decl.
static inline bool bfloat_equal ( const BFloat& a, const BFloat& b );

// PTE 14/9/08
static inline BFloat bfloat_pr_diff (const BFloat& a, const BFloat& b) 
{
  if( b.isZero() ) {
    // subtracting zero does nothing.
    assertBFloatValidity( a );
    return a;
  } else if( a.isZero() ) {
    // a is 0.  Negative!
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff( zero, ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    return a; // 0
  } // End if a is 0
  if( ( a.e == b.e ) && ( a.f <= b.f ) ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if(
      //( (a.e == b.e) && ( ( a.f - b.f ) <= -cBFloatEpsilon ) )
      ( (a.e == b.e) && ( ( a.f - b.f ) <= -(100.0f*cBFloatEpsilon) ) )
    ) {
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff( ";
      bfloat_print( std::cerr, a );
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
      std::cerr << "\t\t(a.e - b.e) is zero" << std::endl;
      std::cerr << "\t\t(a.f - b.f) is " << ( a.f - b.f ) << std::endl;
      return 0.0;
    } // End if it's negative ..
#endif // ALGEBRA_CHECK_NEGATIVE
    return 0.0;
  }
  if( (a.e > b.e) || ( (a.e == b.e) && ( ( a.f - b.f ) > -cBFloatEpsilon ) ) ) {
    // a positive result...
    if( a.e >= b.e + cBFloatConvTableSize ) {
      assertBFloatValidity( a );
      return a;
    } else {
      return BFloat( a.f - b.f * BFloat::aConversionLookup[ a.e - b.e ], a.e, true );
    }
  } else {
    // It could still be positive
    if( a.e > ( b.e - cBFloatConvTableSize ) ) {
      BFMantissa tmp_f = ( a.f * BFloat::aConversionLookup[ b.e - a.e ] ) - b.f;
      if( tmp_f > 0 ) {
        return BFloat( tmp_f, b.e, true );
      } // else negative
    }

    if( bfloat_equal( a, b ) ) {
      return 0.0;
    }

    // a negative result..
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff( ";
      bfloat_print( std::cerr, a );
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
      std::cerr << "\t\t(a.e - b.e) is " << ( a.e - b.e ) << std::endl;
      std::cerr << "\t\t(a.f - b.f) is " << ( a.f - b.f ) << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    return 0.0;

    //if (a.e <= b.e - cBFloatConvTableSize)
    //  return b;
    //else
    //  return BFloat( b.f - a.f * BFloat::aConversionLookup[ b.e - a.e ], b.e );
  } // End if positive .. else negative ..
} // bfloat_pr_diff ( const BFloat &, const BFloat & )
 
// PTE 14/9/08
static inline void bfloat_pr_diff_accum ( BFloat& a, const BFloat& b) 
{
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( b.isZero() ) {
    // subtracting zero does nothing.
    return;
  }
  if( a.isZero() ) {
    // a is 0.  Negative!
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff_accum( zero, ";
      bfloat_print( std::cerr, b );
      std::cerr << " )" << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    return;
  } // End if a is 0
  if( ( a.e == b.e ) && ( a.f <= b.f ) ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if(
      ( ( a.f - b.f ) <= -cBFloatEpsilon )
    ) {
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff_accum( ";
      bfloat_print( std::cerr, a );
      std::cerr << " [ a.f = " << a.f << ", a.e = " << a.e << " ]";
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " [ b.f = " << b.f << ", b.e = " << b.e << " ]";
      std::cerr << " )" << std::endl;
      a.f = 0.0;
      a.e = 0;
      return;
    } // End if it's negative.
#endif // ALGEBRA_CHECK_NEGATIVE
    a.f = 0.0;
    a.e = 0;
    return;
  }
  if( (a.e > b.e) || ( (a.e == b.e) && ( ( a.f - b.f ) > -cBFloatEpsilon ) ) ) {
    // a positive result...
    if( a.e < b.e + cBFloatConvTableSize ) {
      a.f -= b.f * BFloat::aConversionLookup[ a.e - b.e ];
      BFloatNormalise( a );
    }
  } else {
    // It could still be positive
    if( a.e > b.e - cBFloatConvTableSize ) {
      if( a.f * BFloat::aConversionLookup[ b.e - a.e ] > b.f ) {
        a.f = a.f * BFloat::aConversionLookup[ b.e - a.e ] - b.f;
        a.e = b.e;
        BFloatNormalise( a );
        return;
      } // else negative..
    }

    // a negative result..
#ifdef ALGEBRA_CHECK_NEGATIVE
      std::cerr << "BFloat: Negative number" << std::endl;
      std::cerr << "\tin bfloat_pr_diff_accum( ";
      bfloat_print( std::cerr, a );
      std::cerr << " [ a.f = " << a.f << ", a.e = " << a.e << " ]";
      std::cerr << ", ";
      bfloat_print( std::cerr, b );
      std::cerr << " [ b.f = " << b.f << ", b.e = " << b.e << " ]";
      std::cerr << " )" << std::endl;
#endif // ALGEBRA_CHECK_NEGATIVE
    a.f = 0.0;
    a.e = 0;
  } // End if positive .. else negative ..
} // bfloat_pr_diff_accum ( BFloat &, const BFloat & )

static inline bool bfloat_less ( const BFloat& a, const BFloat& b) 
{
  // ASSUMPTION: a and b are non-negative.
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( a.isZero() ) {
    return ( !b.isZero() );
  }
  if( b.isZero() ) {
    return false;
  }
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize)
      return false;
    else
      return a.f < b.f * BFloat::aConversionLookup[ a.e - b.e ];
  }
  if (a.e <= b.e - cBFloatConvTableSize)
    return true;
  else
    return a.f * BFloat::aConversionLookup[ b.e - a.e ] < b.f;
} // static bfloat_less ( BFloat const &, BFloat const & )
  
static inline bool bfloat_equal ( const BFloat& a, const BFloat& b )
{
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( a.isZero() ) {
    return ( b.isZero() );
  }
  if( b.isZero() ) {
    return false;
  }
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize)
      return false;
    else
      return a.f == b.f * BFloat::aConversionLookup[ a.e - b.e ];
  }
  if (a.e <= b.e - cBFloatConvTableSize)
    return false;
  else
    return a.f * BFloat::aConversionLookup[ b.e - a.e ] == b.f;
} // static bfloat_equal( BFloat const &, BFloat const & )

static inline bool bfloat_lessequal ( const BFloat& a, const BFloat& b )
{
  assertBFloatValidity( a );
  assertBFloatValidity( b );
  if( a.isZero() ) {
    return true;
  }
  if( b.isZero() ) {
    return false;
  }
  if (a.e > b.e) {
    if (a.e >= b.e + cBFloatConvTableSize)
      return false;
    else
      return a.f <= b.f * BFloat::aConversionLookup[ a.e - b.e ];
  }
  if (a.e <= b.e - cBFloatConvTableSize)
    return true;
  else
    return a.f * BFloat::aConversionLookup[ b.e - a.e ] <= b.f;
} // static bfloat_lessequal( BFloat const &, BFloat const & )


//
// Wrapper to allow BFloats to be used by Algebra template
//
struct BFloatMethods
{
  typedef BFloat Value;
  static inline double to_double (BFloat const & iX) { return bfloat2double(iX); }
  static inline BFloat from_double (double iP) { return double2bfloat(iP); }
  static inline BFloat pmul ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_product(iX,iY); }
  static inline BFloat pmuldouble ( BFloat const & iX, double iY) { return bfloat_pr_double_product(iX,iY); }
  // PTE 14/9/08
  static inline BFloat ppow ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_power(iX,iY); }
  // PTE 14/9/08
  static inline BFloat ppowdouble ( BFloat const & iX, double iY) { return bfloat_pr_double_power(iX,iY); }
  static inline BFloat pdiv ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_quotient(iX,iY); }
  static inline BFloat psum ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_sum(iX,iY); }
  static inline BFloat pdiff ( BFloat const & iX, BFloat const & iY) { return bfloat_pr_diff(iX,iY); }
  static inline BFloat doubleexp ( double iX ) { return bfloat_doubleexp(iX); }
  static inline double doublelog ( BFloat const & iX ) { return bfloat_doublelog(iX); }
  static inline void pmulacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_product_accum(iX,iY); }
  static inline void pmulaccdouble ( BFloat& iX, double iY ) { bfloat_pr_double_product_accum(iX,iY); }
  static inline void pdivacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_quotient_accum(iX,iY); }
  static inline void psumacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_sum_accum(iX,iY); }
  static inline void pdiffacc ( BFloat& iX, BFloat const & iY ) { bfloat_pr_diff_accum(iX,iY); }
  static inline bool less ( BFloat const & iX, BFloat const & iY ) { return bfloat_less(iX,iY); }
  static inline bool equal ( BFloat const & iX, BFloat const & iY ) { return bfloat_equal(iX,iY); }
  static inline bool lessequal ( BFloat const & iX, BFloat const & iY ) { return bfloat_lessequal(iX,iY); }
  static inline std::ostream& print ( std::ostream& iOut, BFloat const & iX ) { return bfloat_print( iOut, iX ); }
};


//
// Wrapper to use Algebra with a double or a float.
//
// PTE 15/9/08
template <typename RealspaceType>
struct RealspaceMethods
{
  typedef RealspaceType Value;
  static inline double to_double ( Value const & iX ) { return (double)iX.x; }
  static inline Value from_double ( double const & iP ) { return Value(iP); }
  static inline Value pmul ( Value const & iX, Value const & iY ) { return iX.x*iY.x; }
  static inline Value pmuldouble ( Value const & iX, double const & iY ) { return iX.x*iY; }
  // PTE 14/9/08
  static inline Value ppow ( Value const & iX, Value const & iY ) { return ::pow(iX.x,iY.x); }
  // PTE 14/9/08
  static inline Value ppowdouble ( Value const & iX, double const & iY ) { return ::pow(iX.x,iY); }
  static inline Value pdiv ( Value const & iX, Value const & iY ) { return iX.x/iY.x; }
  static inline Value psum ( Value const & iX, Value const & iY ) { return iX.x+iY.x; }
  static inline Value pdiff ( Value const & iX, Value const & iY ) { return iX.x-iY.x; }
  static inline Value doubleexp ( double const & iX ) { return Value(exp(iX)); }
  static inline double doublelog ( Value const & iX ) { return (double)log(iX.x); }
  static inline void pmulacc ( Value& iX, Value const & iY ) { iX.x*=iY.x; }
  static inline void pmulaccdouble ( Value& iX, double const & iY ) { iX.x*=static_cast<double>(Value(iY)); }
  static inline void pdivacc ( Value& iX, Value const & iY ) { iX.x /= iY.x; }
  static inline void psumacc ( Value& iX, Value const & iY ) { iX.x += iY.x; }
  static inline void pdiffacc ( Value& iX, Value const & iY ) { iX.x -= iY.x; }
  static inline bool less ( Value const & iX, Value const & iY ) { return ( !equal( iX, iY ) && (iX.x<iY.x) ); } //{ return iX.x<iY.x; }
  // TODO: Implement a numeric_limits<DoubleRealspace> that inherits from numeric_limits<double>, and likewise for FloatRealspace, etc.
  //static inline bool equal ( Value const & iX, Value const & iY ) { return ( ( iX.x == iY.x ) || ( ( iX.x > iY.x ) ? ( ( iX.x-iY.x ) <= std::numeric_limits<RealspaceType>::epsilon() ) : ( ( iY.x-iX.x ) <= std::numeric_limits<RealspaceType>::epsilon() ) ) ); } //{ return iX.x==iY.x; }
  //static inline bool equal ( Value const & iX, Value const & iY ) { return ( ( iX.x == iY.x ) || ( ( iX.x > iY.x ) ? ( ( iX.x-iY.x ) <= std::numeric_limits<float>::epsilon() ) : ( ( iY.x-iX.x ) <= std::numeric_limits<float>::epsilon() ) ) ); } //{ return iX.x==iY.x; }
  static inline bool equal ( Value const & iX, Value const & iY ) { return ( ( iX.x == iY.x ) || ( ( iX.x > iY.x ) ? ( ( iX.x-iY.x ) <= std::numeric_limits<typename Value::XType>::epsilon() ) : ( ( iY.x-iX.x ) <= std::numeric_limits<typename Value::XType>::epsilon() ) ) ); } //{ return iX.x==iY.x; }
  static inline bool lessequal ( Value const & iX, Value const & iY ) { return iX.x<=iY.x; }
  static inline std::ostream& print ( std::ostream& iOut, Value const & iX ) { if( iX.x < 0 ) { iOut << '-'; return bfloat_print( iOut, double2bfloat( -iX.x ) ); } return bfloat_print( iOut, double2bfloat(iX.x) ); }
}; // End class RealspaceMethods

// PTE 15/9/08
class DoubleRealspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  typedef double XType;
  double x;

 public:
  DoubleRealspace () {}
  DoubleRealspace ( double x ) : x (x) {}
  template <typename other_type>
  DoubleRealspace ( other_type const & from_val )
  {
    x = toDouble( from_val );
  }
  //operator double& () { return x; }
  operator double () const { return x; }
  void clear () {
    x = 0.0;
  }
  DoubleRealspace &
  operator= ( double const & d )
  {
    x = d;
    return *this;
  } // operator=( double )
  template <typename other_type>
  DoubleRealspace &
  operator= ( other_type const & other_val )
  {
    x = toDouble( other_val );

    return *this;
  } // operator= ( other_type const & )
  
  inline double Value () const {
    return x;
  } // double Value()

  bool isZero () const
  {
    return ( x == 0.0 );
  } // isZero() const

  bool isOne () const
  {
    return ( x == 1.0 );
  } // isOne() const

  bool isInfinity () const
  {
    return std::isinf( x );
  } // isInfinity() const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN() const

}; // End class DoubleRealspace

// Added by <NAME> 15/9/08
inline static double
toDouble ( DoubleRealspace const & v )
{
  return v.x;
} // toDouble( DoubleRealspace const & )

inline static long double
toLongDouble ( DoubleRealspace const & v )
{
  return v.x;
} // toLongDouble( DoubleRealspace const & )

inline static double
toLogDouble ( DoubleRealspace const & v )
{
  return log( v.x );
} // toLogDouble( DoubleReaspace const & )

///// mark
// PTE 15/9/08
class LongDoubleRealspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  typedef long double XType;
  long double x;

 public:
  LongDoubleRealspace () {}
  LongDoubleRealspace ( double x ) : x (x) {}
  LongDoubleRealspace ( long double x ) : x (x) {}
  LongDoubleRealspace ( BFloat const & from_val )
  {
    x = toLongDouble( from_val );
  }
  template <typename other_type>
  LongDoubleRealspace ( other_type const & from_val )
  {
    x = toLongDouble( from_val );
  }
  //operator double& () { return x; }
  operator double () const { return static_cast<double>( x ); }
  operator long double () const { return x; }
  void clear () {
    x = 0.0;
  }
  LongDoubleRealspace &
  operator= ( long double const & d )
  {
    x = d;
    return *this;
  } // operator=( long double )
  template <typename other_type>
  LongDoubleRealspace &
  operator= ( other_type const & other_val )
  {
    x = toLongDouble( other_val );

    return *this;
  } // operator= ( other_type const & )
  
  inline long double Value () const {
    return x;
  } // long double Value()

  bool isZero () const
  {
    return ( x == static_cast<long double>( 0.0 ) );
  } // isZero() const

  bool isOne () const
  {
    return ( x == static_cast<long double>( 1.0 ) );
  } // isOne() const

  bool isInfinity () const
  {
    return std::isinf( x );
  } // isInfinity() const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN() const

}; // End class LongDoubleRealspace

// Added by <NAME> 15/9/08
inline static double
toDouble ( LongDoubleRealspace const & v )
{
  return static_cast<double>( v.x );
} // toDouble( LongDoubleRealspace const & )

inline static double
toLongDouble ( LongDoubleRealspace const & v )
{
  return v.x;
} // toDouble( LongDoubleRealspace const & )

inline static double
toLogDouble ( LongDoubleRealspace const & v )
{
  return log( v.x );
} // toLogDouble( DoubleReaspace const & )
///// endmark

///////////////////////
class FloatRealspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

 public:
  typedef float XType;
  float x;

 public:
  FloatRealspace () {}
  FloatRealspace ( float x ) : x (x) {}
  FloatRealspace ( double x ) : x (x) {}
  template <typename other_type>
  FloatRealspace ( other_type const & from_val )
  {
    x = toDouble( from_val );
  }
  //operator double& () { return x; }
  operator double () const { return x; }
  void clear () {
    x = 0.0;
  }
  FloatRealspace &
  operator= ( float const & d )
  {
    x = d;
    return *this;
  } // operator=( float )
  FloatRealspace &
  operator= ( double const & d )
  {
    x = d;
    return *this;
  } // operator=( double )
  template <typename other_type>
  FloatRealspace &
  operator= ( other_type const & other_val )
  {
    x = toDouble( other_val );

    return *this;
  } // operator= ( other_type const & )
  
  inline double Value () const {
    return x;
  } // double Value()

  bool isZero () const
  {
    return ( x == 0.0f );
  } // isZero() const

  bool isOne () const
  {
    return ( x == 1.0f );
  } // isOne() const

  bool isInfinity () const
  {
    return std::isinf( x );
  } // isInfinity() const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN() const

}; // End class FloatRealspace

// Added by <NAME> 15/9/08
inline static double
toDouble ( FloatRealspace const & v )
{
  return v.x;
} // toDouble( FloatRealspace const & )

inline static long double
toLongDouble ( FloatRealspace const & v )
{
  return v.x;
} // toLongDouble( FloatRealspace const & )

inline static double
toLogDouble ( FloatRealspace const & v )
{
  return log( v.x );
} // toLogDouble( FloatReaspace const & )

///////////////////////
//
// Simple log-space numbers - don't use, except possibly for Viterbi
//
template <typename LogspaceType>
struct LogspaceMethods
{
  typedef LogspaceType Value;
  static inline double to_double ( Value const & iX ) {
    if( iX.isInfinity() ) {
      return std::numeric_limits<double>::infinity();
    } else if( iX.isNaN() ) {
      return std::numeric_limits<double>::quiet_NaN();
    } else if( iX.isZero() ) {
      // log( 0 ) = -inf
      return 0.0;
    } // else {
    return ::exp(iX.x);
    //}
  } // to_double( Value const & )

  static inline Value from_double ( double const & iP ) {
    if( iP == std::numeric_limits<double>::infinity() ) {
      return Value( iP );
    } else if( iP != iP /*NaN*/ ) {
      // ?
      assert( iP == iP );
      return Value( iP );
    } else if( iP == 0.0 ) {
      // log( 0 ) = -inf
      return Value( -std::numeric_limits<double>::infinity() );
    } else {
      return Value( ::log( iP ) );
    }
  } // from_double( double const & )

  static inline Value pmul ( Value const & iX, Value const & iY ) { return iX.x+iY.x; }
  static inline Value pmuldouble ( Value const & iX, double const & iY) { return iX.x+::log(iY); }
  // PTE 14/9/08
  static inline Value ppow ( Value const & iX, Value const & iY) { return iX.x*iY.x; }
  // PTE 14/9/08
  static inline Value ppowdouble ( Value const & iX, double const & iY) { return iX.x*::log(iY); }
  static inline Value pdiv ( Value const & iX, Value const & iY) { return iX.x-iY.x; }
  static inline Value psum ( Value const & iX, Value const & iY) { return logspace_add(iX.x,iY.x); }
  static inline Value pdiff ( Value const & iX, Value const & iY) { return logspace_diff(iX.x,iY.x); }
  static inline Value doubleexp ( double const & iX ) { return iX; }
  static inline double doublelog ( Value const & iX ) { return iX.x; }
  static inline void pmulacc ( Value& iX, Value const & iY) { iX.x+=iY.x; }
  static inline void pmulaccdouble ( Value& iX, double const & iY) { iX.x+=::log(iY); }
  static inline void pdivacc ( Value& iX, Value const & iY) { iX.x -= iY.x; }
  static inline void psumacc ( Value& iX, Value const & iY) { iX.x = logspace_add(iX.x,iY.x); }
  static inline void pdiffacc ( Value& iX, Value const & iY) { iX.x = logspace_diff(iX.x,iY.x); }
  static inline bool less ( Value const & iX, Value const & iY) { return iX.x<iY.x; }
  static inline bool equal ( Value const & iX, Value const & iY) { return iX.x==iY.x; }
  static inline bool lessequal ( Value const & iX, Value const & iY) { return iX.x<=iY.x; }
  static inline std::ostream& print ( std::ostream& iOut, Value const & iX ) {
    if( ( iX.x != iX.x ) /* NaN */ ) {
      iOut << iX.x;
      return iOut;
    } else if( iX.x == std::numeric_limits<double>::infinity() ) {
      iOut << iX.x;
      return iOut;
    } else if( iX.x == -std::numeric_limits<double>::infinity() ) {
      // -inf is log( 0 )
      iOut << "0.0";
      return iOut;
    }
    //iOut << iX.x;
    //return iOut;
    return bfloat_print( iOut, bfloat_doubleexp(iX.x) );
  }
};

class Logspace {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize ( Archive & ar, const unsigned int /* file_version */ )
  {
    ar & BOOST_SERIALIZATION_NVP( x );
  } // serialize( Archive &, const unsigned int )
#endif // NBOOST_SERIALIZATION

  // PTE 14/9/08
 public:
  double x;

 public:
  Logspace () {}
  // TODO: ERE I AM.  Can't do this.  Need to expect arg double is already logged.
  Logspace ( double const & other_val ) {
    //operator=( other_val );
    x = other_val;
  }
  Logspace ( Logspace const & other_logspace )
  {
    x = other_logspace.x;
  }
  Logspace ( LongDoubleRealspace const & other_val )
  {
    operator=( other_val );
  } // <init>( LongDoubleRealspace const & )
  Logspace ( FloatRealspace const & other_val )
  {
    operator=( other_val );
  } // <init>( FloatRealspace const & )
  Logspace ( BFloat const & other_val )
  {
    operator=( other_val );
  } // <init>( BFloat const & )
  //operator double&(){ return x; }
  // PTE 15/9/08
  //operator double () const { return x; } //return ::exp(x); }
  void clear() {
    // Set to log( 0 ).
    x = -std::numeric_limits<double>::infinity();
  }

  Logspace &
  operator= ( Logspace const & other_logspace )
  {
    x = other_logspace.x;

    return *this;
  } // operator= ( Logspace const & )
  Logspace &
  operator= ( LongDoubleRealspace const & other_val )
  {
    return operator=( toDouble( other_val ) );
  } // operator= ( LongDoubleRealspace const & )
  Logspace &
  operator= ( FloatRealspace const & other_val )
  {
    return operator=( toDouble( other_val ) );
  } // operator= ( FloatRealspace const & )
  Logspace &
  operator= ( BFloat const & other_val )
  {
    if( other_val.isZero() ) {
      x = -std::numeric_limits<double>::infinity();
    } else if( other_val.isInfinity() ) {
      x = std::numeric_limits<double>::infinity();
    } else {
      x = toLogDouble( other_val );
    }
    return *this;
  } // operator= ( BFloat const & )
  Logspace &
  operator= ( double const & other_val )
  {
    if( other_val == std::numeric_limits<double>::infinity() ) {
      x = other_val;
    } else if( other_val != other_val /* NaN */ ) {
      // ?
      assert( other_val == other_val );
      x = other_val;
    } else if( other_val == 0.0 ) {
      // log( 0 ) = -inf
      x = -std::numeric_limits<double>::infinity();
    } else {
      x = ::log( other_val );
    }

    return *this;
  } // operator= ( double const & )

  bool isZero () const
  {
    return ( x == -std::numeric_limits<double>::infinity() );
  } // isZero () const

  bool isOne () const
  {
    return ( x == 0 );
  } // isOne () const

  bool isInfinity () const
  {
    return ( x == std::numeric_limits<double>::infinity() );
  } // isInfinity () const

  bool isNaN () const
  {
    return ( x != x );
  } // isNaN () const

  //friend double logspace_add( double iX, double iY );
  //friend double logspace_diff( double iX, double iY );
  //friend bool LogspaceMethods<Logspace>::less( Logspace iX, Logspace iY);
}; // End class Logspace

// OLD.
//inline Logspace logspace_addsmall( Logspace iX, Logspace iY ) {
//  if (iX.x - iY.x > 36.7) return iX;
//  return Logspace( iX.x + ::log(1.0+exp(iY.x-iX.x)) );
//}

// Paul changed from "iX>iY" to "iX.x>iY.x" to work around a g++ "i686-apple-darwin8-g++-4.0.1 (GCC) 4.0.1 (Apple Computer, Inc. build 5367)" complaint of an "ambiguous overload for operator>".
// PTE 15/9/08
inline double
logspace_add ( double lhs, double const & rhs )
{
  static const double log2 = ::log( 2.0 );
  if( lhs == rhs ) {
    lhs += log2;
  } else if( lhs > rhs ) {
    lhs += ::log1p(::exp(rhs-lhs));
  } else if( lhs < rhs ) {
    lhs = rhs + ::log1p(::exp(lhs-rhs));
  }
  return lhs;
} // logspace_add( double, double const & )

inline double
logspace_diff ( double lhs, double const & rhs )
{
  if( lhs <= rhs ) {
#ifdef ALGEBRA_CHECK_NEGATIVE
    if(
      ( ( lhs - rhs ) < -cLogspaceEpsilon )
    ) {
      std::cerr << "Logspace: Negative number" << std::endl;
      std::cerr << "\tin logspace_diff( " << lhs << ", " << rhs << " )" << std::endl;
    }
#endif // ALGEBRA_CHECK_NEGATIVE
    return -std::numeric_limits<double>::infinity();
  }
  // assert( lhs > rhs )
  return ( lhs + ::log1p(-exp(rhs-lhs)) );
} // logspace_diff( double, double const & )

// Added by <NAME> 15/9/08
inline static double
toDouble ( Logspace const & v )
{
  return ::exp( v.x );
} // toDouble( Logspace const & )

inline static long double
toLongDouble ( Logspace const & v )
{
  return ::exp( static_cast<long double>( v.x ) );
} // toLongDouble( Logspace const & )

inline static double
toLogDouble ( Logspace const & v )
{
  return v.x;
} // toLogDouble( Logspace const & )

inline
BFloat::
  BFloat ( Logspace const & copy_from ) {
    *this = bfloat_doubleexp( copy_from.x );
  };
inline
BFloat::
  BFloat ( LongDoubleRealspace const & copy_from ) {
    *this = copy_from.x;
  };
inline
BFloat::
  BFloat ( DoubleRealspace const & copy_from ) {
    *this = copy_from.x;
  };
inline
BFloat::
  BFloat ( FloatRealspace const & copy_from ) {
    *this = copy_from.x;
  };

//TAH 11/01/11 Modified next four functions to compile with g++ v>4.5
inline
BFloat &
  BFloat::operator= ( Logspace const & copy_from ) {
    *this = bfloat_doubleexp( copy_from.x );
    return *this;
  };

inline
BFloat &
  BFloat::operator= ( DoubleRealspace const & copy_from ) {
    *this = copy_from.x;
    return *this;
  };
inline
BFloat &
  BFloat::operator= ( LongDoubleRealspace const & copy_from ) {
    *this = static_cast<double>( copy_from.x );
    return *this;
  };
inline
BFloat &
  BFloat::operator= ( FloatRealspace const & copy_from ) {
    *this = copy_from.x;
    return *this;
  };


//
// Algebra class - Wrapper for overloading all arithmetic operators to use a different algebra.
//
// <NAME>, 19/3/03
// Based on logprob.h by by <NAME>.
//

template <class AlgebraMethods>
class Algebra {
// Boost serialization
#ifndef NBOOST_SERIALIZATION
private:
  friend class boost::serialization::access;
  template<class Archive>
{'fraction_non_alphanumeric': 0.14634657501231257, 'fraction_numerical': 0.010353029937860863, 'mean_word_length': 2.6913317603372917, 'pattern_counts': {'":': 0, '<': 698, '<?xml version=': 0, '>': 269, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3128042', 'n_tokens_mistral': 35572, 'n_tokens_neox': 34473, 'n_words': 9470}