#include "rpneval.h"
#include <stdlib.h>
#include "strutilxx.h"
#include "srandgen.h"
#include <iostream>
#include <math.h>

namespace SOPHYA {

/*!
  \class SOPHYA::RPNExpressionEvaluator
  \ingroup SysTools
  Arithmetic expression (double precision float) evaluator 
  in Reverse Polish Notation (RPN). This is an HP calculator 
  like syntax. Space are used for separating the string 
  expression into tokens. 
*/

RPNExpressionEvaluator::RPNExpressionEvaluator(string const & sex)
{
  vector<string> exe;
  FillVStringFrString(sex, exe, ' ');
  int rc = EvalRPNExpr(exe, 0);
  if (rc != 0) throw RPNExprException("RPNExpressionEvaluator() - syntax error");  
}

RPNExpressionEvaluator::RPNExpressionEvaluator(vector<string> & exe, int off)
{
  int rc = EvalRPNExpr(exe, off);
  if (rc != 0) throw RPNExprException("RPNExpressionEvaluator() - syntax error");
}

RPNExpressionEvaluator::~RPNExpressionEvaluator()
{
}

/* Operations sur le stack RPN */
/* --Methode-- */
double RPNExpressionEvaluator::Evaluate() const
{
  double x;
  if ( CheckStack( x) ) 
    throw RPNExprException("RPNExpressionEvaluator::Evaluate() EmptyStack");
  else return x;
}
  
/* --Methode-- */
int RPNExpressionEvaluator::EvalRPNExpr(vector<string> & args, int off)
{
  
  if (args.size() <= off)  return 0;
  double x,y;
  x = y = 0.;
  stack<double> rpnstack_;  // Stack des operations en RPN 
  for(int k=off; k<args.size(); k++) {
    // Les 4 operations de base + - * / 
    if (args[k] == "+") { 
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = y+x;
    }
    else if (args[k] == "-") { 
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = y-x;
    }
    else if (args[k] == "*") { 
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = y*x;
    }
    else if (args[k] == "/") { 
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = y/x;
    }
    else if (args[k] == "%") { 
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = (int)y % (int)x;
    }
    // Les constantes : e , pi
    else if (args[k] == "e") {
      rpnstack_.push(M_E);
    }
    else if (args[k] == "pi") {
      rpnstack_.push(M_PI);
    }
    // Les fonctions usuelles a 1 argument f(x)
    else if (args[k] == "cos") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = cos(x);
    }
    else if (args[k] == "sin") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = sin(x);
    }
    else if (args[k] == "tan") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = tan(x);
    }
    else if (args[k] == "acos") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = acos(x);
    }
    else if (args[k] == "asin") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = asin(x);
    }
    else if (args[k] == "atan") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = atan(x);
    }
    else if (args[k] == "chs") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = -x;
    }
    else if (args[k] == "sqrt") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = sqrt(x);
    }
    else if (args[k] == "sq") {  // x^2 
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = x*x;
    }
    else if (args[k] == "log") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = log(x);
    }
    else if (args[k] == "log10") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = log10(x);
    }
    else if (args[k] == "exp") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = exp(x);
    }
    else if (args[k] == "fabs") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = fabs(x);
    }
    else if (args[k] == "floor") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = floor(x);
    }
    else if (args[k] == "ceil") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = ceil(x);
    }
    // trunc et nint vire - ca ne compile pas sous linux - Reza 01/2003
    else if (args[k] == "deg2rad") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = x*M_PI/180.;
    }
    else if (args[k] == "rad2deg") {
      if ( CheckStack( x) ) return(1);
      rpnstack_.top() = x*180./M_PI;
    }
    // Les fonctions usuelles a 2 argument f(x,y)
    else if (args[k] == "pow") { 
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = pow(y,x);
    }
    else if (args[k] == "atan2") { 
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = atan2(x,y);
    }
    // generateur aleatoire
    else if (args[k] == "rand") { 
      double rnd = drand01();
      rpnstack_.push(rnd);
    }
    else if (args[k] == "norand") { 
      double rnd = GauRnd(0., 1.);
      rpnstack_.push(rnd);
    }
    // Fonction a N arguments  - Somme, produit, etc ...
    else if ((args[k] == "sum") || (args[k] == "mean") || (args[k] == "sigmean") || 
	     (args[k] == "sigma") || (args[k] == "sigma2") ) {
      double sx, sx2;
      int nn = SumStack( sx, sx2);
      if (args[k] == "sum") rpnstack_.push(sx);
      else { 
	if (nn == 0) return 1;
	double fnn = nn;
	if ((args[k] == "sigma") || (args[k] == "sigmean")) 
	  rpnstack_.push(sqrt(sx2/fnn-(x*x/(fnn*fnn))));
	else if ((args[k] == "mean") || (args[k] == "sigmean"))  rpnstack_.push(sx/fnn);
	else rpnstack_.push(sx2/fnn-(x*x/(fnn*fnn)));
      }
    }
    else if (args[k] == "product") { 
      double px;
      int nn = ProductStack( px);
      if (nn == 0) return(1);
      rpnstack_.push(px);
    }
    // Fonctions de manipulation de stack
    else if (args[k] == "print") { 
      PrintStack();
    }
    else if (args[k] == "x<>y") {
      if ( CheckStack( x, y) ) return(1);
      rpnstack_.top() = x;  rpnstack_.push(y);
    }
    else if (args[k] == "pop") {
      rpnstack_.pop();
    }
    else if (args[k] == "push") {
      if (rpnstack_.empty()) rpnstack_.push(0.);
      else rpnstack_.push(rpnstack_.top());
    }
    // On met un nombre sur le stack 
    else {
      char * esptr;
      x = strtod(args[k].c_str(), &esptr);
      //      if (ctof(args[k].c_str(),&x) < 0) {
      if (esptr == args[k].c_str()) return 2; 
      rpnstack_.push(x);
    }

  }
  return(0);
}

inline void RPNExpressionEvaluator::PrintStack()
{
  if (rpnstack_.empty()) 
    cout << "RPNExpressionEvaluator::PrintStack() Empty stack " << endl;
  else {
    stack<double> s;
    s = rpnstack_;
    int k = 0;
    cout << "RPNExpressionEvaluator::PrintStack() Size()= " << s.size() << endl;
    while( !s.empty() ) {
      cout << "    " << k << ":  " << s.top() << "  ";
      if (k == 0)  cout << " (x) " << endl;
      else if (k == 1) cout << " (y) " << endl;
      else if (k == 2) cout << " (z) " << endl;
      else cout << endl;
      s.pop(); k++;
    }
  }
  
} 

int RPNExpressionEvaluator::SumStack(double& sx, double& sx2) 
{
  sx = sx2 = 0.;
  int nn = 0;
  double x = 0.;
  while( !rpnstack_.empty() ) {
    x = rpnstack_.top(); rpnstack_.pop();
    sx += x; sx2 += x*x;
    nn++;
  } 
  return(nn);  
}

int RPNExpressionEvaluator::ProductStack(double& px) 
{
  px = 1.;
  int nn = 0;
  double x = 0.;
  while( !rpnstack_.empty() ) {
    x = rpnstack_.top(); rpnstack_.pop();
    px *= x;  nn++;
  } 
  return(nn);  
}


} // End of namespace SOPHYA