00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "MdcRecoUtil/ComPackExpFloat.h"
00025
00026 #include <math.h>
00027 #include <iostream>
00028 using std::endl;
00029 using std::ostream;
00030
00031 ComPackExpFloat::ComPackExpFloat(unsigned nbits,
00032 unsigned maxexponent,
00033 double start,
00034 double stop,
00035 bool center ) :
00036 ComPackBase<double>(start,stop,nbits),
00037 _maxexp(maxexponent) {
00038 if(_maxexp > 62 ) {
00039 std::cout<<"ErrMsg(fatal) "<< "Can't pack exponents larger than 62" << endl;
00040 abort();
00041 }
00042 static double invln2(1.0/log(2.0));
00043 long long one(1);
00044
00045 _mansft = _maxexp>0?(unsigned)(log(double(_maxexp))*invln2+1.01):0;
00046 _expmsk = (1<<_mansft)-1;
00047 _maxman = (1<<(_bitRange-_mansft));
00048 _manmsk = _maxman-1;
00049 double invmaxman = 1.0/_maxman;
00050 long long maxnorm = (one<<(maxexponent+1))-1;
00051
00052 if(center) {
00053 long long norm = one<<maxexponent;
00054 double alpha = 1.0/((double)maxnorm*_maxman*2.0 - (norm + 1.0));
00055 _minVal -= _valRange*alpha;
00056 _maxVal += _valRange*norm*alpha;
00057 _valRange = _maxVal - _minVal;
00058
00059 }
00060
00061 _invrange = maxnorm/_valRange;
00062 double invmaxnorm = _valRange/maxnorm;
00063 double manoff = (0.5 * invmaxman+1.0) * invmaxnorm;
00064 double valoff = _minVal - invmaxnorm;
00065 double manfac = invmaxnorm * invmaxman;
00066
00067 _expfac = new double[_maxexp+1];
00068 _expoff = new double[_maxexp+1];
00069 for(unsigned iexp=0;iexp<=_maxexp;iexp++){
00070 double expf = (one<<iexp);
00071 _expoff[iexp] = valoff + manoff*expf;
00072 _expfac[iexp] = manfac*expf;
00073 }
00074 }
00075
00076 ComPackExpFloat::~ComPackExpFloat()
00077 {
00078 delete [] _expfac;
00079 delete [] _expoff;
00080 }
00081
00082 ComPackBase<double>::StatusCode
00083 ComPackExpFloat::pack (const double value, d_ULong & packdata) const {
00084 static double invln2(1.0/log(2.0));
00085 static long long one(1);
00086 ComPackBase<double>::StatusCode retval(TAG_OK);
00087
00088 double renorm = 1.0+(value-_minVal)*_invrange;
00089 if(renorm<1.0){
00090 renorm=1.0;
00091 retval = TAG_VAL_ROUND_UP;
00092 }
00093 unsigned iexp = unsigned(log(renorm)*invln2);
00094 unsigned iman(0);
00095 if(iexp<=_maxexp){
00096 long long expon = one<<iexp;
00097 iman = (unsigned)( _maxman*(renorm/expon - 1.0) );
00098
00099 if(iman==_maxman){
00100 if(iexp==_maxexp)
00101 iman=_maxman-1;
00102 else {
00103 iman=0;
00104 iexp++;
00105 }
00106 }
00107 } else {
00108 iexp=_maxexp;
00109 iman = _maxman-1;
00110 retval = TAG_VAL_ROUND_DOWN;
00111 }
00112 packdata = (iexp&_expmsk) | (iman&_manmsk)<<_mansft;
00113 return retval;
00114 }
00115
00116 ComPackBase<double>::StatusCode
00117 ComPackExpFloat::unpack (const d_ULong packdata, double & value) const {
00118 size_t iman = (packdata>>_mansft)&_manmsk;
00119 size_t iexp = packdata&_expmsk;
00120 value = _expoff[iexp] + iman * _expfac[iexp];
00121 return TAG_OK;
00122 }
00123
00124 void
00125 ComPackExpFloat::print(ostream& os) const {
00126 os << "Exponential packer for range " << _minVal << " to " << _maxVal << endl;
00127 os << "Maximum exponent = " << _maxexp << endl;
00128 os << "Maximum mantissa = " << _maxman << endl;
00129 os << "Mantissa storage shift, mask = " << _mansft << " , " << _manmsk << endl;
00130 }
00131