ComPackExpFloat.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: ComPackExpFloat.cxx,v 1.2 2009/12/23 02:59:56 zhangy Exp $
//
// Description:
//      class ComPackExpFloat
//
// Environment:
//      Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//      Dave Brown               10/12/00
//
// Copyright Information:
//      Copyright (C) 2000 LBNL
//
// History:
//  Migration for BESIII MDC
//
//
//--------------------------------------------------------------------------
 
// #include "BaBar/BaBar.hh"
#include "MdcRecoUtil/ComPackExpFloat.h"
// #include "ErrLogger/ErrLog.hh"
#include <iostream>
#include <math.h>
using std::endl;
using std::ostream;
 
ComPackExpFloat::ComPackExpFloat( unsigned nbits, unsigned maxexponent, double start,
                                  double stop, bool center )
    : ComPackBase<double>( start, stop, nbits ), _maxexp( maxexponent ) {
  if ( _maxexp > 62 )
  {
    std::cout << "ErrMsg(fatal) "
              << "Can't pack exponents larger than 62" << endl;
    abort();
  }
  static double invln2( 1.0 / log( 2.0 ) );
  long long     one( 1 );
  // compute how many bits used by exponent
  _mansft = _maxexp > 0 ? (unsigned)( log( double( _maxexp ) ) * invln2 + 1.01 ) : 0;
  _expmsk = ( 1 << _mansft ) - 1;
  _maxman = ( 1 << ( _bitRange - _mansft ) );
  _manmsk = _maxman - 1;
  double    invmaxman = 1.0 / _maxman;
  long long maxnorm   = ( one << ( maxexponent + 1 ) ) - 1;
  // correct the range if we're centering the bins
  if ( center )
  {
    long long norm  = one << maxexponent;
    double    alpha = 1.0 / ( (double)maxnorm * _maxman * 2.0 - ( norm + 1.0 ) );
    _minVal -= _valRange * alpha;
    _maxVal += _valRange * norm * alpha;
    _valRange = _maxVal - _minVal;
    //    cout << "alpha = " << alpha << " norm = " << norm << endl;
  }
  // finish the calculations
  _invrange         = maxnorm / _valRange;
  double invmaxnorm = _valRange / maxnorm;
  double manoff     = ( 0.5 * invmaxman + 1.0 ) * invmaxnorm;
  double valoff     = _minVal - invmaxnorm;
  double manfac     = invmaxnorm * invmaxman;
  // pre-compute some factors and store them in an array
  _expfac = new double[_maxexp + 1];
  _expoff = new double[_maxexp + 1];
  for ( unsigned iexp = 0; iexp <= _maxexp; iexp++ )
  {
    double expf   = ( one << iexp );
    _expoff[iexp] = valoff + manoff * expf;
    _expfac[iexp] = manfac * expf;
  }
}
 
ComPackExpFloat::~ComPackExpFloat() {
  delete[] _expfac;
  delete[] _expoff;
}
 
ComPackBase<double>::StatusCode ComPackExpFloat::pack( const double value,
                                                       d_ULong&     packdata ) const {
  static double                   invln2( 1.0 / log( 2.0 ) );
  static long long                one( 1 );
  ComPackBase<double>::StatusCode retval( TAG_OK );
  // first, figure out the exponent
  double renorm = 1.0 + ( value - _minVal ) * _invrange;
  if ( renorm < 1.0 )
  {
    renorm = 1.0;
    retval = TAG_VAL_ROUND_UP;
  }
  unsigned iexp = unsigned( log( renorm ) * invln2 );
  unsigned iman( 0 );
  if ( iexp <= _maxexp )
  {
    long long expon = one << iexp;
    iman            = (unsigned)( _maxman * ( renorm / expon - 1.0 ) );
    // deal with floating-point rounding, which could give a 'non-physical' result
    if ( iman == _maxman )
    {
      if ( iexp == _maxexp ) iman = _maxman - 1;
      else
      {
        iman = 0;
        iexp++;
      }
    }
  }
  else
  {
    iexp   = _maxexp;
    iman   = _maxman - 1;
    retval = TAG_VAL_ROUND_DOWN;
  }
  packdata = ( iexp & _expmsk ) | ( iman & _manmsk ) << _mansft;
  return retval;
}
 
ComPackBase<double>::StatusCode ComPackExpFloat::unpack( const d_ULong packdata,
                                                         double&       value ) const {
  size_t iman = ( packdata >> _mansft ) & _manmsk;
  size_t iexp = packdata & _expmsk;
  value       = _expoff[iexp] + iman * _expfac[iexp];
  return TAG_OK;
}
 
void ComPackExpFloat::print( ostream& os ) const {
  os << "Exponential packer for range " << _minVal << " to " << _maxVal << endl;
  os << "Maximum exponent = " << _maxexp << endl;
  os << "Maximum mantissa = " << _maxman << endl;
  os << "Mantissa storage shift, mask = " << _mansft << " , " << _manmsk << endl;
}