/* This file is part of the KDE project
   Copyright (C) 2005 Tomas Mecir <mecirt@gmail.com>

   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Library General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later version.

   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Library General Public License for more details.

   You should have received a copy of the GNU Library General Public License
   along with this library; see the file COPYING.LIB.  If not, write to
   the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 * Boston, MA 02110-1301, USA.
*/

#include "valuecalc.h"

#include "valueconverter.h"

#include <kdebug.h>
#include <errno.h>
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

using namespace KSpread;


// Array-walk functions registered on ValueCalc object

void awSum (ValueCalc *c, Value &res, Value val, Value)
{
  if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
    res = c->add (res, val);
}

void awSumA (ValueCalc *c, Value &res, Value val, Value)
{
  if (!val.isEmpty())
    res = c->add (res, val);
}

void awSumSq (ValueCalc *c, Value &res, Value val, Value)
{
  if (!val.isEmpty())
    res = c->add (res, c->sqr (val));
}

void awSumSqA (ValueCalc *c, Value &res, Value val, Value)
{
  if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
    res = c->add (res, c->sqr (val));
}

void awCount (ValueCalc *c, Value &res, Value val, Value)
{
  if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
    res = c->add (res, 1);
}

void awCountA (ValueCalc *c, Value &res, Value val, Value)
{
  if (!val.isEmpty())
    res = c->add (res, 1);
}

void awMax (ValueCalc *c, Value &res, Value val, Value)
{
  if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
    if (res.isEmpty())
      res = val;
    else
      if (c->greater (val, res)) res = val;
}

void awMaxA (ValueCalc *c, Value &res, Value val, Value)
{
  if (!val.isEmpty())
    if (res.isEmpty())
      // convert to number, so that we don't return string/bool
      res = c->conv()->asNumeric (val);
    else
      if (c->greater (val, res))
        // convert to number, so that we don't return string/bool
        res = c->conv()->asNumeric (val);
}

void awMin (ValueCalc *c, Value &res, Value val, Value)
{
  if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
    if (res.isEmpty())
      res = val;
    else
      if (c->lower (val, res)) res = val;
}

void awMinA (ValueCalc *c, Value &res, Value val, Value)
{
  if (!val.isEmpty())
    if (res.isEmpty())
      // convert to number, so that we don't return string/bool
      res = c->conv()->asNumeric (val);
    else
      if (c->lower (val, res))
        // convert to number, so that we don't return string/bool
        res = c->conv()->asNumeric (val);
}

void awProd (ValueCalc *c, Value &res, Value val, Value)
{
  if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
    res = c->mul (res, val);
}

void awProdA (ValueCalc *c, Value &res, Value val, Value)
{
  if (!val.isEmpty())
    res = c->mul (res, val);
}

// sum of squares of deviations, used to compute standard deviation
void awDevSq (ValueCalc *c, Value &res, Value val,
    Value avg)
{
  if (!val.isEmpty())
    res = c->add (res, c->sqr (c->sub (val, avg)));
}

// sum of squares of deviations, used to compute standard deviation
void awDevSqA (ValueCalc *c, Value &res, Value val,
    Value avg)
{
  if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
    res = c->add (res, c->sqr (c->sub (val, avg)));
}


bool isDate (Value val) {
  Value::Format fmt = val.format();
  if ((fmt == Value::fmt_Date) || (fmt == Value::fmt_DateTime))
    return true;
  return false;
}

// ***********************
// ****** ValueCalc ******
// ***********************

ValueCalc::ValueCalc (ValueConverter* c): converter( c )
{
  // initialize the random number generator
  srand (time (0));
  
  // register array-walk functions
  registerAwFunc ("sum", awSum);
  registerAwFunc ("suma", awSumA);
  registerAwFunc ("sumsq", awSumSq);
  registerAwFunc ("sumsqa", awSumSqA);
  registerAwFunc ("count", awCount);
  registerAwFunc ("counta", awCountA);
  registerAwFunc ("max", awMax);
  registerAwFunc ("maxa", awMaxA);
  registerAwFunc ("min", awMin);
  registerAwFunc ("mina", awMinA);
  registerAwFunc ("prod", awProd);
  registerAwFunc ("proda", awProdA);
  registerAwFunc ("devsq", awDevSq);
  registerAwFunc ("devsqa", awDevSq);
}

Value ValueCalc::add (const Value &a, const Value &b)
{
  if (a.isError()) return a;
  if (b.isError()) return b;
  Value res;
  if (a.isInteger() && b.isEmpty() || a.isEmpty() && b.isInteger() 
      || a.isInteger() && b.isInteger())
  {
    int aa, bb;
    aa = converter->asInteger (a).asInteger();
    bb = converter->asInteger (b).asInteger();
    res = Value (aa + bb);
  }
  else
  {
    double aa, bb;
    aa = converter->asFloat (a).asFloat();
    bb = converter->asFloat (b).asFloat();
    res = Value (aa + bb);
  }

  if (a.isNumber() || a.isEmpty())
    res.setFormat (format (a.format(), b.format()));
  // operation on two dates should produce a number
  if (isDate(a) && isDate(b))
    res.setFormat (Value::fmt_Number);

  return res;
}

Value ValueCalc::sub (const Value &a, const Value &b)
{
  if (a.isError()) return a;
  if (b.isError()) return b;
  double aa, bb;
  aa = converter->asFloat (a).asFloat();
  bb = converter->asFloat (b).asFloat();
  Value res = Value (aa - bb);

  if (a.isNumber() || a.isEmpty())
    res.setFormat (format (a.format(), b.format()));
  // operation on two dates should produce a number
  if (isDate(a) && isDate(b))
    res.setFormat (Value::fmt_Number);

  return res;
}

Value ValueCalc::mul (const Value &a, const Value &b)
{
  if (a.isError()) return a;
  if (b.isError()) return b;
  double aa, bb;
  aa = converter->asFloat (a).asFloat();
  bb = converter->asFloat (b).asFloat();
  Value res = Value (aa * bb);

  if (a.isNumber() || a.isEmpty())
    res.setFormat (format (a.format(), b.format()));
  // operation on two dates should produce a number
  if (isDate(a) && isDate(b))
    res.setFormat (Value::fmt_Number);

  return res;
}

Value ValueCalc::div (const Value &a, const Value &b)
{
  if (a.isError()) return a;
  if (b.isError()) return b;
  double aa, bb;
  aa = converter->asFloat (a).asFloat();
  bb = converter->asFloat (b).asFloat();
  Value res;
  if (bb == 0.0)
    return Value::errorDIV0();
  else
    res = Value (aa / bb);

  if (a.isNumber() || a.isEmpty())
    res.setFormat (format (a.format(), b.format()));
  // operation on two dates should produce a number
  if (isDate(a) && isDate(b))
    res.setFormat (Value::fmt_Number);

  return res;
}

Value ValueCalc::mod (const Value &a, const Value &b)
{
  if (a.isError()) return a;
  if (b.isError()) return b;
  double aa, bb;
  aa = converter->asFloat (a).asFloat();
  bb = converter->asFloat (b).asFloat();
  Value res;
  if (bb == 0.0)
    return Value::errorDIV0();
  else
  {
    double m = fmod (aa, bb);
    // the following adjustments are needed by OpenFormula:
    // can't simply use fixed increases/decreases, because the implementation
    // of fmod may differ on various platforms, and we should always return
    // the same results ...
    if ((bb > 0) && (aa < 0))  // result must be positive here
      while (m < 0) m += bb;
    if (bb < 0)  // result must be negative here, but not lower than bb
    {
      // bb is negative, hence the following two are correct
      while (m < bb) m -= bb;  // same as m+=fabs(bb)
      while (m > 0) m += bb;   // same as m-=fabs(bb)
    }

    res = Value (m);
  }

  if (a.isNumber() || a.isEmpty())
    res.setFormat (format (a.format(), b.format()));
  if (isDate(a) && isDate(b))
    res.setFormat (Value::fmt_Number);

  return res;
}

Value ValueCalc::pow (const Value &a, const Value &b)
{
  if (a.isError()) return a;
  if (b.isError()) return b;
  double aa, bb;
  aa = converter->asFloat (a).asFloat();
  bb = converter->asFloat (b).asFloat();
  Value res = Value (::pow (aa, bb));

  if (a.isNumber() || a.isEmpty())
    res.setFormat (format (a.format(), b.format()));
  // operation on date(s) should produce a number
  if (isDate(a) || isDate(b))
    res.setFormat (Value::fmt_Number);

  return res;
}

Value ValueCalc::sqr (const Value &a)
{
  if (a.isError()) return a;
  return mul (a, a);
}

Value ValueCalc::sqrt (const Value &a)
{
  if (a.isError()) return a;
  Value res = Value (::sqrt (converter->asFloat(a).asFloat()));
  if (a.isNumber() || a.isEmpty())
    res.setFormat (a.format());
  // operation on date(s) should produce a number
  if (isDate(a))
    res.setFormat (Value::fmt_Number);

  return res;
}

Value ValueCalc::add (const Value &a, double b)
{
  if (a.isError()) return a;
  Value res = Value (converter->asFloat(a).asFloat() + b);

  if (a.isNumber() || a.isEmpty())
    res.setFormat (a.format());

  return res;
}

Value ValueCalc::sub (const Value &a, double b)
{
  if (a.isError()) return a;
  Value res = Value (converter->asFloat(a).asFloat() - b);

  if (a.isNumber() || a.isEmpty())
    res.setFormat (a.format());

  return res;
}

Value ValueCalc::mul (const Value &a, double b)
{
  if (a.isError()) return a;
  Value res = Value (converter->asFloat(a).asFloat() * b);

  if (a.isNumber() || a.isEmpty())
    res.setFormat (a.format());

  return res;
}

Value ValueCalc::div (const Value &a, double b)
{
  if (a.isError()) return a;
  Value res;
  if (b == 0.0)
    return Value::errorDIV0();

  res = Value (converter->asFloat(a).asFloat() / b);

  if (a.isNumber() || a.isEmpty())
    res.setFormat (a.format());

  return res;
}

Value ValueCalc::pow (const Value &a, double b)
{
  if (a.isError()) return a;
  Value res = Value (::pow (converter->asFloat(a).asFloat(), b));

  if (a.isNumber() || a.isEmpty())
    res.setFormat (a.format());

  return res;
}

Value ValueCalc::abs (const Value &a)
{
  if (a.isError()) return a;
  return Value (fabs (converter->asFloat (a).asFloat()));
}

bool ValueCalc::isZero (const Value &a)
{
  if (a.isError()) return false;
  return (converter->asFloat (a).asFloat() == 0.0);
}

bool ValueCalc::isEven (const Value &a)
{
  if (a.isError()) return false;
  return ((converter->asInteger (a).asInteger() % 2) == 0);
}

bool ValueCalc::equal (const Value &a, const Value &b)
{
  return (converter->asFloat (a).asFloat() == converter->asFloat (b).asFloat());
}

/*********************************************************************
 *
 * Helper function to avoid problems with rounding floating point
 * values. Idea for this kind of solution taken from Openoffice.
 *
 *********************************************************************/
bool ValueCalc::approxEqual (const Value &a, const Value &b)
{
  double aa = converter->asFloat (a).asFloat();
  double bb = converter->asFloat (b).asFloat();
  if (aa == bb)
    return true;
  double x = aa - bb;
  return (x < 0.0 ? -x : x)  <  ((aa < 0.0 ? -aa : aa) * DBL_EPSILON);
}

bool ValueCalc::strEqual (const Value &a, const Value &b)
{
  return (converter->asString (a).asString() == converter->asString (b).asString());
}

bool ValueCalc::greater (const Value &a, const Value &b)
{
  double aa = converter->asFloat (a).asFloat();
  double bb = converter->asFloat (b).asFloat();
  return (aa > bb);
}

bool ValueCalc::gequal (const Value &a, const Value &b)
{
  return (greater (a,b) || approxEqual (a,b));
}

bool ValueCalc::lower (const Value &a, const Value &b)
{
  return greater (b, a);
}

Value ValueCalc::roundDown (const Value &a,
    const Value &digits) {
  return roundDown (a, converter->asInteger (digits).asInteger());
}

Value ValueCalc::roundUp (const Value &a,
    const Value &digits) {
  return roundUp (a, converter->asInteger (digits).asInteger());
}

Value ValueCalc::round (const Value &a,
    const Value &digits) {
  return round (a, converter->asInteger (digits).asInteger());
}

Value ValueCalc::roundDown (const Value &a, int digits)
{
  // shift in one direction, round, shift back
  Value val = a;
  if (digits > 0)
    for (int i = 0; i < digits; ++i)
      val = mul (val, 10);
  if (digits < 0)
    for (int i = 0; i < digits; ++i)
      val = div (val, 10);
  
  val = Value (floor (converter->asFloat (val).asFloat()));
  
  if (digits > 0)
    for (int i = 0; i < digits; ++i)
      val = div (val, 10);
  if (digits < 0)
    for (int i = 0; i < digits; ++i)
      val = mul (val, 10);
  return val;
}

Value ValueCalc::roundUp (const Value &a, int digits)
{
  // shift in one direction, round, shift back
  Value val = a;
  if (digits > 0)
    for (int i = 0; i < digits; ++i)
      val = mul (val, 10);
  if (digits < 0)
    for (int i = 0; i < digits; ++i)
      val = div (val, 10);
  
  val = Value (ceil (converter->asFloat (val).asFloat()));
  
  if (digits > 0)
    for (int i = 0; i < digits; ++i)
      val = div (val, 10);
  if (digits < 0)
    for (int i = 0; i < digits; ++i)
      val = mul (val, 10);
  return val;
}

Value ValueCalc::round (const Value &a, int digits)
{
  // shift in one direction, round, shift back
  Value val = a;
  if (digits > 0)
    for (int i = 0; i < digits; ++i)
      val = mul (val, 10);
  if (digits < 0)
    for (int i = 0; i < digits; ++i)
      val = div (val, 10);

  val = Value (int(converter->asFloat (val).asFloat()+0.5));

  if (digits > 0)
    for (int i = 0; i < digits; ++i)
      val = div (val, 10);
  if (digits < 0)
    for (int i = 0; i < digits; ++i)
      val = mul (val, 10);
  return val;
}

int ValueCalc::sign (const Value &a)
{
  double val = converter->asFloat (a).asFloat ();
  if (val == 0) return 0;
  if (val > 0) return 1;
  return -1;
}


Value ValueCalc::log (const Value &number,
    const Value &base)
{
  double logbase = converter->asFloat (base).asFloat();
  if (logbase == 1.0)
    return Value::errorDIV0();
  if (logbase <= 0.0)
    return Value::errorNA();

  logbase = log10 (logbase);
  Value res = Value (log10 (converter->asFloat (number).asFloat()) / logbase);

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::ln (const Value &number)
{
  Value res = Value (::log (converter->asFloat (number).asFloat()));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::log (const Value &number, double base)
{
  if (base <= 0.0)
    return Value::errorNA();
  if (base == 1.0)
    return Value::errorDIV0();

  double num = converter->asFloat (number).asFloat();
  Value res = Value (log10 (num) / log10 (base));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::exp (const Value &number)
{
  return Value (::exp (converter->asFloat (number).asFloat()));
}

Value ValueCalc::pi ()
{
  // retun PI in double-precision
  // if arbitrary precision gets in, this should be extended to return
  // more if need be
  return Value (M_PI);
}

Value ValueCalc::eps ()
{
  // #### This should adjust according to the actual number system used
  // (float, double, long double, ...)
  return Value (DBL_EPSILON);
}

Value ValueCalc::random (double range)
{
  return Value (range * (double) rand() / (RAND_MAX + 1.0));
}

Value ValueCalc::random (Value range)
{
  return random (converter->asFloat (range).asFloat());
}

Value ValueCalc::fact (const Value &which)
{
  // we can simply use integers - no one is going to compute factorial of
  // anything bigger than 2^32
  return fact (converter->asInteger (which).asInteger());
}

Value ValueCalc::fact (const Value &which,
    const Value &end)
{
  // we can simply use integers - no one is going to compute factorial of
  // anything bigger than 2^32
  return fact (converter->asInteger (which).asInteger(),
      converter->asInteger (end).asInteger ());
}

Value ValueCalc::fact (int which, int end) {
  if (which < 0)
    return Value (-1);
  if (which == 0)
    return Value (1);
  // no multiplication if val==end
  if (which == end)
    return Value (1);
    
  return (mul (fact (which-1, end), which));
}

Value ValueCalc::factDouble (int which)
{
  if (which < 0)
    return Value (-1);
  if ((which == 0) || (which == 1))
    return Value (1);
    
  return (mul (factDouble (which-2), which));
}

Value ValueCalc::factDouble (Value which)
{
  return factDouble (converter->asInteger (which).asInteger());
}

Value ValueCalc::combin (int n, int k)
{
  if (n >= 15)
  {
    double result = ::exp(lgamma (n + 1) - lgamma (k + 1) - lgamma (n-k+1));
    return Value (floor(result + 0.5));
  }
  else
    return div (div (fact (n), fact (k)), fact (n - k));
}

Value ValueCalc::combin (Value n, Value k)
{
  int nn = converter->asInteger (n).asInteger();
  int kk = converter->asInteger (k).asInteger();
  return combin (nn, kk);
}

Value ValueCalc::gcd (const Value &a, const Value &b)
{
  // Euler's GCD algorithm
  Value aa = round (a);
  Value bb = round (b);

  if (approxEqual (aa, bb)) return aa;
  
  if (aa.isZero()) return bb;
  if (bb.isZero()) return aa;
  

  if (greater (aa, bb))
    return gcd (bb, mod (aa, bb));
  else
    return gcd (aa, mod (bb, aa));
}

Value ValueCalc::lcm (const Value &a, const Value &b)
{
  Value aa = round (a);
  Value bb = round (b);

  if (approxEqual (aa, bb)) return aa;
  
  if (aa.isZero()) return bb;
  if (bb.isZero()) return aa;
  
  Value g = gcd (aa, bb);
  if (g.isZero())  // GCD is zero for some weird reason
    return mul (aa, bb);

  return div (mul (aa, bb), g);
}

Value ValueCalc::base (const Value &val, int base, int prec)
{
  if (base == 10) return round (val, prec);
  if (prec < 0) prec = 2;
  if ((base < 2) || (base > 36))
    return Value::errorVALUE();

  double value = converter->asFloat (val).asFloat();
  TQString result = TQString::number ((int)value, base);

  if (prec > 0)
  {
    result += "."; value = value - (int)value;

    int ix;
    for( int i = 0; i < prec; i++ )
    {
      ix = (int) value * base;
      result += "0123456789abcdefghijklmnopqrstuvwxyz"[ix];
      value = base * (value - (double)ix/base);
    }
  }

  return Value (result.upper());
}

Value ValueCalc::fromBase (const Value &val, int base)
{
  TQString str = converter->asString (val).asString();
  bool ok;
  double num = str.toLong (&ok, base);
  if (ok)
    return Value (num);
  return Value::errorVALUE();
}


Value ValueCalc::sin (const Value &number)
{
  Value res = Value (::sin (converter->asFloat (number).asFloat()));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::cos (const Value &number)
{
  Value res = Value (::cos (converter->asFloat (number).asFloat()));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::tg (const Value &number)
{
  Value res = Value (::tan (converter->asFloat (number).asFloat()));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::cotg (const Value &number)
{
  Value res = Value (div (1, ::tan (converter->asFloat (number).asFloat())));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::asin (const Value &number)
{
  errno = 0;
  Value res = Value (::asin (converter->asFloat (number).asFloat()));
  if (errno)
    return Value::errorVALUE();

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::acos (const Value &number)
{
  errno = 0;
  Value res = Value (::acos (converter->asFloat (number).asFloat()));
  if (errno)
    return Value::errorVALUE();

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::atg (const Value &number)
{
  errno = 0;
  Value res = Value (::atan (converter->asFloat (number).asFloat()));
  if (errno)
    return Value::errorVALUE();

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::atan2 (const Value &y, const Value &x)
{
  double yy = converter->asFloat (y).asFloat();
  double xx = converter->asFloat (x).asFloat();
  return Value (::atan2 (yy, xx));
}

Value ValueCalc::sinh (const Value &number)
{
  Value res = Value (::sinh (converter->asFloat (number).asFloat()));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::cosh (const Value &number)
{
  Value res = Value (::cosh (converter->asFloat (number).asFloat()));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::tgh (const Value &number)
{
  Value res = Value (::tanh (converter->asFloat (number).asFloat()));

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::asinh (const Value &number)
{
  errno = 0;
  Value res = Value (::asinh (converter->asFloat (number).asFloat()));
  if (errno)
    return Value::errorVALUE();

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::acosh (const Value &number)
{
  errno = 0;
  Value res = Value (::acosh (converter->asFloat (number).asFloat()));
  if (errno)
    return Value::errorVALUE();

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::atgh (const Value &number)
{
  errno = 0;
  Value res = Value (::atanh (converter->asFloat (number).asFloat()));
  if (errno)
    return Value::errorVALUE();

  if (number.isNumber() || number.isEmpty())
    res.setFormat (number.format());

  return res;
}

Value ValueCalc::phi (Value x)
{
  Value constant (0.39894228040143268);
  
  // constant * exp(-(x * x) / 2.0);
  Value x2neg = mul (sqr (x), -1);
  return mul (constant, exp (div (x2neg, 2.0)));
}

static double taylor_helper (double* pPolynom, uint nMax, double x)
{
  double nVal = pPolynom[nMax];
  for (int i = nMax-1; i >= 0; i--) {
    nVal = pPolynom[i] + (nVal * x);
  }
  return nVal;
}

Value ValueCalc::gauss (Value xx)
// this is a weird function
{
  double x = converter->asFloat (xx).asFloat();
  
  double t0[] =
    { 0.39894228040143268, -0.06649038006690545,  0.00997355701003582,
     -0.00118732821548045,  0.00011543468761616, -0.00000944465625950,
      0.00000066596935163, -0.00000004122667415,  0.00000000227352982,
      0.00000000011301172,  0.00000000000511243, -0.00000000000021218 };
  double t2[] =
    { 0.47724986805182079,  0.05399096651318805, -0.05399096651318805,
      0.02699548325659403, -0.00449924720943234, -0.00224962360471617,
      0.00134977416282970, -0.00011783742691370, -0.00011515930357476,
      0.00003704737285544,  0.00000282690796889, -0.00000354513195524,
      0.00000037669563126,  0.00000019202407921, -0.00000005226908590,
     -0.00000000491799345,  0.00000000366377919, -0.00000000015981997,
     -0.00000000017381238,  0.00000000002624031,  0.00000000000560919,
     -0.00000000000172127, -0.00000000000008634,  0.00000000000007894 };
  double t4[] =
    { 0.49996832875816688,  0.00013383022576489, -0.00026766045152977,
      0.00033457556441221, -0.00028996548915725,  0.00018178605666397,
     -0.00008252863922168,  0.00002551802519049, -0.00000391665839292,
     -0.00000074018205222,  0.00000064422023359, -0.00000017370155340,
      0.00000000909595465,  0.00000000944943118, -0.00000000329957075,
      0.00000000029492075,  0.00000000011874477, -0.00000000004420396,
      0.00000000000361422,  0.00000000000143638, -0.00000000000045848 };
  double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };

  double xAbs = fabs(x);
  uint xShort = static_cast<uint>(floor(xAbs));
  double nVal = 0.0;
  if (xShort == 0)
    nVal = taylor_helper(t0, 11, (xAbs * xAbs)) * xAbs;
  else if ((xShort >= 1) && (xShort <= 2))
    nVal = taylor_helper(t2, 23, (xAbs - 2.0));
  else if ((xShort >= 3) && (xShort <= 4))
    nVal = taylor_helper(t4, 20, (xAbs - 4.0));
  else
  {
    double phiAbs = converter->asFloat (phi (xAbs)).asFloat();
    nVal = 0.5 + phiAbs * taylor_helper(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs;
  }

  if (x < 0.0)
    return Value (-nVal);
  else
    return Value (nVal);
}

Value ValueCalc::gaussinv (Value xx)
// this is a weird function
{
  double x = converter->asFloat (xx).asFloat();
  
  double q,t,z;

  q=x-0.5;

  if(fabs(q)<=.425)
  {
    t=0.180625-q*q;

    z=
    q*
    (
      (
        (
          (
            (
              (
                (
                  t*2509.0809287301226727+33430.575583588128105
                )
                *t+67265.770927008700853
              )
              *t+45921.953931549871457
            )
            *t+13731.693765509461125
          )
          *t+1971.5909503065514427
        )
        *t+133.14166789178437745
      )
      *t+3.387132872796366608
    )
    /
    (
      (
        (
          (
            (
              (
                (
                  t*5226.495278852854561+28729.085735721942674
                )
                *t+39307.89580009271061
              )
              *t+21213.794301586595867
            )
            *t+5394.1960214247511077
          )
          *t+687.1870074920579083
        )
        *t+42.313330701600911252
      )
      *t+1.0
    );

  }
  else
  {
    if(q>0)  t=1-x;
    else    t=x;

    t=::sqrt(-::log(t));

    if(t<=5.0)
    {
      t+=-1.6;

      z=
      (
        (
          (
            (
              (
                (
                  (
                    t*7.7454501427834140764e-4+0.0227238449892691845833
                  )
                  *t+0.24178072517745061177
                )
                *t+1.27045825245236838258
              )
              *t+3.64784832476320460504
            )
            *t+5.7694972214606914055
          )
          *t+4.6303378461565452959
        )
        *t+1.42343711074968357734
      )
      /
      (
        (
          (
            (
              (
                (
                  (
                    t*1.05075007164441684324e-9+5.475938084995344946e-4
                  )
                  *t+0.0151986665636164571966
                )
                *t+0.14810397642748007459
              )
              *t+0.68976733498510000455
            )
            *t+1.6763848301838038494
          )
          *t+2.05319162663775882187
        )
        *t+1.0
      );

    }
    else
    {
      t+=-5.0;

      z=
      (
        (
          (
            (
              (
                (
                  (
                    t*2.01033439929228813265e-7+2.71155556874348757815e-5
                  )
                  *t+0.0012426609473880784386
                )
                *t+0.026532189526576123093
              )
              *t+0.29656057182850489123
            )
            *t+1.7848265399172913358
          )
          *t+5.4637849111641143699
        )
        *t+6.6579046435011037772
      )
      /
      (
        (
          (
            (
              (
                (
                  (
                    t*2.04426310338993978564e-15+1.4215117583164458887e-7
                  )
                  *t+1.8463183175100546818e-5
                )
                *t+7.868691311456132591e-4
              )
              *t+0.0148753612908506148525
            )
            *t+0.13692988092273580531
          )
          *t+0.59983220655588793769
        )
        *t+1.0
      );

    }

    if(q<0.0) z=-z;
  }

  return Value (z);
}

//helper for GetGamma and GetLogGamma
static double GammaHelp(double& x, bool& bReflect)
{
  double c[6] = {76.18009173, -86.50532033, 24.01409822,
                 -1.231739516, 0.120858003E-2, -0.536382E-5};
  if (x >= 1.0)
    {
      bReflect = false;
      x -= 1.0;
    }
  else
    {
      bReflect = true;
      x = 1.0 - x;
    }
  double s, anum;
  s = 1.0;
  anum = x;
  for (uint i = 0; i < 6; i++)
    {
      anum += 1.0;
      s += c[i]/anum;
    }
  s *= 2.506628275;   // sqrt(2*PI)
  return s;
}

Value ValueCalc::GetGamma (Value _x)
{
  double x = converter->asFloat (_x).asFloat();

  bool bReflect;
  double G = GammaHelp(x, bReflect);
  G = ::pow(x+5.5,x+0.5)*G/::exp(x+5.5);
  if (bReflect)
    G = M_PI*x/(G*::sin(M_PI*x));
  return Value (G);
}

Value ValueCalc::GetLogGamma (Value _x)
{
  double x = converter->asFloat (_x).asFloat();

  bool bReflect;
  double G = GammaHelp(x, bReflect);
  G = (x+0.5)*::log(x+5.5)+::log(G)-(x+5.5);
  if (bReflect)
    G = ::log(M_PI*x)-G-::log(::sin(M_PI*x));
  return Value (G);
}

Value ValueCalc::GetGammaDist (Value _x, Value _alpha,
    Value _beta)
{
  double x = converter->asFloat (_x).asFloat();
  double alpha = converter->asFloat (_alpha).asFloat();
  double beta = converter->asFloat (_beta).asFloat();
  
  if (x == 0.0)
    return Value (0.0);

  x /= beta;
  double gamma = alpha;

  double c = 0.918938533204672741;
  double d[10] = {
    0.833333333333333333E-1,
    -0.277777777777777778E-2,
    0.793650793650793651E-3,
    -0.595238095238095238E-3,
    0.841750841750841751E-3,
    -0.191752691752691753E-2,
    0.641025641025641025E-2,
    -0.295506535947712418E-1,
    0.179644372368830573,
    -0.139243221690590111E1
  };

  double dx = x;
  double dgamma = gamma;
  int maxit = 10000;

  double z = dgamma;
  double den = 1.0;
  while ( z < 10.0 ) {
    den *= z;
    z += 1.0;
  }

  double z2 = z*z;
  double z3 = z*z2;
  double z4 = z2*z2;
  double z5 = z2*z3;
  double a = ( z - 0.5 ) * ::log(z) - z + c;
  double b = d[0]/z + d[1]/z3 + d[2]/z5 + d[3]/(z2*z5) + d[4]/(z4*z5) +
    d[5]/(z*z5*z5) + d[6]/(z3*z5*z5) + d[7]/(z5*z5*z5) + d[8]/(z2*z5*z5*z5);
  // double g = exp(a+b) / den;

  double sum = 1.0 / dgamma;
  double term = 1.0 / dgamma;
  double cut1 = dx - dgamma;
  double cut2 = dx * 10000000000.0;

  for ( int i=1; i<=maxit; i++ ) {
    double ai = i;
    term = dx * term / ( dgamma + ai );
    sum += term;
    double cutoff = cut1 + ( cut2 * term / sum );
    if ( ai > cutoff ) {
      double t = sum;
      // return pow( dx, dgamma ) * exp( -dx ) * t / g;
      return Value (::exp( dgamma * ::log(dx) - dx - a - b ) * t * den);
    }
  }

  return Value (1.0);             // should not happen ...
}

Value ValueCalc::GetBeta (Value _x, Value _alpha,
    Value _beta)
{
  if (equal (_beta, 1.0))
    return pow (_x, _alpha);
  else if (equal (_alpha, 1.0))
    // 1.0 - pow (1.0-_x, _beta)
    return sub (1.0, pow (sub (1.0, _x), _beta));

  double x = converter->asFloat (_x).asFloat();
  double alpha = converter->asFloat (_alpha).asFloat();
  double beta = converter->asFloat (_beta).asFloat();
    
  double fEps = 1.0E-8;
  bool bReflect;
  double cf, fA, fB;

  if (x < (alpha+1.0)/(alpha+beta+1.0)) {
    bReflect = false;
    fA = alpha;
    fB = beta;
  }
  else {
    bReflect = true;
    fA = beta;
    fB = alpha;
    x = 1.0 - x;
  }
  if (x < fEps)
    cf = 0.0;
  else {
    double a1, b1, a2, b2, fnorm, rm, apl2m, d2m, d2m1, cfnew;
    a1 = 1.0; b1 = 1.0;
    b2 = 1.0 - (fA+fB)*x/(fA+1.0);
    if (b2 == 0.0) {
      a2 = b2;
      fnorm = 1.0;
      cf = 1.0;
    }
    else {
      a2 = 1.0;
      fnorm = 1.0/b2;
      cf = a2*fnorm;
    }
    cfnew = 1.0;
    for (uint j = 1; j <= 100; j++) {
      rm = (double) j;
      apl2m = fA + 2.0*rm;
      d2m = rm*(fB-rm)*x/((apl2m-1.0)*apl2m);
      d2m1 = -(fA+rm)*(fA+fB+rm)*x/(apl2m*(apl2m+1.0));
      a1 = (a2+d2m*a1)*fnorm;
      b1 = (b2+d2m*b1)*fnorm;
      a2 = a1 + d2m1*a2*fnorm;
      b2 = b1 + d2m1*b2*fnorm;
      if (b2 != 0.0) {
        fnorm = 1.0/b2;
        cfnew = a2*fnorm;
        if (fabs(cf-cfnew)/cf < fEps)
          j = 101;
        else
          cf = cfnew;
      }
    }
    if (fB < fEps)
      b1 = 1.0E30;
    else
      b1 = ::exp(GetLogGamma(fA).asFloat()+GetLogGamma(fB).asFloat()-
          GetLogGamma(fA+fB).asFloat());

    cf *= ::pow(x, fA)*::pow(1.0-x,fB)/(fA*b1);
  }
  if (bReflect)
    return Value (1.0-cf);
  else
    return Value (cf);
}

// ------------------------------------------------------


/*
 *
 * The code for calculating Bessel functions is taken
 * from CCMATH, a mathematics library source.code.
 *
 * Original copyright follows:
 *
 *  Copyright (C)  2000   Daniel A. Atkinson    All rights reserved.
 *  This code may be redistributed under the terms of the GNU library
 *  public license (LGPL).
 */

static double ccmath_gaml(double x)
{ double g,h;
  for(g=1.; x<30. ;g*=x,x+=1.); h=x*x;
  g=(x-.5)*log(x)-x+.918938533204672-log(g);
  g+=(1.-(1./6.-(1./3.-1./(4.*h))/(7.*h))/(5.*h))/(12.*x);
  return g;
}

static double ccmath_psi(int m)
{ double s= -.577215664901533; int k;
  for(k=1; k<m ;++k) s+=1./k;
  return s;
}

static double ccmath_ibes(double v,double x)
{ double y,s,t,tp; int p,m;
  y=x-9.; if(y>0.) y*=y; tp=v*v*.2+25.;
  if(y<tp){ x/=2.; m=(int)x;
    if(x>0.) s=t=exp(v*log(x)-ccmath_gaml(v+1.));
    else{ if(v>0.) return 0.; else if(v==0.) return 1.;}
    for(p=1,x*=x;;++p){ t*=x/(p*(v+=1.)); s+=t;
      if(p>m && t<1.e-13*s) break;
     }
   }
  else{ double u,a0=1.57079632679490;
    s=t=1./sqrt(x*a0); x*=2.; u=0.;
    for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
      t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>=tp) break;
      if(!(p&1)) s+=t; else u-=t;
     }
    x/=2.; s=cosh(x)*s+sinh(x)*u;
   }
  return s;
}

static double ccmath_kbes(double v,double x)
{ double y,s,t,tp,f,a0=1.57079632679490;
  int p,k,m;
  if(x==0.) return HUGE_VAL;
  y=x-10.5; if(y>0.) y*=y; tp=25.+.185*v*v;
  if(y<tp && modf(v+.5,&t)!=0.){ y=1.5+.5*v;
    if(x<y){ x/=2.; m=(int)x; tp=t=exp(v*log(x)-ccmath_gaml(v+1.));
      if(modf(v,&y)==0.){ k=(int)y; tp*=v;
        f=2.*log(x)-ccmath_psi(1)-ccmath_psi(k+1);
        t/=2.; if(!(k&1)) t= -t; s=f*t;
        for(p=1,x*=x;;++p){ f-=1./p+1./(v+=1.);
          t*=x/(p*v); s+=(y=t*f);
          if(p>m && fabs(y)<1.e-14) break; }
        if(k>0){ x= -x; s+=(t=1./(tp*2.));
          for(p=1,--k; k>0 ;++p,--k) s+=(t*=x/(p*k)); }
       }
      else{ f=1./(t*v*2.); t*=a0/sin(2.*a0*v); s=f-t;
        for(p=1,x*=x,tp=v;;++p){
          t*=x/(p*(v+=1.)); f*= -x/(p*(tp-=1.));
          s+=(y=f-t); if(p>m && fabs(y)<1.e-14) break; }
       }
     }
    else{ double tq,h,w,z,r;
      t=12./pow(x,.333); k=(int)(t*t); y=2.*(x+k);
      m=(int)v; v-=m; tp=v*v-.25; v+=1.; tq=v*v-.25;
      for(s=h=1.,r=f=z=w=0.; k>0 ;--k,y-=2.){
        t=(y*h-(k+1)*z)/(k-1-tp/k); z=h; f+=(h=t);
        t=(y*s-(k+1)*w)/(k-1-tq/k); w=s; r+=(s=t);  }
      t=sqrt(a0/x)*exp(-x); s*=t/r; h*=t/f; x/=2.; if(m==0) s=h;
      for(k=1; k<m ;++k){ t=v*s/x+h; h=s; s=t; v+=1.;}
     }
   }
  else{ s=t=sqrt(a0/x); x*=2.;
    for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
      t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>=tp) break; s+=t; }
    s*=exp(-x/2.);
   }
  return s;
}

static double ccmath_jbes(double v,double x)
{ double y,s,t,tp; int p,m;
  y=x-8.5; if(y>0.) y*=y; tp=v*v/4.+13.69;
  if(y<tp){ x/=2.; m=(int)x;
    if(x>0.) s=t=exp(v*log(x)-ccmath_gaml(v+1.));
    else{ if(v>0.) return 0.; else if(v==0.) return 1.;}
    for(p=1,x*= -x;;++p){ t*=x/(p*(v+=1.)); s+=t;
      if(p>m && fabs(t)<1.e-13) break;
     }
   }
  else{ double u,a0=1.57079632679490;
    s=t=1./sqrt(x*a0); x*=2.; u=0.;
    for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
      t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>=tp) break;
      if(!(p&1)){ t= -t; s+=t;} else u-=t;
     }
    y=x/2.-(v+.5)*a0; s=cos(y)*s+sin(y)*u;
   }
  return s;
}

static double ccmath_nbes(double v,double x)
{ double y,s,t,tp,u,f,a0=3.14159265358979;
  int p,k,m;
  y=x-8.5; if(y>0.) y*=y; tp=v*v/4.+13.69;
  if(y<tp){ if(x==0.) return HUGE_VAL;
    x/=2.; m=(int)x; u=t=exp(v*log(x)-ccmath_gaml(v+1.));
    if(modf(v,&y)==0.){ k=(int)y; u*=v;
      f=2.*log(x)-ccmath_psi(1)-ccmath_psi(k+1);
      t/=a0; x*= -x; s=f*t;
      for(p=1;;++p){ f-=1./p+1./(v+=1.);
        t*=x/(p*v); s+=(y=t*f); if(p>m && fabs(y)<1.e-13) break; }
      if(k>0){ x= -x; s-=(t=1./(u*a0));
        for(p=1,--k; k>0 ;++p,--k) s-=(t*=x/(p*k)); }
     }
    else{ f=1./(t*v*a0); t/=tan(a0*v); s=t-f;
      for(p=1,x*=x,u=v;;++p){
        t*= -x/(p*(v+=1.)); f*=x/(p*(u-=1.));
        s+=(y=t-f); if(p>m && fabs(y)<1.e-13) break; }
     }
   }
  else{ x*=2.; s=t=2./sqrt(x*a0); u=0.;
    for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
      t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>tp) break;
      if(!(p&1)){ t= -t; s+=t;} else u+=t;
     }
    y=(x-(v+.5)*a0)/2.; s=sin(y)*s+cos(y)*u;
   }
  return s;
}


/* ---------- end of CCMATH code ---------- */


Value ValueCalc::besseli (Value v, Value x)
{
  double vv = converter->asFloat (v).asFloat ();
  double xx = converter->asFloat (x).asFloat ();
  return Value (ccmath_ibes (vv, xx));
}

Value ValueCalc::besselj (Value v, Value x)
{
  double vv = converter->asFloat (v).asFloat ();
  double xx = converter->asFloat (x).asFloat ();
  return Value (ccmath_jbes (vv, xx));
}

Value ValueCalc::besselk (Value v, Value x)
{
  double vv = converter->asFloat (v).asFloat ();
  double xx = converter->asFloat (x).asFloat ();
  return Value (ccmath_kbes (vv, xx));
}

Value ValueCalc::besseln (Value v, Value x)
{
  double vv = converter->asFloat (v).asFloat ();
  double xx = converter->asFloat (x).asFloat ();
  return Value (ccmath_nbes (vv, xx));
}

// ------------------------------------------------------
  
Value ValueCalc::erf (Value x)
{
  return Value (::erf (converter->asFloat (x).asFloat()));
}

Value ValueCalc::erfc (Value x)
{
  return Value (::erfc (converter->asFloat (x).asFloat()));
}

// ------------------------------------------------------

void ValueCalc::arrayWalk (const Value &range,
    Value &res, arrayWalkFunc func, Value param)
{
  if (res.isError()) return;
  if (!range.isArray ())
  {
    func (this, res, range, param);
    return;
  }

  int rows = range.rows ();
  int cols = range.columns ();
  for (int r = 0; r < rows; r++)
    for (int c = 0; c < cols; c++)
    {
      Value v = range.element (c, r);
      if (v.isArray())
        arrayWalk (v, res, func, param);
      else {
        func (this, res, v, param);
        if (res.format() == Value::fmt_None)
          res.setFormat (v.format());
      }
    }
}

void ValueCalc::arrayWalk (TQValueVector<Value> &range,
    Value &res, arrayWalkFunc func, Value param)
{
  if (res.isError()) return;
  for (unsigned int i = 0; i < range.count(); ++i)
    arrayWalk (range[i], res, func, param);
}

void ValueCalc::twoArrayWalk (const Value &a1, const Value &a2,
    Value &res, arrayWalkFunc func)
{
  if (res.isError()) return;
  if (!a1.isArray ())
  {
    func (this, res, a1, a2);
    return;
  }

  int rows = a1.rows ();
  int cols = a1.columns ();
  int rows2 = a2.rows ();
  int cols2 = a2.columns ();
  if ((rows != rows2) || (cols != cols2)) {
    res = Value::errorVALUE();
    return;
  }
  for (int r = 0; r < rows; r++)
    for (int c = 0; c < cols; c++)
    {
      Value v1 = a1.element (c, r);
      Value v2 = a2.element (c, r);
      if (v1.isArray() && v2.isArray())
        twoArrayWalk (v1, v2, res, func);
      else {
        func (this, res, v1, v2);
        if (res.format() == Value::fmt_None)
          res.setFormat (format (v1.format(), v2.format()));
      }
    }
}

void ValueCalc::twoArrayWalk (TQValueVector<Value> &a1,
    TQValueVector<Value> &a2, Value &res, arrayWalkFunc func)
{
  if (res.isError()) return;
  if (a1.count() != a2.count()) {
    res = Value::errorVALUE();
    return;
  }
  for (unsigned int i = 0; i < a1.count(); ++i)
    twoArrayWalk (a1[i], a2[i], res, func);
}

arrayWalkFunc ValueCalc::awFunc (const TQString &name)
{
  if (awFuncs.count(name))
    return awFuncs[name];
  return 0;
}

void ValueCalc::registerAwFunc (const TQString &name, arrayWalkFunc func)
{
  awFuncs[name] = func;
}

// ------------------------------------------------------

Value ValueCalc::sum (const Value &range, bool full)
{
  Value res;
  arrayWalk (range, res, awFunc (full ? "suma" : "sum"), 0);
  return res;
}

Value ValueCalc::sum (TQValueVector<Value> range, bool full)
{
  Value res;
  arrayWalk (range, res, awFunc (full ? "suma" : "sum"), 0);
  return res;
}

// sum of squares
Value ValueCalc::sumsq (const Value &range, bool full)
{
  Value res;
  arrayWalk (range, res, awFunc (full ? "sumsqa" : "sumsq"), 0);
  return res;
}

Value ValueCalc::sumIf (const Value &range,
    const Value &checkRange, const Condition &cond)
{
  if (!range.isArray())
  {
    if (matches (cond, checkRange.element (0, 0)))
      return converter->asNumeric (range);
    return Value (0.0);
  }

  //if we are here, we have an array
  Value res;

  unsigned int rows = range.rows ();
  unsigned int cols = range.columns ();
  for (unsigned int r = 0; r < rows; r++)
    for (unsigned int c = 0; c < cols; c++)
    {
      Value v = range.element (c, r);
      Value newcheck = v;
      if ((c < checkRange.columns()) && (r < checkRange.rows()))
        newcheck = checkRange.element (c, r);
      
      if (v.isArray())
        res = add (res, sumIf (v, newcheck, cond));
      else
        if (matches (cond, newcheck))
          res = add (res, v);
    }

  return res;
}

int ValueCalc::count (const Value &range, bool full)
{
  Value res = 0;
  arrayWalk (range, res, awFunc (full ? "counta" : "count"), 0);
  return converter->asInteger (res).asInteger ();
}

int ValueCalc::count (TQValueVector<Value> range, bool full)
{
  Value res = 0;
  arrayWalk (range, res, awFunc (full ? "counta" : "count"), 0);
  return converter->asInteger (res).asInteger ();
}

int ValueCalc::countIf (const Value &range, const Condition &cond)
{
  if (!range.isArray())
  {
    if (matches (cond, range))
      return range.isEmpty() ? 0 : 1;
    return 0;
  }

  int res = 0;

  int cols = range.columns ();
  int rows = range.rows ();
  for (int r = 0; r < rows; r++)
    for (int c = 0; c < cols; c++)
    {
      Value v = range.element (c, r);
      
      if (v.isArray())
        res += countIf (v, cond);
      else
        if (matches (cond, v))
          res++;
    }

  return res;
}

Value ValueCalc::avg (const Value &range, bool full)
{
  int cnt = count (range, full);
  if (cnt)
    return div (sum (range, full), cnt);
  return Value (0.0);
}

Value ValueCalc::avg (TQValueVector<Value> range, bool full)
{
  int cnt = count (range, full);
  if (cnt)
    return div (sum (range, full), cnt);
  return Value (0.0);
}

Value ValueCalc::max (const Value &range, bool full)
{
  Value res;
  arrayWalk (range, res, awFunc (full ? "maxa" : "max"), 0);
  return res;
}

Value ValueCalc::max (TQValueVector<Value> range, bool full)
{
  Value res;
  arrayWalk (range, res, awFunc (full ? "maxa" : "max"), 0);
  return res;
}

Value ValueCalc::min (const Value &range, bool full)
{
  Value res;
  arrayWalk (range, res, awFunc (full ? "mina" : "min"), 0);
  return res;
}

Value ValueCalc::min (TQValueVector<Value> range, bool full)
{
  Value res;
  arrayWalk (range, res, awFunc (full ? "mina" : "min"), 0);
  return res;
}

Value ValueCalc::product (const Value &range, Value init,
    bool full)
{
  Value res = init;
  if (isZero (init))  // special handling of a zero, due to excel-compat
  {
    if (count (range, full) == 0)
      return init;
    res = 1.0;
  }  
  arrayWalk (range, res, awFunc (full ? "proda" : "prod"), 0);
  return res;
}

Value ValueCalc::product (TQValueVector<Value> range,
    Value init, bool full)
{
  Value res = init;
  if (isZero (init))  // special handling of a zero, due to excel-compat
  {
    if (count (range, full) == 0)
      return init;
    res = 1.0;
  }  
  arrayWalk (range, res, awFunc (full ? "proda" : "prod"), 0);
  return res;
}

Value ValueCalc::stddev (const Value &range, bool full)
{
  return stddev (range, avg (range, full), full);
}

Value ValueCalc::stddev (const Value &range, Value avg,
    bool full)
{
  Value res;
  int cnt = count (range, full);
  arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  return sqrt (div (res, cnt-1));
}

Value ValueCalc::stddev (TQValueVector<Value> range, bool full)
{
  return stddev (range, avg (range, full), full);
}

Value ValueCalc::stddev (TQValueVector<Value> range,
    Value avg, bool full)
{
  Value res;
  int cnt = count (range, full);
  arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  return sqrt (div (res, cnt-1));
}

Value ValueCalc::stddevP (const Value &range, bool full)
{
  return stddevP (range, avg (range, full), full);
}

Value ValueCalc::stddevP (const Value &range, Value avg,
    bool full)
{
  Value res;
  int cnt = count (range, full);
  arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  return sqrt (div (res, cnt));
}

Value ValueCalc::stddevP (TQValueVector<Value> range, bool full)
{
  return stddevP (range, avg (range, full), full);
}

Value ValueCalc::stddevP (TQValueVector<Value> range,
    Value avg, bool full)
{
  Value res;
  int cnt = count (range, full);
  arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  return sqrt (div (res, cnt));
}

Value::Format ValueCalc::format (Value::Format a,
    Value::Format b)
{
  if ((a == Value::fmt_None) || (a == Value::fmt_Boolean))
    return b;
  return a;
}


// ------------------------------------------------------

void ValueCalc::getCond (Condition &cond, Value val)
{
  // not a string - we simply take it as a numeric value
  // that also handles floats, logical values, date/time and such
  if (!val.isString()) {
    cond.comp = isEqual;
    cond.type = numeric;
    cond.value = converter->asFloat (val).asFloat();
    return;
  }
  TQString text = converter->asString (val).asString();
  cond.comp = isEqual;
  text = text.stripWhiteSpace();

  if ( text.startsWith( "<=" ) )
  {
    cond.comp = lessEqual;
    text = text.remove( 0, 2 );
  }
  else if ( text.startsWith( ">=" ) )
  {
    cond.comp = greaterEqual;
    text = text.remove( 0, 2 );
  }
  else if ( text.startsWith( "!=" ) || text.startsWith( "<>" ) )
  {
    cond.comp = notEqual;
    text = text.remove( 0, 2 );
  }
  else if ( text.startsWith( "==" ) )
  {
    cond.comp = isEqual;
    text = text.remove( 0, 2 );
  }
  else if ( text.startsWith( "<" ) )
  {
    cond.comp = isLess;
    text = text.remove( 0, 1 );
  }
  else if ( text.startsWith( ">" ) )
  {
    cond.comp = isGreater;
    text = text.remove( 0, 1 );
  }
  else if ( text.startsWith( "=" ) )
  {
    cond.comp = isEqual;
    text = text.remove( 0, 1 );
  }

  text = text.stripWhiteSpace();

  bool ok = false;
  double d = text.toDouble( &ok );
  if ( ok )
  {
    cond.type = numeric;
    cond.value = d;
  }
  else
  {
    cond.type = string;
    cond.stringValue = text;
  }
}

bool ValueCalc::matches (const Condition &cond, Value val)
{
  if (val.isEmpty())
	return false;
  if (cond.type == numeric) {
    double d = converter->asFloat (val).asFloat();
    switch ( cond.comp )
    {
      case isEqual:
      if (approxEqual (d, cond.value)) return true;
      break;
  
      case isLess:
      if (d < cond.value) return true;
      break;
  
      case isGreater:
      if (d > cond.value) return true;
      break;
  
      case lessEqual:
      if (d <= cond.value) return true;
      break;
  
      case greaterEqual:
      if (d >= cond.value) return true;
      break;
  
      case notEqual:
      if (d != cond.value) return true;
      break;
    }
  } else {
    TQString d = converter->asString (val).asString();
    switch ( cond.comp )
    {
      case isEqual:
      if (d == cond.stringValue) return true;
      break;
  
      case isLess:
      if (d < cond.stringValue) return true;
      break;
  
      case isGreater:
      if (d > cond.stringValue) return true;
      break;
  
      case lessEqual:
      if (d <= cond.stringValue) return true;
      break;
  
      case greaterEqual:
      if (d >= cond.stringValue) return true;
      break;
  
      case notEqual:
      if (d != cond.stringValue) return true;
      break;
    }
  }
  return false;
}