// -*- mode: c++ -*-
//
// $Id: BrMath.h,v 1.2 2002/02/26 08:58:59 cholm Exp $
// $Author: cholm $
// $Date: 2002/02/26 08:58:59 $
// $Copyright: (C) 2002 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//
///////////////////////////////////////////////////////////////////////////
//                                                                       //
// BrMath                                                                //
//                                                                       //
// Encapsulate math routines. For the time being avoid templates.        //
// Not so happy with ROOT implementation. Edited for more uniform access //
//                                                                       //
///////////////////////////////////////////////////////////////////////////
#ifndef BRAT_BrMath
#define BRAT_BrMath
#include "Rtypes.h"

// This typedef causes trouble for `Red Hat G++-2.96'
// I believe that ISO C++ does not allow for redefinitions of
// typedefs, even if they're the same.  Hence, this should not be
// here.  I've left it in case some compilers have some trouble here. 
// typedef unsigned char Bool_t;

class BrMath {
public:
   static const double Pi, E;

   // Trigo
   static double Sin(double);
   static double Cos(double);
   static double Tan(double);
   static double ASin(double);
   static double ACos(double);
   static double ATan(double);
   static double ATan2(double, double);
   static double Hypot(double x, double y);

   // Misc
   static double Sqrt(double x);
   static double Ceil(double x);
   static double Floor(double x);
   static double Exp(double);
   static double Power(double x, double y);
   static double Log(double x);
   static double Log2(double x);
   static double Log10(double x);
   static int    Nint(float x);
   static int    Nint(double x);
   static int    Int(float x);
   static int    Int(double x);

   // Some integer math
   static long NextPrime(long x);   // least prime number greater than x
   static long Sqrt(long x);
   static long Hypot(long x, long y);     // sqrt(px*px + py*py)

   // Abs
   static short Abs(short d);
   static int   Abs(int d);
   static long  Abs(long d);
   static float Abs(float d);
   static double Abs(double d);

   // Even/Odd
   static Bool_t Even(long a);
   static Bool_t Odd(long a);

   // Signum
   static short Sign(short a);
   static int Sign(int a);
   static long Sign(long a);
   static double Sign(double a);

   // Min
   static short Min(short a, short b);
   static unsigned short Min(unsigned short a, unsigned short b);
   static int Min(int a, int b);
   static unsigned int Min(unsigned int a, unsigned int b);
   static long Min(long a, long b);
   static unsigned long Min(unsigned long a, unsigned long b);
   static float Min(float a, float b);
   static double Min(double a, double b);

   // Max
   static short Max(short a, short b);
   static unsigned short Max(unsigned short a, unsigned short b);
   static int Max(int a, int b);
   static unsigned int Max(unsigned int a, unsigned int b);
   static long Max(long a, long b);
   static unsigned long Max(unsigned long a, unsigned long b);
   static float Max(float a, float b);
   static double Max(double a, double b);

   // Range
   static short Range(short lb, short ub, short x);
   static int Range(int lb, int ub, int x);
   static long Range(long lb, long ub, long x);
   static unsigned long Range(unsigned long lb, unsigned long ub, unsigned long x);
   static double Range(double lb, double ub, double x);

   static Bool_t MrqMin( const Int_t npts, const Float_t *x, const Float_t *y, 
	       const Float_t *dy, const Int_t npars, Double_t *a, 
	       Double_t &alam, Double_t **cov, Double_t &chi );
   static Bool_t MrqCof( const Int_t npts, const Float_t *x, const Float_t *y, 
               const Float_t *dy, const Int_t npars, Double_t *a, 
	       Double_t **alpha, Double_t *beta, Double_t &chi );
   static Bool_t FuncTRF( const Int_t npars, Double_t t, Double_t *par, Double_t &y, 
		Double_t *deriv );
   static Bool_t rsInvert( Int_t N, Double_t **A, Int_t IDIM );
   static Bool_t rsEqn( Int_t N, Double_t **A, Int_t IDIM, Int_t K, Double_t **B );

   static float Erf(float x);
   static double Erf(double x);

   static float FGauss();
   static int MOD(int dividend,int divisor);

   static void vzero(int *arr,int len);

   static void DEQINV(int n,double **a,int idim,float *r,int &ifail,int k,double **b);
   static void DFACT(int n,double **a,int idim,float *r,int &ifail, double &det,int &jfail);
   static void DFEQN(int n,double **a,int idim,float *r,int k,double **b);
   static void DFINV(int n,double **a,int idim,float *r); 

   static Bool_t BTWN(float a,float b,float c);
   static Bool_t BTWN(int a,int b,int c);

};

//---- Misc --------------------------------------------------------------------

inline int BrMath::Nint(float x)
{
   int i = int(x + 0.5);
   if (x + 0.5 == float(i) && i & 1) i--;
   return i;
}

inline int BrMath::Nint(double x)
{
   int i = int(x + 0.5);
   if (x + 0.5 == double(i) && i & 1) i--;
   return i;
}

//__________________________________________
inline int BrMath::Int(float x)
{
   int i;
   if( x > 0 ) {
      i = int(x + 0.5);
      if (x + 0.5 == float(i) && i & 1) i--;
   } 
   else {
      i = int(x - 0.5);
      if (x - 0.5 == float(i) && i & 1) i++;
   }
   return i;
}

inline int BrMath::Int(double x)
{
   int i;
   if( x > 0 ) {
      i = int(x + 0.5);
      if (x + 0.5 == double(i) && i & 1) i--;
   }
   else {
      i = int(x - 0.5);
      if (x - 0.5 == double(i) && i & 1) i++;
   }
   return i;
}

//---- Sign --------------------------------------------------------------------

inline short BrMath::Sign(short a)
   { return (a == 0) ? 0 : ( (a > 0) ? 1 : -1 ); }

inline int BrMath::Sign(int a)
   { return (a == 0) ? 0 : ( (a > 0) ? 1 : -1 ); }

inline long BrMath::Sign(long a)
   { return (a == 0) ? 0 : ( (a > 0) ? 1 : -1 ); }

inline double BrMath::Sign(double a)
   { return (a == 0.0) ? 0 : ( (a > 0.0) ? 1 : -1 ); }

//---- Even/odd ----------------------------------------------------------------

inline Bool_t BrMath::Even(long a)
   { return ! (a & 1); }

inline Bool_t BrMath::Odd(long a)
   { return (a & 1); }

//---- Abs ---------------------------------------------------------------------

inline short BrMath::Abs(short d)
   { return (d > 0) ? d : -d; }

inline int BrMath::Abs(int d)
   { return (d > 0) ? d : -d; }

inline long BrMath::Abs(long d)
   { return (d > 0) ? d : -d; }

inline float BrMath::Abs(float d)
   { return (d > 0) ? d : -d; }

inline double BrMath::Abs(double d)
   { return (d > 0) ? d : -d; }

//---- Min ---------------------------------------------------------------------

inline short BrMath::Min(short a, short b)
   { return a <= b ? a : b; }

inline unsigned short BrMath::Min(unsigned short a, unsigned short b)
   { return a <= b ? a : b; }

inline int BrMath::Min(int a, int b)
   { return a <= b ? a : b; }

inline unsigned int BrMath::Min(unsigned int a, unsigned int b)
   { return a <= b ? a : b; }

inline long BrMath::Min(long a, long b)
   { return a <= b ? a : b; }

inline unsigned long BrMath::Min(unsigned long a, unsigned long b)
   { return a <= b ? a : b; }

inline float BrMath::Min(float a, float b)
   { return a <= b ? a : b; }

inline double BrMath::Min(double a, double b)
   { return a <= b ? a : b; }

//---- Max ---------------------------------------------------------------------

inline short BrMath::Max(short a, short b)
   { return a >= b ? a : b; }

inline unsigned short BrMath::Max(unsigned short a, unsigned short b)
   { return a >= b ? a : b; }

inline int BrMath::Max(int a, int b)
   { return a >= b ? a : b; }

inline unsigned int BrMath::Max(unsigned int a, unsigned int b)
   { return a >= b ? a : b; }

inline long BrMath::Max(long a, long b)
   { return a >= b ? a : b; }

inline unsigned long BrMath::Max(unsigned long a, unsigned long b)
   { return a >= b ? a : b; }

inline float BrMath::Max(float a, float b)
   { return a >= b ? a : b; }

inline double BrMath::Max(double a, double b)
   { return a >= b ? a : b; }

//---- Range -------------------------------------------------------------------

inline short BrMath::Range(short lb, short ub, short x)
   { return x < lb ? lb : (x > ub ? ub : x); }

inline int BrMath::Range(int lb, int ub, int x)
   { return x < lb ? lb : (x > ub ? ub : x); }

inline long BrMath::Range(long lb, long ub, long x)
   { return x < lb ? lb : (x > ub ? ub : x); }

inline unsigned long BrMath::Range(unsigned long lb, unsigned long ub, unsigned long x)
   { return x < lb ? lb : (x > ub ? ub : x); }

inline double BrMath::Range(double lb, double ub, double x)
   { return x < lb ? lb : (x > ub ? ub : x); }

//---- Trig and other functions ------------------------------------------------
// And why not include math.h??? K.Olchanski 19-DEC-2000
#include <math.h>

inline double BrMath::Sin(double x)
   { return sin(x); }

inline double BrMath::Cos(double x)
   { return cos(x); }

inline double BrMath::Tan(double x)
   { return tan(x); }

inline double BrMath::Sqrt(double x)
   { return sqrt(x); }

inline double BrMath::Exp(double x)
   { return exp(x); }

inline double BrMath::Power(double x, double y)
   { return pow(x, y); }

inline double BrMath::Log(double x)
   { return log(x); }

inline double BrMath::Log10(double x)
   { return log10(x); }

#endif
//____________________________________________________________________
//  $Log: BrMath.h,v $
//  Revision 1.2  2002/02/26 08:58:59  cholm
//  Added const'ness to many methods - needed by ISO C++
//
//  Revision 1.1.1.1  2001/06/21 14:55:19  hagel
//  Initial revision of brat2
//
//  Revision 1.8  2000/12/20 01:20:00  brahmlib
//  Removed bogus math library function sincos()
//
//  Revision 1.7  2000/12/20 00:56:05  hagel
//  Replaced bogus sin,cos&co declarations with #include math.h
//
//  Revision 1.6  2000/03/24 21:09:37  videbaek
//  Moved linear equation routines from TpcClusterFinder class.
//
//  Revision 1.5  1999/01/21 17:54:32  hagel
//  Test log change
//
//  Revision 1.4  1999/01/21 17:49:36  hagel
//  Test log change
//
//  Revision 1.3  1999/01/21 17:48:06  hagel
//  Change #ifndef's
//