You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
2054 lines
48 KiB
2054 lines
48 KiB
/* 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, bool CS)
|
|
{
|
|
TQString aa = converter->asString(a).asString();
|
|
TQString bb = converter->asString(b).asString();
|
|
if (!CS) {
|
|
aa = aa.lower();
|
|
bb = bb.lower();
|
|
}
|
|
return (aa == bb);
|
|
}
|
|
|
|
bool ValueCalc::strGreater (const Value &a, const Value &b, bool CS)
|
|
{
|
|
TQString aa = converter->asString(a).asString();
|
|
TQString bb = converter->asString(b).asString();
|
|
if (!CS) {
|
|
aa = aa.lower();
|
|
bb = bb.lower();
|
|
}
|
|
return (aa > bb);
|
|
}
|
|
|
|
bool ValueCalc::strGequal (const Value &a, const Value &b, bool CS)
|
|
{
|
|
TQString aa = converter->asString(a).asString();
|
|
TQString bb = converter->asString(b).asString();
|
|
if (!CS) {
|
|
aa = aa.lower();
|
|
bb = bb.lower();
|
|
}
|
|
return (aa >= bb);
|
|
}
|
|
|
|
bool ValueCalc::strLower (const Value &a, const Value &b, bool CS)
|
|
{
|
|
return strGreater (b, a, CS);
|
|
}
|
|
|
|
bool ValueCalc::naturalEqual (const Value &a, const Value &b, bool CS)
|
|
{
|
|
if (a.isNumber() && b.isNumber())
|
|
return approxEqual(a, b);
|
|
if (a.allowComparison(b))
|
|
return a.equal(b);
|
|
return strEqual(a, b, CS);
|
|
}
|
|
|
|
bool ValueCalc::naturalGreater (const Value &a, const Value &b, bool CS)
|
|
{
|
|
if (a.isNumber() && b.isNumber())
|
|
return greater(a, b);
|
|
if (a.allowComparison(b))
|
|
return a.greater(b);
|
|
return strGreater(a, b, CS);
|
|
}
|
|
|
|
bool ValueCalc::naturalGequal (const Value &a, const Value &b, bool CS)
|
|
{
|
|
return naturalGreater(a, b, CS) || naturalEqual(a, b, CS);
|
|
}
|
|
|
|
bool ValueCalc::naturalLower (const Value &a, const Value &b, bool CS)
|
|
{
|
|
return naturalGreater(b, a, CS);
|
|
}
|
|
|
|
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;
|
|
}
|
|
|