/*-------------------------------------------------------------------------
| The software accompanies the paper                                      |
|                                                                         |
| Classes and Objects of Chemical Thermodynamics in Object-Oriented       |
| Programming. 2. A Class of Chemical Species                             |
|                                                                         |
| E.B. Rudnyi                                                             |
| E-mail: rudnyi@comp.chem.msu.su                                         |
| Homepages: http://www.chem.msu.su/~rudnyi/welcome.html                  |
|                                                                         |
| presented at Second Electronic Computational Chemistry Conference,      |
| November 1995, http://hackberry.chem.niu.edu/ECCC2/                     |
|                                                                         |
| E.B. Rudnyi, Copyright (c) 1995                                         |
|                                                                         |
| Permission to use, copy, modify, distribute and sell this software and  |
| its documentation for any purpose is hereby granted without fee,        |
| provided that the above copyright notice with the name of the paper and |
| the conference appear in all copies and that both that copyright notice |
| and this permission notice appear in supporting documentation.          |
| E.B. Rudnyi makes no representations about the suitability of this      |
| software for any purpose. It is provided "as is" without expressed or   |
| implied warranty.                                                       |
--------------------------------------------------------------------------*/

#include <iostream.h>
#include <strstrea.h>
#include <ctype.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "func_tp.h"

const int BUF_TOK = 200;
const int BUF_DBL = 40;
const int STACK = 20;

double global::T;
double global::p = 1.;
double global::step1 = 10*sqrt(DBL_EPSILON);
double global::step2 = 10*pow(DBL_EPSILON, 1./3.);
int global::digits = 8;

double (*func_Tp::func[NFUNC])(double x) = {acos, fabs, asin,
  atan, cos, exp, log10, log, sin, sqrt, tan};
char *func_Tp::func_name[NFUNC] = {"acos", "abs", "asin", "atan",
  "cos", "exp", "log10", "log", "sin", "sqrt", "tan"};
int func_Tp::func_len[NFUNC] = {4, 3, 4, 4, 3, 3, 5, 3, 3, 4, 3};
char func_Tp::tok_name[NSYMB] = {'+', '-', '*', '/', '^', 'T', 'p', '(', ')',
  '-'};
func_Tp::func_token* func_Tp::curr_tok;
func_Tp::func_token* func_Tp::token;
double* func_Tp::value;
func_Tp::func_err func_Tp::errno;

func_Tp::func_token func_Tp::get_token(const char *&in, double &value) {
  int i;

  while (isspace(*in)) in++;
  for (i = 0; i < NSYMB - 1; i++) {
    if (toupper(*in) == toupper(tok_name[i])) {
      if (tok_name[i] == 'T') {
        if (toupper(in[1]) == 'A' && toupper(in[2]) == 'N') {
          in += 3;
          return TAN;
        }
      }
      in++;
      return (func_token)(PLUS + i);
    }
  }
  switch(toupper(*in)) {
    case '0': case '1': case '2': case '3': case '4':
    case '5': case '6': case '7': case '8': case '9':
    case '.':
      char *end;
      value = strtod(in, &end);
      in = end;
      return NUMBER;
    case ':':
      in++;
      return END;
    case '\0': case ';':
      return END;
    default:
      for (i = 0; i < NFUNC; i++) {
        if (!strnicmp(in, func_name[i], func_len[i])) {
          in += func_len[i];
          return (func_token)i;
        }
      }
      errno = WRTOK;
      return END;
  }
}

void func_Tp::expr(int must) {
  func_token tmp;

  term(must);
  for(;;)
    switch (tmp = *curr_tok) {
      case PLUS: case MINUS:
        curr_tok++;
        term(2);
        *token++ = tmp;
        break;
      case END: case RP:
        return;
      default:
        errno = NOTERM;
        *curr_tok = END;
        return;
    }
}

void func_Tp::term(int must) {
  func_token tmp;

  fact(must);
  for(;;)
    switch (tmp = *curr_tok) {
      case MUL: case DIV:
        curr_tok++;
        fact(2);
        *token++ = tmp;
        break;
      case END: case PLUS: case MINUS: case RP:
	return;
      default:
        errno = NOFACT;
        *curr_tok = END;
        return;
   }
}

void func_Tp::fact(int must) {
  prim(must);
  for(;;)
    switch (*curr_tok) {
      case POW:
        curr_tok++;
        prim(2);
        *token++ = POW;
        break;
      case END: case PLUS: case MINUS: case MUL: case DIV: case RP:
        return;
      default:
        errno = NOPRIM;
        *curr_tok = END;
        return;
    }
}

void func_Tp::prim(int must) {
  func_token tmp;

  switch (*curr_tok) {
    case NUMBER: case NAME_T: case NAME_p:
      *token++ = *curr_tok++;
      return;
    case LP:
      *token++ = *curr_tok++;
      expr(1);
      if (*curr_tok != RP) {
        errno = NORP;
        *curr_tok = END;
        return;
      }
      *token++ = *curr_tok++;
      return;
    case ACOS: case ABS: case ASIN: case ATAN: case COS: case EXP:
    case LOG: case LOG10: case SIN: case SQRT: case TAN:
      tmp = *curr_tok++;
      if (*curr_tok != LP) {
	errno = NOLP;
        *curr_tok = END;
	return;
      }
      *token++ = *curr_tok++;
      expr(1);
      if (*curr_tok != RP) {
	errno = NORP;
        *curr_tok = END;
	return;
      }
      *token++ = *curr_tok++;
      *token++ = tmp;
      return;
    case MINUS:
      if (must < 2) {
        curr_tok++;
        prim(2);
        *token++ = UMIN;
      }
      else {
        *curr_tok = END;
        errno = NOPRIM;
      }
      return;
    case RP:
      errno = NOLP;
      *curr_tok = END;
      return;
    case END:
      if (must) errno = NOPRIM;
      return;
    default:
      *curr_tok = END;
      errno = NOPRIM;
      return;
  }
}

void func_Tp::exprpr(int print) const {
  func_token *tmp;

  termpr(print);
  for (;;) {
    if (token == start_tok + ntokens || *token == RP)
      return;
    if (token > start_tok + ntokens) {
      errno = BAD;
      return;
    }
    tmp = token;
    termpr(0);
    switch (*token) {
      case PLUS: case MINUS:
        if (print) {
          *curr_tok++ = *token;
          token = tmp;
          termpr(1);
        }
        token++;
        break;
      default:
        token = start_tok + ntokens + 1;
        break;
    }
  }
}

void func_Tp::termpr(int print) const {
  func_token *tmp;

  factpr(print);
  for(;;) {
    if (token >= start_tok + ntokens ||
             *token == PLUS || *token == MINUS || *token == RP)
      return;
    tmp = token;
    factpr(0);
    switch (*token) {
      case MUL: case DIV:
        if (print) {
          *curr_tok++ = *token;
          token = tmp;
          factpr(1);
        }
        token++;
        break;
      default:
        token = tmp;
        return;
    }
  }
}

void func_Tp::factpr(int print) const {
  func_token *tmp;

  primpr(print);
  for(;;) {
    if (token >= start_tok + ntokens ||
          *token == PLUS || *token == MINUS
          || *token == MUL || *token == DIV || *token == RP )
      return;
    tmp = token;
    primpr(0);
    switch (*token) {
      case POW:
        if (print) {
          *curr_tok++ = *token;
          token = tmp;
          primpr(1);
        }
        token++;
        break;
      default:
        token = tmp;
        return;
    }
  }
}

void func_Tp::primpr(int print) const {
  func_token *tmp;

  switch (*token) {
    case NUMBER: case NAME_T: case NAME_p:
      if (token[1] == UMIN) {
        if (print) {
          *curr_tok++ = UMIN;
          *curr_tok++ = *token;
        }
        token += 2;
      }
      else {
        if (print)
          *curr_tok++ = *token;
        token++;
      }
      return;
    case LP:
      tmp = ++token;
      exprpr(0);
      if (*token != RP) {
        token = start_tok + ntokens + 1;
        return;
      }
      token++;
      if (*token == UMIN || *token < PLUS) {
        if (print) {
          *curr_tok++ = *token;
          token = tmp;
          *curr_tok++ = LP;
          exprpr(1);
          *curr_tok++ = RP;
          token += 2;
        }
        else
          token++;
      }
      else if (print) {
        token = tmp;
        *curr_tok++ = LP;
        exprpr(1);
        *curr_tok++ = *token++;
      }
      return;
    default:
      token = start_tok + ntokens + 1;
      return;
  }
}

void func_Tp::copy(const func_Tp &old) {
  errno = OKAY;
  clear();
  if (!old.ntokens)
    return;
  if ((start_tok = new func_token[ntokens = old.ntokens]) == NULL) {
    init();
    return;
  }
  if (old.nvalues) {
    if ((start_val = new double[nvalues = old.nvalues]) == NULL) {
      delete [] start_tok;
      init();
      return;
    }
  }
  memcpy(start_tok, old.start_tok, sizeof(func_token)*ntokens);
  if (nvalues) memcpy(start_val, old.start_val, sizeof(double)*nvalues);
}

const char* func_Tp::anal(const char *in) {
  errno = OKAY;
  clear();
  if (*in == '\0')
    return in;

  func_token buf_tok[BUF_TOK];
  double buf_dbl[BUF_DBL];

  while ((buf_tok[ntokens] = get_token(in, buf_dbl[nvalues])) != END) {
    if (buf_tok[ntokens++] == NUMBER)
      nvalues++;
    if (ntokens == BUF_TOK || nvalues == BUF_DBL) {
      errno = NOBUF;
      break;
    }
  }
  if (errno != OKAY || ntokens == 0) {
    init();
    return in;
  }
  if ((start_tok = new func_token[ntokens]) == NULL) {
    init();
    return in;
  }
  if (nvalues) {
    if ((start_val = new double[nvalues]) == NULL) {
      delete [] start_tok;
      init();
      return in;
    }
  }
  else
    start_val = NULL;
  token = start_tok;
  curr_tok = buf_tok;
  if (nvalues) memcpy(start_val, buf_dbl, sizeof(double)*nvalues);
  expr(0);
  if (errno != OKAY)
    clear();
  return in;
}

char* func_Tp::to_str(char *buf) const {
  char *str = buf;

  errno = OKAY;
  if (!ntokens) {
    *buf = '\0';
    return buf;
  }

/*
// to get postfix notation
  value = start_val;
  for(token = start_tok; token < start_tok + ntokens; token++) {
    if (*token < PLUS)
      str = stpcpy(str, func_name[*token]);
    else if (*token == NUMBER) {
      gcvt(*value++, global::digits, str);
      str = str + strlen(str);
    }
    else
      *str++ = tok_name[*token - PLUS];
  }
  *str++ = '\n';
*/

  token = start_tok;
  value = start_val;
  func_token buffer[BUF_TOK];
  curr_tok = buffer;
  exprpr(1);

  for(token = buffer; token < buffer + ntokens; token++) {
    if (*token < PLUS)
      str = stpcpy(str, func_name[*token]);
    else if (*token == NUMBER) {
      gcvt(*value++, global::digits, str);
      str += strlen(str);
    }
    else
      *str++ = tok_name[*token - PLUS];
  }
  *str = '\0';
  return buf;
}

#define IFTOP   if (tp == STACK) {errno = NOSTK; return 0.;}
#define IFEMPTY if (!(--tp)) {errno = BAD; return 0.;}

double func_Tp::est(const double &T, const double &p) {
  errno = OKAY;
  if (!ntokens) return 0.;

  double stack[STACK];
  int tp = 0;

  value = start_val;
  for(token = start_tok; token < start_tok + ntokens; token++) {
    switch (*token) {
      case ACOS: case ABS: case ASIN: case ATAN: case COS: case EXP:
      case LOG: case LOG10: case SIN: case SQRT: case TAN:
        stack[tp - 1] = (*func[*token])(stack[tp - 1]);
        break;
      case UMIN:
        stack[tp - 1] = -stack[tp - 1];
        break;
      case PLUS:
        IFEMPTY
        stack[tp - 1] += stack[tp];
        break;
      case MINUS:
        IFEMPTY
        stack[tp - 1] -= stack[tp];
        break;
      case MUL:
        IFEMPTY
        stack[tp - 1] *= stack[tp];
        break;
      case DIV:
        IFEMPTY
        stack[tp - 1] /= stack[tp];
        break;
      case POW:
        IFEMPTY
        stack[tp - 1] = pow(stack[tp - 1], stack[tp]);
        break;
      case NUMBER:
        IFTOP
        stack[tp++] = *value++;
        break;
      case NAME_T:
        IFTOP
        stack[tp++] = T;
        break;
      case NAME_p:
        IFTOP
        stack[tp++] = p;
        break;
      case LP:
        break;
      case RP:
        break;
      default:
        errno = WRTOK;
        return 0.;
    }
  }
  if (--tp) {
    errno = BAD;
    return HUGE_VAL;
  }
  return stack[tp];
}

double func_Tp::df(global::derivative d, const double &T, const double &p) {
  double f1, f2, f3, f4;
  double stepT, stepp;
  switch (d) {
    case global::dT:
      stepT = global::step1*(T + global::step1);
      f1 = est(T - stepT, p);
      f2 = est(T + stepT, p);
      return (f2 - f1)/stepT/2.;
    case global::dp:
      stepp = global::step1*(p + global::step1);
      f1 = est(T, p - stepp);
      f2 = est(T, p + stepp);
      return (f2 - f1)/stepp/2.;
    case global::dT2:
      stepT = global::step2*(T + global::step2);
      f1 = est(T - stepT, p);
      f2 = est(T, p);
      f3 = est(T + stepT, p);
      return ((f3 - f2) - (f2 - f1))/stepT/stepT;
    case global::dp2:
      stepp = global::step2*(p + global::step2);
      f1 = est(T, p - stepp);
      f2 = est(T, p);
      f3 = est(T, p + stepp);
      return ((f3 - f2) - (f2 - f1))/stepp/stepp;
    case global::dTp: case global::dpT:
      stepT = global::step2*(T + global::step2);
      stepp = global::step2*(p + global::step2);
      f1 = est(T - stepT, p - stepp);
      f2 = est(T + stepT, p - stepp);
      f3 = est(T - stepT, p + stepp);
      f4 = est(T + stepT, p + stepp);
      return ((f4 - f3) - (f2 - f1))/stepT/stepp/4.;
    default:
      return HUGE_VAL;
  }
}

void oper(func_Tp &res, const func_Tp &left,
                      const func_Tp &right, const char ope) {
  func_Tp::func_token op;
  func_Tp::errno = func_Tp::OKAY;
  res.clear();
  switch (ope) {
    case '+':
      op = func_Tp::PLUS;
      break;
    case '-':
      op = func_Tp::MINUS;
      break;
    case '*':
      op = func_Tp::MUL;
      break;
    case '/':
      op = func_Tp::DIV;
      break;
    default:
      return;
  }
  if (left) {
    if (right) {
      if ((res.start_tok = new func_Tp::func_token[res.ntokens =
                           left.ntokens + right.ntokens + 5]) == NULL) {
        res.init();
        return;
      }
      if ((res.nvalues = left.nvalues + right.nvalues) != 0) {
        if ((res.start_val = new double[res.nvalues]) == NULL) {
          delete [] res.start_tok;
          res.init();
          return;
        }
      }
      else
        res.start_val = NULL;
      res.token = res.start_tok;
      *res.token++ = func_Tp::LP;
      memcpy(res.token, left.start_tok, sizeof(func_Tp::func_token)*left.ntokens);
      res.token += left.ntokens;
      *res.token++ = func_Tp::RP;
      *res.token++ = func_Tp::LP;
      memcpy(res.token, right.start_tok, sizeof(func_Tp::func_token)*right.ntokens);
      res.token += right.ntokens;
      *res.token++ = func_Tp::RP;
      *res.token++ = op;
      if (res.nvalues) {
        if (left.nvalues)
          memcpy(res.start_val, left.start_val, sizeof(double)*left.nvalues);
        if (right.nvalues)
          memcpy(res.start_val + left.nvalues,
                     right.start_val, sizeof(double)*right.nvalues);
      }
    }
    else
      res = left;
  }
  else if (right)
      res = right;
}

int func_Tp::len() const {
  int len = 0;
  char *buf = new char[global::digits+10];
  value = start_val;
  if (!ntokens) return len;
  for(token = start_tok; token < start_tok + ntokens; token++) {
    if (*token < PLUS)
      len = len + strlen(func_name[*token]);
    else if (*token == NUMBER)
      len = len + strlen(gcvt(*value++, global::digits, buf));
    else
      len++;
    }
  delete [] buf;
  return len;
}

ostream& operator<<(ostream &out, const func_Tp &old) {
  char *buf = new char[old.len() + 10];
  if (!buf) return out;
  out << old.to_str(buf);
  delete [] buf;
  return out;
}

istream& operator>>(istream &in, func_Tp &to) {
  ostrstream out;
  char ch;
  while (in >> ws, ch = in.peek(),
         (ch != EOF && ch != ';')) {
    in >> ch;
    if (ch == ':')
      break;
    out << ch;
  }
  out << ends;
  to.anal(out.str());
  delete out.str();
  return in;
}

func_Tp operator+(const func_Tp &left, const func_Tp &right) {
  func_Tp res;
  oper(res, left, right, '+');
  return res;
}

func_Tp operator-(const func_Tp &left, const func_Tp &right) {
  func_Tp res;
  oper(res, left, right, '-');
  return res;
}

func_Tp operator*(const func_Tp &left, const func_Tp &right) {
  func_Tp res;
  oper(res, left, right, '*');
  return res;
}

func_Tp operator/(const func_Tp &left, const func_Tp &right) {
  func_Tp res;
  oper(res, left, right, '/');
  return res;
}

lim_Tp::lim_err lim_Tp::errno;

lim_Tp::lim_token lim_Tp::get_token(const char *&in, double &value) {
  while (isspace(*in)) in++;
  switch(toupper(*in)) {
    case ':':
     in++;
     return END;
    case ',':
     in++;
     return COMMA;
    case '0': case '1': case '2': case '3': case '4':
    case '5': case '6': case '7': case '8': case '9':
    case '.':
      char *end;
      value = strtod(in, &end);
      in = end;
      return NUMBER;
    case '\0': case ';':
      return END;
    default:
      errno = WRTOK;
      return END;
  }
}

const char* lim_Tp::anal(const char *in) {
  double value;
  double *par = &Ti;
  lim_token curr_tok;

  errno = OKAY;
  init();
  for (int i = 0; i < 4; i++) {
    curr_tok = get_token(in, value);
    if (curr_tok == END) break;
    if (curr_tok == NUMBER) {
      par[i] = value;
      curr_tok = get_token(in, value);
      if (curr_tok == END) break;
      if (curr_tok != COMMA || i == 3) {
        errno = WRTOK;
        break;
      }
    }
    else if (curr_tok == COMMA && i != 3)
      continue;
    else {
      errno = WRTOK;
      break;
    }
  }
  if (errno != OKAY)
    init();
  if (Tf < Ti || pf < pi) {
    errno = WRORD;
    init();
  }
  return in;
}

char* lim_Tp::to_str(char *buf) const {
  char *str = buf;

  if (Ti != 0.) {
    gcvt(Ti, global::digits, str);
    str += strlen(str);
  }
  *str++ = ',';
  if (Tf != HUGE_VAL) {
    gcvt(Tf, global::digits, str);
    str += strlen(str);
  }
  *str++ = ',';
  if (pi != 0.) {
    gcvt(pi, global::digits, str);
    str += strlen(str);
  }
  *str++ = ',';
  if (pf != HUGE_VAL) {
    gcvt(pf, global::digits, str);
    str += strlen(str);
  }
  *str = '\0';
  return buf;
}

lim_Tp::lim_cmp lim_Tp::compare(const lim_Tp &second) const {
  if (Tf <= second.Ti)
    return LT;
  else if (Tf == second.Tf && Ti == second.Ti) {
      if (pf <= second.pi)
        return LT;
      else if (pf == second.pf && pi == second.pi)
        return EQ;
      else if (pi >= second.pf)
        return GT;
  }
  else if (Ti >= second.Tf)
    return GT;
  if (max(pi, second.pi) < min(pf, second.pf))
    return CROSS;
  else
    return CROSS_T;
}

int lim_Tp::len() const{
  char *buf = new char[global::digits + 10];
  int len = 3;
  if (Ti != 0.)
    len += strlen(gcvt(Ti, global::digits, buf));
  if (Tf != HUGE_VAL)
    len += strlen(gcvt(Tf, global::digits, buf));
  if (pi != 0.)
    len += strlen(gcvt(pi, global::digits, buf));
  if (pf != HUGE_VAL)
    len += strlen(gcvt(pf, global::digits, buf));
  delete [] buf;
  return len;
}

ostream& operator<<(ostream &out, const lim_Tp &old) {
  char *buf = new char[old.len() + 10];
  if (!buf) return out;
  out << old.to_str(buf);
  delete [] buf;
  return out;
}

istream& operator>>(istream &in, lim_Tp &to) {
  ostrstream out;
  char ch;
  while (in >> ws, ch = in.peek(),
         (ch != EOF && ch != ';')) {
    in >> ch;
    if (ch == ':')
      break;
    out << ch;
  }
  out << ends;
  to.anal(out.str());
  delete out.str();
  return in;
}

lim_Tp cross(const lim_Tp &first, const lim_Tp &second) {
  lim_Tp res;
  res.Ti = max(first.Ti, second.Ti);
  res.Tf = min(first.Tf, second.Tf);
  res.pi = max(first.pi, second.pi);
  res.pf = min(first.pf, second.pf);
  if (res.Tf < res.Ti || res.pf < res.pi) {
    res.init();
    lim_Tp::errno = lim_Tp::WRORD;
  }
  return res;
}