source: branches/nameserver/GDE/MUSCLE/src/intmath.h

// IntMath.h: Header for doing fractional math with integers for speed.

#ifndef IntMath_h
#define IntMath_h

typedef float BASETYPE;
//typedef double BASETYPE;

// Scaling factor used to store certain floating point
// values as integers to a few significant figures.
//const int INTSCALE = 1000;
const int INTSCALE = 1;

// Type for a probability in range 0.0 to 1.0.
typedef BASETYPE PROB;

// Type for an log-odds integer score.
// Stored as log2(PROB)*INTSCALE.
//typedef int   SCORE;
typedef BASETYPE SCORE;

// Type for a weight.
// Stored as w*INTSCALE where w is in range 0.0 to 1.0.
//typedef unsigned WEIGHT;
typedef BASETYPE WEIGHT;

// Type for a fractional weighted count stored as n*WEIGHT/N
// where n=measured count (integer >= 0) and N is total for
// the distribution (e.g., n=number of residues of a given
// type in a column, N=number of residues in the column).
// Hence values in an FCOUNT variable range from 0..INTSCALE
// as an integer, representing "true" values 0.0 to 1.0.
//typedef unsigned FCOUNT;
typedef BASETYPE FCOUNT;

// Representation of -infinity. Value should
// be large and negative, but not so large
// that adding a few of them overflows.
// TODO: Multiplied by 10 to work around bug
// when aligning Bali 1ckaA in ref4, which is
// so long that B->Mmax got to -infinity, causing
// traceback to fail.
//const int MINUS_INFINITY = -10000000;
const BASETYPE MINUS_INFINITY = (BASETYPE) -1e37;
const BASETYPE PLUS_INFINITY = (BASETYPE) 1e37;

// Probability relative to a null model
typedef double RPROB;

PROB ScoreToProb(SCORE Score);
SCORE ProbToScore(PROB Prob);
SCORE DoubleToScore(double d);
WEIGHT DoubleToWeight(double d);
double WeightToDouble(WEIGHT w);
SCORE MulScoreWeight(SCORE Score, WEIGHT Weight);
bool ScoreEq(SCORE s1, SCORE s2);
bool BTEq(double b1, double b2);

static double ScoreToDouble(SCORE Score)
        {
        return (double) Score / (double) INTSCALE;
        }

#if     0
// In-line assembler for Result = (x*y)/z
// Note that imul and idiv will do 64-bit arithmetic
// on 32-bit operands, so this shouldn't overflow
// Can't write this efficiently in C/C++ (would
// often overlow 32 bits).
#define MulDivAssign(Result, x, y, z)   \
        {                                                                       \
        int X = (x);                                            \
        int Y = (y);                                            \
        int Z = (z);                                            \
        _asm mov        eax,X                                   \
        _asm imul       Y                                               \
        _asm mov        ecx,Z                                   \
        _asm idiv       ecx                                             \
        _asm mov        Result,eax                              \
        }
#else
#define MulDivAssign(Result, x, y, z)   Result = (((x)*(y))/(z))
#endif

#define MulScoreWeight(r, s, w)         MulDivAssign(r, s, w, INTSCALE)
#define MulWeightWCount(r, wt, wc)      MulDivAssign(r, wt, wc, INTSCALE)
#define MulFCountScore(r, fc, sc)       MulDivAssign(r, fc, sc, INTSCALE)

#if     _DEBUG

static inline SCORE Add2(SCORE a, SCORE b)
        {
        if (MINUS_INFINITY == a)
                return MINUS_INFINITY;
        if (MINUS_INFINITY == b)
                return MINUS_INFINITY;
        SCORE sum = a + b;
        if (sum < MINUS_INFINITY)
                return MINUS_INFINITY;
//      assert(sum < OVERFLOW_WARN);
        return sum;
        }

static inline SCORE Add3(SCORE a, SCORE b, SCORE c)
        {
        return Add2(Add2(a, b), c);
        }

static inline SCORE Add4(SCORE a, SCORE b, SCORE c, SCORE d)
        {
        return Add2(Add2(a, b), Add2(c, d));
        }

static inline SCORE Add5(SCORE a, SCORE b, SCORE c, SCORE d, SCORE e)
        {
        return Add3(Add2(a, b), Add2(c, d), e);
        }

static inline SCORE Add6(SCORE a, SCORE b, SCORE c, SCORE d, SCORE e, SCORE f)
        {
        return Add3(Add2(a, b), Add2(c, d), Add2(e, f));
        }

static inline SCORE Add7(SCORE a, SCORE b, SCORE c, SCORE d, SCORE e, SCORE f, SCORE g)
        {
        return Add4(Add2(a, b), Add2(c, d), Add2(e, f), g);
        }

static inline SCORE Mul2(SCORE a, SCORE b)
        {
        if (MINUS_INFINITY == a)
                return MINUS_INFINITY;
        if (MINUS_INFINITY == b)
                return MINUS_INFINITY;
        //__int64 prod = (__int64) a * (__int64) b;
        //assert((SCORE) prod == prod);
        //return (SCORE) prod;
        return a*b;
        }

static inline SCORE Sub2(SCORE a, SCORE b)
        {
        if (MINUS_INFINITY == a)
                return MINUS_INFINITY;
        if (MINUS_INFINITY == b)
                return MINUS_INFINITY;
        SCORE diff = a - b;
        if (diff < MINUS_INFINITY)
                return MINUS_INFINITY;
//      assert(diff < OVERFLOW_WARN);
        return diff;
        }

static inline SCORE Div2(SCORE a, int b)
        {
        if (MINUS_INFINITY == a)
                return MINUS_INFINITY;
        return a/b;
        }

//static inline SCORE MulScoreWeight(SCORE s, WEIGHT w)
//      {
//      SCORE Prod = s*(SCORE) w;
//      assert(Prod < OVERFLOW_WARN);
//      extern void Log(const char Format[], ...);
//      if (Prod/(SCORE) w != s)
//              Log("**WARRNING MulScoreWeight Prod=%d w=%d Prod/w=%d s=%d\n",
//                Prod,
//                w,
//                Prod/(SCORE) w,
//                s);
//      assert(Prod/ (SCORE) w == s);
//      return Prod/INTSCALE;
//      }
//
//static inline WCOUNT MulWeightWCount(WEIGHT wt, WCOUNT wc)
//      {
//      return (wt*wc)/INTSCALE;
//      }

#else
#define Add2(a, b)                                      ((a) + (b))
#define Sub2(a, b)                                      ((MINUS_INFINITY == (a)) ? MINUS_INFINITY : ((a) - (b)))
#define Div2(a, b)                                      ((MINUS_INFINITY == (a)) ? MINUS_INFINITY : ((a) / (b)))
#define Add3(a, b, c)                           ((a) + (b) + (c))
#define Add4(a, b, c, d)                        ((a) + (b) + (c) + (d))
#define Add5(a, b, c, d, e)                     ((a) + (b) + (c) + (d) + (e))
#define Add6(a, b, c, d, e, f)          ((a) + (b) + (c) + (d) + (e) + (f))
#define Add7(a, b, c, d, e, f, g)       ((a) + (b) + (c) + (d) + (e) + (f) + (g))
//#define       MulScoreWeight(s, w)            (((s)*(SCORE) (w))/INTSCALE)
#define Mul2(a, b)                                      ((a)*(b))
#endif

//static inline SCORE MulFCountScore(FCOUNT fc, SCORE sc)
//      {
//// Fast way to say "if (fc >= 2^15 || sc >= 2^15)":
//      if ((fc | sc) & 0xffff1000)
//              {
//              SCORE Score = ((fc+5)/10)*sc;
//              assert(Score < assert);
//              OVERFLOW_WARN(Score > MINUS_INFINITY);
//              return Score/(INTSCALE/10);
//              }
//      SCORE Score = fc*sc;
//      assert(Score < OVERFLOW_WARN);
//      assert(Score > MINUS_INFINITY);
//      return Score/INTSCALE;
//      }

#endif  // IntMath_h