// @(#)root/hist:$Id: TFormulaPrimitive.cxx 29672 2009-08-03 15:37:43Z pcanal $
// Author: Marian Ivanov, 2005

/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
* All rights reserved.                                                  *
*                                                                       *
* For the licensing terms see $ROOTSYS/LICENSE.                         *
* For the list of contributors see $ROOTSYS/README/CREDITS.             *
*************************************************************************/

#include <math.h>

#include "TFormulaPrimitive.h"
#include "TMath.h"
#ifdef WIN32
#pragma optimize("",off)
#endif

void TMath_GenerInterface();

ClassImp(TFormulaPrimitive)

//______________________________________________________________________________
// The Formula Primitive class
//
//    Helper class for TFormula to speed up TFormula evaluation  
//    TFormula can use all functions registered in the list of TFormulaPrimitives 
//    User can add new function to the list of primitives
//    if FormulaPrimitive with given name is already defined new primitive is ignored  
//    Example:
//      TFormulaPrimitive::AddFormula(new TFormulaPrimitive("Pow2","Pow2",TFastFun::Pow2));
//      TF1 f1("f1","Pow2(x)");
//      
//    
//
//    TFormulaPrimitive is used to get direct acces to the function pointers
//    GenFunc     -  pointers  to the static function
//    TFunc       -  pointers  to the data member functions
//
//    The following sufixes are currently used, to describe function arguments:
//    ------------------------------------------------------------------------
//    G     - generic layout - pointer to double (arguments), pointer to double (parameters)
//    10    - double
//    110   - double, double
//    1110  - double, double, double


//______________________________________________________________________________
// TFormula primitive
//
TObjArray * TFormulaPrimitive::fgListOfFunction = 0;


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive() : TNamed()
{
   // Default constructor.

   fType=0;
   fNArguments=0;
   fNParameters=0;
   fIsStatic=kTRUE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula,
                                     GenFunc0 fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType= 0;
   fFunc0 = fpointer; 
   fNArguments=0;
   fNParameters=0;
   fIsStatic=kTRUE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula, 
                                     GenFunc10 fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType= 10;
   fFunc10 = fpointer; 
   fNArguments=1;
   fNParameters=0;
   fIsStatic=kTRUE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula,
                                     GenFunc110 fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType= 110;
   fFunc110 = fpointer; 
   fNArguments=2;
   fNParameters=0;
   fIsStatic=kTRUE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula, 
                                     GenFunc1110 fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType= 1110;
   fFunc1110 = fpointer; 
   fNArguments=3;
   fNParameters=0;
   fIsStatic=kTRUE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula, 
                                     GenFuncG fpointer,Int_t npar) : TNamed(name,formula)
{
   // Constructor.

   fType= -1;
   fFuncG = fpointer; 
   fNArguments=2; 
   fNParameters=npar;
   fIsStatic=kTRUE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula,
                                     TFuncG fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType= 0;
   fTFuncG=fpointer;  
   fNArguments=0;
   fNParameters=0;
   fIsStatic =kFALSE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula,
                                     TFunc0 fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType  = 0;
   fTFunc0= fpointer;
   fNArguments=0;
   fNParameters=0;
   fIsStatic =kFALSE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula,
                                     TFunc10 fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType=-10;
   fTFunc10=fpointer;
   fNArguments=1;
   fNParameters=0;
   fIsStatic =kFALSE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula,
                                     TFunc110 fpointer) : TNamed(name,formula)
{
   // Constructor.

   fType=-110;
   fTFunc110=fpointer;
   fNArguments=2;
   fNParameters=0;
   fIsStatic =kFALSE;
}


//______________________________________________________________________________
TFormulaPrimitive::TFormulaPrimitive(const char *name,const char *formula,
                                     TFunc1110 fpointer) :TNamed(name,formula)
{
   // Constructor.

   fType=-1110;
   fTFunc1110=fpointer;
   fNArguments=3;
   fNParameters=0;
   fIsStatic =kFALSE;
}


//______________________________________________________________________________
Int_t TFormulaPrimitive::AddFormula(TFormulaPrimitive * formula)
{
   // Add formula to the list of primitive formulas.
   // If primitive formula already defined do nothing.

   if (fgListOfFunction == 0) BuildBasicFormulas();
   if (FindFormula(formula->GetName(),formula->fNArguments)){
      delete formula;
      return 0;
   }
   fgListOfFunction->AddLast(formula);
   return 1;
}


//______________________________________________________________________________
Double_t TFormulaPrimitive::Eval(Double_t* x)
{
   // Eval primitive function at point x.

   if (fIsStatic ==kFALSE) return 0;

   if (fType==0) return  fFunc0();
   if (fType==10) {
      return fFunc10(x[0]);
   }
   if (fType==110) {
      return fFunc110(x[0],x[1]);
   }
   if (fType==1110) {
      return fFunc1110(x[0],x[1],x[2]);
   }
   return 0;
} 


//______________________________________________________________________________
Double_t  TFormulaPrimitive::Eval(TObject *o, Double_t *x)
{
   // Eval member function of object o at point x.

   if (fIsStatic ==kTRUE) return 0;
   if (fType== 0)    return (*o.*fTFunc0)();
   if (fType==-10)   return (*o.*fTFunc10)(*x);
   if (fType==-110)  return (*o.*fTFunc110)(x[0],x[1]);
   if (fType==-1110) return (*o.*fTFunc1110)(x[0],x[1],x[2]);
   return 0;
}


//______________________________________________________________________________
Double_t TFormulaPrimitive::Eval(Double_t *x, Double_t *param)
{
   // Eval primitive parametric function.

   return fFuncG(x,param);
}


#define RTFastFun__POLY(var)                                          \
{                                                                     \
   Double_t res= param[var-1]+param[var]*x[0];                        \
   for (Int_t j=var-1 ;j>0;j--) res      = param[j-1]+x[0]*res;       \
   return res;                                                        \
}

namespace TFastFun {
   //
   // Namespace with basic primitive functions registered by TFormulaPrimitive
   // all function registerd by TFormulaPrimitive can be used in TFormula
   //
   Double_t Pow2(Double_t x){return x*x;}
   Double_t Pow3(Double_t x){return x*x*x;}
   Double_t Pow4(Double_t x){return x*x*x*x;}
   Double_t Pow5(Double_t x){return x*x*x*x*x;}
   inline   Double_t FPoln(Double_t *x, Double_t *param, Int_t npar);
   Double_t FPol0(Double_t * /*x*/, Double_t *param){ return param[0];}
   Double_t FPol1(Double_t *x, Double_t *param){ return param[0]+param[1]*x[0];}
   Double_t FPol2(Double_t *x, Double_t *param){ return param[0]+x[0]*(param[1]+param[2]*x[0]);} 
   Double_t FPol3(Double_t *x, Double_t *param){ return param[0]+x[0]*(param[1]+x[0]*(param[2]+param[3]*x[0]));}
   Double_t FPol4(Double_t *x, Double_t *param){ RTFastFun__POLY(4)}
   Double_t FPol5(Double_t *x, Double_t *param){ RTFastFun__POLY(5)}
   Double_t FPol6(Double_t *x, Double_t *param){ RTFastFun__POLY(6)}
   Double_t FPol7(Double_t *x, Double_t *param){ RTFastFun__POLY(7)}
   Double_t FPol8(Double_t *x, Double_t *param){ RTFastFun__POLY(8)}
   Double_t FPol9(Double_t *x, Double_t *param){ RTFastFun__POLY(9)}
   Double_t FPol10(Double_t *x, Double_t *param){ RTFastFun__POLY(10)}
   //
   //
   Double_t PlusXY(Double_t x,Double_t y){return x+y;}
   Double_t MinusXY(Double_t x,Double_t y){return x-y;}
   Double_t MultXY(Double_t x,Double_t y){return x*y;}
   Double_t DivXY(Double_t x, Double_t y){return TMath::Abs(y)>0 ? x/y:0;} 
   Double_t XpYpZ(Double_t x, Double_t y, Double_t z){ return x+y+z;}
   Double_t XxYxZ(Double_t x, Double_t y, Double_t z){ return x*y*z;}
   Double_t XxYpZ(Double_t x, Double_t y, Double_t z){ return x*(y+z);}
   Double_t XpYxZ(Double_t x, Double_t y, Double_t z){ return x+(y*z);}
   Double_t Gaus(Double_t x, Double_t mean, Double_t sigma);
   Double_t Gausn(Double_t x, Double_t mean, Double_t sigma);
   Double_t Landau(Double_t x, Double_t mean, Double_t sigma){return TMath::Landau(x,mean,sigma,kFALSE);}
   Double_t Landaun(Double_t x, Double_t mean, Double_t sigma){return TMath::Landau(x,mean,sigma,kTRUE);}
   Double_t Sqrt(Double_t x) {return x>0?sqrt(x):0;}
   //  
   Double_t Sign(Double_t x){return (x<0)? -1:1;}
   Double_t Nint(Double_t x){return TMath::Nint(x);}
   Double_t Abs(Double_t x){return TMath::Abs(x);}
   //logical
   Double_t XandY(Double_t x, Double_t y){ return (x*y>0.1);}
   Double_t XorY(Double_t x, Double_t y) { return (x+y>0.1);}
   Double_t XgY(Double_t x, Double_t y) {return (x>y);}
   Double_t XgeY(Double_t x, Double_t y) {return (x>=y);}
   Double_t XlY(Double_t x, Double_t y) {return (x<y);}
   Double_t XleY(Double_t x, Double_t y) {return (x<=y);}
   Double_t XeY(Double_t x,Double_t y) {return (x==y);}
   Double_t XneY(Double_t x,Double_t y) {return (x!=y);}
   Double_t XNot(Double_t x){ return (x<0.1);}
};


//______________________________________________________________________________
TFormulaPrimitive* TFormulaPrimitive::FindFormula(const char* name)
{
   // Find the formula in the list of formulas.
   if (!fgListOfFunction) {
      BuildBasicFormulas();
   }      
   Int_t nobjects = fgListOfFunction->GetEntries();
   for (Int_t i = 0; i < nobjects; ++i) {
      TFormulaPrimitive *formula = (TFormulaPrimitive*)fgListOfFunction->At(i);
      if (formula && 0==strcmp(name, formula->GetName())) return formula;
   }
   return 0;
}


//______________________________________________________________________________
TFormulaPrimitive* TFormulaPrimitive::FindFormula(const char* name, UInt_t nargs)
{
   // Find the formula in the list of formulas.
   
   if (!fgListOfFunction) {
      BuildBasicFormulas();
   }
   Int_t nobjects = fgListOfFunction->GetEntries();
   for (Int_t i = 0; i < nobjects; ++i) {
      TFormulaPrimitive *prim = (TFormulaPrimitive*)fgListOfFunction->At(i);
      if (prim) {
         bool match = ( ((UInt_t)prim->fNArguments) == nargs );
         if (match && 0==strcmp(name, prim->GetName())) return prim;
      }
   }
   return 0;
}


//______________________________________________________________________________
TFormulaPrimitive* TFormulaPrimitive::FindFormula(const char* name, const char *args)
{
   // Find the formula in the list of formulas.
   
   // let's count the argument(s)
   if (args) {
      Int_t nargs = 0;
      if (args[0]!=')') {
         nargs = 1;
         int nest = 0;
         for(UInt_t c = 0; c < strlen(args); ++c ) {
            switch (args[c]) {
               case '(': ++nest; break;
               case ')': --nest; break;
               case '<': ++nest; break;
               case '>': --nest; break;
               case ',': nargs += (nest==0); break;
            }
         }
      }
      return FindFormula(name,nargs);
   } else {
      return FindFormula(name);
   }
   return 0;
}


//______________________________________________________________________________
Double_t TFastFun::FPoln(Double_t *x, Double_t *param, Int_t npar)
{
   // FPoln.

   Double_t res = 0; Double_t temp=1;                                  
   for (Int_t j=npar ;j>=0;j--) {                                            
      res  += temp*param[j];                                           
      temp *= *x;                                                  
   }                                                                   
   return res;
}


//______________________________________________________________________________
Double_t TFastFun::Gaus(Double_t x, Double_t mean, Double_t sigma)
{
   // Gauss.

   if (sigma == 0) return 1.e30;
   Double_t arg = (x-mean)/sigma;
   return TMath::Exp(-0.5*arg*arg);
}


//______________________________________________________________________________
Double_t TFastFun::Gausn(Double_t x, Double_t mean, Double_t sigma)
{ 
   // Normalize gauss.

   if (sigma == 0)  return 0;
   Double_t arg = (x-mean)/sigma;
   return TMath::Exp(-0.5*arg*arg)/(2.50662827463100024*sigma);  //sqrt(2*Pi)=2.50662827463100024
}


//______________________________________________________________________________
Int_t TFormulaPrimitive::BuildBasicFormulas()
{
   // Built-in functions.

   if (fgListOfFunction==0) fgListOfFunction = new TObjArray(1000);
   //
   // logical
   //
   AddFormula(new TFormulaPrimitive("XandY","XandY",TFastFun::XandY));
   AddFormula(new TFormulaPrimitive("XorY","XorY",TFastFun::XorY));
   AddFormula(new TFormulaPrimitive("XNot","XNot",TFastFun::XNot));
   AddFormula(new TFormulaPrimitive("XlY","XlY",TFastFun::XlY));
   AddFormula(new TFormulaPrimitive("XleY","XleY",TFastFun::XleY));
   AddFormula(new TFormulaPrimitive("XgY","XgY",TFastFun::XgY));
   AddFormula(new TFormulaPrimitive("XgeY","XgeY",TFastFun::XgeY));
   AddFormula(new TFormulaPrimitive("XeY","XeY",TFastFun::XeY));
   AddFormula(new TFormulaPrimitive("XneY","XneY",TFastFun::XneY));
   // addition  + multiplication
   AddFormula(new TFormulaPrimitive("PlusXY","PlusXY",TFastFun::PlusXY));
   AddFormula(new TFormulaPrimitive("MinusXY","MinusXY",TFastFun::MinusXY));
   AddFormula(new TFormulaPrimitive("MultXY","MultXY",TFastFun::MultXY));
   AddFormula(new TFormulaPrimitive("DivXY","DivXY",TFastFun::DivXY));
   AddFormula(new TFormulaPrimitive("XpYpZ","XpYpZ",TFastFun::XpYpZ));
   AddFormula(new TFormulaPrimitive("XxYxZ","XxYxZ",TFastFun::XxYxZ));
   AddFormula(new TFormulaPrimitive("XxYpZ","XxYpZ",TFastFun::XxYpZ));
   AddFormula(new TFormulaPrimitive("XpYxZ","XpYxZ",TFastFun::XpYxZ));
   //
   //
   AddFormula(new TFormulaPrimitive("Gaus","Gaus",TFastFun::Gaus));
   AddFormula(new TFormulaPrimitive("Gausn","Gausn",TFastFun::Gausn));
   AddFormula(new TFormulaPrimitive("Landau","Landau",TFastFun::Landau));
   AddFormula(new TFormulaPrimitive("Landaun","Landaun",TFastFun::Landaun));
   //
   //
   // polynoms
   //
   //
   AddFormula(new TFormulaPrimitive("Pol0","Pol0",(GenFuncG)TFastFun::FPol0,1)); 
   AddFormula(new TFormulaPrimitive("Pol1","Pol1",(GenFuncG)TFastFun::FPol1,2)); 
   AddFormula(new TFormulaPrimitive("Pol2","Pol2",(GenFuncG)TFastFun::FPol2,3));                   
   AddFormula(new TFormulaPrimitive("Pol3","Pol3",(GenFuncG)TFastFun::FPol3,4));                       
   AddFormula(new TFormulaPrimitive("Pol4","Pol4",(GenFuncG)TFastFun::FPol4,5));   
   AddFormula(new TFormulaPrimitive("Pol5","Pol5",(GenFuncG)TFastFun::FPol5,6));
   AddFormula(new TFormulaPrimitive("Pol6","Pol6",(GenFuncG)TFastFun::FPol6,7));
   AddFormula(new TFormulaPrimitive("Pol7","Pol7",(GenFuncG)TFastFun::FPol7,8));
   AddFormula(new TFormulaPrimitive("Pol8","Pol8",(GenFuncG)TFastFun::FPol8,9));
   AddFormula(new TFormulaPrimitive("Pol9","Pol9",(GenFuncG)TFastFun::FPol9,10));
   AddFormula(new TFormulaPrimitive("Pol10","Pol10",(GenFuncG)TFastFun::FPol10,11));
   //
   // pows
   AddFormula(new TFormulaPrimitive("Pow2","Pow2",TFastFun::Pow2));                        
   AddFormula(new TFormulaPrimitive("Pow3","Pow3",TFastFun::Pow3));                   
   AddFormula(new TFormulaPrimitive("Pow4","Pow4",TFastFun::Pow4));                        
   AddFormula(new TFormulaPrimitive("Pow5","Pow5",TFastFun::Pow5));                         
   //
   //
   AddFormula(new TFormulaPrimitive("TMath::Cos","TMath::Cos",cos));              // 10
   AddFormula(new TFormulaPrimitive("cos","cos",cos));                            // 10
   AddFormula(new TFormulaPrimitive("TMath::Sin","TMath::Sin",sin));              // 11
   AddFormula(new TFormulaPrimitive("sin","sin",sin));                            // 11
   AddFormula(new TFormulaPrimitive("TMath::Tan","TMath::Tan",tan));              // 12
   AddFormula(new TFormulaPrimitive("tan","tan",tan));                            // 12
   AddFormula(new TFormulaPrimitive("TMath::ACos","TMath::ACos",acos));           // 13
   AddFormula(new TFormulaPrimitive("acos","acos",acos));                         // 13
   AddFormula(new TFormulaPrimitive("TMath::ASin","TMath::ASin",asin));           // 14
   AddFormula(new TFormulaPrimitive("asin","asin",asin));                         // 14
   AddFormula(new TFormulaPrimitive("TMath::ATan","TMath::ATan",atan));           // 15
   AddFormula(new TFormulaPrimitive("atan","atan",atan));                         // 15
   AddFormula(new TFormulaPrimitive("TMath::ATan2","TMath::ATan2",atan2));        // 16
   AddFormula(new TFormulaPrimitive("atan2","atan2",atan2));                      // 16
   //   kpow      = 20, ksq = 21, ksqrt     = 22, 
   AddFormula(new TFormulaPrimitive("pow","pow",TMath::Power));                 //20
   AddFormula(new TFormulaPrimitive("sq","sq",TFastFun::Pow2));                 //21
   AddFormula(new TFormulaPrimitive("sqrt","sqrt",TFastFun::Sqrt));             //22
   // kmin      = 24, kmax = 25,
   AddFormula(new TFormulaPrimitive("min","min",(GenFunc110)TMath::Min));       //24
   AddFormula(new TFormulaPrimitive("max","max",(GenFunc110)TMath::Max));       //25
   // klog      = 30, kexp = 31, klog10 = 32,
   AddFormula(new TFormulaPrimitive("log","log",TMath::Log));                           //30
   AddFormula(new TFormulaPrimitive("exp","exp",TMath::Exp));                           //31
   AddFormula(new TFormulaPrimitive("log10","log10",TMath::Log10));                     //32
   //  
   //    cosh        70                  acosh        73
   //    sinh        71                  asinh        74
   //    tanh        72                  atanh        75
   //
   AddFormula(new TFormulaPrimitive("TMath::CosH","TMath::Cosh",cosh));                     // 70
   AddFormula(new TFormulaPrimitive("cosh","cosh",cosh));                                   // 70
   AddFormula(new TFormulaPrimitive("TMath::SinH","TMath::SinH",sinh));                     // 71
   AddFormula(new TFormulaPrimitive("sinh","sinh",sinh));                                   // 71
   AddFormula(new TFormulaPrimitive("TMath::TanH","TMath::Tanh",tanh));                     // 72
   AddFormula(new TFormulaPrimitive("tanh","tanh",tanh));                                   // 72
   AddFormula(new TFormulaPrimitive("TMath::ACosH","TMath::ACosh",TMath::ACosH));           // 73
   AddFormula(new TFormulaPrimitive("acosh","acosH",TMath::ACosH));                         // 73
   AddFormula(new TFormulaPrimitive("TMath::ASinH","TMath::ASinh",TMath::ASinH));           // 74
   AddFormula(new TFormulaPrimitive("acosh","acosH",TMath::ASinH));                         // 74
   AddFormula(new TFormulaPrimitive("TMath::ATanH","TMath::ATanh",TMath::ATanH));           // 75
   AddFormula(new TFormulaPrimitive("atanh","atanh",TMath::ATanH));                         // 75
   //  
   AddFormula(new TFormulaPrimitive("TMath::Abs","TMath::Abs",TMath::Abs));
   AddFormula(new TFormulaPrimitive("TMath::BreitWigner","TMath::BreitWigner",TMath::BreitWigner));

   //Disable direct access to TMath::Landau for now because of the default parameter.
   //AddFormula(new TFormulaPrimitive("TMath::Landau","TMath::Landau",(TFormulaPrimitive::GenFunc1110)TMath::Landau));

   TMath_GenerInterface();
   return 1;
}