plProbability32.h

/*=============================================================================
 * Product        : OpenPL 
 * File           : plProbability32.h
 * Author         : Kamel Mekhnacha
 * Creation       : September 2001
 * 
 *=============================================================================
 *     (c) Copyright 2000, Centre National de la Recherche Scientifique,
 *                              all rights reserved
 *=============================================================================
 *
 *------------------------- Description ---------------------------------------
 *   Implements a new arithmetics to handle probabilities.
 *   
 *   
 *-----------------------------------------------------------------------------
 */

#ifndef _plProbability_h_
#define _plProbability_h_


#include <plError.h>
#include <iostream>
#include <math.h>
#include <limits.h>
#include <errno.h>



#define PL_F_SIZE 31

#define PL_MASK_32 0x80000000
#define PL_MASK_31 0x40000000
#define PL_MASK_30 0x20000000
#define PL_MASK_29 0x10000000

#define PL_MAX_31_BITS 0x3FFFFFFF
#define PL_MIN_31_BITS (-0x3FFFFFFF)

#define PL_BIT_32(v) ( (v) & PL_MASK_32)
#define PL_BIT_31(v) ( (v) & PL_MASK_31)

#define PL_LOG_BASE2(v) (log(v)/0.69314718055994530942)

#ifndef PL_LOG_OF_10_BASE2
#define PL_LOG_OF_10_BASE2 3.321928094887362660
#endif

#ifndef PL_LOG_OF_E_BASE2
#define PL_LOG_OF_E_BASE2 1.442695040888963
#endif



/*****************************************************************************
 *                           Classes Definitions
 *****************************************************************************/



class plProbability
{
protected:
  unsigned long f;
  long e;

protected:

public:
  inline plProbability();

  inline plProbability(const long double &f);

  inline plProbability(const double &f);
  
  inline plProbability(const float &f);

  inline plProbability(const int &i);

  inline plProbability(const unsigned int &i);

  inline plProbability( unsigned long f, long e);

  inline operator float() const;

  inline operator double() const;

  inline operator long double() const;

  inline  plProbability operator+(const plProbability& pc) const;
  inline  plProbability operator+(const long double& fl) const;
  inline  plProbability operator+(const double& fl) const;
  inline  plProbability operator+(const float& fl) const;

  inline  plProbability operator-(const plProbability& pc) const;
  inline  plProbability operator-(const long double& fl) const;
  inline  plProbability operator-(const double& fl) const;
  inline  plProbability operator-(const float& fl) const;

  inline  plProbability& operator+=(const plProbability& pc);
  inline  plProbability& operator+=(const long double& fl);
  inline  plProbability& operator+=(const double& fl);
  inline  plProbability& operator+=(const float& fl);

  inline  plProbability& operator-=(const plProbability& pc);
  inline  plProbability& operator-=(const long double& fl);
  inline  plProbability& operator-=(const double& fl);
  inline  plProbability& operator-=(const float& fl);
  
  inline  plProbability operator*(const plProbability& pc) const;
  inline  plProbability operator*(const long double& fl) const;
  inline  plProbability operator*(const double& fl) const;
  inline  plProbability operator*(const float& fl) const;
  
  inline  plProbability& operator*=(const plProbability& pc);
  inline  plProbability& operator*=(const long double& fl);
  inline  plProbability& operator*=(const double& fl);
  inline  plProbability& operator*=(const float& fl);

  inline  plProbability operator/(const plProbability& pc) const;
  inline  plProbability operator/(const long double& fl) const;
  inline  plProbability operator/(const double& fl) const;
  inline  plProbability operator/(const float& fl) const;

  inline  plProbability& operator/=(const plProbability& pc);
  inline  plProbability& operator/=(const long double& fl);
  inline  plProbability& operator/=(const double& fl);
  inline  plProbability& operator/=(const float& fl);
    
  inline bool operator==(const plProbability& p) const; 
  inline bool operator==(const long double&) const;
  inline bool operator==(const double&) const;
  inline bool operator==(const float&) const;
  
  inline bool operator!=(const plProbability& p) const; 
  inline bool operator!=(const long double&) const;
  inline bool operator!=(const double&) const;
  inline bool operator!=(const float&) const;
  

  inline bool operator<=(const plProbability& p) const;
  inline bool operator<=(const long double&) const;
  inline bool operator<=(const double&) const;
  inline bool operator<=(const float&) const;
  
  inline bool operator>=(const plProbability& p) const;
  inline bool operator>=(const long double&) const;
  inline bool operator>=(const double&) const;
  inline bool operator>=(const float&) const;

  inline bool operator<(const plProbability& p) const;
  inline bool operator<(const long double&) const;
  inline bool operator<(const double&) const;
  inline bool operator<(const float&) const;
  
  inline bool operator>(const plProbability& p) const;
  inline bool operator>(const long double&) const;
  inline bool operator>(const double&) const;
  inline bool operator>(const float&) const;

  inline bool isNaN() const;

  inline bool isInfinity() const;
  
  friend inline double log2(const plProbability &pr);
  
  friend inline double log(const plProbability &pr);
  
  friend inline double log10(const plProbability &pr);
  
  friend inline plProbability ppow( const plProbability &pr, long exp );
  
 
  // Friend oprerators
  
  friend inline std::ostream& operator<<(std::ostream& out, const plProbability& p);
  friend inline std::istream& operator>>(std::istream& out, plProbability& p);
  
  
  friend inline long double operator*(const long double& theFloat, const plProbability& p);
  friend inline double operator*(const double& theFloat, const plProbability& p);
  friend inline float operator*(const float& theFloat, const plProbability& p);

  friend inline long double operator/(const long double& theFloat, const plProbability& p);
  friend inline double operator/(const double& theFloat, const plProbability& p);
  friend inline float operator/(const float& theFloat, const plProbability& p);
  

  // Defining some constants
  static const plProbability plZeroProbaConstant;
  static const plProbability plOneProbaConstant;
  static const plProbability plTwoProbaConstant;
  static const plProbability plHalfProbaConstant;

  static const plProbability plInfinityProbaConstant;
  static const plProbability plNaNProbaConstant;
};

//=============================================================================
// Returns "true" is the value is "NaN" and "false" otherwise.
inline bool plProbability::isNaN() const
{
    return *this == plNaNProbaConstant;
}

//=============================================================================
// Returns "true" is the value is "Infinity" and "false" otherwise.
bool plProbability::isInfinity() const
{
    return (e == LONG_MAX) && (f != 0);
}

#ifdef PL_NOT_IEEE754
// We use a generic (but slow) method!!!
//=============================================================================
inline plProbability::plProbability(const long double &fl)
{
  if(fl == 0.0L){
    *this = plZeroProbaConstant;
  }
  else{
    e = (long)(PL_LOG_BASE2(fl))+1;
    f = (unsigned long) (fl*pow(2.0, PL_F_SIZE - (double)e ));

    if(PL_BIT_32(f)){
      f >>=1; e++;
    }
    else{
      if( (f!=0) && (!PL_BIT_31(f)) ) {
    while(!PL_BIT_31(f)) {
      f <<=1; e--;
    }
      }
    }
  }
}

//=============================================================================
inline plProbability::plProbability(const double &fl)
{
  if(fl == 0.0){
    *this = plZeroProbaConstant;
  }
  else{
    e = (long)(PL_LOG_BASE2(fl))+1;
    f = (unsigned long) (fl*pow(2.0, PL_F_SIZE - (double)e ));

    if(PL_BIT_32(f)){
      f >>=1; e++;
    }
    else{
      if( (f!=0) && (!PL_BIT_31(f)) ) {
    while(!PL_BIT_31(f)) {
      f <<=1; e--;
    }
      }
    }
  }
}

//=============================================================================
inline plProbability::plProbability(const float &fl)
{
  if(fl == 0.0f){
    *this = plZeroProbaConstant;
  }
  else{
    e = (long)(PL_LOG_BASE2(fl))+1;
    f = (unsigned long) (fl*pow(2.0, PL_F_SIZE - (double)e ));

    if(PL_BIT_32(f)){
      f >>=1; e++;
    }
    else{
      if( (f!=0) && (!PL_BIT_31(f)) ) {
    while(!PL_BIT_31(f)) {
      f <<=1; e--;
    }
      }
    }
  }
}


// Conversion to float
//=============================================================================
inline plProbability::operator float() const
{
  if (*this == plZeroProbaConstant) return 0.0f;

  return  f * pow(2.0, ((double)e - PL_F_SIZE));
}

// Conversion to double
//=============================================================================
inline plProbability::operator double() const
{
  if (*this == plZeroProbaConstant) return 0.0;

  return  f * pow(2.0, ((double)e - PL_F_SIZE));
}

// Conversion to long double
//=============================================================================
inline plProbability::operator long double() const
{
  if (*this == plZeroProbaConstant) return 0.0L;

  return  f * pow(2.0, ((double)e - PL_F_SIZE));
}

#else // PL_NOT_IEEE754 is defined
// We use a precific method for IEEE754 machines 


// Masks and macros for long double conversion
#define PL_MASK_E_LD  0x7FFF
#define PL_GET_E_LD(f) (f & PL_MASK_E_LD)

// Masks and macros for double conversion
#define PL_MASK_E_D  0x7FF00000
#define PL_MASK_FH_D 0x000FFFFF
#define PL_GET_E_D(f) ((f & PL_MASK_E_D) >> 20)
#define PL_GET_FH_D(f)(f & PL_MASK_FH_D)

// Masks and macros for float conversion
#define PL_MASK_E_S  0x7F800000
#define PL_MASK_FH_S  0x007FFFFF
#define PL_GET_E_S(f) ((f & PL_MASK_E_S) >> 23)
#define PL_GET_FH_S(f)(f & PL_MASK_FH_S)

typedef struct {unsigned long w1, w2;} TwoWords;
typedef union {TwoWords longPart; double doublePart;} DoubleUnion;


typedef union { unsigned long longPart; float floatPart;} FloatUnion;

// Converting a "long double" to a plProbability
//=============================================================================
inline plProbability::plProbability(const long double &ld)
{
  if(ld == 0.0L){
    *this = plZeroProbaConstant;
  }
  else{
    
    DoubleUnion dCopy;
    dCopy.doublePart = double(ld);
    unsigned long dh, dl;
#ifdef __BIG_ENDIAN__ //macintosh
    dh = dCopy.longPart.w1;
    dl = dCopy.longPart.w2;
#else
    dh = dCopy.longPart.w2;
    dl = dCopy.longPart.w1;
#endif
    f = ( (PL_GET_FH_D(dh)  | 0x00100000) << 10) | (dl >> 22) ;
    e = PL_GET_E_D(dh) - 1022;
  }
}

// Converting a "double" to a plProbability
// This is a hidden-bit representation!!!
//=============================================================================
inline plProbability::plProbability(const double &fl)
{
  if(fl == 0.0){
    *this = plZeroProbaConstant;
  }
  else{
    DoubleUnion dCopy;
    dCopy.doublePart = fl;
    unsigned long dh, dl;
#ifdef __BIG_ENDIAN__ //macintosh
    dh = dCopy.longPart.w1;
    dl = dCopy.longPart.w2;
#else
    dh = dCopy.longPart.w2;
    dl = dCopy.longPart.w1;
#endif
    f = ( (PL_GET_FH_D(dh)  | 0x00100000) << 10) | (dl >> 22) ;
    e = PL_GET_E_D(dh) - 1022;
  }
}

// Converting a "float" to a plProbability
// This is a hidden-bit representation!!!
//=============================================================================
inline plProbability::plProbability(const float &fl)
{
  if(fl == 0.0f){
    *this = plZeroProbaConstant;
  }
  else{
    FloatUnion flCopy;
    flCopy.floatPart = fl;
    unsigned long af = flCopy.longPart;
    f =  ((PL_GET_FH_S(af) | 0x00800000) << 7);
    e = PL_GET_E_S(af)- 126;
  }
}



// Conversion to float
//=============================================================================
inline plProbability::operator float() const
{
  if(f == 0) return 0.0f;
  if(e < -126) return 0.0f;
  if(e > 0x000000FF) return HUGE_VAL;
  
  FloatUnion flCopy;
  unsigned long fe = ( (e+126) << 23 );
  unsigned long ff = ( (f & 0x3FFFFFFF) >> 7);
  unsigned long fint = ( (fe | ff) & 0x7FFFFFFF);
  flCopy.longPart = fint;
  
  return flCopy.floatPart;
}


// Conversion to double
//=============================================================================
inline plProbability::operator double() const
{
  if(f == 0) return 0.0;
  if(e < -1022) return 0.0;
  if(e > 0x000007FF) return HUGE_VAL;
  
  DoubleUnion dCopy;
  unsigned long dh =  ( ((e + 1022) << 20) | ( (f & 0x3FFFFC00) >> 10) ) & 0x7FFFFFFF;
  unsigned long dl =  (f & 0x000003FF) << 22;
  
  
#ifdef __BIG_ENDIAN__ //macintosh
  dCopy.longPart.w1 = dh;
  dCopy.longPart.w2 = dl;
#else
  dCopy.longPart.w2 = dh;
  dCopy.longPart.w1 = dl;
#endif
  
  return dCopy.doublePart;
}

// Conversion to long double
//=============================================================================
inline plProbability::operator long double() const
{
  return (long double) (double(*this));
}


#endif // PL_NOT_IEEE754

//=============================================================================
inline plProbability::plProbability(const int &i)
{
  *this = double(i);
}

//=============================================================================
inline plProbability::plProbability(const unsigned int &i)
{
  *this = double(i);
}

//=============================================================================
inline plProbability::plProbability()
{
}

//=============================================================================
inline plProbability::plProbability(unsigned long _f,  long _e)
{
  f = _f; 
  e = _e; 
}




// Comparison operators
//=============================================================================
inline bool plProbability::operator==(const plProbability& p) const
{
  return ( (f==p.f) && (e==p.e) );
}

//=============================================================================
inline bool plProbability::operator==(const long double& fl) const 
{
  return ( (*this) == (plProbability)fl ); 
} 


//=============================================================================
inline bool plProbability::operator==(const double& fl) const 
{
  return ( (*this) == (plProbability)fl ); 
} 


//=============================================================================
inline bool plProbability::operator==(const float& fl) const 
{
  return ( (*this) == (plProbability)fl ); 
} 


//=============================================================================
inline bool plProbability::operator!=(const plProbability& p) const
{
  return (!(*this == p));
}

//=============================================================================
inline bool plProbability::operator!=(const long double& fl) const 
{ 
  return ( (*this) != (plProbability)fl ); 
} 

//=============================================================================
inline bool plProbability::operator!=(const double& fl) const 
{ 
  return ( (*this) != (plProbability)fl ); 
} 

//=============================================================================
inline bool plProbability::operator!=(const float& fl) const 
{ 
  return ( (*this) != (plProbability)fl ); 
} 




//=============================================================================
inline bool plProbability::operator<=(const plProbability& p) const
{
  if(f == 0) return true;
  if(p.f == 0) return false;

  if(e < p.e) return true;
  if(e == p.e) return (f <= p.f);
  return false;
}

//=============================================================================
inline bool plProbability::operator<=(const long double& fl) const 
{
  return ( (*this) <= (plProbability)fl ); 
}


//=============================================================================
inline bool plProbability::operator<=(const double& fl) const 
{
  return ( (*this) <= (plProbability)fl ); 
}

//=============================================================================
inline bool plProbability::operator<=(const float& fl) const 
{
  return ( (*this) <= (plProbability)fl ); 
}


//=============================================================================
inline bool plProbability::operator<(const plProbability& p) const
{
  if(f == 0) return (p.f != 0);
  if(p.f == 0) return false;
  if(e < p.e) return true;
  if(e == p.e) return (f < p.f);
  return false;
}

//=============================================================================
inline bool plProbability::operator<(const long double& fl) const 
{
  return ( (*this) < (plProbability)fl );
} 

//=============================================================================
inline bool plProbability::operator<(const double& fl) const 
{
  return ( (*this) < (plProbability)fl );
} 

//=============================================================================
inline bool plProbability::operator<(const float& fl) const 
{
  return ( (*this) < (plProbability)fl );
} 



//=============================================================================
inline bool plProbability::operator>=(const plProbability& p) const 
{
  if(p.f == 0) return true;
  if(f == 0) return false;
  if(e > p.e) return true;
  if(e == p.e) return (f >= p.f);
  return false;
}

//=============================================================================
inline bool plProbability::operator>=(const long double& fl) const
{ 
  return ( (*this) >= (plProbability)fl );
}


//=============================================================================
inline bool plProbability::operator>=(const double& fl) const
{ 
  return ( (*this) >= (plProbability)fl );
}


//=============================================================================
inline bool plProbability::operator>=(const float& fl) const
{ 
  return ( (*this) >= (plProbability)fl );
}


//=============================================================================
inline bool plProbability::operator>(const plProbability& p) const 
{ 
  if(p.f == 0) return (f != 0);
  if(f == 0) return false;
  if(e > p.e) return true;
  if(e == p.e) return (f > p.f);
  return false;
}

//=============================================================================
inline bool plProbability::operator>(const long double& fl) const 
{
  return ( (*this) > (plProbability)fl ); 
}

//=============================================================================
inline bool plProbability::operator>(const double& fl) const 
{
  return ( (*this) > (plProbability)fl ); 
}

//=============================================================================
inline bool plProbability::operator>(const float& fl) const 
{
  return ( (*this) > (plProbability)fl ); 
}


//=============================================================================
inline  plProbability plProbability::operator+(const plProbability& pc) const
{ 
  if(isNaN()) return plNaNProbaConstant;
  if(pc.isNaN()) return plNaNProbaConstant;
  
  if(isInfinity()) return plInfinityProbaConstant;
  if(pc.isInfinity()) return plInfinityProbaConstant;
  
  long first_exponent_minus_second = e - pc.e;
  if(first_exponent_minus_second >= PL_F_SIZE) return *this;
  if(first_exponent_minus_second <= -PL_F_SIZE) return pc;
   
   plProbability pa;
  
  if (first_exponent_minus_second < 0){
    pa.f = (f >> (-first_exponent_minus_second )) + pc.f;
    pa.e = pc.e;
  }
  else{
    pa.f = f + (pc.f >> first_exponent_minus_second);
    pa.e = e;
  }
  
  if(PL_BIT_32(pa.f)){
    pa.f >>= 1;
    pa.e++;
    if( (pa.e > PL_MAX_31_BITS)) return plInfinityProbaConstant;
  }
  
  return pa;
}

//=============================================================================
inline  plProbability plProbability::operator+(const long double& fl) const
{
  return (*this) + (plProbability)fl;
}


//=============================================================================
inline  plProbability plProbability::operator+(const double& fl) const
{
  return (*this) + (plProbability)fl;
}

//=============================================================================
inline  plProbability plProbability::operator+(const float& fl) const
{
  return (*this) + (plProbability)fl;
}

//=============================================================================
inline  plProbability& plProbability::operator+=(const plProbability& pc)
{
  if(isNaN()) {*this = plNaNProbaConstant; return *this;}
  if(pc.isNaN()) {*this = plNaNProbaConstant; return *this;}
  
  if(isInfinity()) {*this = plInfinityProbaConstant; return *this;}
  if(pc.isInfinity()) {*this = plInfinityProbaConstant; return *this;}
  
  long first_exponent_minus_second = e - pc.e;

  if(first_exponent_minus_second >= PL_F_SIZE) return *this;

  if(first_exponent_minus_second <= -PL_F_SIZE){
    f = pc.f; 
    e = pc.e;
    return *this;
  }
  
  if (first_exponent_minus_second < 0){ 
    f = (f >> (-first_exponent_minus_second )) + pc.f;
    e = pc.e;
  }
  else{
    f += (pc.f >> first_exponent_minus_second);
  }
  
  if(PL_BIT_32(f)){
    f >>= 1;
    e++;
    if( (e > PL_MAX_31_BITS)) *this = plInfinityProbaConstant;
  }
  return *this;
}

//=============================================================================
inline  plProbability& plProbability::operator+=(const long double& fl)
{
  return (*this)+= ((plProbability)fl); 
}

//=============================================================================
inline  plProbability& plProbability::operator+=(const double& fl)
{
  return (*this)+= ((plProbability)fl); 
}

//=============================================================================
inline  plProbability& plProbability::operator+=(const float& fl)
{
  return (*this)+= ((plProbability)fl); 
}

//=============================================================================
inline  plProbability plProbability::operator-(const plProbability& pc) const
{
  if(isNaN()) return plNaNProbaConstant;
  if(pc.isNaN()) return plNaNProbaConstant;
  
  if(isInfinity() && !pc.isInfinity()) return plInfinityProbaConstant;
  if(pc.isInfinity()) return plNaNProbaConstant;
  
  if(*this == pc) return plZeroProbaConstant;
  
  if( *this < pc){
    plWarning(34);
    return plNaNProbaConstant;
  }

  
  long first_exponent_minus_second = e - pc.e;
  
  if(first_exponent_minus_second >= PL_F_SIZE) return *this;
  
  plProbability pa;
  pa.f = f - (pc.f >> first_exponent_minus_second);
  pa.e = e;
  
  while(!PL_BIT_31(pa.f)){
    pa.f <<= 1;
    pa.e--;
  }
  
  return pa;
}

//=============================================================================
inline  plProbability plProbability::operator-(const long double& fl) const 
{
  return (*this) - (plProbability)fl;
}

//=============================================================================
inline  plProbability plProbability::operator-(const double& fl) const 
{
  return (*this) - (plProbability)fl;
}

//=============================================================================
inline  plProbability plProbability::operator-(const float& fl) const 
{
  return (*this) - (plProbability)fl;
}

//=============================================================================
inline  plProbability& plProbability::operator-=(const plProbability& pc)
{
  if(isNaN()) return *this;
  if(pc.isNaN()) {*this = pc; return *this;}
   
  if(isInfinity() && !pc.isInfinity() ) {*this = plInfinityProbaConstant; return *this;}
  if(pc.isInfinity()) {*this = plNaNProbaConstant; return *this;}
  
  if(*this == pc) {*this = plZeroProbaConstant; return *this;}
  
  if( *this < pc){
    plWarning(34);
    *this = plNaNProbaConstant;
  }

  
  long first_exponent_minus_second = e - pc.e;
  
  if(first_exponent_minus_second >= PL_F_SIZE) return *this;
  
  f -= (pc.f >> first_exponent_minus_second);
  
  while(!PL_BIT_31(f)){
    f <<= 1;
    e--;
  }
  
  return *this;
}

//=============================================================================
inline  plProbability& plProbability::operator-=(const long double& fl)
{
  return (*this)-= ((plProbability)fl); 
}

//=============================================================================
inline  plProbability& plProbability::operator-=(const double& fl)
{
  return (*this)-= ((plProbability)fl); 
}

//=============================================================================
inline  plProbability& plProbability::operator-=(const float& fl)
{
  return (*this)-= ((plProbability)fl); 
}


//=============================================================================
inline  plProbability plProbability::operator*(const plProbability& pc) const
{
  if(isNaN()) return plNaNProbaConstant;
  if(pc.isNaN()) return plNaNProbaConstant;
  
  if(isInfinity()){
    if(pc.f == 0) return plNaNProbaConstant;
    else return plInfinityProbaConstant;
  }
  
  if(pc.isInfinity()){
    if(f == 0) return plNaNProbaConstant;
    else return plInfinityProbaConstant;
  }
  
  plProbability pa;
  //pa.f = (f >> 15) * (pc.f >> 15);
  pa.f = ((unsigned long long)f * (unsigned long long)pc.f) >> 30;

  pa.e = e + pc.e - 1 ;

  if(PL_BIT_32(pa.f)){
    pa.f >>= 1;
    pa.e++;
  }

  if(pa.e > PL_MAX_31_BITS) return plInfinityProbaConstant;
  
  if(pa.e < PL_MIN_31_BITS) return plZeroProbaConstant;
 
  
  return pa;
}

//=============================================================================
inline  plProbability plProbability::operator*(const long double& fl) const 
{
  return (*this) * (plProbability)fl;
}

//=============================================================================
inline  plProbability plProbability::operator*(const double& fl) const 
{
  return (*this) * (plProbability)fl;
}

//=============================================================================
inline  plProbability plProbability::operator*(const float& fl) const 
{
  return (*this) * (plProbability)fl;
}


//=============================================================================
inline long double operator*(const long double& theFloat, const plProbability& p)
{
  return (theFloat * (long double)p);
}


//=============================================================================
inline double operator*(const double& theFloat, const plProbability& p)
{
  return (theFloat * double(p));
}

//=============================================================================
inline float operator*(const float& theFloat, const plProbability& p)
{
  return (theFloat * float(p));
}

//=============================================================================
inline  plProbability& plProbability::operator*=(const plProbability& pc)
{
  if(isNaN()) {*this = plNaNProbaConstant; return *this;}
  if(pc.isNaN()) {*this = plNaNProbaConstant; return *this;}
  
  if(isInfinity()){
    if(pc.f == 0) {*this = plNaNProbaConstant; return *this;}
    else {*this = plInfinityProbaConstant; return *this;}
  }
  
  if(pc.isInfinity()){
    if(f == 0) {*this = plNaNProbaConstant; return *this;}
    else {*this = plInfinityProbaConstant; return *this;}
  }
  
  
  //f>>=15;
  //f *= (pc.f>>15);
  f = ((unsigned long long)f * (unsigned long long)pc.f) >> 30;

  e += (pc.e - 1) ;
  
  if(PL_BIT_32(f)){
    f >>= 1;
    e ++;
  }

  if(e > PL_MAX_31_BITS){
    *this = plInfinityProbaConstant;
  }
  else{
    if(e < PL_MIN_31_BITS){
      *this = plZeroProbaConstant;
    }
  }
  
  return *this;
}

//=============================================================================
inline  plProbability& plProbability::operator*=(const long double& fl)
{ 
  return (*this)*= ((plProbability)fl);
}


//=============================================================================
inline  plProbability& plProbability::operator*=(const double& fl)
{ 
  return (*this)*= ((plProbability)fl);
}

//=============================================================================
inline  plProbability& plProbability::operator*=(const float& fl)
{ 
  return (*this)*= ((plProbability)fl);
}


//=============================================================================
inline  plProbability plProbability::operator/(const plProbability& pc) const
{
  if(isNaN()) return plNaNProbaConstant;
  if(pc.isNaN()) return plNaNProbaConstant;
  
  if(isInfinity()) {
    if(pc.isInfinity()) return plNaNProbaConstant;
    else return plInfinityProbaConstant;
  }
  
  if(pc.isInfinity()) return plZeroProbaConstant;

  if(pc.f == 0){
    plWarning(33);
    if(f == 0) return plNaNProbaConstant;
    return plInfinityProbaConstant;
  }
  if(f == 0) return plZeroProbaConstant;
  
  plProbability pa;
  pa.e = e - pc.e + 3 ;
  long double f_over_pcf = (long double)f/(long double)pc.f;
  unsigned long int_f_over_pcf = (unsigned long) (PL_MASK_29 * f_over_pcf);
  pa.f = int_f_over_pcf;
  while(pa.f < PL_MASK_31){
    pa.f <<= 1;
    pa.e--;
    if(pa.e < PL_MIN_31_BITS){
      return plZeroProbaConstant;
    }
  }
 
  return pa;
}

//=============================================================================
inline long double operator/(const long double& theFloat, const plProbability& p)
{
  return theFloat/(long double)p;
}

//=============================================================================
inline double operator/(const double& theFloat, const plProbability& p)
{
  return theFloat/double(p);
}

//=============================================================================
inline float operator/(const float& theFloat, const plProbability& p)
{
  return theFloat/float(p);
}

//=============================================================================
inline  plProbability plProbability::operator/(const long double& fl) const 
{ 
  return (*this)/(plProbability)fl;
}

//=============================================================================
inline  plProbability plProbability::operator/(const double& fl) const 
{ 
  return (*this)/(plProbability)fl;
}

//=============================================================================
inline  plProbability plProbability::operator/(const float& fl) const 
{ 
  return (*this)/(plProbability)fl;
}

//=============================================================================
inline  plProbability& plProbability::operator/=(const plProbability& pc)
{
  
  *this =  *this / pc;
  
  return *this;
}

//=============================================================================
inline  plProbability& plProbability::operator/=(const long double& fl)
{
  return (*this) /=((plProbability)fl);
}

//=============================================================================
inline  plProbability& plProbability::operator/=(const double& fl)
{
  return (*this) /=((plProbability)fl);
}

//=============================================================================
inline  plProbability& plProbability::operator/=(const float& fl)
{
  return (*this) /=((plProbability)fl);
}

//=============================================================================
inline std::istream& operator>>(std::istream& in, plProbability& p)
{
    string st;
    in >> st;
    
    if(strcasecmp(st.c_str(), "NAN") == 0) { p = plProbability::plNaNProbaConstant; return in;}
    if(strcasecmp(st.c_str(), "INF") == 0) { p = plProbability::plInfinityProbaConstant; return in;}
    
    long df = 0;
    long ide = 0;
    string dfStr, ideStr;
    
    unsigned int i = 0;
    unsigned int ndec = 0;
    bool afterdot = false;

    while( (st[i] != 'e') && (st[i] != 'E') && (i < st.size()) ) {
        if(afterdot) ndec++;
        
        if(st[i] == '.') afterdot = true;
        else dfStr += st[i];
        
        i++;
    }
    i++;
    while(i < st.size() ) ideStr += st[i++];

    sscanf(dfStr.c_str(), "%ld", &df);
    sscanf(ideStr.c_str(), "%ld", &ide);
    
    if(df == 0){ p = plProbability::plZeroProbaConstant; return in;}
    if(df < 0){ p = plProbability::plNaNProbaConstant; return in;}

    ide -= ndec;
    
    long double fideBin = ide*3.321928094887362660L;
    long ideBin = long(fideBin);
    long double diff = fideBin - ideBin;
    
    long double fdf = df * pow(2.0L, diff);
    
    while(fdf < PL_MASK_31) {
        fdf *= 2.0L;
        ideBin--;
    }
    
    p.e = ideBin + PL_F_SIZE;
    p.f = (unsigned long)fdf;
    
    return in;
}

//=============================================================================
inline std::ostream& operator<<(std::ostream& out,  const plProbability& p)
{
  if(p.isNaN()) {out << "NaN"; return out;}
  if(p.isInfinity()) {out << "Inf"; return out;}
  
  if(p.f ==  0){
    out << "0";
    return out;
  }
  
  long double de = p.e/3.321928094887362660L;
  long ide = long(de);
  long double diff =  de - ide;
  long double df = p.f/(0x80000000 * pow(10.0L, -diff));
  out << df;
  if(ide){
    out << "e";
    if( ide > 0) out << "+";
    out << ide;
  }
  
  //out << "\t---->\tf=" << p.f << " e=" << p.e << hex << "\t---->\tf=" << p.f << " e=" << p.e << endl << dec;
  return out;
}

//=============================================================================
inline double log2(const plProbability &pr)
{
    if(pr.f == 0){
        errno = ERANGE;
        return -HUGE_VAL;
    }
    
    if(pr.isNaN()) return NAN;
    if(pr.isInfinity()) return HUGE_VAL;
    
    return PL_LOG_BASE2((double)pr.f) + pr.e - PL_F_SIZE;
}

//=============================================================================
inline double log(const plProbability &pr)
{
    return log2(pr)/PL_LOG_OF_E_BASE2;
}
  
//=============================================================================
inline double log10(const plProbability &pr)
{
    return log2(pr)/PL_LOG_OF_10_BASE2;
}

//=============================================================================
inline plProbability ppow( const plProbability &pr, long exp )
{
  if( pr.isNaN() )
    {
      return pr;
    }
  else
    {
      if( 0 == exp )
    {
      return plProbability::plOneProbaConstant;
    }
      if( 1 == exp )
    {
      return pr;
    }

      plProbability pr_aux;
      if( exp < 0 )
    {
      pr_aux = plProbability(1)/pr;
      exp = -exp;
    }
      else
    pr_aux = pr;

      double mantisse = pow( (double)pr_aux.f, (double)exp );
      plProbability interm( mantisse );
      if( interm.isNaN() 
      || interm.isInfinity() )
    {
      return plProbability::plNaNProbaConstant;
    }

      long exposant = interm.e + (pr_aux.e - PL_F_SIZE)*exp;
      if( !(LONG_MIN + PL_F_SIZE <= pr_aux.e ) 
      || !( (pr_aux.e - PL_F_SIZE)<= LONG_MAX/exp  )
      || !( LONG_MIN/exp <= (pr_aux.e - PL_F_SIZE) )
      )
    return plProbability::plNaNProbaConstant;

      if( interm.e > 0 )
    {
          if( !( ( pr_aux.e - PL_F_SIZE )*exp <= LONG_MAX - interm.e ) )
        {
          return plProbability::plNaNProbaConstant;
        }
    }
      else
    if( interm.e < 0 )
      {
        if( !(LONG_MIN - interm.e <= ( pr_aux.e - PL_F_SIZE )*exp ) )
          {
        return plProbability::plNaNProbaConstant;
          }
      }

      return plProbability( interm.f, exposant );
    }
}

#endif

---