1 | /* This file is part of the KDE project |
2 | Copyright (C) 2008-2010 Marijn Kruisselbrink <mkruisselbrink@kde.org> |
3 | Copyright (C) 2005-2007 Tomas Mecir <mecirt@gmail.com> |
4 | |
5 | This library is free software; you can redistribute it and/or |
6 | modify it under the terms of the GNU Library General Public |
7 | License as published by the Free Software Foundation; either |
8 | version 2 of the License, or (at your option) any later version. |
9 | |
10 | This library is distributed in the hope that it will be useful, |
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
13 | Library General Public License for more details. |
14 | |
15 | You should have received a copy of the GNU Library General Public License |
16 | along with this library; see the file COPYING.LIB. If not, write to |
17 | the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, |
18 | Boston, MA 02110-1301, USA. |
19 | */ |
20 | |
21 | #include "ValueCalc.h" |
22 | |
23 | #include "Cell.h" |
24 | #include "Number.h" |
25 | #include "ValueConverter.h" |
26 | #include "CalculationSettings.h" |
27 | |
28 | #include <QRegExp> |
29 | |
30 | #include <kdebug.h> |
31 | #include <errno.h> |
32 | #include <float.h> |
33 | #include <math.h> |
34 | #include <stdlib.h> |
35 | #include <time.h> |
36 | |
37 | using namespace Calligra::Sheets; |
38 | |
39 | |
40 | // |
41 | // helper: gammaHelper |
42 | // |
43 | // args[0] = value |
44 | // args[1] = & reflect |
45 | // |
46 | static double GammaHelp(double& x, bool& reflect) |
47 | { |
48 | double c[6] = { 76.18009173, -86.50532033 , 24.01409822, |
49 | -1.231739516, 0.120858003E-2, -0.536382E-5 |
50 | }; |
51 | if (x >= 1.0) { |
52 | reflect = false; |
53 | x -= 1.0; |
54 | } else { |
55 | reflect = true; |
56 | x = 1.0 - x; |
57 | } |
58 | double res, anum; |
59 | res = 1.0; |
60 | anum = x; |
61 | for (int i = 0; i < 6; ++i) { |
62 | anum += 1.0; |
63 | res += c[i] / anum; |
64 | } |
65 | res *= 2.506628275; // sqrt(2*PI) |
66 | |
67 | return res; |
68 | } |
69 | |
70 | // Array-walk functions registered on ValueCalc object |
71 | |
72 | void awSum(ValueCalc *c, Value &res, Value val, Value) |
73 | { |
74 | if (val.isError()) |
75 | res = val; |
76 | else if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString())) |
77 | res = c->add(res, val); |
78 | } |
79 | |
80 | void awSumA(ValueCalc *c, Value &res, Value val, Value) |
81 | { |
82 | if (!val.isEmpty()) |
83 | res = c->add(res, val); |
84 | } |
85 | |
86 | void awSumSq(ValueCalc *c, Value &res, Value val, Value) |
87 | { |
88 | // removed (!val.isBoolean()) to allow conversion from BOOL to int |
89 | if ((!val.isEmpty()) && (!val.isString()) && (!val.isError())) |
90 | res = c->add(res, c->sqr(val)); |
91 | } |
92 | |
93 | void awSumSqA(ValueCalc *c, Value &res, Value val, Value) |
94 | { |
95 | if (!val.isEmpty()) |
96 | res = c->add(res, c->sqr(val)); |
97 | } |
98 | |
99 | void awCount(ValueCalc *c, Value &res, Value val, Value) |
100 | { |
101 | if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()) && (!val.isError())) |
102 | res = c->add(res, 1); |
103 | } |
104 | |
105 | void awCountA(ValueCalc *c, Value &res, Value val, Value) |
106 | { |
107 | if (!val.isEmpty()) |
108 | res = c->add(res, 1); |
109 | } |
110 | |
111 | void awMax(ValueCalc *c, Value &res, Value val, Value) |
112 | { |
113 | // propagate error values |
114 | if (res.isError()) |
115 | return; |
116 | if (val.isError()) |
117 | res = val; |
118 | else if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString())) { |
119 | if (res.isEmpty()) { |
120 | res = val; |
121 | } else { |
122 | if (c->greater(val, res)) res = val; |
123 | } |
124 | } |
125 | } |
126 | |
127 | void awMaxA(ValueCalc *c, Value &res, Value val, Value) |
128 | { |
129 | // propagate error values |
130 | if (res.isError()) |
131 | return; |
132 | if (val.isError()) { |
133 | res = val; |
134 | } else if (!val.isEmpty()) { |
135 | if (res.isEmpty()) |
136 | // convert to number, so that we don't return string/bool |
137 | res = c->conv()->asNumeric(val); |
138 | else if (c->greater(val, res)) |
139 | // convert to number, so that we don't return string/bool |
140 | res = c->conv()->asNumeric(val); |
141 | } |
142 | } |
143 | |
144 | void awMin(ValueCalc *c, Value &res, Value val, Value) |
145 | { |
146 | if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString())) { |
147 | if (res.isEmpty()) |
148 | res = val; |
149 | else if (c->lower(val, res)) res = val; |
150 | } |
151 | } |
152 | |
153 | void awMinA(ValueCalc *c, Value &res, Value val, Value) |
154 | { |
155 | if (!val.isEmpty()) { |
156 | if (res.isEmpty()) |
157 | // convert to number, so that we don't return string/bool |
158 | res = c->conv()->asNumeric(val); |
159 | else if (c->lower(val, res)) |
160 | // convert to number, so that we don't return string/bool |
161 | res = c->conv()->asNumeric(val); |
162 | } |
163 | } |
164 | |
165 | void awProd(ValueCalc *c, Value &res, Value val, Value) |
166 | { |
167 | if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()) && (!val.isError())) |
168 | res = c->mul(res, val); |
169 | } |
170 | |
171 | void awProdA(ValueCalc *c, Value &res, Value val, Value) |
172 | { |
173 | if (!val.isEmpty()) |
174 | res = c->mul(res, val); |
175 | } |
176 | |
177 | // sum of squares of deviations, used to compute standard deviation |
178 | void awDevSq(ValueCalc *c, Value &res, Value val, |
179 | Value avg) |
180 | { |
181 | if ((!val.isEmpty()) && (!val.isBoolean()) && (!val.isString()) && (!val.isError())) |
182 | res = c->add(res, c->sqr(c->sub(val, avg))); |
183 | } |
184 | |
185 | // sum of squares of deviations, used to compute standard deviation |
186 | void awDevSqA(ValueCalc *c, Value &res, Value val, |
187 | Value avg) |
188 | { |
189 | if (!val.isEmpty()) |
190 | res = c->add(res, c->sqr(c->sub(val, avg))); |
191 | } |
192 | |
193 | |
194 | // *********************** |
195 | // ****** ValueCalc ****** |
196 | // *********************** |
197 | |
198 | ValueCalc::ValueCalc(ValueConverter* c): converter(c) |
199 | { |
200 | // initialize the random number generator |
201 | srand(time(0)); |
202 | |
203 | // register array-walk functions |
204 | registerAwFunc("sum" , awSum); |
205 | registerAwFunc("suma" , awSumA); |
206 | registerAwFunc("sumsq" , awSumSq); |
207 | registerAwFunc("sumsqa" , awSumSqA); |
208 | registerAwFunc("count" , awCount); |
209 | registerAwFunc("counta" , awCountA); |
210 | registerAwFunc("max" , awMax); |
211 | registerAwFunc("maxa" , awMaxA); |
212 | registerAwFunc("min" , awMin); |
213 | registerAwFunc("mina" , awMinA); |
214 | registerAwFunc("prod" , awProd); |
215 | registerAwFunc("proda" , awProdA); |
216 | registerAwFunc("devsq" , awDevSq); |
217 | registerAwFunc("devsqa" , awDevSqA); |
218 | } |
219 | |
220 | const CalculationSettings* ValueCalc::settings() const |
221 | { |
222 | return converter->settings(); |
223 | } |
224 | |
225 | Value ValueCalc::add(const Value &a, const Value &b) |
226 | { |
227 | if (a.isError()) return a; |
228 | if (b.isError()) return b; |
229 | if (a.isArray() || b.isArray()) |
230 | return twoArrayMap(a, &ValueCalc::add, b); |
231 | |
232 | Number aa, bb; |
233 | aa = converter->toFloat(a); |
234 | bb = converter->toFloat(b); |
235 | Value res = Value(aa + bb); |
236 | |
237 | if (a.isNumber() || a.isEmpty()) |
238 | res.setFormat(format(a, b)); |
239 | |
240 | return res; |
241 | } |
242 | |
243 | Value ValueCalc::sub(const Value &a, const Value &b) |
244 | { |
245 | if (a.isError()) return a; |
246 | if (b.isError()) return b; |
247 | if (a.isArray() || b.isArray()) |
248 | return twoArrayMap(a, &ValueCalc::sub, b); |
249 | |
250 | Number aa, bb; |
251 | aa = converter->toFloat(a); |
252 | bb = converter->toFloat(b); |
253 | Value res = Value(aa - bb); |
254 | |
255 | if (a.isNumber() || a.isEmpty()) |
256 | res.setFormat(format(a, b)); |
257 | |
258 | return res; |
259 | } |
260 | |
261 | Value ValueCalc::mul(const Value &a, const Value &b) |
262 | { |
263 | if (a.isError()) return a; |
264 | if (b.isError()) return b; |
265 | // This operation is only defined for an array if it is multiplied |
266 | // with a number, that however is commutative, thus swap parameters |
267 | // if necessary. |
268 | if (a.isArray() && !b.isArray()) |
269 | return arrayMap(a, &ValueCalc::mul, b); |
270 | if (b.isArray() && !a.isArray()) |
271 | return arrayMap(b, &ValueCalc::mul, a); |
272 | |
273 | Number aa, bb; |
274 | aa = converter->toFloat(a); |
275 | bb = converter->toFloat(b); |
276 | Value res = Value(aa * bb); |
277 | |
278 | if (a.isNumber() || a.isEmpty()) |
279 | res.setFormat(format(a, b)); |
280 | |
281 | return res; |
282 | } |
283 | |
284 | Value ValueCalc::div(const Value &a, const Value &b) |
285 | { |
286 | if (a.isError()) return a; |
287 | if (b.isError()) return b; |
288 | if (a.isArray() && !b.isArray()) |
289 | return arrayMap(a, &ValueCalc::div, b); |
290 | |
291 | Number aa, bb; |
292 | aa = converter->toFloat(a); |
293 | bb = converter->toFloat(b); |
294 | Value res; |
295 | if (bb == 0.0) |
296 | return Value::errorDIV0(); |
297 | else |
298 | res = Value(aa / bb); |
299 | |
300 | if (a.isNumber() || a.isEmpty()) |
301 | res.setFormat(format(a, b)); |
302 | |
303 | return res; |
304 | } |
305 | |
306 | Value ValueCalc::mod(const Value &a, const Value &b) |
307 | { |
308 | if (a.isError()) return a; |
309 | if (b.isError()) return b; |
310 | if (a.isArray() && !b.isArray()) |
311 | return arrayMap(a, &ValueCalc::mod, b); |
312 | |
313 | Number aa, bb; |
314 | aa = converter->toFloat(a); |
315 | bb = converter->toFloat(b); |
316 | Value res; |
317 | if (bb == 0.0) |
318 | return Value::errorDIV0(); |
319 | else { |
320 | Number m = fmod(aa, bb); |
321 | // the following adjustments are needed by OpenFormula: |
322 | // can't simply use fixed increases/decreases, because the implementation |
323 | // of fmod may differ on various platforms, and we should always return |
324 | // the same results ... |
325 | if ((bb > 0) && (aa < 0)) // result must be positive here |
326 | while (m < 0) m += bb; |
327 | if (bb < 0) { // result must be negative here, but not lower than bb |
328 | // bb is negative, hence the following two are correct |
329 | while (m < bb) m -= bb; // same as m+=fabs(bb) |
330 | while (m > 0) m += bb; // same as m-=fabs(bb) |
331 | } |
332 | |
333 | res = Value(m); |
334 | } |
335 | |
336 | if (a.isNumber() || a.isEmpty()) |
337 | res.setFormat(format(a, b)); |
338 | |
339 | return res; |
340 | } |
341 | |
342 | Value ValueCalc::pow(const Value &a, const Value &b) |
343 | { |
344 | if (a.isError()) return a; |
345 | if (b.isError()) return b; |
346 | if (a.isArray() && !b.isArray()) |
347 | return arrayMap(a, &ValueCalc::pow, b); |
348 | |
349 | Number aa, bb; |
350 | aa = converter->toFloat(a); |
351 | bb = converter->toFloat(b); |
352 | Value res(::pow(aa, bb)); |
353 | |
354 | if (a.isNumber() || a.isEmpty()) |
355 | res.setFormat(format(a, b)); |
356 | |
357 | return res; |
358 | } |
359 | |
360 | Value ValueCalc::sqr(const Value &a) |
361 | { |
362 | if (a.isError()) return a; |
363 | return mul(a, a); |
364 | } |
365 | |
366 | Value ValueCalc::sqrt(const Value &a) |
367 | { |
368 | if (a.isError()) return a; |
369 | Value res = Value(::pow((qreal)converter->toFloat(a), 0.5)); |
370 | if (a.isNumber() || a.isEmpty()) |
371 | res.setFormat(a.format()); |
372 | |
373 | return res; |
374 | } |
375 | |
376 | Value ValueCalc::add(const Value &a, Number b) |
377 | { |
378 | if (a.isError()) return a; |
379 | Value res = Value(converter->toFloat(a) + b); |
380 | |
381 | if (a.isNumber() || a.isEmpty()) |
382 | res.setFormat(a.format()); |
383 | |
384 | return res; |
385 | } |
386 | |
387 | Value ValueCalc::sub(const Value &a, Number b) |
388 | { |
389 | if (a.isError()) return a; |
390 | Value res = Value(converter->toFloat(a) - b); |
391 | |
392 | if (a.isNumber() || a.isEmpty()) |
393 | res.setFormat(a.format()); |
394 | |
395 | return res; |
396 | } |
397 | |
398 | Value ValueCalc::mul(const Value &a, Number b) |
399 | { |
400 | if (a.isError()) return a; |
401 | Value res = Value(converter->toFloat(a) * b); |
402 | |
403 | if (a.isNumber() || a.isEmpty()) |
404 | res.setFormat(a.format()); |
405 | |
406 | return res; |
407 | } |
408 | |
409 | Value ValueCalc::div(const Value &a, Number b) |
410 | { |
411 | if (a.isError()) return a; |
412 | Value res; |
413 | if (b == 0.0) |
414 | return Value::errorDIV0(); |
415 | |
416 | res = Value(converter->toFloat(a) / b); |
417 | |
418 | if (a.isNumber() || a.isEmpty()) |
419 | res.setFormat(a.format()); |
420 | |
421 | return res; |
422 | } |
423 | |
424 | Value ValueCalc::pow(const Value &a, Number b) |
425 | { |
426 | if (a.isError()) return a; |
427 | Value res = Value(::pow(converter->toFloat(a), b)); |
428 | |
429 | if (a.isNumber() || a.isEmpty()) |
430 | res.setFormat(a.format()); |
431 | |
432 | return res; |
433 | } |
434 | |
435 | Value ValueCalc::abs(const Value &a) |
436 | { |
437 | if (a.isError()) return a; |
438 | return Value(fabs(converter->toFloat(a))); |
439 | } |
440 | |
441 | bool ValueCalc::isZero(const Value &a) |
442 | { |
443 | if (a.isError()) return false; |
444 | return (converter->toFloat(a) == 0.0); |
445 | } |
446 | |
447 | bool ValueCalc::isEven(const Value &a) |
448 | { |
449 | if (a.isError()) return false; |
450 | if (gequal(a, Value(0))) { |
451 | return ((converter->toInteger(roundDown(a)) % 2) == 0); |
452 | } else { |
453 | return ((converter->toInteger(roundUp(a)) % 2) == 0); |
454 | } |
455 | } |
456 | |
457 | bool ValueCalc::equal(const Value &a, const Value &b) |
458 | { |
459 | return (converter->toFloat(a) == converter->toFloat(b)); |
460 | } |
461 | |
462 | /********************************************************************* |
463 | * |
464 | * Helper function to avoid problems with rounding floating point |
465 | * values. Idea for this kind of solution taken from Openoffice. |
466 | * |
467 | *********************************************************************/ |
468 | bool ValueCalc::approxEqual(const Value &a, const Value &b) |
469 | { |
470 | Number aa = converter->toFloat(a); |
471 | Number bb = converter->toFloat(b); |
472 | if (aa == bb) |
473 | return true; |
474 | Number x = aa - bb; |
475 | return (x < 0.0 ? -x : x) < ((aa < 0.0 ? -aa : aa) * DBL_EPSILON); |
476 | } |
477 | |
478 | bool ValueCalc::greater(const Value &a, const Value &b) |
479 | { |
480 | Number aa = converter->toFloat(a); |
481 | Number bb = converter->toFloat(b); |
482 | return (aa > bb); |
483 | } |
484 | |
485 | bool ValueCalc::gequal(const Value &a, const Value &b) |
486 | { |
487 | return (greater(a, b) || approxEqual(a, b)); |
488 | } |
489 | |
490 | bool ValueCalc::lower(const Value &a, const Value &b) |
491 | { |
492 | return greater(b, a); |
493 | } |
494 | |
495 | bool ValueCalc::strEqual(const Value &a, const Value &b, bool CalcS) |
496 | { |
497 | QString aa = converter->asString(a).asString(); |
498 | QString bb = converter->asString(b).asString(); |
499 | if (!CalcS) { |
500 | aa = aa.toLower(); |
501 | bb = bb.toLower(); |
502 | } |
503 | return (aa == bb); |
504 | } |
505 | |
506 | bool ValueCalc::strGreater(const Value &a, const Value &b, bool CalcS) |
507 | { |
508 | QString aa = converter->asString(a).asString(); |
509 | QString bb = converter->asString(b).asString(); |
510 | if (!CalcS) { |
511 | aa = aa.toLower(); |
512 | bb = bb.toLower(); |
513 | } |
514 | return (aa > bb); |
515 | } |
516 | |
517 | bool ValueCalc::strGequal(const Value &a, const Value &b, bool CalcS) |
518 | { |
519 | QString aa = converter->asString(a).asString(); |
520 | QString bb = converter->asString(b).asString(); |
521 | if (!CalcS) { |
522 | aa = aa.toLower(); |
523 | bb = bb.toLower(); |
524 | } |
525 | return (aa >= bb); |
526 | } |
527 | |
528 | bool ValueCalc::strLower(const Value &a, const Value &b, bool CalcS) |
529 | { |
530 | return strGreater(b, a, CalcS); |
531 | } |
532 | |
533 | bool ValueCalc::naturalEqual(const Value &a, const Value &b, bool CalcS) |
534 | { |
535 | Value aa = a; |
536 | Value bb = b; |
537 | if (!CalcS) { |
538 | // not case sensitive -> convert strings to lowercase |
539 | if (aa.isString()) aa = Value(aa.asString().toLower()); |
540 | if (bb.isString()) bb = Value(bb.asString().toLower()); |
541 | } |
542 | if (aa.allowComparison(bb)) return aa.equal(bb); |
543 | return strEqual(aa, bb, CalcS); |
544 | } |
545 | |
546 | bool ValueCalc::naturalGreater(const Value &a, const Value &b, bool CalcS) |
547 | { |
548 | Value aa = a; |
549 | Value bb = b; |
550 | if (!CalcS) { |
551 | // not case sensitive -> convert strings to lowercase |
552 | if (aa.isString()) aa = Value(aa.asString().toLower()); |
553 | if (bb.isString()) bb = Value(bb.asString().toLower()); |
554 | } |
555 | if (aa.allowComparison(bb)) return aa.greater(bb); |
556 | return strEqual(aa, bb, CalcS); |
557 | } |
558 | |
559 | bool ValueCalc::naturalGequal(const Value &a, const Value &b, bool CalcS) |
560 | { |
561 | return (naturalGreater(a, b, CalcS) || naturalEqual(a, b, CalcS)); |
562 | } |
563 | |
564 | bool ValueCalc::naturalLower(const Value &a, const Value &b, bool CalcS) |
565 | { |
566 | return naturalGreater(b, a, CalcS); |
567 | } |
568 | |
569 | bool ValueCalc::naturalLequal(const Value &a, const Value &b, bool CalcS) |
570 | { |
571 | return (naturalLower(a, b, CalcS) || naturalEqual(a, b, CalcS)); |
572 | } |
573 | |
574 | Value ValueCalc::roundDown(const Value &a, |
575 | const Value &digits) |
576 | { |
577 | return roundDown(a, converter->asInteger(digits).asInteger()); |
578 | } |
579 | |
580 | Value ValueCalc::roundUp(const Value &a, |
581 | const Value &digits) |
582 | { |
583 | return roundUp(a, converter->asInteger(digits).asInteger()); |
584 | } |
585 | |
586 | Value ValueCalc::round(const Value &a, |
587 | const Value &digits) |
588 | { |
589 | return round(a, converter->asInteger(digits).asInteger()); |
590 | } |
591 | |
592 | Value ValueCalc::roundDown(const Value &a, int digits) |
593 | { |
594 | // shift in one direction, round, shift back |
595 | Value val = a; |
596 | if (digits > 0) |
597 | for (int i = 0; i < digits; ++i) |
598 | val = mul(val, 10); |
599 | if (digits < 0) |
600 | for (int i = 0; i > digits; --i) |
601 | val = div(val, 10); |
602 | |
603 | val = Value(floor(numToDouble(converter->toFloat(val)))); |
604 | |
605 | if (digits > 0) |
606 | for (int i = 0; i < digits; ++i) |
607 | val = div(val, 10); |
608 | if (digits < 0) |
609 | for (int i = 0; i > digits; --i) |
610 | val = mul(val, 10); |
611 | return val; |
612 | } |
613 | |
614 | Value ValueCalc::roundUp(const Value &a, int digits) |
615 | { |
616 | // shift in one direction, round, shift back |
617 | Value val = a; |
618 | if (digits > 0) |
619 | for (int i = 0; i < digits; ++i) |
620 | val = mul(val, 10); |
621 | if (digits < 0) |
622 | for (int i = 0; i > digits; --i) |
623 | val = div(val, 10); |
624 | |
625 | val = Value(ceil(numToDouble(converter->toFloat(val)))); |
626 | |
627 | if (digits > 0) |
628 | for (int i = 0; i < digits; ++i) |
629 | val = div(val, 10); |
630 | if (digits < 0) |
631 | for (int i = 0; i > digits; --i) |
632 | val = mul(val, 10); |
633 | return val; |
634 | } |
635 | |
636 | Value ValueCalc::round(const Value &a, int digits) |
637 | { |
638 | // shift in one direction, round, shift back |
639 | Value val = a; |
640 | if (digits > 0) |
641 | for (int i = 0; i < digits; ++i) |
642 | val = mul(val, 10); |
643 | if (digits < 0) |
644 | for (int i = 0; i > digits; --i) |
645 | val = div(val, 10); |
646 | |
647 | val = Value((lower(val, 0) ? -1 : 1) * qRound(qAbs(converter->toFloat(val)))); |
648 | |
649 | if (digits > 0) |
650 | for (int i = 0; i < digits; ++i) |
651 | val = div(val, 10); |
652 | if (digits < 0) |
653 | for (int i = 0; i > digits; --i) |
654 | val = mul(val, 10); |
655 | return val; |
656 | } |
657 | |
658 | int ValueCalc::sign(const Value &a) |
659 | { |
660 | Number val = converter->toFloat(a); |
661 | if (val == 0) return 0; |
662 | if (val > 0) return 1; |
663 | return -1; |
664 | } |
665 | |
666 | |
667 | Value ValueCalc::log(const Value &number, |
668 | const Value &base) |
669 | { |
670 | Number logbase = converter->toFloat(base); |
671 | if (logbase == 1.0) |
672 | return Value::errorDIV0(); |
673 | if (logbase <= 0.0) |
674 | return Value::errorNA(); |
675 | |
676 | logbase = Calligra::Sheets::log(logbase, 10); |
677 | Value res = Value(Calligra::Sheets::log(converter->toFloat(number), 10) / logbase); |
678 | |
679 | if (number.isNumber() || number.isEmpty()) |
680 | res.setFormat(number.format()); |
681 | |
682 | return res; |
683 | } |
684 | |
685 | Value ValueCalc::ln(const Value &number) |
686 | { |
687 | Value res = Value(::ln(converter->toFloat(number))); |
688 | |
689 | if (number.isNumber() || number.isEmpty()) |
690 | res.setFormat(number.format()); |
691 | |
692 | return res; |
693 | } |
694 | |
695 | Value ValueCalc::log(const Value &number, Number base) |
696 | { |
697 | if (base <= 0.0) |
698 | return Value::errorNA(); |
699 | if (base == 1.0) |
700 | return Value::errorDIV0(); |
701 | |
702 | Number num = converter->toFloat(number); |
703 | Value res = Value(Calligra::Sheets::log(num, base)); |
704 | |
705 | if (number.isNumber() || number.isEmpty()) |
706 | res.setFormat(number.format()); |
707 | |
708 | return res; |
709 | } |
710 | |
711 | Value ValueCalc::exp(const Value &number) |
712 | { |
713 | return Value(::exp(converter->toFloat(number))); |
714 | } |
715 | |
716 | Value ValueCalc::pi() |
717 | { |
718 | // retun PI in double-precision |
719 | // if arbitrary precision gets in, this should be extended to return |
720 | // more if need be |
721 | return Value(M_PI); |
722 | } |
723 | |
724 | Value ValueCalc::eps() |
725 | { |
726 | // #### This should adjust according to the actual number system used |
727 | // (float, double, long double, ...) |
728 | return Value(DBL_EPSILON); |
729 | } |
730 | |
731 | Value ValueCalc::random(Number range) |
732 | { |
733 | return Value(range *(double) rand() / (RAND_MAX + 1.0)); |
734 | } |
735 | |
736 | Value ValueCalc::random(Value range) |
737 | { |
738 | return random(converter->toFloat(range)); |
739 | } |
740 | |
741 | Value ValueCalc::fact(const Value &which) |
742 | { |
743 | // we can simply use integers - no one is going to compute factorial of |
744 | // anything bigger than 2^64 |
745 | return fact(converter->asInteger(which).asInteger()); |
746 | } |
747 | |
748 | Value ValueCalc::fact(const Value &which, |
749 | const Value &end) |
750 | { |
751 | // we can simply use integers - no one is going to compute factorial of |
752 | // anything bigger than 2^64 |
753 | return fact(converter->asInteger(which).asInteger(), |
754 | converter->asInteger(end).asInteger()); |
755 | } |
756 | |
757 | Value ValueCalc::fact(int which, int end) |
758 | { |
759 | if (which < 0) |
760 | return Value(-1); |
761 | if (which == 0) |
762 | return Value(1); |
763 | Value res = Value(1); |
764 | while (which > end) { |
765 | res = mul(res, which); |
766 | --which; |
767 | } |
768 | return res; |
769 | } |
770 | |
771 | Value ValueCalc::factDouble(int which) |
772 | { |
773 | if (which < 0) |
774 | return Value(-1); |
775 | if ((which == 0) || (which == 1)) |
776 | return Value(1); |
777 | |
778 | Value res = Value(1); |
779 | while (which > 1) { |
780 | res = mul(res, which); |
781 | which -= 2; |
782 | } |
783 | return res; |
784 | } |
785 | |
786 | Value ValueCalc::factDouble(Value which) |
787 | { |
788 | return factDouble(converter->asInteger(which).asInteger()); |
789 | } |
790 | |
791 | Value ValueCalc::combin(int n, int k) |
792 | { |
793 | if (n >= 15) { |
794 | Number result = ::exp(Number(lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1))); |
795 | return Value(floor(numToDouble(result + 0.5))); |
796 | } else |
797 | return div(div(fact(n), fact(k)), fact(n - k)); |
798 | } |
799 | |
800 | Value ValueCalc::combin(Value n, Value k) |
801 | { |
802 | int nn = converter->toInteger(n); |
803 | int kk = converter->toInteger(k); |
804 | return combin(nn, kk); |
805 | } |
806 | |
807 | Value ValueCalc::gcd(const Value &a, const Value &b) |
808 | { |
809 | // Euler's GCD algorithm |
810 | Value aa = round(a); |
811 | Value bb = round(b); |
812 | |
813 | if (approxEqual(aa, bb)) return aa; |
814 | |
815 | if (aa.isZero()) return bb; |
816 | if (bb.isZero()) return aa; |
817 | |
818 | |
819 | if (greater(aa, bb)) |
820 | return gcd(bb, mod(aa, bb)); |
821 | else |
822 | return gcd(aa, mod(bb, aa)); |
823 | } |
824 | |
825 | Value ValueCalc::lcm(const Value &a, const Value &b) |
826 | { |
827 | Value aa = round(a); |
828 | Value bb = round(b); |
829 | |
830 | if (approxEqual(aa, bb)) return aa; |
831 | |
832 | if (aa.isZero()) return bb; |
833 | if (bb.isZero()) return aa; |
834 | |
835 | Value g = gcd(aa, bb); |
836 | if (g.isZero()) // GCD is zero for some weird reason |
837 | return mul(aa, bb); |
838 | |
839 | return div(mul(aa, bb), g); |
840 | } |
841 | |
842 | Value ValueCalc::base(const Value &val, int base, int prec, int minLength) |
843 | { |
844 | if (prec < 0) prec = 2; |
845 | if ((base < 2) || (base > 36)) |
846 | return Value::errorVALUE(); |
847 | |
848 | Number value = converter->toFloat(val); |
849 | QString result = QString::number((int)numToDouble(value), base); |
850 | if (result.length() < minLength) |
851 | result = result.rightJustified(minLength, QChar('0')); |
852 | |
853 | if (prec > 0) { |
854 | result += '.'; value = value - (int)numToDouble(value); |
855 | |
856 | int ix; |
857 | for (int i = 0; i < prec; ++i) { |
858 | ix = (int) numToDouble(value * base); |
859 | result += "0123456789abcdefghijklmnopqrstuvwxyz" [ix]; |
860 | value = base * (value - (double)ix / base); |
861 | } |
862 | } |
863 | |
864 | return Value(result.toUpper()); |
865 | } |
866 | |
867 | Value ValueCalc::fromBase(const Value &val, int base) |
868 | { |
869 | QString str = converter->asString(val).asString(); |
870 | bool ok; |
871 | qint64 num = str.toLongLong(&ok, base); |
872 | if (ok) |
873 | return Value(num); |
874 | return Value::errorVALUE(); |
875 | } |
876 | |
877 | |
878 | Value ValueCalc::sin(const Value &number) |
879 | { |
880 | bool ok = true; |
881 | Number n = converter->asFloat(number, &ok).asFloat(); |
882 | if (!ok) |
883 | return Value::errorVALUE(); |
884 | |
885 | Value res = Value(::sin(n)); |
886 | |
887 | if (number.isNumber() || number.isEmpty()) |
888 | res.setFormat(number.format()); |
889 | |
890 | return res; |
891 | } |
892 | |
893 | Value ValueCalc::cos(const Value &number) |
894 | { |
895 | bool ok = true; |
896 | Number n = converter->asFloat(number, &ok).asFloat(); |
897 | if (!ok) |
898 | return Value::errorVALUE(); |
899 | |
900 | Value res = Value(::cos(n)); |
901 | |
902 | if (number.isNumber() || number.isEmpty()) |
903 | res.setFormat(number.format()); |
904 | |
905 | return res; |
906 | } |
907 | |
908 | Value ValueCalc::tg(const Value &number) |
909 | { |
910 | Value res = Value(::tg(converter->toFloat(number))); |
911 | |
912 | if (number.isNumber() || number.isEmpty()) |
913 | res.setFormat(number.format()); |
914 | |
915 | return res; |
916 | } |
917 | |
918 | Value ValueCalc::cotg(const Value &number) |
919 | { |
920 | Value res = div(1, Value(::tg(converter->toFloat(number)))); |
921 | |
922 | if (number.isNumber() || number.isEmpty()) |
923 | res.setFormat(number.format()); |
924 | |
925 | return res; |
926 | } |
927 | |
928 | Value ValueCalc::asin(const Value &number) |
929 | { |
930 | bool ok = true; |
931 | Number n = converter->asFloat(number, &ok).asFloat(); |
932 | if (!ok) |
933 | return Value::errorVALUE(); |
934 | const double d = numToDouble(n); |
935 | if (d < -1.0 || d > 1.0) |
936 | return Value::errorVALUE(); |
937 | |
938 | errno = 0; |
939 | Value res = Value(::asin(n)); |
940 | if (errno) |
941 | return Value::errorVALUE(); |
942 | |
943 | if (number.isNumber() || number.isEmpty()) |
944 | res.setFormat(number.format()); |
945 | return res; |
946 | } |
947 | |
948 | Value ValueCalc::acos(const Value &number) |
949 | { |
950 | bool ok = true; |
951 | Number n = converter->asFloat(number, &ok).asFloat(); |
952 | if (!ok) |
953 | return Value::errorVALUE(); |
954 | const double d = numToDouble(n); |
955 | if (d < -1.0 || d > 1.0) |
956 | return Value::errorVALUE(); |
957 | |
958 | errno = 0; |
959 | Value res = Value(::acos(n)); |
960 | if (errno) |
961 | return Value::errorVALUE(); |
962 | |
963 | if (number.isNumber() || number.isEmpty()) |
964 | res.setFormat(number.format()); |
965 | return res; |
966 | } |
967 | |
968 | Value ValueCalc::atg(const Value &number) |
969 | { |
970 | errno = 0; |
971 | Value res = Value(::atg(converter->toFloat(number))); |
972 | if (errno) |
973 | return Value::errorVALUE(); |
974 | |
975 | if (number.isNumber() || number.isEmpty()) |
976 | res.setFormat(number.format()); |
977 | |
978 | return res; |
979 | } |
980 | |
981 | Value ValueCalc::atan2(const Value &y, const Value &x) |
982 | { |
983 | Number yy = converter->toFloat(y); |
984 | Number xx = converter->toFloat(x); |
985 | return Value(::atan2(yy, xx)); |
986 | } |
987 | |
988 | Value ValueCalc::sinh(const Value &number) |
989 | { |
990 | Value res = Value(::sinh(converter->toFloat(number))); |
991 | |
992 | if (number.isNumber() || number.isEmpty()) |
993 | res.setFormat(number.format()); |
994 | |
995 | return res; |
996 | } |
997 | |
998 | Value ValueCalc::cosh(const Value &number) |
999 | { |
1000 | Value res = Value(::cosh(converter->toFloat(number))); |
1001 | |
1002 | if (number.isNumber() || number.isEmpty()) |
1003 | res.setFormat(number.format()); |
1004 | |
1005 | return res; |
1006 | } |
1007 | |
1008 | Value ValueCalc::tgh(const Value &number) |
1009 | { |
1010 | Value res = Value(::tgh(converter->toFloat(number))); |
1011 | |
1012 | if (number.isNumber() || number.isEmpty()) |
1013 | res.setFormat(number.format()); |
1014 | |
1015 | return res; |
1016 | } |
1017 | |
1018 | Value ValueCalc::asinh(const Value &number) |
1019 | { |
1020 | errno = 0; |
1021 | Value res = Value(::asinh(converter->toFloat(number))); |
1022 | if (errno) |
1023 | return Value::errorVALUE(); |
1024 | |
1025 | if (number.isNumber() || number.isEmpty()) |
1026 | res.setFormat(number.format()); |
1027 | |
1028 | return res; |
1029 | } |
1030 | |
1031 | Value ValueCalc::acosh(const Value &number) |
1032 | { |
1033 | errno = 0; |
1034 | Value res = Value(::acosh(converter->toFloat(number))); |
1035 | if (errno) |
1036 | return Value::errorVALUE(); |
1037 | |
1038 | if (number.isNumber() || number.isEmpty()) |
1039 | res.setFormat(number.format()); |
1040 | |
1041 | return res; |
1042 | } |
1043 | |
1044 | Value ValueCalc::atgh(const Value &number) |
1045 | { |
1046 | errno = 0; |
1047 | Value res = Value(::atgh(converter->toFloat(number))); |
1048 | if (errno) |
1049 | return Value::errorVALUE(); |
1050 | |
1051 | if (number.isNumber() || number.isEmpty()) |
1052 | res.setFormat(number.format()); |
1053 | |
1054 | return res; |
1055 | } |
1056 | |
1057 | Value ValueCalc::phi(Value x) |
1058 | { |
1059 | Value constant(0.39894228040143268); |
1060 | |
1061 | // constant * exp(-(x * x) / 2.0); |
1062 | Value x2neg = mul(sqr(x), -1); |
1063 | return mul(constant, exp(div(x2neg, 2.0))); |
1064 | } |
1065 | |
1066 | static double taylor_helper(double* pPolynom, uint nMax, double x) |
1067 | { |
1068 | double nVal = pPolynom[nMax]; |
1069 | for (int i = nMax - 1; i >= 0; --i) { |
1070 | nVal = pPolynom[i] + (nVal * x); |
1071 | } |
1072 | return nVal; |
1073 | } |
1074 | |
1075 | Value ValueCalc::gauss(Value xx) |
1076 | // this is a weird function |
1077 | { |
1078 | double x = converter->toFloat(xx); |
1079 | |
1080 | double t0[] = { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582, |
1081 | -0.00118732821548045, 0.00011543468761616, -0.00000944465625950, |
1082 | 0.00000066596935163, -0.00000004122667415, 0.00000000227352982, |
1083 | 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 |
1084 | }; |
1085 | double t2[] = { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805, |
1086 | 0.02699548325659403, -0.00449924720943234, -0.00224962360471617, |
1087 | 0.00134977416282970, -0.00011783742691370, -0.00011515930357476, |
1088 | 0.00003704737285544, 0.00000282690796889, -0.00000354513195524, |
1089 | 0.00000037669563126, 0.00000019202407921, -0.00000005226908590, |
1090 | -0.00000000491799345, 0.00000000366377919, -0.00000000015981997, |
1091 | -0.00000000017381238, 0.00000000002624031, 0.00000000000560919, |
1092 | -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 |
1093 | }; |
1094 | double t4[] = { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977, |
1095 | 0.00033457556441221, -0.00028996548915725, 0.00018178605666397, |
1096 | -0.00008252863922168, 0.00002551802519049, -0.00000391665839292, |
1097 | -0.00000074018205222, 0.00000064422023359, -0.00000017370155340, |
1098 | 0.00000000909595465, 0.00000000944943118, -0.00000000329957075, |
1099 | 0.00000000029492075, 0.00000000011874477, -0.00000000004420396, |
1100 | 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 |
1101 | }; |
1102 | double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 }; |
1103 | |
1104 | double xAbs = fabs(x); |
1105 | uint xShort = static_cast<uint>(approxFloor(xAbs)); // approxFloor taken from OOo |
1106 | double nVal = 0.0; |
1107 | if (xShort == 0) |
1108 | nVal = taylor_helper(t0, 11, (xAbs * xAbs)) * xAbs; |
1109 | else if ((xShort >= 1) && (xShort <= 2)) |
1110 | nVal = taylor_helper(t2, 23, (xAbs - 2.0)); |
1111 | else if ((xShort >= 3) && (xShort <= 4)) |
1112 | nVal = taylor_helper(t4, 20, (xAbs - 4.0)); |
1113 | else { |
1114 | double phiAbs = numToDouble(converter->toFloat(phi(Value(xAbs)))); |
1115 | nVal = 0.5 + phiAbs * taylor_helper(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs; |
1116 | } |
1117 | |
1118 | if (x < 0.0) |
1119 | return Value(-nVal); |
1120 | else |
1121 | return Value(nVal); |
1122 | } |
1123 | |
1124 | Value ValueCalc::gaussinv(Value xx) |
1125 | // this is a weird function |
1126 | { |
1127 | double x = numToDouble(converter->toFloat(xx)); |
1128 | |
1129 | double q, t, z; |
1130 | |
1131 | q = x - 0.5; |
1132 | |
1133 | if (fabs(q) <= .425) { |
1134 | t = 0.180625 - q * q; |
1135 | |
1136 | z = |
1137 | q * |
1138 | ( |
1139 | ( |
1140 | ( |
1141 | ( |
1142 | ( |
1143 | ( |
1144 | ( |
1145 | t * 2509.0809287301226727 + 33430.575583588128105 |
1146 | ) |
1147 | * t + 67265.770927008700853 |
1148 | ) |
1149 | * t + 45921.953931549871457 |
1150 | ) |
1151 | * t + 13731.693765509461125 |
1152 | ) |
1153 | * t + 1971.5909503065514427 |
1154 | ) |
1155 | * t + 133.14166789178437745 |
1156 | ) |
1157 | * t + 3.387132872796366608 |
1158 | ) |
1159 | / |
1160 | ( |
1161 | ( |
1162 | ( |
1163 | ( |
1164 | ( |
1165 | ( |
1166 | ( |
1167 | t * 5226.495278852854561 + 28729.085735721942674 |
1168 | ) |
1169 | * t + 39307.89580009271061 |
1170 | ) |
1171 | * t + 21213.794301586595867 |
1172 | ) |
1173 | * t + 5394.1960214247511077 |
1174 | ) |
1175 | * t + 687.1870074920579083 |
1176 | ) |
1177 | * t + 42.313330701600911252 |
1178 | ) |
1179 | * t + 1.0 |
1180 | ); |
1181 | |
1182 | } else { |
1183 | if (q > 0) t = 1 - x; |
1184 | else t = x; |
1185 | |
1186 | t =::sqrt(-::log(t)); |
1187 | |
1188 | if (t <= 5.0) { |
1189 | t += -1.6; |
1190 | |
1191 | z = |
1192 | ( |
1193 | ( |
1194 | ( |
1195 | ( |
1196 | ( |
1197 | ( |
1198 | ( |
1199 | t * 7.7454501427834140764e-4 + 0.0227238449892691845833 |
1200 | ) |
1201 | * t + 0.24178072517745061177 |
1202 | ) |
1203 | * t + 1.27045825245236838258 |
1204 | ) |
1205 | * t + 3.64784832476320460504 |
1206 | ) |
1207 | * t + 5.7694972214606914055 |
1208 | ) |
1209 | * t + 4.6303378461565452959 |
1210 | ) |
1211 | * t + 1.42343711074968357734 |
1212 | ) |
1213 | / |
1214 | ( |
1215 | ( |
1216 | ( |
1217 | ( |
1218 | ( |
1219 | ( |
1220 | ( |
1221 | t * 1.05075007164441684324e-9 + 5.475938084995344946e-4 |
1222 | ) |
1223 | * t + 0.0151986665636164571966 |
1224 | ) |
1225 | * t + 0.14810397642748007459 |
1226 | ) |
1227 | * t + 0.68976733498510000455 |
1228 | ) |
1229 | * t + 1.6763848301838038494 |
1230 | ) |
1231 | * t + 2.05319162663775882187 |
1232 | ) |
1233 | * t + 1.0 |
1234 | ); |
1235 | |
1236 | } else { |
1237 | t += -5.0; |
1238 | |
1239 | z = |
1240 | ( |
1241 | ( |
1242 | ( |
1243 | ( |
1244 | ( |
1245 | ( |
1246 | ( |
1247 | t * 2.01033439929228813265e-7 + 2.71155556874348757815e-5 |
1248 | ) |
1249 | * t + 0.0012426609473880784386 |
1250 | ) |
1251 | * t + 0.026532189526576123093 |
1252 | ) |
1253 | * t + 0.29656057182850489123 |
1254 | ) |
1255 | * t + 1.7848265399172913358 |
1256 | ) |
1257 | * t + 5.4637849111641143699 |
1258 | ) |
1259 | * t + 6.6579046435011037772 |
1260 | ) |
1261 | / |
1262 | ( |
1263 | ( |
1264 | ( |
1265 | ( |
1266 | ( |
1267 | ( |
1268 | ( |
1269 | t * 2.04426310338993978564e-15 + 1.4215117583164458887e-7 |
1270 | ) |
1271 | * t + 1.8463183175100546818e-5 |
1272 | ) |
1273 | * t + 7.868691311456132591e-4 |
1274 | ) |
1275 | * t + 0.0148753612908506148525 |
1276 | ) |
1277 | * t + 0.13692988092273580531 |
1278 | ) |
1279 | * t + 0.59983220655588793769 |
1280 | ) |
1281 | * t + 1.0 |
1282 | ); |
1283 | |
1284 | } |
1285 | |
1286 | if (q < 0.0) z = -z; |
1287 | } |
1288 | |
1289 | return Value(z); |
1290 | } |
1291 | |
1292 | // |
1293 | // GetGamma |
1294 | // |
1295 | Value ValueCalc::GetGamma(Value value) |
1296 | { |
1297 | double val = conv()->asFloat(value).asFloat(); |
1298 | |
1299 | bool reflect; |
1300 | |
1301 | double gamma = GammaHelp(val, reflect); |
1302 | |
1303 | gamma = ::pow(val + 5.5, val + 0.5) * gamma /::exp(val + 5.5); |
1304 | |
1305 | if (reflect) |
1306 | gamma = M_PI * val / (gamma*::sin(M_PI * val)); |
1307 | |
1308 | return Value(gamma); |
1309 | } |
1310 | |
1311 | Value ValueCalc::GetLogGamma(Value _x) |
1312 | { |
1313 | double x = numToDouble(converter->toFloat(_x)); |
1314 | |
1315 | bool bReflect; |
1316 | double G = GammaHelp(x, bReflect); |
1317 | G = (x + 0.5)*::log(x + 5.5) +::log(G) - (x + 5.5); |
1318 | if (bReflect) |
1319 | G = ::log(M_PI * x) - G -::log(::sin(M_PI * x)); |
1320 | return Value(G); |
1321 | } |
1322 | |
1323 | // Value ValueCalc::GetGammaDist (Value _x, Value _alpha, |
1324 | // Value _beta) |
1325 | // { |
1326 | // double x = numToDouble (converter->toFloat (_x)); |
1327 | // double alpha = numToDouble (converter->toFloat (_alpha)); |
1328 | // double beta = numToDouble (converter->toFloat (_beta)); |
1329 | // |
1330 | // if (x == 0.0) |
1331 | // return Value (0.0); |
1332 | // |
1333 | // x /= beta; |
1334 | // double gamma = alpha; |
1335 | // |
1336 | // double c = 0.918938533204672741; |
1337 | // double d[10] = { |
1338 | // 0.833333333333333333E-1, |
1339 | // -0.277777777777777778E-2, |
1340 | // 0.793650793650793651E-3, |
1341 | // -0.595238095238095238E-3, |
1342 | // 0.841750841750841751E-3, |
1343 | // -0.191752691752691753E-2, |
1344 | // 0.641025641025641025E-2, |
1345 | // -0.295506535947712418E-1, |
1346 | // 0.179644372368830573, |
1347 | // -0.139243221690590111E1 |
1348 | // }; |
1349 | // |
1350 | // double dx = x; |
1351 | // double dgamma = gamma; |
1352 | // int maxit = 10000; |
1353 | // |
1354 | // double z = dgamma; |
1355 | // double den = 1.0; |
1356 | // while ( z < 10.0 ) { |
1357 | // den *= z; |
1358 | // z += 1.0; |
1359 | // } |
1360 | // |
1361 | // double z2 = z*z; |
1362 | // double z3 = z*z2; |
1363 | // double z4 = z2*z2; |
1364 | // double z5 = z2*z3; |
1365 | // double a = ( z - 0.5 ) * ::log10(z) - z + c; |
1366 | // double b = d[0]/z + d[1]/z3 + d[2]/z5 + d[3]/(z2*z5) + d[4]/(z4*z5) + |
1367 | // d[5]/(z*z5*z5) + d[6]/(z3*z5*z5) + d[7]/(z5*z5*z5) + d[8]/(z2*z5*z5*z5); |
1368 | // // double g = exp(a+b) / den; |
1369 | // |
1370 | // double sum = 1.0 / dgamma; |
1371 | // double term = 1.0 / dgamma; |
1372 | // double cut1 = dx - dgamma; |
1373 | // double cut2 = dx * 10000000000.0; |
1374 | // |
1375 | // for ( int i=1; i<=maxit; ++i ) { |
1376 | // double ai = i; |
1377 | // term = dx * term / ( dgamma + ai ); |
1378 | // sum += term; |
1379 | // double cutoff = cut1 + ( cut2 * term / sum ); |
1380 | // if ( ai > cutoff ) { |
1381 | // double t = sum; |
1382 | // // return pow( dx, dgamma ) * exp( -dx ) * t / g; |
1383 | // return Value (::exp( dgamma * ::log(dx) - dx - a - b ) * t * den); |
1384 | // } |
1385 | // } |
1386 | // |
1387 | // return Value (1.0); // should not happen ... |
1388 | // } |
1389 | |
1390 | Value ValueCalc::GetGammaDist(Value _x, Value _alpha, Value _beta) |
1391 | { |
1392 | // Algorithm adopted from StatLib (http://lib.stat.cmu.edu/apstat/239): |
1393 | // Algorithm AS 239: Chi-Squared and Incomplete Gamma Integral |
1394 | // B. L. Shea |
1395 | // Applied Statistics, Vol. 37, No. 3 (1988), pp. 466-473 |
1396 | |
1397 | // TODO - normal approx |
1398 | |
1399 | double x = numToDouble(converter->toFloat(_x)); // x |
1400 | double alpha = numToDouble(converter->toFloat(_alpha)); // alpha |
1401 | double beta = numToDouble(converter->toFloat(_beta)); // beta |
1402 | |
1403 | // debug info |
1404 | // kDebug()<<"GetGammaDist( x="<<x<<", alpha="<<alpha<<", beta="<<beta<<" )"; |
1405 | |
1406 | int lower_tail = 1; // |
1407 | int pearson; // flag is set if pearson was used |
1408 | |
1409 | const static double xlarge = 1.0e+37; |
1410 | |
1411 | double res; // result |
1412 | double pn1, pn2, pn3, pn4, pn5, pn6, a, b, c, an, osum, sum; |
1413 | long n; |
1414 | |
1415 | x /= beta; |
1416 | // kDebug()<<"-> x=x/beta ="<<x; |
1417 | |
1418 | // check constraints |
1419 | if (x <= 0.0) |
1420 | return Value(0.0); |
1421 | |
1422 | if (x <= 1.0 || x < alpha) { |
1423 | // |
1424 | // Result = e^log(alpha * log(x) - x - log( gamma( alpha+1 ) ) ) |
1425 | // |
1426 | |
1427 | pearson = 1; // set flag -> use pearson's series expansion. |
1428 | res = alpha * ::log(x) - x - ::log(GetGamma(Value(alpha + 1.0)).asFloat()); |
1429 | // kDebug()<<"Pearson res="<<res; |
1430 | |
1431 | // x x x |
1432 | // sum = 1.0 + --------- + --------- * --------- + ... |
1433 | // alpha+1 alpha+1 alpha+2 |
1434 | c = 1.0; |
1435 | sum = 1.0; |
1436 | a = alpha; |
1437 | do { |
1438 | a += 1.0; |
1439 | c *= x / a; |
1440 | sum += c; |
1441 | } while (c > DBL_EPSILON * sum); |
1442 | } else { |
1443 | // x >= max( 1, alpha) |
1444 | pearson = 0; // clear flag -> use a continued fraction expansion |
1445 | |
1446 | // TODO use GetLogGamma? |
1447 | res = alpha * ::log(x) - x - ::log(GetGamma(Value(alpha)).asFloat()); |
1448 | |
1449 | // kDebug()<<"Continued fraction expression res="<<res; |
1450 | |
1451 | // |
1452 | // |
1453 | // |
1454 | |
1455 | a = 1.0 - alpha; |
1456 | b = a + x + 1.0; |
1457 | pn1 = 1.0; |
1458 | pn2 = x; |
1459 | pn3 = x + 1.0; |
1460 | pn4 = x * b; |
1461 | sum = pn3 / pn4; |
1462 | for (n = 1; ; ++n) { |
1463 | // kDebug()<<"n="<<n<<" sum="<< sum; |
1464 | a += 1.0; // = n+1 -alpha |
1465 | b += 2.0; // = 2(n+1)-alph+x |
1466 | an = a * n; |
1467 | pn5 = b * pn3 - an * pn1; |
1468 | pn6 = b * pn4 - an * pn2; |
1469 | // kDebug()<<"a ="<<a<<" an="<<an<<" b="<<b<<" pn5="<<pn5<<" pn6="<<pn6; |
1470 | if (fabs(pn6) > 0.0) { |
1471 | osum = sum; |
1472 | sum = pn5 / pn6; |
1473 | // kDebug()<<"sum ="<<sum<<" osum="<<osum; |
1474 | if (fabs(osum - sum) <= DBL_EPSILON * fmin(1.0, sum)) |
1475 | break; |
1476 | } |
1477 | pn1 = pn3; |
1478 | pn2 = pn4; |
1479 | pn3 = pn5; |
1480 | pn4 = pn6; |
1481 | if (fabs(pn5) >= xlarge) { |
1482 | // re-scale the terms in continued fraction if they are large |
1483 | kDebug() << "the terms are to large -> rescaleling by " << xlarge; |
1484 | pn1 /= xlarge; |
1485 | pn2 /= xlarge; |
1486 | pn3 /= xlarge; |
1487 | pn4 /= xlarge; |
1488 | } |
1489 | } |
1490 | } |
1491 | |
1492 | res += ::log(sum); |
1493 | lower_tail = (lower_tail == pearson); |
1494 | |
1495 | if (lower_tail) |
1496 | return Value(::exp(res)); |
1497 | else |
1498 | return Value(-1 * ::expm1(res)); |
1499 | } |
1500 | |
1501 | Value ValueCalc::GetBeta(Value _x, Value _alpha, |
1502 | Value _beta) |
1503 | { |
1504 | // kDebug()<<"GetBeta: x= " << _x << " alpha= " << _alpha << " beta=" << _beta; |
1505 | if (equal(_beta, Value(1.0))) |
1506 | return pow(_x, _alpha); |
1507 | else if (equal(_alpha, Value(1.0))) |
1508 | // 1.0 - pow (1.0-_x, _beta) |
1509 | return sub(Value(1.0), pow(sub(Value(1.0), _x), _beta)); |
1510 | |
1511 | double x = numToDouble(converter->toFloat(_x)); |
1512 | double alpha = numToDouble(converter->toFloat(_alpha)); |
1513 | double beta = numToDouble(converter->toFloat(_beta)); |
1514 | |
1515 | double fEps = 1.0E-8; |
1516 | bool bReflect; |
1517 | double cf, fA, fB; |
1518 | |
1519 | if (x < (alpha + 1.0) / (alpha + beta + 1.0)) { |
1520 | bReflect = false; |
1521 | fA = alpha; |
1522 | fB = beta; |
1523 | } else { |
1524 | bReflect = true; |
1525 | fA = beta; |
1526 | fB = alpha; |
1527 | x = 1.0 - x; |
1528 | } |
1529 | if (x < fEps) |
1530 | cf = 0.0; |
1531 | else { |
1532 | double a1, b1, a2, b2, fnorm, rm, apl2m, d2m, d2m1, cfnew; |
1533 | a1 = 1.0; b1 = 1.0; |
1534 | b2 = 1.0 - (fA + fB) * x / (fA + 1.0); |
1535 | if (b2 == 0.0) { |
1536 | a2 = b2; |
1537 | fnorm = 1.0; |
1538 | cf = 1.0; |
1539 | } else { |
1540 | a2 = 1.0; |
1541 | fnorm = 1.0 / b2; |
1542 | cf = a2 * fnorm; |
1543 | } |
1544 | cfnew = 1.0; |
1545 | for (uint j = 1; j <= 100; ++j) { |
1546 | rm = (double) j; |
1547 | apl2m = fA + 2.0 * rm; |
1548 | d2m = rm * (fB - rm) * x / ((apl2m - 1.0) * apl2m); |
1549 | d2m1 = -(fA + rm) * (fA + fB + rm) * x / (apl2m * (apl2m + 1.0)); |
1550 | a1 = (a2 + d2m * a1) * fnorm; |
1551 | b1 = (b2 + d2m * b1) * fnorm; |
1552 | a2 = a1 + d2m1 * a2 * fnorm; |
1553 | b2 = b1 + d2m1 * b2 * fnorm; |
1554 | if (b2 != 0.0) { |
1555 | fnorm = 1.0 / b2; |
1556 | cfnew = a2 * fnorm; |
1557 | if (fabs(cf - cfnew) / cf < fEps) |
1558 | j = 101; |
1559 | else |
1560 | cf = cfnew; |
1561 | } |
1562 | } |
1563 | if (fB < fEps) |
1564 | b1 = 1.0E30; |
1565 | else |
1566 | b1 = ::exp(numToDouble(GetLogGamma(Value(fA)).asFloat() + GetLogGamma(Value(fB)).asFloat() - |
1567 | GetLogGamma(Value(fA + fB)).asFloat())); |
1568 | |
1569 | cf *= ::pow(x, fA)*::pow(1.0 - x, fB) / (fA * b1); |
1570 | } |
1571 | if (bReflect) |
1572 | return Value(1.0 - cf); |
1573 | else |
1574 | return Value(cf); |
1575 | } |
1576 | |
1577 | // ------------------------------------------------------ |
1578 | |
1579 | |
1580 | /* |
1581 | * |
1582 | * The code for calculating Bessel functions is taken |
1583 | * from CCMATH, a mathematics library source.code. |
1584 | * |
1585 | * Original copyright follows: |
1586 | * |
1587 | * Copyright (C) 2000 Daniel A. Atkinson All rights reserved. |
1588 | * This code may be redistributed under the terms of the GNU library |
1589 | * public license (LGPL). |
1590 | */ |
1591 | |
1592 | static double ccmath_gaml(double x) |
1593 | { |
1594 | double g, h = 0; /* NB must be called with 0<=x<29 */ |
1595 | for (g = 1.; x < 30. ; g *= x, x += 1.) h = x * x; |
1596 | g = (x - .5) * log(x) - x + .918938533204672 - log(g); |
1597 | g += (1. - (1. / 6. - (1. / 3. - 1. / (4.*h)) / (7.*h)) / (5.*h)) / (12.*x); |
1598 | return g; |
1599 | } |
1600 | |
1601 | static double ccmath_psi(int m) |
1602 | { |
1603 | double s = -.577215664901533; int k; |
1604 | for (k = 1; k < m ; ++k) s += 1. / k; |
1605 | return s; |
1606 | } |
1607 | |
1608 | static double ccmath_ibes(double v, double x) |
1609 | { |
1610 | double y, s = 0., t = 0., tp; int p, m; |
1611 | y = x - 9.; if (y > 0.) y *= y; tp = v * v * .2 + 25.; |
1612 | if (y < tp) { |
1613 | x /= 2.; m = (int)x; |
1614 | if (x > 0.) s = t = exp(v * log(x) - ccmath_gaml(v + 1.)); |
1615 | else { |
1616 | if (v > 0.) return 0.; else if (v == 0.) return 1.; |
1617 | } |
1618 | for (p = 1, x *= x;; ++p) { |
1619 | t *= x / (p * (v += 1.)); s += t; |
1620 | if (p > m && t < 1.e-13*s) break; |
1621 | } |
1622 | } else { |
1623 | double u, a0 = 1.57079632679490; |
1624 | s = t = 1. / sqrt(x * a0); x *= 2.; u = 0.; |
1625 | for (p = 1, y = .5; (tp = fabs(t)) > 1.e-14 ; ++p, y += 1.) { |
1626 | t *= (v + y) * (v - y) / (p * x); if (y > v && fabs(t) >= tp) break; |
1627 | if (!(p&1)) s += t; else u -= t; |
1628 | } |
1629 | x /= 2.; s = cosh(x) * s + sinh(x) * u; |
1630 | } |
1631 | return s; |
1632 | } |
1633 | |
1634 | static double ccmath_kbes(double v, double x) |
1635 | { |
1636 | double y, s, t, tp, f, a0 = 1.57079632679490; |
1637 | int p, k, m; |
1638 | if (x == 0.) return HUGE_VAL; |
1639 | y = x - 10.5; if (y > 0.) y *= y; tp = 25. + .185 * v * v; |
1640 | if (y < tp && modf(v + .5, &t) != 0.) { |
1641 | y = 1.5 + .5 * v; |
1642 | if (x < y) { |
1643 | x /= 2.; m = (int)x; tp = t = exp(v * log(x) - ccmath_gaml(v + 1.)); |
1644 | if (modf(v, &y) == 0.) { |
1645 | k = (int)y; tp *= v; |
1646 | f = 2.*log(x) - ccmath_psi(1) - ccmath_psi(k + 1); |
1647 | t /= 2.; if (!(k&1)) t = -t; s = f * t; |
1648 | for (p = 1, x *= x;; ++p) { |
1649 | f -= 1. / p + 1. / (v += 1.); |
1650 | t *= x / (p * v); s += (y = t * f); |
1651 | if (p > m && fabs(y) < 1.e-14) break; |
1652 | } |
1653 | if (k > 0) { |
1654 | x = -x; s += (t = 1. / (tp * 2.)); |
1655 | for (p = 1, --k; k > 0 ; ++p, --k) s += (t *= x / (p * k)); |
1656 | } |
1657 | } else { |
1658 | f = 1. / (t * v * 2.); t *= a0 / sin(2.*a0 * v); s = f - t; |
1659 | for (p = 1, x *= x, tp = v;; ++p) { |
1660 | t *= x / (p * (v += 1.)); f *= -x / (p * (tp -= 1.)); |
1661 | s += (y = f - t); if (p > m && fabs(y) < 1.e-14) break; |
1662 | } |
1663 | } |
1664 | } else { |
1665 | double tq, h, w, z, r; |
1666 | t = 12. / pow(x, .333); k = (int)(t * t); y = 2.*(x + k); |
1667 | m = (int)v; v -= m; tp = v * v - .25; v += 1.; tq = v * v - .25; |
1668 | for (s = h = 1., r = f = z = w = 0.; k > 0 ; --k, y -= 2.) { |
1669 | t = (y * h - (k + 1) * z) / (k - 1 - tp / k); z = h; f += (h = t); |
1670 | t = (y * s - (k + 1) * w) / (k - 1 - tq / k); w = s; r += (s = t); |
1671 | } |
1672 | t = sqrt(a0 / x) * exp(-x); s *= t / r; h *= t / f; x /= 2.; if (m == 0) s = h; |
1673 | for (k = 1; k < m ; ++k) { |
1674 | t = v * s / x + h; h = s; s = t; v += 1.; |
1675 | } |
1676 | } |
1677 | } else { |
1678 | s = t = sqrt(a0 / x); x *= 2.; |
1679 | for (p = 1, y = .5; (tp = fabs(t)) > 1.e-14 ; ++p, y += 1.) { |
1680 | t *= (v + y) * (v - y) / (p * x); if (y > v && fabs(t) >= tp) break; s += t; |
1681 | } |
1682 | s *= exp(-x / 2.); |
1683 | } |
1684 | return s; |
1685 | } |
1686 | |
1687 | static double ccmath_jbes(double v, double x) |
1688 | { |
1689 | double y, s = 0., t = 0., tp; int p, m; |
1690 | y = x - 8.5; if (y > 0.) y *= y; tp = v * v / 4. + 13.69; |
1691 | if (y < tp) { |
1692 | x /= 2.; m = (int)x; |
1693 | if (x > 0.) s = t = exp(v * log(x) - ccmath_gaml(v + 1.)); |
1694 | else { |
1695 | if (v > 0.) return 0.; else if (v == 0.) return 1.; |
1696 | } |
1697 | for (p = 1, x *= -x;; ++p) { |
1698 | t *= x / (p * (v += 1.)); s += t; |
1699 | if (p > m && fabs(t) < 1.e-13) break; |
1700 | } |
1701 | } else { |
1702 | double u, a0 = 1.57079632679490; |
1703 | s = t = 1. / sqrt(x * a0); x *= 2.; u = 0.; |
1704 | for (p = 1, y = .5; (tp = fabs(t)) > 1.e-14 ; ++p, y += 1.) { |
1705 | t *= (v + y) * (v - y) / (p * x); if (y > v && fabs(t) >= tp) break; |
1706 | if (!(p&1)) { |
1707 | t = -t; s += t; |
1708 | } else u -= t; |
1709 | } |
1710 | y = x / 2. - (v + .5) * a0; s = cos(y) * s + sin(y) * u; |
1711 | } |
1712 | return s; |
1713 | } |
1714 | |
1715 | static double ccmath_nbes(double v, double x) |
1716 | { |
1717 | double y, s, t, tp, u, f, a0 = 3.14159265358979; |
1718 | int p, k, m; |
1719 | y = x - 8.5; if (y > 0.) y *= y; tp = v * v / 4. + 13.69; |
1720 | if (y < tp) { |
1721 | if (x == 0.) return HUGE_VAL; |
1722 | x /= 2.; m = (int)x; u = t = exp(v * log(x) - ccmath_gaml(v + 1.)); |
1723 | if (modf(v, &y) == 0.) { |
1724 | k = (int)y; u *= v; |
1725 | f = 2.*log(x) - ccmath_psi(1) - ccmath_psi(k + 1); |
1726 | t /= a0; x *= -x; s = f * t; |
1727 | for (p = 1;; ++p) { |
1728 | f -= 1. / p + 1. / (v += 1.); |
1729 | t *= x / (p * v); s += (y = t * f); if (p > m && fabs(y) < 1.e-13) break; |
1730 | } |
1731 | if (k > 0) { |
1732 | x = -x; s -= (t = 1. / (u * a0)); |
1733 | for (p = 1, --k; k > 0 ; ++p, --k) s -= (t *= x / (p * k)); |
1734 | } |
1735 | } else { |
1736 | f = 1. / (t * v * a0); t /= tan(a0 * v); s = t - f; |
1737 | for (p = 1, x *= x, u = v;; ++p) { |
1738 | t *= -x / (p * (v += 1.)); f *= x / (p * (u -= 1.)); |
1739 | s += (y = t - f); if (p > m && fabs(y) < 1.e-13) break; |
1740 | } |
1741 | } |
1742 | } else { |
1743 | x *= 2.; s = t = 2. / sqrt(x * a0); u = 0.; |
1744 | for (p = 1, y = .5; (tp = fabs(t)) > 1.e-14 ; ++p, y += 1.) { |
1745 | t *= (v + y) * (v - y) / (p * x); if (y > v && fabs(t) > tp) break; |
1746 | if (!(p&1)) { |
1747 | t = -t; s += t; |
1748 | } else u += t; |
1749 | } |
1750 | y = (x - (v + .5) * a0) / 2.; s = sin(y) * s + cos(y) * u; |
1751 | } |
1752 | return s; |
1753 | } |
1754 | |
1755 | |
1756 | /* ---------- end of CCMATH code ---------- */ |
1757 | |
1758 | template <typename func_ptr> |
1759 | Value CalcBessel( |
1760 | func_ptr *func, |
1761 | ValueConverter *converter, |
1762 | Value v, Value x) |
1763 | { |
1764 | double vv = numToDouble(converter->toFloat(v)); |
1765 | double xx = numToDouble(converter->toFloat(x)); |
1766 | // vv must be a non-negative integer and <29 for implementation reasons |
1767 | // xx must be non-negative |
1768 | if (xx >= 0 && vv >= 0 && vv < 29 && vv == floor(vv)) |
1769 | return Value(func(vv, xx)); |
1770 | return Value::errorVALUE(); |
1771 | } |
1772 | |
1773 | Value ValueCalc::besseli(Value v, Value x) |
1774 | { |
1775 | return CalcBessel(ccmath_ibes, converter, v, x); |
1776 | } |
1777 | |
1778 | Value ValueCalc::besselj(Value v, Value x) |
1779 | { |
1780 | return CalcBessel(ccmath_jbes, converter, v, x); |
1781 | } |
1782 | |
1783 | Value ValueCalc::besselk(Value v, Value x) |
1784 | { |
1785 | return CalcBessel(ccmath_kbes, converter, v, x); |
1786 | } |
1787 | |
1788 | Value ValueCalc::besseln(Value v, Value x) |
1789 | { |
1790 | return CalcBessel(ccmath_nbes, converter, v, x); |
1791 | } |
1792 | |
1793 | |
1794 | // ------------------------------------------------------ |
1795 | |
1796 | Value ValueCalc::erf(Value x) |
1797 | { |
1798 | return Value(::erf(numToDouble(converter->toFloat(x)))); |
1799 | } |
1800 | |
1801 | Value ValueCalc::erfc(Value x) |
1802 | { |
1803 | return Value(::erfc(numToDouble(converter->toFloat(x)))); |
1804 | } |
1805 | |
1806 | // ------------------------------------------------------ |
1807 | |
1808 | void ValueCalc::arrayWalk(const Value &range, |
1809 | Value &res, arrayWalkFunc func, Value param) |
1810 | { |
1811 | if (res.isError()) return; |
1812 | if (!range.isArray()) { |
1813 | func(this, res, range, param); |
1814 | return; |
1815 | } |
1816 | |
1817 | // iterate over the non-empty entries |
1818 | for (uint i = 0; i < range.count(); ++i) { |
1819 | Value v = range.element(i); |
1820 | if (v.isArray()) |
1821 | arrayWalk(v, res, func, param); |
1822 | else { |
1823 | func(this, res, v, param); |
1824 | if (res.format() == Value::fmt_None) |
1825 | res.setFormat(v.format()); |
1826 | } |
1827 | } |
1828 | } |
1829 | |
1830 | void ValueCalc::arrayWalk(QVector<Value> &range, |
1831 | Value &res, arrayWalkFunc func, Value param) |
1832 | { |
1833 | if (res.isError()) return; |
1834 | for (int i = 0; i < range.count(); ++i) |
1835 | arrayWalk(range[i], res, func, param); |
1836 | } |
1837 | |
1838 | Value ValueCalc::arrayMap(const Value &array, arrayMapFunc func, const Value ¶m) |
1839 | { |
1840 | Value res( Value::Array ); |
1841 | for (unsigned row = 0; row < array.rows(); ++row) { |
1842 | for (unsigned col = 0; col < array.columns(); ++col) { |
1843 | Value element = array.element( col, row ); |
1844 | Value _res = (this->*func)( element, param ); |
1845 | res.setElement( col, row, _res ); |
1846 | } |
1847 | } |
1848 | return res; |
1849 | } |
1850 | |
1851 | Value ValueCalc::twoArrayMap(const Value &array1, arrayMapFunc func, const Value &array2) |
1852 | { |
1853 | Value res( Value::Array ); |
1854 | // Map each element in one array with the respective element in the other array |
1855 | unsigned rows = qMax(array1.rows(), array2.rows()); |
1856 | unsigned columns = qMax(array1.columns(), array2.columns()); |
1857 | for (unsigned row = 0; row < rows; ++row) { |
1858 | for (unsigned col = 0; col < columns; ++col) { |
1859 | // Value::element() will return an empty value if element(col, row) does not exist. |
1860 | Value element1 = array1.element( col, row ); |
1861 | Value element2 = array2.element( col, row ); |
1862 | Value _res = (this->*func)( element1, element2 ); |
1863 | res.setElement( col, row, _res ); |
1864 | } |
1865 | } |
1866 | return res; |
1867 | } |
1868 | |
1869 | void ValueCalc::twoArrayWalk(const Value &a1, const Value &a2, |
1870 | Value &res, arrayWalkFunc func) |
1871 | { |
1872 | if (res.isError()) return; |
1873 | if (!a1.isArray()) { |
1874 | func(this, res, a1, a2); |
1875 | return; |
1876 | } |
1877 | |
1878 | int rows = a1.rows(); |
1879 | int cols = a1.columns(); |
1880 | int rows2 = a2.rows(); |
1881 | int cols2 = a2.columns(); |
1882 | if ((rows != rows2) || (cols != cols2)) { |
1883 | res = Value::errorVALUE(); |
1884 | return; |
1885 | } |
1886 | for (int r = 0; r < rows; ++r) |
1887 | for (int c = 0; c < cols; ++c) { |
1888 | Value v1 = a1.element(c, r); |
1889 | Value v2 = a2.element(c, r); |
1890 | if (v1.isArray() && v2.isArray()) |
1891 | twoArrayWalk(v1, v2, res, func); |
1892 | else { |
1893 | func(this, res, v1, v2); |
1894 | if (res.format() == Value::fmt_None) |
1895 | res.setFormat(format(v1, v2)); |
1896 | } |
1897 | } |
1898 | } |
1899 | |
1900 | void ValueCalc::twoArrayWalk(QVector<Value> &a1, |
1901 | QVector<Value> &a2, Value &res, arrayWalkFunc func) |
1902 | { |
1903 | if (res.isError()) return; |
1904 | if (a1.count() != a2.count()) { |
1905 | res = Value::errorVALUE(); |
1906 | return; |
1907 | } |
1908 | for (int i = 0; i < a1.count(); ++i) |
1909 | twoArrayWalk(a1[i], a2[i], res, func); |
1910 | } |
1911 | |
1912 | arrayWalkFunc ValueCalc::awFunc(const QString &name) |
1913 | { |
1914 | if (awFuncs.count(name)) |
1915 | return awFuncs[name]; |
1916 | return 0; |
1917 | } |
1918 | |
1919 | void ValueCalc::registerAwFunc(const QString &name, arrayWalkFunc func) |
1920 | { |
1921 | awFuncs[name] = func; |
1922 | } |
1923 | |
1924 | // ------------------------------------------------------ |
1925 | |
1926 | Value ValueCalc::sum(const Value &range, bool full) |
1927 | { |
1928 | Value res(0); |
1929 | arrayWalk(range, res, full ? awSumA : awSum, Value(0)); |
1930 | return res; |
1931 | } |
1932 | |
1933 | Value ValueCalc::sum(QVector<Value> range, bool full) |
1934 | { |
1935 | Value res(0); |
1936 | arrayWalk(range, res, full ? awSumA : awSum, Value(0)); |
1937 | return res; |
1938 | } |
1939 | |
1940 | // sum of squares |
1941 | Value ValueCalc::sumsq(const Value &range, bool full) |
1942 | { |
1943 | Value res(0); |
1944 | arrayWalk(range, res, full ? awSumSqA : awSumSq, Value(0)); |
1945 | return res; |
1946 | } |
1947 | |
1948 | Value ValueCalc::sumIf(const Value &range, const Condition &cond) |
1949 | { |
1950 | if(range.isError()) |
1951 | return range; |
1952 | |
1953 | if (!range.isArray()) { |
1954 | if (matches(cond, range.element(0, 0))) { |
1955 | //kDebug()<<"return non array value "<<range; |
1956 | return range; |
1957 | } |
1958 | return Value(0.0); |
1959 | } |
1960 | |
1961 | //if we are here, we have an array |
1962 | Value res(0); |
1963 | Value tmp; |
1964 | |
1965 | unsigned int rows = range.rows(); |
1966 | unsigned int cols = range.columns(); |
1967 | for (unsigned int r = 0; r < rows; ++r) |
1968 | for (unsigned int c = 0; c < cols; ++c) { |
1969 | Value v = range.element(c, r); |
1970 | |
1971 | if (v.isArray()) |
1972 | tmp = sumIf(v, cond); |
1973 | if (tmp.isNumber()) {// only add numbers, no conversion from string allowed |
1974 | res = add(res, tmp); |
1975 | } else if (matches(cond, v)) { |
1976 | if (v.isNumber()) {// only add numbers, no conversion from string allowed |
1977 | //kDebug()<<"add "<<v; |
1978 | res = add(res, v); |
1979 | } |
1980 | } |
1981 | } |
1982 | |
1983 | return res; |
1984 | } |
1985 | |
1986 | Value ValueCalc::sumIf(const Cell &sumRangeStart, const Value &range, const Condition &cond) |
1987 | { |
1988 | if(range.isError()) |
1989 | return range; |
1990 | |
1991 | if (!range.isArray()) { |
1992 | if (matches(cond, range.element(0, 0))) { |
1993 | //kDebug()<<"return non array value "<<range; |
1994 | return sumRangeStart.value(); |
1995 | } |
1996 | return Value(0.0); |
1997 | } |
1998 | |
1999 | //if we are here, we have an array |
2000 | Value res(0); |
2001 | Value tmp; |
2002 | |
2003 | unsigned int rows = range.rows(); |
2004 | unsigned int cols = range.columns(); |
2005 | for (unsigned int r = 0; r < rows; ++r) |
2006 | for (unsigned int c = 0; c < cols; ++c) { |
2007 | Value v = range.element(c, r); |
2008 | if (v.isArray()) |
2009 | return Value::errorVALUE(); |
2010 | |
2011 | if (matches(cond, v)) { |
2012 | Value val = Cell(sumRangeStart.sheet(), sumRangeStart.column() + c, sumRangeStart.row() + r).value(); |
2013 | if (val.isNumber()) {// only add numbers, no conversion from string allowed |
2014 | res = add(res, val); |
2015 | } |
2016 | } |
2017 | } |
2018 | |
2019 | return res; |
2020 | } |
2021 | |
2022 | Value ValueCalc::sumIfs(const Cell &sumRangeStart, QList<Value> range, QList<Condition> cond, const float limit) |
2023 | { |
2024 | if(range[0].isError()) |
2025 | return range[0]; |
2026 | |
2027 | Value res(0); |
2028 | Value val; |
2029 | |
2030 | unsigned int rows = range[0].rows(); |
2031 | unsigned int cols = range[0].columns(); |
2032 | for (unsigned int r = 0; r < rows; ++r) { |
2033 | for (unsigned int c = 0; c < cols; ++c) { |
2034 | for (unsigned int i = 1; i <= limit ; ++i) { |
2035 | |
2036 | if(range[i].isError()) |
2037 | return range[0]; |
2038 | |
2039 | if (!range[i].isArray()) { |
2040 | if (matches(cond[i-1], range[i].element(0, 0))) { |
2041 | return sumRangeStart.value(); |
2042 | } |
2043 | return Value(0.0); |
2044 | } |
2045 | |
2046 | Value v = range[i].element(c, r); |
2047 | if (v.isArray()) { |
2048 | return Value::errorVALUE(); |
2049 | } |
2050 | |
2051 | if (!matches(cond[i-1], v)) { |
2052 | val = Value(0.0); |
2053 | break; |
2054 | } |
2055 | val = range[0].element(c, r); |
2056 | } |
2057 | if (val.isNumber()) {// only add numbers, no conversion from string allowed |
2058 | res = add(res, val); |
2059 | } |
2060 | } |
2061 | } |
2062 | return res; |
2063 | } |
2064 | |
2065 | Value ValueCalc::averageIf(const Value &range, const Condition &cond) |
2066 | { |
2067 | if(range.isError()) |
2068 | return range; |
2069 | |
2070 | if (!range.isArray()) { |
2071 | if (matches(cond, range.element(0, 0))) { |
2072 | return range; |
2073 | } |
2074 | return Value(0.0); |
2075 | } |
2076 | |
2077 | Value res(0); |
2078 | Value tmp; |
2079 | |
2080 | unsigned int rows = range.rows(); |
2081 | unsigned int cols = range.columns(); |
2082 | unsigned int cnt = 0; |
2083 | for (unsigned int r = 0; r < rows; ++r) |
2084 | for (unsigned int c = 0; c < cols; ++c) { |
2085 | Value v = range.element(c, r); |
2086 | |
2087 | if (v.isArray()) |
2088 | tmp = averageIf(v, cond); |
2089 | if (tmp.isNumber()) {// only add numbers, no conversion from string allowed |
2090 | res = add(res, tmp); |
2091 | } else if (matches(cond, v)) { |
2092 | if (v.isNumber()) {// only add numbers, no conversion from string allowed |
2093 | //kDebug()<<"add "<<v; |
2094 | res = add(res, v); |
2095 | ++cnt; |
2096 | } |
2097 | } |
2098 | } |
2099 | |
2100 | res = div(res, cnt); |
2101 | return res; |
2102 | } |
2103 | |
2104 | Value ValueCalc::averageIf(const Cell &avgRangeStart, const Value &range, const Condition &cond) |
2105 | { |
2106 | if(range.isError()) |
2107 | return range; |
2108 | |
2109 | if (!range.isArray()) { |
2110 | if (matches(cond, range.element(0, 0))) { |
2111 | return avgRangeStart.value(); |
2112 | } |
2113 | return Value(0.0); |
2114 | } |
2115 | |
2116 | Value res(0); |
2117 | |
2118 | unsigned int rows = range.rows(); |
2119 | unsigned int cols = range.columns(); |
2120 | unsigned int cnt = 0; |
2121 | for (unsigned int r = 0; r < rows; ++r) |
2122 | for (unsigned int c = 0; c < cols; ++c) { |
2123 | Value v = range.element(c, r); |
2124 | |
2125 | if (v.isArray()) |
2126 | return Value::errorVALUE(); |
2127 | |
2128 | if (matches(cond, v)) { |
2129 | Value val = Cell(avgRangeStart.sheet(), avgRangeStart.column() + c, avgRangeStart.row() + r).value(); |
2130 | if (val.isNumber()) {// only add numbers, no conversion from string allowed |
2131 | //kDebug()<<"add "<<val; |
2132 | res = add(res, val); |
2133 | ++cnt; |
2134 | } |
2135 | } |
2136 | } |
2137 | |
2138 | res = div(res, cnt); |
2139 | return res; |
2140 | } |
2141 | |
2142 | Value ValueCalc::averageIfs(const Cell &avgRangeStart, QList<Value> range, QList<Condition> cond, const float limit) |
2143 | { |
2144 | if(range[0].isError()) |
2145 | return range[0]; |
2146 | |
2147 | Value res(0); |
2148 | Value val; |
2149 | |
2150 | unsigned int rows = range[0].rows(); |
2151 | unsigned int cols = range[0].columns(); |
2152 | unsigned int cnt = 0; |
2153 | for (unsigned int r = 0; r < rows; ++r) { |
2154 | for (unsigned int c = 0; c < cols; ++c) { |
2155 | bool flag = true; |
2156 | for (unsigned int i = 1; i <= limit ; ++i) { |
2157 | |
2158 | if(range[i].isError()) |
2159 | return range[0]; |
2160 | |
2161 | if (!range[i].isArray()) { |
2162 | if (matches(cond[i-1], range[i].element(0, 0))) { |
2163 | return avgRangeStart.value(); |
2164 | } |
2165 | return Value(0.0); |
2166 | } |
2167 | |
2168 | Value v = range[i].element(c, r); |
2169 | if (v.isArray()) { |
2170 | return Value::errorVALUE(); |
2171 | } |
2172 | |
2173 | if (!matches(cond[i-1], v)) { |
2174 | flag = false; |
2175 | val = Value(0.0); |
2176 | break; |
2177 | } |
2178 | val = range[0].element(c, r); |
2179 | } |
2180 | if (flag) { |
2181 | ++cnt; |
2182 | flag = true; |
2183 | } |
2184 | if (val.isNumber()) {// only add numbers, no conversion from string allowed |
2185 | res = add(res, val); |
2186 | } |
2187 | } |
2188 | } |
2189 | res = div(res, cnt); |
2190 | return res; |
2191 | } |
2192 | |
2193 | int ValueCalc::count(const Value &range, bool full) |
2194 | { |
2195 | Value res(0); |
2196 | arrayWalk(range, res, full ? awCountA : awCount, Value(0)); |
2197 | return converter->asInteger(res).asInteger(); |
2198 | } |
2199 | |
2200 | int ValueCalc::count(QVector<Value> range, bool full) |
2201 | { |
2202 | Value res(0); |
2203 | arrayWalk(range, res, full ? awCountA : awCount, Value(0)); |
2204 | return converter->asInteger(res).asInteger(); |
2205 | } |
2206 | |
2207 | int ValueCalc::countIf(const Value &range, const Condition &cond) |
2208 | { |
2209 | if (!range.isArray()) { |
2210 | if (matches(cond, range)) |
2211 | return range.isEmpty() ? 0 : 1; |
2212 | return 0; |
2213 | } |
2214 | |
2215 | int res = 0; |
2216 | |
2217 | // iterate over the non-empty entries |
2218 | for (uint i = 0; i < range.count(); ++i) { |
2219 | Value v = range.element(i); |
2220 | |
2221 | if (v.isArray()) |
2222 | res += countIf(v, cond); |
2223 | else if (matches(cond, v)) |
2224 | ++res; |
2225 | } |
2226 | |
2227 | return res; |
2228 | } |
2229 | |
2230 | Value ValueCalc::countIfs(const Cell &cntRangeStart, QList<Value> range, QList<Condition> cond, const float limit) |
2231 | { |
2232 | if (!range[0].isArray()) |
2233 | return Value(0.0); |
2234 | |
2235 | if(range[0].isError()) |
2236 | return range[0]; |
2237 | |
2238 | Value res(0); |
2239 | |
2240 | unsigned int rows = range[0].rows(); |
2241 | unsigned int cols = range[0].columns(); |
2242 | for (unsigned int r = 0; r < rows; ++r) { |
2243 | for (unsigned int c = 0; c < cols; ++c) { |
2244 | bool flag = true; |
2245 | for (unsigned int i = 0; i <= limit ; ++i) { |
2246 | |
2247 | if(range[i].isError()) |
2248 | return range[0]; |
2249 | |
2250 | if (!range[i].isArray()) { |
2251 | if (matches(cond[i], range[i].element(0, 0))) { |
2252 | return cntRangeStart.value(); |
2253 | } |
2254 | return Value(0.0); |
2255 | } |
2256 | |
2257 | Value v = range[i].element(c, r); |
2258 | if (v.isArray()) { |
2259 | return Value::errorVALUE(); |
2260 | } |
2261 | |
2262 | if (!matches(cond[i], v)) { |
2263 | flag = false; |
2264 | break; |
2265 | } |
2266 | } |
2267 | if (flag) { |
2268 | res = add(res, 1); |
2269 | flag = true; |
2270 | } |
2271 | } |
2272 | } |
2273 | return res; |
2274 | } |
2275 | |
2276 | Value ValueCalc::avg(const Value &range, bool full) |
2277 | { |
2278 | int cnt = count(range, full); |
2279 | if (cnt) |
2280 | return div(sum(range, full), cnt); |
2281 | return Value(0.0); |
2282 | } |
2283 | |
2284 | Value ValueCalc::avg(QVector<Value> range, bool full) |
2285 | { |
2286 | int cnt = count(range, full); |
2287 | if (cnt) |
2288 | return div(sum(range, full), cnt); |
2289 | return Value(0.0); |
2290 | } |
2291 | |
2292 | Value ValueCalc::max(const Value &range, bool full) |
2293 | { |
2294 | Value res; |
2295 | arrayWalk(range, res, full ? awMaxA : awMax, Value(0)); |
2296 | return res; |
2297 | } |
2298 | |
2299 | Value ValueCalc::max(QVector<Value> range, bool full) |
2300 | { |
2301 | Value res; |
2302 | arrayWalk(range, res, full ? awMaxA: awMax, Value(0)); |
2303 | return res; |
2304 | } |
2305 | |
2306 | Value ValueCalc::min(const Value &range, bool full) |
2307 | { |
2308 | Value res; |
2309 | arrayWalk(range, res, full ? awMinA : awMin, Value(0)); |
2310 | return res; |
2311 | } |
2312 | |
2313 | Value ValueCalc::min(QVector<Value> range, bool full) |
2314 | { |
2315 | Value res; |
2316 | arrayWalk(range, res, full ? awMinA : awMin, Value(0)); |
2317 | return res; |
2318 | } |
2319 | |
2320 | Value ValueCalc::product(const Value &range, Value init, |
2321 | bool full) |
2322 | { |
2323 | Value res = init; |
2324 | if (isZero(init)) { // special handling of a zero, due to excel-compat |
2325 | if (count(range, full) == 0) |
2326 | return init; |
2327 | res = Value(1.0); |
2328 | } |
2329 | arrayWalk(range, res, full ? awProdA : awProd, Value(0)); |
2330 | return res; |
2331 | } |
2332 | |
2333 | Value ValueCalc::product(QVector<Value> range, |
2334 | Value init, bool full) |
2335 | { |
2336 | Value res = init; |
2337 | if (isZero(init)) { // special handling of a zero, due to excel-compat |
2338 | if (count(range, full) == 0) |
2339 | return init; |
2340 | res = Value(1.0); |
2341 | } |
2342 | arrayWalk(range, res, full ? awProdA : awProd, Value(0)); |
2343 | return res; |
2344 | } |
2345 | |
2346 | Value ValueCalc::stddev(const Value &range, bool full) |
2347 | { |
2348 | return stddev(range, avg(range, full), full); |
2349 | } |
2350 | |
2351 | Value ValueCalc::stddev(const Value &range, Value avg, |
2352 | bool full) |
2353 | { |
2354 | Value res; |
2355 | int cnt = count(range, full); |
2356 | arrayWalk(range, res, full ? awDevSqA : awDevSq, avg); |
2357 | return sqrt(div(res, cnt - 1)); |
2358 | } |
2359 | |
2360 | Value ValueCalc::stddev(QVector<Value> range, bool full) |
2361 | { |
2362 | return stddev(range, avg(range, full), full); |
2363 | } |
2364 | |
2365 | Value ValueCalc::stddev(QVector<Value> range, |
2366 | Value avg, bool full) |
2367 | { |
2368 | Value res; |
2369 | int cnt = count(range, full); |
2370 | arrayWalk(range, res, full ? awDevSqA : awDevSq, avg); |
2371 | return sqrt(div(res, cnt - 1)); |
2372 | } |
2373 | |
2374 | Value ValueCalc::stddevP(const Value &range, bool full) |
2375 | { |
2376 | return stddevP(range, avg(range, full), full); |
2377 | } |
2378 | |
2379 | Value ValueCalc::stddevP(const Value &range, Value avg, |
2380 | bool full) |
2381 | { |
2382 | Value res; |
2383 | int cnt = count(range, full); |
2384 | arrayWalk(range, res, full ? awDevSqA : awDevSq, avg); |
2385 | return sqrt(div(res, cnt)); |
2386 | } |
2387 | |
2388 | Value ValueCalc::stddevP(QVector<Value> range, bool full) |
2389 | { |
2390 | return stddevP(range, avg(range, full), full); |
2391 | } |
2392 | |
2393 | Value ValueCalc::stddevP(QVector<Value> range, |
2394 | Value avg, bool full) |
2395 | { |
2396 | Value res; |
2397 | int cnt = count(range, full); |
2398 | arrayWalk(range, res, full ? awDevSqA : awDevSq, avg); |
2399 | return sqrt(div(res, cnt)); |
2400 | } |
2401 | |
2402 | bool isDate(Value::Format fmt) |
2403 | { |
2404 | if ((fmt == Value::fmt_Date) || (fmt == Value::fmt_DateTime)) |
2405 | return true; |
2406 | return false; |
2407 | } |
2408 | |
2409 | Value::Format ValueCalc::format(Value a, Value b) |
2410 | { |
2411 | Value::Format af = a.format(); |
2412 | Value::Format bf = b.format(); |
2413 | |
2414 | // operation on two dates should produce a number |
2415 | if (isDate(af) && isDate(bf)) |
2416 | return Value::fmt_Number; |
2417 | |
2418 | if ((af == Value::fmt_None) || (af == Value::fmt_Boolean)) |
2419 | return bf; |
2420 | return af; |
2421 | } |
2422 | |
2423 | // ------------------------------------------------------ |
2424 | |
2425 | void ValueCalc::getCond(Condition &cond, Value val) |
2426 | { |
2427 | // not a string - we simply take it as a numeric value |
2428 | // that also handles floats, logical values, date/time and such |
2429 | if (!val.isString()) { |
2430 | cond.comp = isEqual; |
2431 | cond.type = numeric; |
2432 | cond.value = converter->toFloat(val); |
2433 | return; |
2434 | } |
2435 | QString text = converter->asString(val).asString(); |
2436 | cond.comp = isEqual; |
2437 | text = text.trimmed(); |
2438 | |
2439 | if (text.startsWith(QLatin1String("<=" ))) { |
2440 | cond.comp = lessEqual; |
2441 | text.remove(0, 2); |
2442 | } else if (text.startsWith(QLatin1String(">=" ))) { |
2443 | cond.comp = greaterEqual; |
2444 | text.remove(0, 2); |
2445 | } else if (text.startsWith(QLatin1String("!=" )) || text.startsWith(QLatin1String("<>" ))) { |
2446 | cond.comp = notEqual; |
2447 | text.remove(0, 2); |
2448 | } else if (text.startsWith(QLatin1String("==" ))) { |
2449 | cond.comp = isEqual; |
2450 | text.remove(0, 2); |
2451 | } else if (text.startsWith(QLatin1Char('<'))) { |
2452 | cond.comp = isLess; |
2453 | text.remove(0, 1); |
2454 | } else if (text.startsWith(QLatin1Char('>'))) { |
2455 | cond.comp = isGreater; |
2456 | text.remove(0, 1); |
2457 | } else if (text.startsWith(QLatin1Char('='))) { |
2458 | cond.comp = isEqual; |
2459 | text.remove(0, 1); |
2460 | } else { // character comparision |
2461 | cond.type = string; |
2462 | cond.stringValue = text; |
2463 | if (settings()->useWildcards()) { // HOST-USE-WILDCARDS Excel like wildcard matching |
2464 | cond.comp = wildcardMatch; |
2465 | } else if (settings()->useRegularExpressions()) { // HOST-USE-REGULAR-EXPRESSION ODF like regex matching |
2466 | cond.comp = regexMatch; |
2467 | } else { // Simple string matching |
2468 | cond.comp = stringMatch; |
2469 | } |
2470 | return; |
2471 | } |
2472 | |
2473 | text = text.trimmed(); |
2474 | |
2475 | bool ok = false; |
2476 | double d = text.toDouble(&ok); |
2477 | if (ok) { |
2478 | cond.type = numeric; |
2479 | cond.value = d; |
2480 | } else { |
2481 | cond.type = string; |
2482 | cond.stringValue = text; |
2483 | } |
2484 | //TODO: date values |
2485 | } |
2486 | |
2487 | bool ValueCalc::matches(const Condition &cond, Value val) |
2488 | { |
2489 | if (val.isEmpty()) |
2490 | return false; |
2491 | if (cond.type == numeric) { |
2492 | Number d = converter->toFloat(val); |
2493 | switch (cond.comp) { |
2494 | case isEqual: |
2495 | if (approxEqual(Value(d), Value(cond.value))) return true; |
2496 | break; |
2497 | |
2498 | case isLess: |
2499 | if (d < cond.value) return true; |
2500 | break; |
2501 | |
2502 | case isGreater: |
2503 | if (d > cond.value) return true; |
2504 | break; |
2505 | |
2506 | case lessEqual: |
2507 | if (d <= cond.value) return true; |
2508 | break; |
2509 | |
2510 | case greaterEqual: |
2511 | if (d >= cond.value) return true; |
2512 | break; |
2513 | |
2514 | case notEqual: |
2515 | if (d != cond.value) return true; |
2516 | break; |
2517 | |
2518 | default: |
2519 | break; |
2520 | } |
2521 | } else { |
2522 | QString d = converter->asString(val).asString(); |
2523 | switch (cond.comp) { |
2524 | case isEqual: |
2525 | if (d == cond.stringValue) return true; |
2526 | break; |
2527 | |
2528 | case isLess: |
2529 | if (d < cond.stringValue) return true; |
2530 | break; |
2531 | |
2532 | case isGreater: |
2533 | if (d > cond.stringValue) return true; |
2534 | break; |
2535 | |
2536 | case lessEqual: |
2537 | if (d <= cond.stringValue) return true; |
2538 | break; |
2539 | |
2540 | case greaterEqual: |
2541 | if (d >= cond.stringValue) return true; |
2542 | break; |
2543 | |
2544 | case notEqual: |
2545 | if (d != cond.stringValue) return true; |
2546 | break; |
2547 | |
2548 | case stringMatch: |
2549 | if (d.toLower() == cond.stringValue.toLower()) return true; |
2550 | break; |
2551 | |
2552 | case regexMatch: { |
2553 | QRegExp rx; |
2554 | rx.setPattern(cond.stringValue); |
2555 | rx.setPatternSyntax(QRegExp::RegExp); |
2556 | rx.setCaseSensitivity(Qt::CaseInsensitive); |
2557 | if (rx.exactMatch(d)) return true; |
2558 | } break; |
2559 | |
2560 | case wildcardMatch: { |
2561 | QRegExp rx; |
2562 | rx.setPattern(cond.stringValue); |
2563 | rx.setPatternSyntax(QRegExp::Wildcard); |
2564 | rx.setCaseSensitivity(Qt::CaseInsensitive); |
2565 | if (rx.exactMatch(d)) return true; |
2566 | } break; |
2567 | |
2568 | } |
2569 | } |
2570 | return false; |
2571 | } |
2572 | |
2573 | |