KOffice – TDE office suite
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.

1993 lines
47KB

  1. /* This file is part of the KDE project
  2. Copyright (C) 2005 Tomas Mecir <mecirt@gmail.com>
  3. This library is free software; you can redistribute it and/or
  4. modify it under the terms of the GNU Library General Public
  5. License as published by the Free Software Foundation; either
  6. version 2 of the License, or (at your option) any later version.
  7. This library is distributed in the hope that it will be useful,
  8. but WITHOUT ANY WARRANTY; without even the implied warranty of
  9. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  10. Library General Public License for more details.
  11. You should have received a copy of the GNU Library General Public License
  12. along with this library; see the file COPYING.LIB. If not, write to
  13. the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  14. * Boston, MA 02110-1301, USA.
  15. */
  16. #include "valuecalc.h"
  17. #include "valueconverter.h"
  18. #include <kdebug.h>
  19. #include <errno.h>
  20. #include <float.h>
  21. #include <math.h>
  22. #include <stdlib.h>
  23. #include <time.h>
  24. using namespace KSpread;
  25. // Array-walk functions registered on ValueCalc object
  26. void awSum (ValueCalc *c, Value &res, Value val, Value)
  27. {
  28. if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
  29. res = c->add (res, val);
  30. }
  31. void awSumA (ValueCalc *c, Value &res, Value val, Value)
  32. {
  33. if (!val.isEmpty())
  34. res = c->add (res, val);
  35. }
  36. void awSumSq (ValueCalc *c, Value &res, Value val, Value)
  37. {
  38. if (!val.isEmpty())
  39. res = c->add (res, c->sqr (val));
  40. }
  41. void awSumSqA (ValueCalc *c, Value &res, Value val, Value)
  42. {
  43. if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
  44. res = c->add (res, c->sqr (val));
  45. }
  46. void awCount (ValueCalc *c, Value &res, Value val, Value)
  47. {
  48. if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
  49. res = c->add (res, 1);
  50. }
  51. void awCountA (ValueCalc *c, Value &res, Value val, Value)
  52. {
  53. if (!val.isEmpty())
  54. res = c->add (res, 1);
  55. }
  56. void awMax (ValueCalc *c, Value &res, Value val, Value)
  57. {
  58. if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
  59. if (res.isEmpty())
  60. res = val;
  61. else
  62. if (c->greater (val, res)) res = val;
  63. }
  64. void awMaxA (ValueCalc *c, Value &res, Value val, Value)
  65. {
  66. if (!val.isEmpty())
  67. if (res.isEmpty())
  68. // convert to number, so that we don't return string/bool
  69. res = c->conv()->asNumeric (val);
  70. else
  71. if (c->greater (val, res))
  72. // convert to number, so that we don't return string/bool
  73. res = c->conv()->asNumeric (val);
  74. }
  75. void awMin (ValueCalc *c, Value &res, Value val, Value)
  76. {
  77. if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
  78. if (res.isEmpty())
  79. res = val;
  80. else
  81. if (c->lower (val, res)) res = val;
  82. }
  83. void awMinA (ValueCalc *c, Value &res, Value val, Value)
  84. {
  85. if (!val.isEmpty())
  86. if (res.isEmpty())
  87. // convert to number, so that we don't return string/bool
  88. res = c->conv()->asNumeric (val);
  89. else
  90. if (c->lower (val, res))
  91. // convert to number, so that we don't return string/bool
  92. res = c->conv()->asNumeric (val);
  93. }
  94. void awProd (ValueCalc *c, Value &res, Value val, Value)
  95. {
  96. if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
  97. res = c->mul (res, val);
  98. }
  99. void awProdA (ValueCalc *c, Value &res, Value val, Value)
  100. {
  101. if (!val.isEmpty())
  102. res = c->mul (res, val);
  103. }
  104. // sum of squares of deviations, used to compute standard deviation
  105. void awDevSq (ValueCalc *c, Value &res, Value val,
  106. Value avg)
  107. {
  108. if (!val.isEmpty())
  109. res = c->add (res, c->sqr (c->sub (val, avg)));
  110. }
  111. // sum of squares of deviations, used to compute standard deviation
  112. void awDevSqA (ValueCalc *c, Value &res, Value val,
  113. Value avg)
  114. {
  115. if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()))
  116. res = c->add (res, c->sqr (c->sub (val, avg)));
  117. }
  118. bool isDate (Value val) {
  119. Value::Format fmt = val.format();
  120. if ((fmt == Value::fmt_Date) || (fmt == Value::fmt_DateTime))
  121. return true;
  122. return false;
  123. }
  124. // ***********************
  125. // ****** ValueCalc ******
  126. // ***********************
  127. ValueCalc::ValueCalc (ValueConverter* c): converter( c )
  128. {
  129. // initialize the random number generator
  130. srand (time (0));
  131. // register array-walk functions
  132. registerAwFunc ("sum", awSum);
  133. registerAwFunc ("suma", awSumA);
  134. registerAwFunc ("sumsq", awSumSq);
  135. registerAwFunc ("sumsqa", awSumSqA);
  136. registerAwFunc ("count", awCount);
  137. registerAwFunc ("counta", awCountA);
  138. registerAwFunc ("max", awMax);
  139. registerAwFunc ("maxa", awMaxA);
  140. registerAwFunc ("min", awMin);
  141. registerAwFunc ("mina", awMinA);
  142. registerAwFunc ("prod", awProd);
  143. registerAwFunc ("proda", awProdA);
  144. registerAwFunc ("devsq", awDevSq);
  145. registerAwFunc ("devsqa", awDevSq);
  146. }
  147. Value ValueCalc::add (const Value &a, const Value &b)
  148. {
  149. if (a.isError()) return a;
  150. if (b.isError()) return b;
  151. Value res;
  152. if (a.isInteger() && b.isEmpty() || a.isEmpty() && b.isInteger()
  153. || a.isInteger() && b.isInteger())
  154. {
  155. int aa, bb;
  156. aa = converter->asInteger (a).asInteger();
  157. bb = converter->asInteger (b).asInteger();
  158. res = Value (aa + bb);
  159. }
  160. else
  161. {
  162. double aa, bb;
  163. aa = converter->asFloat (a).asFloat();
  164. bb = converter->asFloat (b).asFloat();
  165. res = Value (aa + bb);
  166. }
  167. if (a.isNumber() || a.isEmpty())
  168. res.setFormat (format (a.format(), b.format()));
  169. // operation on two dates should produce a number
  170. if (isDate(a) && isDate(b))
  171. res.setFormat (Value::fmt_Number);
  172. return res;
  173. }
  174. Value ValueCalc::sub (const Value &a, const Value &b)
  175. {
  176. if (a.isError()) return a;
  177. if (b.isError()) return b;
  178. double aa, bb;
  179. aa = converter->asFloat (a).asFloat();
  180. bb = converter->asFloat (b).asFloat();
  181. Value res = Value (aa - bb);
  182. if (a.isNumber() || a.isEmpty())
  183. res.setFormat (format (a.format(), b.format()));
  184. // operation on two dates should produce a number
  185. if (isDate(a) && isDate(b))
  186. res.setFormat (Value::fmt_Number);
  187. return res;
  188. }
  189. Value ValueCalc::mul (const Value &a, const Value &b)
  190. {
  191. if (a.isError()) return a;
  192. if (b.isError()) return b;
  193. double aa, bb;
  194. aa = converter->asFloat (a).asFloat();
  195. bb = converter->asFloat (b).asFloat();
  196. Value res = Value (aa * bb);
  197. if (a.isNumber() || a.isEmpty())
  198. res.setFormat (format (a.format(), b.format()));
  199. // operation on two dates should produce a number
  200. if (isDate(a) && isDate(b))
  201. res.setFormat (Value::fmt_Number);
  202. return res;
  203. }
  204. Value ValueCalc::div (const Value &a, const Value &b)
  205. {
  206. if (a.isError()) return a;
  207. if (b.isError()) return b;
  208. double aa, bb;
  209. aa = converter->asFloat (a).asFloat();
  210. bb = converter->asFloat (b).asFloat();
  211. Value res;
  212. if (bb == 0.0)
  213. return Value::errorDIV0();
  214. else
  215. res = Value (aa / bb);
  216. if (a.isNumber() || a.isEmpty())
  217. res.setFormat (format (a.format(), b.format()));
  218. // operation on two dates should produce a number
  219. if (isDate(a) && isDate(b))
  220. res.setFormat (Value::fmt_Number);
  221. return res;
  222. }
  223. Value ValueCalc::mod (const Value &a, const Value &b)
  224. {
  225. if (a.isError()) return a;
  226. if (b.isError()) return b;
  227. double aa, bb;
  228. aa = converter->asFloat (a).asFloat();
  229. bb = converter->asFloat (b).asFloat();
  230. Value res;
  231. if (bb == 0.0)
  232. return Value::errorDIV0();
  233. else
  234. {
  235. double m = fmod (aa, bb);
  236. // the following adjustments are needed by OpenFormula:
  237. // can't simply use fixed increases/decreases, because the implementation
  238. // of fmod may differ on various platforms, and we should always return
  239. // the same results ...
  240. if ((bb > 0) && (aa < 0)) // result must be positive here
  241. while (m < 0) m += bb;
  242. if (bb < 0) // result must be negative here, but not lower than bb
  243. {
  244. // bb is negative, hence the following two are correct
  245. while (m < bb) m -= bb; // same as m+=fabs(bb)
  246. while (m > 0) m += bb; // same as m-=fabs(bb)
  247. }
  248. res = Value (m);
  249. }
  250. if (a.isNumber() || a.isEmpty())
  251. res.setFormat (format (a.format(), b.format()));
  252. if (isDate(a) && isDate(b))
  253. res.setFormat (Value::fmt_Number);
  254. return res;
  255. }
  256. Value ValueCalc::pow (const Value &a, const Value &b)
  257. {
  258. if (a.isError()) return a;
  259. if (b.isError()) return b;
  260. double aa, bb;
  261. aa = converter->asFloat (a).asFloat();
  262. bb = converter->asFloat (b).asFloat();
  263. Value res = Value (::pow (aa, bb));
  264. if (a.isNumber() || a.isEmpty())
  265. res.setFormat (format (a.format(), b.format()));
  266. // operation on date(s) should produce a number
  267. if (isDate(a) || isDate(b))
  268. res.setFormat (Value::fmt_Number);
  269. return res;
  270. }
  271. Value ValueCalc::sqr (const Value &a)
  272. {
  273. if (a.isError()) return a;
  274. return mul (a, a);
  275. }
  276. Value ValueCalc::sqrt (const Value &a)
  277. {
  278. if (a.isError()) return a;
  279. Value res = Value (::sqrt (converter->asFloat(a).asFloat()));
  280. if (a.isNumber() || a.isEmpty())
  281. res.setFormat (a.format());
  282. // operation on date(s) should produce a number
  283. if (isDate(a))
  284. res.setFormat (Value::fmt_Number);
  285. return res;
  286. }
  287. Value ValueCalc::add (const Value &a, double b)
  288. {
  289. if (a.isError()) return a;
  290. Value res = Value (converter->asFloat(a).asFloat() + b);
  291. if (a.isNumber() || a.isEmpty())
  292. res.setFormat (a.format());
  293. return res;
  294. }
  295. Value ValueCalc::sub (const Value &a, double b)
  296. {
  297. if (a.isError()) return a;
  298. Value res = Value (converter->asFloat(a).asFloat() - b);
  299. if (a.isNumber() || a.isEmpty())
  300. res.setFormat (a.format());
  301. return res;
  302. }
  303. Value ValueCalc::mul (const Value &a, double b)
  304. {
  305. if (a.isError()) return a;
  306. Value res = Value (converter->asFloat(a).asFloat() * b);
  307. if (a.isNumber() || a.isEmpty())
  308. res.setFormat (a.format());
  309. return res;
  310. }
  311. Value ValueCalc::div (const Value &a, double b)
  312. {
  313. if (a.isError()) return a;
  314. Value res;
  315. if (b == 0.0)
  316. return Value::errorDIV0();
  317. res = Value (converter->asFloat(a).asFloat() / b);
  318. if (a.isNumber() || a.isEmpty())
  319. res.setFormat (a.format());
  320. return res;
  321. }
  322. Value ValueCalc::pow (const Value &a, double b)
  323. {
  324. if (a.isError()) return a;
  325. Value res = Value (::pow (converter->asFloat(a).asFloat(), b));
  326. if (a.isNumber() || a.isEmpty())
  327. res.setFormat (a.format());
  328. return res;
  329. }
  330. Value ValueCalc::abs (const Value &a)
  331. {
  332. if (a.isError()) return a;
  333. return Value (fabs (converter->asFloat (a).asFloat()));
  334. }
  335. bool ValueCalc::isZero (const Value &a)
  336. {
  337. if (a.isError()) return false;
  338. return (converter->asFloat (a).asFloat() == 0.0);
  339. }
  340. bool ValueCalc::isEven (const Value &a)
  341. {
  342. if (a.isError()) return false;
  343. return ((converter->asInteger (a).asInteger() % 2) == 0);
  344. }
  345. bool ValueCalc::equal (const Value &a, const Value &b)
  346. {
  347. return (converter->asFloat (a).asFloat() == converter->asFloat (b).asFloat());
  348. }
  349. /*********************************************************************
  350. *
  351. * Helper function to avoid problems with rounding floating point
  352. * values. Idea for this kind of solution taken from Openoffice.
  353. *
  354. *********************************************************************/
  355. bool ValueCalc::approxEqual (const Value &a, const Value &b)
  356. {
  357. double aa = converter->asFloat (a).asFloat();
  358. double bb = converter->asFloat (b).asFloat();
  359. if (aa == bb)
  360. return true;
  361. double x = aa - bb;
  362. return (x < 0.0 ? -x : x) < ((aa < 0.0 ? -aa : aa) * DBL_EPSILON);
  363. }
  364. bool ValueCalc::strEqual (const Value &a, const Value &b)
  365. {
  366. return (converter->asString (a).asString() == converter->asString (b).asString());
  367. }
  368. bool ValueCalc::greater (const Value &a, const Value &b)
  369. {
  370. double aa = converter->asFloat (a).asFloat();
  371. double bb = converter->asFloat (b).asFloat();
  372. return (aa > bb);
  373. }
  374. bool ValueCalc::gequal (const Value &a, const Value &b)
  375. {
  376. return (greater (a,b) || approxEqual (a,b));
  377. }
  378. bool ValueCalc::lower (const Value &a, const Value &b)
  379. {
  380. return greater (b, a);
  381. }
  382. Value ValueCalc::roundDown (const Value &a,
  383. const Value &digits) {
  384. return roundDown (a, converter->asInteger (digits).asInteger());
  385. }
  386. Value ValueCalc::roundUp (const Value &a,
  387. const Value &digits) {
  388. return roundUp (a, converter->asInteger (digits).asInteger());
  389. }
  390. Value ValueCalc::round (const Value &a,
  391. const Value &digits) {
  392. return round (a, converter->asInteger (digits).asInteger());
  393. }
  394. Value ValueCalc::roundDown (const Value &a, int digits)
  395. {
  396. // shift in one direction, round, shift back
  397. Value val = a;
  398. if (digits > 0)
  399. for (int i = 0; i < digits; ++i)
  400. val = mul (val, 10);
  401. if (digits < 0)
  402. for (int i = 0; i < digits; ++i)
  403. val = div (val, 10);
  404. val = Value (floor (converter->asFloat (val).asFloat()));
  405. if (digits > 0)
  406. for (int i = 0; i < digits; ++i)
  407. val = div (val, 10);
  408. if (digits < 0)
  409. for (int i = 0; i < digits; ++i)
  410. val = mul (val, 10);
  411. return val;
  412. }
  413. Value ValueCalc::roundUp (const Value &a, int digits)
  414. {
  415. // shift in one direction, round, shift back
  416. Value val = a;
  417. if (digits > 0)
  418. for (int i = 0; i < digits; ++i)
  419. val = mul (val, 10);
  420. if (digits < 0)
  421. for (int i = 0; i < digits; ++i)
  422. val = div (val, 10);
  423. val = Value (ceil (converter->asFloat (val).asFloat()));
  424. if (digits > 0)
  425. for (int i = 0; i < digits; ++i)
  426. val = div (val, 10);
  427. if (digits < 0)
  428. for (int i = 0; i < digits; ++i)
  429. val = mul (val, 10);
  430. return val;
  431. }
  432. Value ValueCalc::round (const Value &a, int digits)
  433. {
  434. // shift in one direction, round, shift back
  435. Value val = a;
  436. if (digits > 0)
  437. for (int i = 0; i < digits; ++i)
  438. val = mul (val, 10);
  439. if (digits < 0)
  440. for (int i = 0; i < digits; ++i)
  441. val = div (val, 10);
  442. val = Value (int(converter->asFloat (val).asFloat()+0.5));
  443. if (digits > 0)
  444. for (int i = 0; i < digits; ++i)
  445. val = div (val, 10);
  446. if (digits < 0)
  447. for (int i = 0; i < digits; ++i)
  448. val = mul (val, 10);
  449. return val;
  450. }
  451. int ValueCalc::sign (const Value &a)
  452. {
  453. double val = converter->asFloat (a).asFloat ();
  454. if (val == 0) return 0;
  455. if (val > 0) return 1;
  456. return -1;
  457. }
  458. Value ValueCalc::log (const Value &number,
  459. const Value &base)
  460. {
  461. double logbase = converter->asFloat (base).asFloat();
  462. if (logbase == 1.0)
  463. return Value::errorDIV0();
  464. if (logbase <= 0.0)
  465. return Value::errorNA();
  466. logbase = log10 (logbase);
  467. Value res = Value (log10 (converter->asFloat (number).asFloat()) / logbase);
  468. if (number.isNumber() || number.isEmpty())
  469. res.setFormat (number.format());
  470. return res;
  471. }
  472. Value ValueCalc::ln (const Value &number)
  473. {
  474. Value res = Value (::log (converter->asFloat (number).asFloat()));
  475. if (number.isNumber() || number.isEmpty())
  476. res.setFormat (number.format());
  477. return res;
  478. }
  479. Value ValueCalc::log (const Value &number, double base)
  480. {
  481. if (base <= 0.0)
  482. return Value::errorNA();
  483. if (base == 1.0)
  484. return Value::errorDIV0();
  485. double num = converter->asFloat (number).asFloat();
  486. Value res = Value (log10 (num) / log10 (base));
  487. if (number.isNumber() || number.isEmpty())
  488. res.setFormat (number.format());
  489. return res;
  490. }
  491. Value ValueCalc::exp (const Value &number)
  492. {
  493. return Value (::exp (converter->asFloat (number).asFloat()));
  494. }
  495. Value ValueCalc::pi ()
  496. {
  497. // retun PI in double-precision
  498. // if arbitrary precision gets in, this should be extended to return
  499. // more if need be
  500. return Value (M_PI);
  501. }
  502. Value ValueCalc::eps ()
  503. {
  504. // #### This should adjust according to the actual number system used
  505. // (float, double, long double, ...)
  506. return Value (DBL_EPSILON);
  507. }
  508. Value ValueCalc::random (double range)
  509. {
  510. return Value (range * (double) rand() / (RAND_MAX + 1.0));
  511. }
  512. Value ValueCalc::random (Value range)
  513. {
  514. return random (converter->asFloat (range).asFloat());
  515. }
  516. Value ValueCalc::fact (const Value &which)
  517. {
  518. // we can simply use integers - no one is going to compute factorial of
  519. // anything bigger than 2^32
  520. return fact (converter->asInteger (which).asInteger());
  521. }
  522. Value ValueCalc::fact (const Value &which,
  523. const Value &end)
  524. {
  525. // we can simply use integers - no one is going to compute factorial of
  526. // anything bigger than 2^32
  527. return fact (converter->asInteger (which).asInteger(),
  528. converter->asInteger (end).asInteger ());
  529. }
  530. Value ValueCalc::fact (int which, int end) {
  531. if (which < 0)
  532. return Value (-1);
  533. if (which == 0)
  534. return Value (1);
  535. // no multiplication if val==end
  536. if (which == end)
  537. return Value (1);
  538. return (mul (fact (which-1, end), which));
  539. }
  540. Value ValueCalc::factDouble (int which)
  541. {
  542. if (which < 0)
  543. return Value (-1);
  544. if ((which == 0) || (which == 1))
  545. return Value (1);
  546. return (mul (factDouble (which-2), which));
  547. }
  548. Value ValueCalc::factDouble (Value which)
  549. {
  550. return factDouble (converter->asInteger (which).asInteger());
  551. }
  552. Value ValueCalc::combin (int n, int k)
  553. {
  554. if (n >= 15)
  555. {
  556. double result = ::exp(lgamma (n + 1) - lgamma (k + 1) - lgamma (n-k+1));
  557. return Value (floor(result + 0.5));
  558. }
  559. else
  560. return div (div (fact (n), fact (k)), fact (n - k));
  561. }
  562. Value ValueCalc::combin (Value n, Value k)
  563. {
  564. int nn = converter->asInteger (n).asInteger();
  565. int kk = converter->asInteger (k).asInteger();
  566. return combin (nn, kk);
  567. }
  568. Value ValueCalc::gcd (const Value &a, const Value &b)
  569. {
  570. // Euler's GCD algorithm
  571. Value aa = round (a);
  572. Value bb = round (b);
  573. if (approxEqual (aa, bb)) return aa;
  574. if (aa.isZero()) return bb;
  575. if (bb.isZero()) return aa;
  576. if (greater (aa, bb))
  577. return gcd (bb, mod (aa, bb));
  578. else
  579. return gcd (aa, mod (bb, aa));
  580. }
  581. Value ValueCalc::lcm (const Value &a, const Value &b)
  582. {
  583. Value aa = round (a);
  584. Value bb = round (b);
  585. if (approxEqual (aa, bb)) return aa;
  586. if (aa.isZero()) return bb;
  587. if (bb.isZero()) return aa;
  588. Value g = gcd (aa, bb);
  589. if (g.isZero()) // GCD is zero for some weird reason
  590. return mul (aa, bb);
  591. return div (mul (aa, bb), g);
  592. }
  593. Value ValueCalc::base (const Value &val, int base, int prec)
  594. {
  595. if (base == 10) return round (val, prec);
  596. if (prec < 0) prec = 2;
  597. if ((base < 2) || (base > 36))
  598. return Value::errorVALUE();
  599. double value = converter->asFloat (val).asFloat();
  600. TQString result = TQString::number ((int)value, base);
  601. if (prec > 0)
  602. {
  603. result += "."; value = value - (int)value;
  604. int ix;
  605. for( int i = 0; i < prec; i++ )
  606. {
  607. ix = (int) value * base;
  608. result += "0123456789abcdefghijklmnopqrstuvwxyz"[ix];
  609. value = base * (value - (double)ix/base);
  610. }
  611. }
  612. return Value (result.upper());
  613. }
  614. Value ValueCalc::fromBase (const Value &val, int base)
  615. {
  616. TQString str = converter->asString (val).asString();
  617. bool ok;
  618. double num = str.toLong (&ok, base);
  619. if (ok)
  620. return Value (num);
  621. return Value::errorVALUE();
  622. }
  623. Value ValueCalc::sin (const Value &number)
  624. {
  625. Value res = Value (::sin (converter->asFloat (number).asFloat()));
  626. if (number.isNumber() || number.isEmpty())
  627. res.setFormat (number.format());
  628. return res;
  629. }
  630. Value ValueCalc::cos (const Value &number)
  631. {
  632. Value res = Value (::cos (converter->asFloat (number).asFloat()));
  633. if (number.isNumber() || number.isEmpty())
  634. res.setFormat (number.format());
  635. return res;
  636. }
  637. Value ValueCalc::tg (const Value &number)
  638. {
  639. Value res = Value (::tan (converter->asFloat (number).asFloat()));
  640. if (number.isNumber() || number.isEmpty())
  641. res.setFormat (number.format());
  642. return res;
  643. }
  644. Value ValueCalc::cotg (const Value &number)
  645. {
  646. Value res = Value (div (1, ::tan (converter->asFloat (number).asFloat())));
  647. if (number.isNumber() || number.isEmpty())
  648. res.setFormat (number.format());
  649. return res;
  650. }
  651. Value ValueCalc::asin (const Value &number)
  652. {
  653. errno = 0;
  654. Value res = Value (::asin (converter->asFloat (number).asFloat()));
  655. if (errno)
  656. return Value::errorVALUE();
  657. if (number.isNumber() || number.isEmpty())
  658. res.setFormat (number.format());
  659. return res;
  660. }
  661. Value ValueCalc::acos (const Value &number)
  662. {
  663. errno = 0;
  664. Value res = Value (::acos (converter->asFloat (number).asFloat()));
  665. if (errno)
  666. return Value::errorVALUE();
  667. if (number.isNumber() || number.isEmpty())
  668. res.setFormat (number.format());
  669. return res;
  670. }
  671. Value ValueCalc::atg (const Value &number)
  672. {
  673. errno = 0;
  674. Value res = Value (::atan (converter->asFloat (number).asFloat()));
  675. if (errno)
  676. return Value::errorVALUE();
  677. if (number.isNumber() || number.isEmpty())
  678. res.setFormat (number.format());
  679. return res;
  680. }
  681. Value ValueCalc::atan2 (const Value &y, const Value &x)
  682. {
  683. double yy = converter->asFloat (y).asFloat();
  684. double xx = converter->asFloat (x).asFloat();
  685. return Value (::atan2 (yy, xx));
  686. }
  687. Value ValueCalc::sinh (const Value &number)
  688. {
  689. Value res = Value (::sinh (converter->asFloat (number).asFloat()));
  690. if (number.isNumber() || number.isEmpty())
  691. res.setFormat (number.format());
  692. return res;
  693. }
  694. Value ValueCalc::cosh (const Value &number)
  695. {
  696. Value res = Value (::cosh (converter->asFloat (number).asFloat()));
  697. if (number.isNumber() || number.isEmpty())
  698. res.setFormat (number.format());
  699. return res;
  700. }
  701. Value ValueCalc::tgh (const Value &number)
  702. {
  703. Value res = Value (::tanh (converter->asFloat (number).asFloat()));
  704. if (number.isNumber() || number.isEmpty())
  705. res.setFormat (number.format());
  706. return res;
  707. }
  708. Value ValueCalc::asinh (const Value &number)
  709. {
  710. errno = 0;
  711. Value res = Value (::asinh (converter->asFloat (number).asFloat()));
  712. if (errno)
  713. return Value::errorVALUE();
  714. if (number.isNumber() || number.isEmpty())
  715. res.setFormat (number.format());
  716. return res;
  717. }
  718. Value ValueCalc::acosh (const Value &number)
  719. {
  720. errno = 0;
  721. Value res = Value (::acosh (converter->asFloat (number).asFloat()));
  722. if (errno)
  723. return Value::errorVALUE();
  724. if (number.isNumber() || number.isEmpty())
  725. res.setFormat (number.format());
  726. return res;
  727. }
  728. Value ValueCalc::atgh (const Value &number)
  729. {
  730. errno = 0;
  731. Value res = Value (::atanh (converter->asFloat (number).asFloat()));
  732. if (errno)
  733. return Value::errorVALUE();
  734. if (number.isNumber() || number.isEmpty())
  735. res.setFormat (number.format());
  736. return res;
  737. }
  738. Value ValueCalc::phi (Value x)
  739. {
  740. Value constant (0.39894228040143268);
  741. // constant * exp(-(x * x) / 2.0);
  742. Value x2neg = mul (sqr (x), -1);
  743. return mul (constant, exp (div (x2neg, 2.0)));
  744. }
  745. static double taylor_helper (double* pPolynom, uint nMax, double x)
  746. {
  747. double nVal = pPolynom[nMax];
  748. for (int i = nMax-1; i >= 0; i--) {
  749. nVal = pPolynom[i] + (nVal * x);
  750. }
  751. return nVal;
  752. }
  753. Value ValueCalc::gauss (Value xx)
  754. // this is a weird function
  755. {
  756. double x = converter->asFloat (xx).asFloat();
  757. double t0[] =
  758. { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,
  759. -0.00118732821548045, 0.00011543468761616, -0.00000944465625950,
  760. 0.00000066596935163, -0.00000004122667415, 0.00000000227352982,
  761. 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };
  762. double t2[] =
  763. { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,
  764. 0.02699548325659403, -0.00449924720943234, -0.00224962360471617,
  765. 0.00134977416282970, -0.00011783742691370, -0.00011515930357476,
  766. 0.00003704737285544, 0.00000282690796889, -0.00000354513195524,
  767. 0.00000037669563126, 0.00000019202407921, -0.00000005226908590,
  768. -0.00000000491799345, 0.00000000366377919, -0.00000000015981997,
  769. -0.00000000017381238, 0.00000000002624031, 0.00000000000560919,
  770. -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };
  771. double t4[] =
  772. { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,
  773. 0.00033457556441221, -0.00028996548915725, 0.00018178605666397,
  774. -0.00008252863922168, 0.00002551802519049, -0.00000391665839292,
  775. -0.00000074018205222, 0.00000064422023359, -0.00000017370155340,
  776. 0.00000000909595465, 0.00000000944943118, -0.00000000329957075,
  777. 0.00000000029492075, 0.00000000011874477, -0.00000000004420396,
  778. 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };
  779. double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };
  780. double xAbs = fabs(x);
  781. uint xShort = static_cast<uint>(floor(xAbs));
  782. double nVal = 0.0;
  783. if (xShort == 0)
  784. nVal = taylor_helper(t0, 11, (xAbs * xAbs)) * xAbs;
  785. else if ((xShort >= 1) && (xShort <= 2))
  786. nVal = taylor_helper(t2, 23, (xAbs - 2.0));
  787. else if ((xShort >= 3) && (xShort <= 4))
  788. nVal = taylor_helper(t4, 20, (xAbs - 4.0));
  789. else
  790. {
  791. double phiAbs = converter->asFloat (phi (xAbs)).asFloat();
  792. nVal = 0.5 + phiAbs * taylor_helper(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs;
  793. }
  794. if (x < 0.0)
  795. return Value (-nVal);
  796. else
  797. return Value (nVal);
  798. }
  799. Value ValueCalc::gaussinv (Value xx)
  800. // this is a weird function
  801. {
  802. double x = converter->asFloat (xx).asFloat();
  803. double q,t,z;
  804. q=x-0.5;
  805. if(fabs(q)<=.425)
  806. {
  807. t=0.180625-q*q;
  808. z=
  809. q*
  810. (
  811. (
  812. (
  813. (
  814. (
  815. (
  816. (
  817. t*2509.0809287301226727+33430.575583588128105
  818. )
  819. *t+67265.770927008700853
  820. )
  821. *t+45921.953931549871457
  822. )
  823. *t+13731.693765509461125
  824. )
  825. *t+1971.5909503065514427
  826. )
  827. *t+133.14166789178437745
  828. )
  829. *t+3.387132872796366608
  830. )
  831. /
  832. (
  833. (
  834. (
  835. (
  836. (
  837. (
  838. (
  839. t*5226.495278852854561+28729.085735721942674
  840. )
  841. *t+39307.89580009271061
  842. )
  843. *t+21213.794301586595867
  844. )
  845. *t+5394.1960214247511077
  846. )
  847. *t+687.1870074920579083
  848. )
  849. *t+42.313330701600911252
  850. )
  851. *t+1.0
  852. );
  853. }
  854. else
  855. {
  856. if(q>0) t=1-x;
  857. else t=x;
  858. t=::sqrt(-::log(t));
  859. if(t<=5.0)
  860. {
  861. t+=-1.6;
  862. z=
  863. (
  864. (
  865. (
  866. (
  867. (
  868. (
  869. (
  870. t*7.7454501427834140764e-4+0.0227238449892691845833
  871. )
  872. *t+0.24178072517745061177
  873. )
  874. *t+1.27045825245236838258
  875. )
  876. *t+3.64784832476320460504
  877. )
  878. *t+5.7694972214606914055
  879. )
  880. *t+4.6303378461565452959
  881. )
  882. *t+1.42343711074968357734
  883. )
  884. /
  885. (
  886. (
  887. (
  888. (
  889. (
  890. (
  891. (
  892. t*1.05075007164441684324e-9+5.475938084995344946e-4
  893. )
  894. *t+0.0151986665636164571966
  895. )
  896. *t+0.14810397642748007459
  897. )
  898. *t+0.68976733498510000455
  899. )
  900. *t+1.6763848301838038494
  901. )
  902. *t+2.05319162663775882187
  903. )
  904. *t+1.0
  905. );
  906. }
  907. else
  908. {
  909. t+=-5.0;
  910. z=
  911. (
  912. (
  913. (
  914. (
  915. (
  916. (
  917. (
  918. t*2.01033439929228813265e-7+2.71155556874348757815e-5
  919. )
  920. *t+0.0012426609473880784386
  921. )
  922. *t+0.026532189526576123093
  923. )
  924. *t+0.29656057182850489123
  925. )
  926. *t+1.7848265399172913358
  927. )
  928. *t+5.4637849111641143699
  929. )
  930. *t+6.6579046435011037772
  931. )
  932. /
  933. (
  934. (
  935. (
  936. (
  937. (
  938. (
  939. (
  940. t*2.04426310338993978564e-15+1.4215117583164458887e-7
  941. )
  942. *t+1.8463183175100546818e-5
  943. )
  944. *t+7.868691311456132591e-4
  945. )
  946. *t+0.0148753612908506148525
  947. )
  948. *t+0.13692988092273580531
  949. )
  950. *t+0.59983220655588793769
  951. )
  952. *t+1.0
  953. );
  954. }
  955. if(q<0.0) z=-z;
  956. }
  957. return Value (z);
  958. }
  959. //helper for GetGamma and GetLogGamma
  960. static double GammaHelp(double& x, bool& bReflect)
  961. {
  962. double c[6] = {76.18009173, -86.50532033, 24.01409822,
  963. -1.231739516, 0.120858003E-2, -0.536382E-5};
  964. if (x >= 1.0)
  965. {
  966. bReflect = false;
  967. x -= 1.0;
  968. }
  969. else
  970. {
  971. bReflect = true;
  972. x = 1.0 - x;
  973. }
  974. double s, anum;
  975. s = 1.0;
  976. anum = x;
  977. for (uint i = 0; i < 6; i++)
  978. {
  979. anum += 1.0;
  980. s += c[i]/anum;
  981. }
  982. s *= 2.506628275; // sqrt(2*PI)
  983. return s;
  984. }
  985. Value ValueCalc::GetGamma (Value _x)
  986. {
  987. double x = converter->asFloat (_x).asFloat();
  988. bool bReflect;
  989. double G = GammaHelp(x, bReflect);
  990. G = ::pow(x+5.5,x+0.5)*G/::exp(x+5.5);
  991. if (bReflect)
  992. G = M_PI*x/(G*::sin(M_PI*x));
  993. return Value (G);
  994. }
  995. Value ValueCalc::GetLogGamma (Value _x)
  996. {
  997. double x = converter->asFloat (_x).asFloat();
  998. bool bReflect;
  999. double G = GammaHelp(x, bReflect);
  1000. G = (x+0.5)*::log(x+5.5)+::log(G)-(x+5.5);
  1001. if (bReflect)
  1002. G = ::log(M_PI*x)-G-::log(::sin(M_PI*x));
  1003. return Value (G);
  1004. }
  1005. Value ValueCalc::GetGammaDist (Value _x, Value _alpha,
  1006. Value _beta)
  1007. {
  1008. double x = converter->asFloat (_x).asFloat();
  1009. double alpha = converter->asFloat (_alpha).asFloat();
  1010. double beta = converter->asFloat (_beta).asFloat();
  1011. if (x == 0.0)
  1012. return Value (0.0);
  1013. x /= beta;
  1014. double gamma = alpha;
  1015. double c = 0.918938533204672741;
  1016. double d[10] = {
  1017. 0.833333333333333333E-1,
  1018. -0.277777777777777778E-2,
  1019. 0.793650793650793651E-3,
  1020. -0.595238095238095238E-3,
  1021. 0.841750841750841751E-3,
  1022. -0.191752691752691753E-2,
  1023. 0.641025641025641025E-2,
  1024. -0.295506535947712418E-1,
  1025. 0.179644372368830573,
  1026. -0.139243221690590111E1
  1027. };
  1028. double dx = x;
  1029. double dgamma = gamma;
  1030. int maxit = 10000;
  1031. double z = dgamma;
  1032. double den = 1.0;
  1033. while ( z < 10.0 ) {
  1034. den *= z;
  1035. z += 1.0;
  1036. }
  1037. double z2 = z*z;
  1038. double z3 = z*z2;
  1039. double z4 = z2*z2;
  1040. double z5 = z2*z3;
  1041. double a = ( z - 0.5 ) * ::log(z) - z + c;
  1042. double b = d[0]/z + d[1]/z3 + d[2]/z5 + d[3]/(z2*z5) + d[4]/(z4*z5) +
  1043. d[5]/(z*z5*z5) + d[6]/(z3*z5*z5) + d[7]/(z5*z5*z5) + d[8]/(z2*z5*z5*z5);
  1044. // double g = exp(a+b) / den;
  1045. double sum = 1.0 / dgamma;
  1046. double term = 1.0 / dgamma;
  1047. double cut1 = dx - dgamma;
  1048. double cut2 = dx * 10000000000.0;
  1049. for ( int i=1; i<=maxit; i++ ) {
  1050. double ai = i;
  1051. term = dx * term / ( dgamma + ai );
  1052. sum += term;
  1053. double cutoff = cut1 + ( cut2 * term / sum );
  1054. if ( ai > cutoff ) {
  1055. double t = sum;
  1056. // return pow( dx, dgamma ) * exp( -dx ) * t / g;
  1057. return Value (::exp( dgamma * ::log(dx) - dx - a - b ) * t * den);
  1058. }
  1059. }
  1060. return Value (1.0); // should not happen ...
  1061. }
  1062. Value ValueCalc::GetBeta (Value _x, Value _alpha,
  1063. Value _beta)
  1064. {
  1065. if (equal (_beta, 1.0))
  1066. return pow (_x, _alpha);
  1067. else if (equal (_alpha, 1.0))
  1068. // 1.0 - pow (1.0-_x, _beta)
  1069. return sub (1.0, pow (sub (1.0, _x), _beta));
  1070. double x = converter->asFloat (_x).asFloat();
  1071. double alpha = converter->asFloat (_alpha).asFloat();
  1072. double beta = converter->asFloat (_beta).asFloat();
  1073. double fEps = 1.0E-8;
  1074. bool bReflect;
  1075. double cf, fA, fB;
  1076. if (x < (alpha+1.0)/(alpha+beta+1.0)) {
  1077. bReflect = false;
  1078. fA = alpha;
  1079. fB = beta;
  1080. }
  1081. else {
  1082. bReflect = true;
  1083. fA = beta;
  1084. fB = alpha;
  1085. x = 1.0 - x;
  1086. }
  1087. if (x < fEps)
  1088. cf = 0.0;
  1089. else {
  1090. double a1, b1, a2, b2, fnorm, rm, apl2m, d2m, d2m1, cfnew;
  1091. a1 = 1.0; b1 = 1.0;
  1092. b2 = 1.0 - (fA+fB)*x/(fA+1.0);
  1093. if (b2 == 0.0) {
  1094. a2 = b2;
  1095. fnorm = 1.0;
  1096. cf = 1.0;
  1097. }
  1098. else {
  1099. a2 = 1.0;
  1100. fnorm = 1.0/b2;
  1101. cf = a2*fnorm;
  1102. }
  1103. cfnew = 1.0;
  1104. for (uint j = 1; j <= 100; j++) {
  1105. rm = (double) j;
  1106. apl2m = fA + 2.0*rm;
  1107. d2m = rm*(fB-rm)*x/((apl2m-1.0)*apl2m);
  1108. d2m1 = -(fA+rm)*(fA+fB+rm)*x/(apl2m*(apl2m+1.0));
  1109. a1 = (a2+d2m*a1)*fnorm;
  1110. b1 = (b2+d2m*b1)*fnorm;
  1111. a2 = a1 + d2m1*a2*fnorm;
  1112. b2 = b1 + d2m1*b2*fnorm;
  1113. if (b2 != 0.0) {
  1114. fnorm = 1.0/b2;
  1115. cfnew = a2*fnorm;
  1116. if (fabs(cf-cfnew)/cf < fEps)
  1117. j = 101;
  1118. else
  1119. cf = cfnew;
  1120. }
  1121. }
  1122. if (fB < fEps)
  1123. b1 = 1.0E30;
  1124. else
  1125. b1 = ::exp(GetLogGamma(fA).asFloat()+GetLogGamma(fB).asFloat()-
  1126. GetLogGamma(fA+fB).asFloat());
  1127. cf *= ::pow(x, fA)*::pow(1.0-x,fB)/(fA*b1);
  1128. }
  1129. if (bReflect)
  1130. return Value (1.0-cf);
  1131. else
  1132. return Value (cf);
  1133. }
  1134. // ------------------------------------------------------
  1135. /*
  1136. *
  1137. * The code for calculating Bessel functions is taken
  1138. * from CCMATH, a mathematics library source.code.
  1139. *
  1140. * Original copyright follows:
  1141. *
  1142. * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
  1143. * This code may be redistributed under the terms of the GNU library
  1144. * public license (LGPL).
  1145. */
  1146. static double ccmath_gaml(double x)
  1147. { double g,h;
  1148. for(g=1.; x<30. ;g*=x,x+=1.); h=x*x;
  1149. g=(x-.5)*log(x)-x+.918938533204672-log(g);
  1150. g+=(1.-(1./6.-(1./3.-1./(4.*h))/(7.*h))/(5.*h))/(12.*x);
  1151. return g;
  1152. }
  1153. static double ccmath_psi(int m)
  1154. { double s= -.577215664901533; int k;
  1155. for(k=1; k<m ;++k) s+=1./k;
  1156. return s;
  1157. }
  1158. static double ccmath_ibes(double v,double x)
  1159. { double y,s,t,tp; int p,m;
  1160. y=x-9.; if(y>0.) y*=y; tp=v*v*.2+25.;
  1161. if(y<tp){ x/=2.; m=(int)x;
  1162. if(x>0.) s=t=exp(v*log(x)-ccmath_gaml(v+1.));
  1163. else{ if(v>0.) return 0.; else if(v==0.) return 1.;}
  1164. for(p=1,x*=x;;++p){ t*=x/(p*(v+=1.)); s+=t;
  1165. if(p>m && t<1.e-13*s) break;
  1166. }
  1167. }
  1168. else{ double u,a0=1.57079632679490;
  1169. s=t=1./sqrt(x*a0); x*=2.; u=0.;
  1170. for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
  1171. t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>=tp) break;
  1172. if(!(p&1)) s+=t; else u-=t;
  1173. }
  1174. x/=2.; s=cosh(x)*s+sinh(x)*u;
  1175. }
  1176. return s;
  1177. }
  1178. static double ccmath_kbes(double v,double x)
  1179. { double y,s,t,tp,f,a0=1.57079632679490;
  1180. int p,k,m;
  1181. if(x==0.) return HUGE_VAL;
  1182. y=x-10.5; if(y>0.) y*=y; tp=25.+.185*v*v;
  1183. if(y<tp && modf(v+.5,&t)!=0.){ y=1.5+.5*v;
  1184. if(x<y){ x/=2.; m=(int)x; tp=t=exp(v*log(x)-ccmath_gaml(v+1.));
  1185. if(modf(v,&y)==0.){ k=(int)y; tp*=v;
  1186. f=2.*log(x)-ccmath_psi(1)-ccmath_psi(k+1);
  1187. t/=2.; if(!(k&1)) t= -t; s=f*t;
  1188. for(p=1,x*=x;;++p){ f-=1./p+1./(v+=1.);
  1189. t*=x/(p*v); s+=(y=t*f);
  1190. if(p>m && fabs(y)<1.e-14) break; }
  1191. if(k>0){ x= -x; s+=(t=1./(tp*2.));
  1192. for(p=1,--k; k>0 ;++p,--k) s+=(t*=x/(p*k)); }
  1193. }
  1194. else{ f=1./(t*v*2.); t*=a0/sin(2.*a0*v); s=f-t;
  1195. for(p=1,x*=x,tp=v;;++p){
  1196. t*=x/(p*(v+=1.)); f*= -x/(p*(tp-=1.));
  1197. s+=(y=f-t); if(p>m && fabs(y)<1.e-14) break; }
  1198. }
  1199. }
  1200. else{ double tq,h,w,z,r;
  1201. t=12./pow(x,.333); k=(int)(t*t); y=2.*(x+k);
  1202. m=(int)v; v-=m; tp=v*v-.25; v+=1.; tq=v*v-.25;
  1203. for(s=h=1.,r=f=z=w=0.; k>0 ;--k,y-=2.){
  1204. t=(y*h-(k+1)*z)/(k-1-tp/k); z=h; f+=(h=t);
  1205. t=(y*s-(k+1)*w)/(k-1-tq/k); w=s; r+=(s=t); }
  1206. t=sqrt(a0/x)*exp(-x); s*=t/r; h*=t/f; x/=2.; if(m==0) s=h;
  1207. for(k=1; k<m ;++k){ t=v*s/x+h; h=s; s=t; v+=1.;}
  1208. }
  1209. }
  1210. else{ s=t=sqrt(a0/x); x*=2.;
  1211. for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
  1212. t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>=tp) break; s+=t; }
  1213. s*=exp(-x/2.);
  1214. }
  1215. return s;
  1216. }
  1217. static double ccmath_jbes(double v,double x)
  1218. { double y,s,t,tp; int p,m;
  1219. y=x-8.5; if(y>0.) y*=y; tp=v*v/4.+13.69;
  1220. if(y<tp){ x/=2.; m=(int)x;
  1221. if(x>0.) s=t=exp(v*log(x)-ccmath_gaml(v+1.));
  1222. else{ if(v>0.) return 0.; else if(v==0.) return 1.;}
  1223. for(p=1,x*= -x;;++p){ t*=x/(p*(v+=1.)); s+=t;
  1224. if(p>m && fabs(t)<1.e-13) break;
  1225. }
  1226. }
  1227. else{ double u,a0=1.57079632679490;
  1228. s=t=1./sqrt(x*a0); x*=2.; u=0.;
  1229. for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
  1230. t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>=tp) break;
  1231. if(!(p&1)){ t= -t; s+=t;} else u-=t;
  1232. }
  1233. y=x/2.-(v+.5)*a0; s=cos(y)*s+sin(y)*u;
  1234. }
  1235. return s;
  1236. }
  1237. static double ccmath_nbes(double v,double x)
  1238. { double y,s,t,tp,u,f,a0=3.14159265358979;
  1239. int p,k,m;
  1240. y=x-8.5; if(y>0.) y*=y; tp=v*v/4.+13.69;
  1241. if(y<tp){ if(x==0.) return HUGE_VAL;
  1242. x/=2.; m=(int)x; u=t=exp(v*log(x)-ccmath_gaml(v+1.));
  1243. if(modf(v,&y)==0.){ k=(int)y; u*=v;
  1244. f=2.*log(x)-ccmath_psi(1)-ccmath_psi(k+1);
  1245. t/=a0; x*= -x; s=f*t;
  1246. for(p=1;;++p){ f-=1./p+1./(v+=1.);
  1247. t*=x/(p*v); s+=(y=t*f); if(p>m && fabs(y)<1.e-13) break; }
  1248. if(k>0){ x= -x; s-=(t=1./(u*a0));
  1249. for(p=1,--k; k>0 ;++p,--k) s-=(t*=x/(p*k)); }
  1250. }
  1251. else{ f=1./(t*v*a0); t/=tan(a0*v); s=t-f;
  1252. for(p=1,x*=x,u=v;;++p){
  1253. t*= -x/(p*(v+=1.)); f*=x/(p*(u-=1.));
  1254. s+=(y=t-f); if(p>m && fabs(y)<1.e-13) break; }
  1255. }
  1256. }
  1257. else{ x*=2.; s=t=2./sqrt(x*a0); u=0.;
  1258. for(p=1,y=.5; (tp=fabs(t))>1.e-14 ;++p,y+=1.){
  1259. t*=(v+y)*(v-y)/(p*x); if(y>v && fabs(t)>tp) break;
  1260. if(!(p&1)){ t= -t; s+=t;} else u+=t;
  1261. }
  1262. y=(x-(v+.5)*a0)/2.; s=sin(y)*s+cos(y)*u;
  1263. }
  1264. return s;
  1265. }
  1266. /* ---------- end of CCMATH code ---------- */
  1267. Value ValueCalc::besseli (Value v, Value x)
  1268. {
  1269. double vv = converter->asFloat (v).asFloat ();
  1270. double xx = converter->asFloat (x).asFloat ();
  1271. return Value (ccmath_ibes (vv, xx));
  1272. }
  1273. Value ValueCalc::besselj (Value v, Value x)
  1274. {
  1275. double vv = converter->asFloat (v).asFloat ();
  1276. double xx = converter->asFloat (x).asFloat ();
  1277. return Value (ccmath_jbes (vv, xx));
  1278. }
  1279. Value ValueCalc::besselk (Value v, Value x)
  1280. {
  1281. double vv = converter->asFloat (v).asFloat ();
  1282. double xx = converter->asFloat (x).asFloat ();
  1283. return Value (ccmath_kbes (vv, xx));
  1284. }
  1285. Value ValueCalc::besseln (Value v, Value x)
  1286. {
  1287. double vv = converter->asFloat (v).asFloat ();
  1288. double xx = converter->asFloat (x).asFloat ();
  1289. return Value (ccmath_nbes (vv, xx));
  1290. }
  1291. // ------------------------------------------------------
  1292. Value ValueCalc::erf (Value x)
  1293. {
  1294. return Value (::erf (converter->asFloat (x).asFloat()));
  1295. }
  1296. Value ValueCalc::erfc (Value x)
  1297. {
  1298. return Value (::erfc (converter->asFloat (x).asFloat()));
  1299. }
  1300. // ------------------------------------------------------
  1301. void ValueCalc::arrayWalk (const Value &range,
  1302. Value &res, arrayWalkFunc func, Value param)
  1303. {
  1304. if (res.isError()) return;
  1305. if (!range.isArray ())
  1306. {
  1307. func (this, res, range, param);
  1308. return;
  1309. }
  1310. int rows = range.rows ();
  1311. int cols = range.columns ();
  1312. for (int r = 0; r < rows; r++)
  1313. for (int c = 0; c < cols; c++)
  1314. {
  1315. Value v = range.element (c, r);
  1316. if (v.isArray())
  1317. arrayWalk (v, res, func, param);
  1318. else {
  1319. func (this, res, v, param);
  1320. if (res.format() == Value::fmt_None)
  1321. res.setFormat (v.format());
  1322. }
  1323. }
  1324. }
  1325. void ValueCalc::arrayWalk (TQValueVector<Value> &range,
  1326. Value &res, arrayWalkFunc func, Value param)
  1327. {
  1328. if (res.isError()) return;
  1329. for (unsigned int i = 0; i < range.count(); ++i)
  1330. arrayWalk (range[i], res, func, param);
  1331. }
  1332. void ValueCalc::twoArrayWalk (const Value &a1, const Value &a2,
  1333. Value &res, arrayWalkFunc func)
  1334. {
  1335. if (res.isError()) return;
  1336. if (!a1.isArray ())
  1337. {
  1338. func (this, res, a1, a2);
  1339. return;
  1340. }
  1341. int rows = a1.rows ();
  1342. int cols = a1.columns ();
  1343. int rows2 = a2.rows ();
  1344. int cols2 = a2.columns ();
  1345. if ((rows != rows2) || (cols != cols2)) {
  1346. res = Value::errorVALUE();
  1347. return;
  1348. }
  1349. for (int r = 0; r < rows; r++)
  1350. for (int c = 0; c < cols; c++)
  1351. {
  1352. Value v1 = a1.element (c, r);
  1353. Value v2 = a2.element (c, r);
  1354. if (v1.isArray() && v2.isArray())
  1355. twoArrayWalk (v1, v2, res, func);
  1356. else {
  1357. func (this, res, v1, v2);
  1358. if (res.format() == Value::fmt_None)
  1359. res.setFormat (format (v1.format(), v2.format()));
  1360. }
  1361. }
  1362. }
  1363. void ValueCalc::twoArrayWalk (TQValueVector<Value> &a1,
  1364. TQValueVector<Value> &a2, Value &res, arrayWalkFunc func)
  1365. {
  1366. if (res.isError()) return;
  1367. if (a1.count() != a2.count()) {
  1368. res = Value::errorVALUE();
  1369. return;
  1370. }
  1371. for (unsigned int i = 0; i < a1.count(); ++i)
  1372. twoArrayWalk (a1[i], a2[i], res, func);
  1373. }
  1374. arrayWalkFunc ValueCalc::awFunc (const TQString &name)
  1375. {
  1376. if (awFuncs.count(name))
  1377. return awFuncs[name];
  1378. return 0;
  1379. }
  1380. void ValueCalc::registerAwFunc (const TQString &name, arrayWalkFunc func)
  1381. {
  1382. awFuncs[name] = func;
  1383. }
  1384. // ------------------------------------------------------
  1385. Value ValueCalc::sum (const Value &range, bool full)
  1386. {
  1387. Value res;
  1388. arrayWalk (range, res, awFunc (full ? "suma" : "sum"), 0);
  1389. return res;
  1390. }
  1391. Value ValueCalc::sum (TQValueVector<Value> range, bool full)
  1392. {
  1393. Value res;
  1394. arrayWalk (range, res, awFunc (full ? "suma" : "sum"), 0);
  1395. return res;
  1396. }
  1397. // sum of squares
  1398. Value ValueCalc::sumsq (const Value &range, bool full)
  1399. {
  1400. Value res;
  1401. arrayWalk (range, res, awFunc (full ? "sumsqa" : "sumsq"), 0);
  1402. return res;
  1403. }
  1404. Value ValueCalc::sumIf (const Value &range,
  1405. const Value &checkRange, const Condition &cond)
  1406. {
  1407. if (!range.isArray())
  1408. {
  1409. if (matches (cond, checkRange.element (0, 0)))
  1410. return converter->asNumeric (range);
  1411. return Value (0.0);
  1412. }
  1413. //if we are here, we have an array
  1414. Value res;
  1415. unsigned int rows = range.rows ();
  1416. unsigned int cols = range.columns ();
  1417. for (unsigned int r = 0; r < rows; r++)
  1418. for (unsigned int c = 0; c < cols; c++)
  1419. {
  1420. Value v = range.element (c, r);
  1421. Value newcheck = v;
  1422. if ((c < checkRange.columns()) && (r < checkRange.rows()))
  1423. newcheck = checkRange.element (c, r);
  1424. if (v.isArray())
  1425. res = add (res, sumIf (v, newcheck, cond));
  1426. else
  1427. if (matches (cond, newcheck))
  1428. res = add (res, v);
  1429. }
  1430. return res;
  1431. }
  1432. int ValueCalc::count (const Value &range, bool full)
  1433. {
  1434. Value res = 0;
  1435. arrayWalk (range, res, awFunc (full ? "counta" : "count"), 0);
  1436. return converter->asInteger (res).asInteger ();
  1437. }
  1438. int ValueCalc::count (TQValueVector<Value> range, bool full)
  1439. {
  1440. Value res = 0;
  1441. arrayWalk (range, res, awFunc (full ? "counta" : "count"), 0);
  1442. return converter->asInteger (res).asInteger ();
  1443. }
  1444. int ValueCalc::countIf (const Value &range, const Condition &cond)
  1445. {
  1446. if (!range.isArray())
  1447. {
  1448. if (matches (cond, range))
  1449. return range.isEmpty() ? 0 : 1;
  1450. return 0;
  1451. }
  1452. int res = 0;
  1453. int cols = range.columns ();
  1454. int rows = range.rows ();
  1455. for (int r = 0; r < rows; r++)
  1456. for (int c = 0; c < cols; c++)
  1457. {
  1458. Value v = range.element (c, r);
  1459. if (v.isArray())
  1460. res += countIf (v, cond);
  1461. else
  1462. if (matches (cond, v))
  1463. res++;
  1464. }
  1465. return res;
  1466. }
  1467. Value ValueCalc::avg (const Value &range, bool full)
  1468. {
  1469. int cnt = count (range, full);
  1470. if (cnt)
  1471. return div (sum (range, full), cnt);
  1472. return Value (0.0);
  1473. }
  1474. Value ValueCalc::avg (TQValueVector<Value> range, bool full)
  1475. {
  1476. int cnt = count (range, full);
  1477. if (cnt)
  1478. return div (sum (range, full), cnt);
  1479. return Value (0.0);
  1480. }
  1481. Value ValueCalc::max (const Value &range, bool full)
  1482. {
  1483. Value res;
  1484. arrayWalk (range, res, awFunc (full ? "maxa" : "max"), 0);
  1485. return res;
  1486. }
  1487. Value ValueCalc::max (TQValueVector<Value> range, bool full)
  1488. {
  1489. Value res;
  1490. arrayWalk (range, res, awFunc (full ? "maxa" : "max"), 0);
  1491. return res;
  1492. }
  1493. Value ValueCalc::min (const Value &range, bool full)
  1494. {
  1495. Value res;
  1496. arrayWalk (range, res, awFunc (full ? "mina" : "min"), 0);
  1497. return res;
  1498. }
  1499. Value ValueCalc::min (TQValueVector<Value> range, bool full)
  1500. {
  1501. Value res;
  1502. arrayWalk (range, res, awFunc (full ? "mina" : "min"), 0);
  1503. return res;
  1504. }
  1505. Value ValueCalc::product (const Value &range, Value init,
  1506. bool full)
  1507. {
  1508. Value res = init;
  1509. if (isZero (init)) // special handling of a zero, due to excel-compat
  1510. {
  1511. if (count (range, full) == 0)
  1512. return init;
  1513. res = 1.0;
  1514. }
  1515. arrayWalk (range, res, awFunc (full ? "proda" : "prod"), 0);
  1516. return res;
  1517. }
  1518. Value ValueCalc::product (TQValueVector<Value> range,
  1519. Value init, bool full)
  1520. {
  1521. Value res = init;
  1522. if (isZero (init)) // special handling of a zero, due to excel-compat
  1523. {
  1524. if (count (range, full) == 0)
  1525. return init;
  1526. res = 1.0;
  1527. }
  1528. arrayWalk (range, res, awFunc (full ? "proda" : "prod"), 0);
  1529. return res;
  1530. }
  1531. Value ValueCalc::stddev (const Value &range, bool full)
  1532. {
  1533. return stddev (range, avg (range, full), full);
  1534. }
  1535. Value ValueCalc::stddev (const Value &range, Value avg,
  1536. bool full)
  1537. {
  1538. Value res;
  1539. int cnt = count (range, full);
  1540. arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  1541. return sqrt (div (res, cnt-1));
  1542. }
  1543. Value ValueCalc::stddev (TQValueVector<Value> range, bool full)
  1544. {
  1545. return stddev (range, avg (range, full), full);
  1546. }
  1547. Value ValueCalc::stddev (TQValueVector<Value> range,
  1548. Value avg, bool full)
  1549. {
  1550. Value res;
  1551. int cnt = count (range, full);
  1552. arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  1553. return sqrt (div (res, cnt-1));
  1554. }
  1555. Value ValueCalc::stddevP (const Value &range, bool full)
  1556. {
  1557. return stddevP (range, avg (range, full), full);
  1558. }
  1559. Value ValueCalc::stddevP (const Value &range, Value avg,
  1560. bool full)
  1561. {
  1562. Value res;
  1563. int cnt = count (range, full);
  1564. arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  1565. return sqrt (div (res, cnt));
  1566. }
  1567. Value ValueCalc::stddevP (TQValueVector<Value> range, bool full)
  1568. {
  1569. return stddevP (range, avg (range, full), full);
  1570. }
  1571. Value ValueCalc::stddevP (TQValueVector<Value> range,
  1572. Value avg, bool full)
  1573. {
  1574. Value res;
  1575. int cnt = count (range, full);
  1576. arrayWalk (range, res, awFunc (full ? "devsqa" : "devsq"), avg);
  1577. return sqrt (div (res, cnt));
  1578. }
  1579. Value::Format ValueCalc::format (Value::Format a,
  1580. Value::Format b)
  1581. {
  1582. if ((a == Value::fmt_None) || (a == Value::fmt_Boolean))
  1583. return b;
  1584. return a;
  1585. }
  1586. // ------------------------------------------------------
  1587. void ValueCalc::getCond (Condition &cond, Value val)
  1588. {
  1589. // not a string - we simply take it as a numeric value
  1590. // that also handles floats, logical values, date/time and such
  1591. if (!val.isString()) {
  1592. cond.comp = isEqual;
  1593. cond.type = numeric;
  1594. cond.value = converter->asFloat (val).asFloat();
  1595. return;
  1596. }
  1597. TQString text = converter->asString (val).asString();
  1598. cond.comp = isEqual;
  1599. text = text.stripWhiteSpace();
  1600. if ( text.startsWith( "<=" ) )
  1601. {
  1602. cond.comp = lessEqual;
  1603. text = text.remove( 0, 2 );
  1604. }
  1605. else if ( text.startsWith( ">=" ) )
  1606. {
  1607. cond.comp = greaterEqual;
  1608. text = text.remove( 0, 2 );
  1609. }
  1610. else if ( text.startsWith( "!=" ) || text.startsWith( "<>" ) )
  1611. {
  1612. cond.comp = notEqual;
  1613. text = text.remove( 0, 2 );
  1614. }
  1615. else if ( text.startsWith( "==" ) )
  1616. {
  1617. cond.comp = isEqual;
  1618. text = text.remove( 0, 2 );
  1619. }
  1620. else if ( text.startsWith( "<" ) )
  1621. {
  1622. cond.comp = isLess;
  1623. text = text.remove( 0, 1 );
  1624. }
  1625. else if ( text.startsWith( ">" ) )
  1626. {
  1627. cond.comp = isGreater;
  1628. text = text.remove( 0, 1 );
  1629. }
  1630. else if ( text.startsWith( "=" ) )
  1631. {
  1632. cond.comp = isEqual;
  1633. text = text.remove( 0, 1 );
  1634. }
  1635. text = text.stripWhiteSpace();
  1636. bool ok = false;
  1637. double d = text.toDouble( &ok );
  1638. if ( ok )
  1639. {
  1640. cond.type = numeric;
  1641. cond.value = d;
  1642. }
  1643. else
  1644. {
  1645. cond.type = string;
  1646. cond.stringValue = text;
  1647. }
  1648. }
  1649. bool ValueCalc::matches (const Condition &cond, Value val)
  1650. {
  1651. if (val.isEmpty())
  1652. return false;
  1653. if (cond.type == numeric) {
  1654. double d = converter->asFloat (val).asFloat();
  1655. switch ( cond.comp )
  1656. {
  1657. case isEqual:
  1658. if (approxEqual (d, cond.value)) return true;
  1659. break;
  1660. case isLess:
  1661. if (d < cond.value) return true;
  1662. break;
  1663. case isGreater:
  1664. if (d > cond.value) return true;
  1665. break;
  1666. case lessEqual:
  1667. if (d <= cond.value) return true;
  1668. break;
  1669. case greaterEqual:
  1670. if (d >= cond.value) return true;
  1671. break;
  1672. case notEqual:
  1673. if (d != cond.value) return true;
  1674. break;
  1675. }
  1676. } else {
  1677. TQString d = converter->asString (val).asString();
  1678. switch ( cond.comp )
  1679. {
  1680. case isEqual:
  1681. if (d == cond.stringValue) return true;
  1682. break;
  1683. case isLess:
  1684. if (d < cond.stringValue) return true;
  1685. break;
  1686. case isGreater:
  1687. if (d > cond.stringValue) return true;
  1688. break;
  1689. case lessEqual:
  1690. if (d <= cond.stringValue) return true;
  1691. break;
  1692. case greaterEqual:
  1693. if (d >= cond.stringValue) return true;
  1694. break;
  1695. case notEqual:
  1696. if (d != cond.stringValue) return true;
  1697. break;
  1698. }
  1699. }
  1700. return false;
  1701. }