Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /*
3 : * This file is part of the LibreOffice project.
4 : *
5 : * This Source Code Form is subject to the terms of the Mozilla Public
6 : * License, v. 2.0. If a copy of the MPL was not distributed with this
7 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 : *
9 : * This file incorporates work covered by the following license notice:
10 : *
11 : * Licensed to the Apache Software Foundation (ASF) under one or more
12 : * contributor license agreements. See the NOTICE file distributed
13 : * with this work for additional information regarding copyright
14 : * ownership. The ASF licenses this file to you under the Apache
15 : * License, Version 2.0 (the "License"); you may not use this file
16 : * except in compliance with the License. You may obtain a copy of
17 : * the License at http://www.apache.org/licenses/LICENSE-2.0 .
18 : */
19 :
20 : #include <tools/solar.h>
21 : #include <stdlib.h>
22 : #include <string.h>
23 :
24 : #include "interpre.hxx"
25 : #include "global.hxx"
26 : #include "compiler.hxx"
27 : #include "formulacell.hxx"
28 : #include "document.hxx"
29 : #include "dociter.hxx"
30 : #include "scmatrix.hxx"
31 : #include "globstr.hrc"
32 :
33 : #include <math.h>
34 : #include <vector>
35 : #include <algorithm>
36 : #include <boost/math/special_functions/log1p.hpp>
37 : #include <comphelper/random.hxx>
38 :
39 : using ::std::vector;
40 : using namespace formula;
41 :
42 : // STATIC DATA
43 : #define MAX_ANZ_DOUBLE_FOR_SORT 100000
44 :
45 : const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental
46 : const double fMachEps = ::std::numeric_limits<double>::epsilon();
47 :
48 70 : class ScDistFunc
49 : {
50 : public:
51 : virtual double GetValue(double x) const = 0;
52 :
53 : protected:
54 70 : ~ScDistFunc() {}
55 : };
56 :
57 : // iteration for inverse distributions
58 :
59 : //template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, bool& rConvError )
60 :
61 : /** u*w<0.0 fails for values near zero */
62 1368 : static inline bool lcl_HasChangeOfSign( double u, double w )
63 : {
64 1368 : return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);
65 : }
66 :
67 70 : static double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError )
68 : {
69 70 : rConvError = false;
70 70 : const double fYEps = 1.0E-307;
71 70 : const double fXEps = ::std::numeric_limits<double>::epsilon();
72 :
73 : OSL_ENSURE(fAx<fBx, "IterateInverse: wrong interval");
74 :
75 : // find enclosing interval
76 :
77 70 : double fAy = rFunction.GetValue(fAx);
78 70 : double fBy = rFunction.GetValue(fBx);
79 : double fTemp;
80 : unsigned short nCount;
81 138 : for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++)
82 : {
83 68 : if (fabs(fAy) <= fabs(fBy))
84 : {
85 39 : fTemp = fAx;
86 39 : fAx += 2.0 * (fAx - fBx);
87 39 : if (fAx < 0.0)
88 39 : fAx = 0.0;
89 39 : fBx = fTemp;
90 39 : fBy = fAy;
91 39 : fAy = rFunction.GetValue(fAx);
92 : }
93 : else
94 : {
95 29 : fTemp = fBx;
96 29 : fBx += 2.0 * (fBx - fAx);
97 29 : fAx = fTemp;
98 29 : fAy = fBy;
99 29 : fBy = rFunction.GetValue(fBx);
100 : }
101 : }
102 :
103 70 : if (fAy == 0.0)
104 0 : return fAx;
105 70 : if (fBy == 0.0)
106 0 : return fBx;
107 70 : if (!lcl_HasChangeOfSign( fAy, fBy))
108 : {
109 0 : rConvError = true;
110 0 : return 0.0;
111 : }
112 : // inverse quadric interpolation with additional brackets
113 : // set three points
114 70 : double fPx = fAx;
115 70 : double fPy = fAy;
116 70 : double fQx = fBx;
117 70 : double fQy = fBy;
118 70 : double fRx = fAx;
119 70 : double fRy = fAy;
120 70 : double fSx = 0.5 * (fAx + fBx); // potential next point
121 70 : bool bHasToInterpolate = true;
122 70 : nCount = 0;
123 6206 : while ( nCount < 500 && fabs(fRy) > fYEps &&
124 6122 : (fBx-fAx) > ::std::max( fabs(fAx), fabs(fBx)) * fXEps )
125 : {
126 1160 : if (bHasToInterpolate)
127 : {
128 731 : if (fPy!=fQy && fQy!=fRy && fRy!=fPy)
129 : {
130 460 : fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)
131 460 : + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)
132 460 : + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);
133 460 : bHasToInterpolate = (fAx < fSx) && (fSx < fBx); // inside the brackets?
134 : }
135 : else
136 271 : bHasToInterpolate = false;
137 : }
138 1160 : if(!bHasToInterpolate)
139 : {
140 754 : fSx = 0.5 * (fAx + fBx);
141 : // reset points
142 754 : fQx = fBx; fQy = fBy;
143 754 : bHasToInterpolate = true;
144 : }
145 : // shift points for next interpolation
146 1160 : fPx = fQx; fQx = fRx; fRx = fSx;
147 1160 : fPy = fQy; fQy = fRy; fRy = rFunction.GetValue(fSx);
148 : // update brackets
149 1160 : if (lcl_HasChangeOfSign( fAy, fRy))
150 : {
151 717 : fBx = fRx; fBy = fRy;
152 : }
153 : else
154 : {
155 443 : fAx = fRx; fAy = fRy;
156 : }
157 : // if last interration brought to small advance, then do bisection next
158 : // time, for safety
159 1160 : bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));
160 1160 : ++nCount;
161 : }
162 70 : return fRx;
163 : }
164 :
165 : // Allgemeine Funktionen
166 :
167 0 : void ScInterpreter::ScNoName()
168 : {
169 0 : PushError(errNoName);
170 0 : }
171 :
172 13 : void ScInterpreter::ScBadName()
173 : {
174 13 : short nParamCount = GetByte();
175 26 : while (nParamCount-- > 0)
176 : {
177 0 : PopError();
178 : }
179 13 : PushError( errNoName);
180 13 : }
181 :
182 21 : double ScInterpreter::phi(double x)
183 : {
184 21 : return 0.39894228040143268 * exp(-(x * x) / 2.0);
185 : }
186 :
187 27 : double ScInterpreter::integralPhi(double x)
188 : { // Using gauss(x)+0.5 has severe cancellation errors for x<-4
189 27 : return 0.5 * ::rtl::math::erfc(-x * 0.7071067811865475); // * 1/sqrt(2)
190 : }
191 :
192 9 : double ScInterpreter::taylor(const double* pPolynom, sal_uInt16 nMax, double x)
193 : {
194 9 : double nVal = pPolynom[nMax];
195 144 : for (short i = nMax-1; i >= 0; i--)
196 : {
197 135 : nVal = pPolynom[i] + (nVal * x);
198 : }
199 9 : return nVal;
200 : }
201 :
202 9 : double ScInterpreter::gauss(double x)
203 : {
204 :
205 9 : double xAbs = fabs(x);
206 9 : sal_uInt16 xShort = (sal_uInt16)::rtl::math::approxFloor(xAbs);
207 9 : double nVal = 0.0;
208 9 : if (xShort == 0)
209 : {
210 : static const double t0[] =
211 : { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,
212 : -0.00118732821548045, 0.00011543468761616, -0.00000944465625950,
213 : 0.00000066596935163, -0.00000004122667415, 0.00000000227352982,
214 : 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };
215 6 : nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;
216 : }
217 3 : else if ((xShort >= 1) && (xShort <= 2))
218 : {
219 : static const double t2[] =
220 : { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,
221 : 0.02699548325659403, -0.00449924720943234, -0.00224962360471617,
222 : 0.00134977416282970, -0.00011783742691370, -0.00011515930357476,
223 : 0.00003704737285544, 0.00000282690796889, -0.00000354513195524,
224 : 0.00000037669563126, 0.00000019202407921, -0.00000005226908590,
225 : -0.00000000491799345, 0.00000000366377919, -0.00000000015981997,
226 : -0.00000000017381238, 0.00000000002624031, 0.00000000000560919,
227 : -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };
228 3 : nVal = taylor(t2, 23, (xAbs - 2.0));
229 : }
230 0 : else if ((xShort >= 3) && (xShort <= 4))
231 : {
232 : static const double t4[] =
233 : { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,
234 : 0.00033457556441221, -0.00028996548915725, 0.00018178605666397,
235 : -0.00008252863922168, 0.00002551802519049, -0.00000391665839292,
236 : -0.00000074018205222, 0.00000064422023359, -0.00000017370155340,
237 : 0.00000000909595465, 0.00000000944943118, -0.00000000329957075,
238 : 0.00000000029492075, 0.00000000011874477, -0.00000000004420396,
239 : 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };
240 0 : nVal = taylor(t4, 20, (xAbs - 4.0));
241 : }
242 : else
243 : {
244 : static const double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };
245 0 : nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs;
246 : }
247 9 : if (x < 0.0)
248 3 : return -nVal;
249 : else
250 6 : return nVal;
251 : }
252 :
253 : // #i26836# new gaussinv implementation by Martin Eitzenberger <m.eitzenberger@unix.net>
254 :
255 36 : double ScInterpreter::gaussinv(double x)
256 : {
257 : double q,t,z;
258 :
259 36 : q=x-0.5;
260 :
261 36 : if(fabs(q)<=.425)
262 : {
263 21 : t=0.180625-q*q;
264 :
265 : z=
266 21 : q*
267 : (
268 : (
269 : (
270 : (
271 : (
272 : (
273 : (
274 21 : t*2509.0809287301226727+33430.575583588128105
275 : )
276 21 : *t+67265.770927008700853
277 : )
278 21 : *t+45921.953931549871457
279 : )
280 21 : *t+13731.693765509461125
281 : )
282 21 : *t+1971.5909503065514427
283 : )
284 21 : *t+133.14166789178437745
285 : )
286 21 : *t+3.387132872796366608
287 : )
288 : /
289 : (
290 : (
291 : (
292 : (
293 : (
294 : (
295 : (
296 21 : t*5226.495278852854561+28729.085735721942674
297 : )
298 21 : *t+39307.89580009271061
299 : )
300 21 : *t+21213.794301586595867
301 : )
302 21 : *t+5394.1960214247511077
303 : )
304 21 : *t+687.1870074920579083
305 : )
306 21 : *t+42.313330701600911252
307 : )
308 21 : *t+1.0
309 21 : );
310 :
311 : }
312 : else
313 : {
314 15 : if(q>0) t=1-x;
315 6 : else t=x;
316 :
317 15 : t=sqrt(-log(t));
318 :
319 15 : if(t<=5.0)
320 : {
321 15 : t+=-1.6;
322 :
323 : z=
324 : (
325 : (
326 : (
327 : (
328 : (
329 : (
330 : (
331 15 : t*7.7454501427834140764e-4+0.0227238449892691845833
332 : )
333 15 : *t+0.24178072517745061177
334 : )
335 15 : *t+1.27045825245236838258
336 : )
337 15 : *t+3.64784832476320460504
338 : )
339 15 : *t+5.7694972214606914055
340 : )
341 15 : *t+4.6303378461565452959
342 : )
343 15 : *t+1.42343711074968357734
344 : )
345 : /
346 : (
347 : (
348 : (
349 : (
350 : (
351 : (
352 : (
353 15 : t*1.05075007164441684324e-9+5.475938084995344946e-4
354 : )
355 15 : *t+0.0151986665636164571966
356 : )
357 15 : *t+0.14810397642748007459
358 : )
359 15 : *t+0.68976733498510000455
360 : )
361 15 : *t+1.6763848301838038494
362 : )
363 15 : *t+2.05319162663775882187
364 : )
365 15 : *t+1.0
366 15 : );
367 :
368 : }
369 : else
370 : {
371 0 : t+=-5.0;
372 :
373 : z=
374 : (
375 : (
376 : (
377 : (
378 : (
379 : (
380 : (
381 0 : t*2.01033439929228813265e-7+2.71155556874348757815e-5
382 : )
383 0 : *t+0.0012426609473880784386
384 : )
385 0 : *t+0.026532189526576123093
386 : )
387 0 : *t+0.29656057182850489123
388 : )
389 0 : *t+1.7848265399172913358
390 : )
391 0 : *t+5.4637849111641143699
392 : )
393 0 : *t+6.6579046435011037772
394 : )
395 : /
396 : (
397 : (
398 : (
399 : (
400 : (
401 : (
402 : (
403 0 : t*2.04426310338993978564e-15+1.4215117583164458887e-7
404 : )
405 0 : *t+1.8463183175100546818e-5
406 : )
407 0 : *t+7.868691311456132591e-4
408 : )
409 0 : *t+0.0148753612908506148525
410 : )
411 0 : *t+0.13692988092273580531
412 : )
413 0 : *t+0.59983220655588793769
414 : )
415 0 : *t+1.0
416 0 : );
417 :
418 : }
419 :
420 15 : if(q<0.0) z=-z;
421 : }
422 :
423 36 : return z;
424 : }
425 :
426 6 : double ScInterpreter::Fakultaet(double x)
427 : {
428 6 : x = ::rtl::math::approxFloor(x);
429 6 : if (x < 0.0)
430 0 : return 0.0;
431 6 : else if (x == 0.0)
432 3 : return 1.0;
433 3 : else if (x <= 170.0)
434 : {
435 3 : double fTemp = x;
436 9 : while (fTemp > 2.0)
437 : {
438 3 : fTemp--;
439 3 : x *= fTemp;
440 : }
441 : }
442 : else
443 0 : SetError(errNoValue);
444 3 : return x;
445 : }
446 :
447 6 : double ScInterpreter::BinomKoeff(double n, double k)
448 : {
449 : // this method has been duplicated as BinomialCoefficient()
450 : // in scaddins/source/analysis/analysishelper.cxx
451 :
452 6 : double nVal = 0.0;
453 6 : k = ::rtl::math::approxFloor(k);
454 6 : if (n < k)
455 0 : nVal = 0.0;
456 6 : else if (k == 0.0)
457 0 : nVal = 1.0;
458 : else
459 : {
460 6 : nVal = n/k;
461 6 : n--;
462 6 : k--;
463 18 : while (k > 0.0)
464 : {
465 6 : nVal *= n/k;
466 6 : k--;
467 6 : n--;
468 : }
469 :
470 : }
471 6 : return nVal;
472 : }
473 :
474 : // The algorithm is based on lanczos13m53 in lanczos.hpp
475 : // in math library from http://www.boost.org
476 : /** you must ensure fZ>0
477 : Uses a variant of the Lanczos sum with a rational function. */
478 4275 : static double lcl_getLanczosSum(double fZ)
479 : {
480 : static const double fNum[13] ={
481 : 23531376880.41075968857200767445163675473,
482 : 42919803642.64909876895789904700198885093,
483 : 35711959237.35566804944018545154716670596,
484 : 17921034426.03720969991975575445893111267,
485 : 6039542586.35202800506429164430729792107,
486 : 1439720407.311721673663223072794912393972,
487 : 248874557.8620541565114603864132294232163,
488 : 31426415.58540019438061423162831820536287,
489 : 2876370.628935372441225409051620849613599,
490 : 186056.2653952234950402949897160456992822,
491 : 8071.672002365816210638002902272250613822,
492 : 210.8242777515793458725097339207133627117,
493 : 2.506628274631000270164908177133837338626
494 : };
495 : static const double fDenom[13] = {
496 : 0,
497 : 39916800,
498 : 120543840,
499 : 150917976,
500 : 105258076,
501 : 45995730,
502 : 13339535,
503 : 2637558,
504 : 357423,
505 : 32670,
506 : 1925,
507 : 66,
508 : 1
509 : };
510 : // Horner scheme
511 : double fSumNum;
512 : double fSumDenom;
513 : int nI;
514 4275 : if (fZ<=1.0)
515 : {
516 810 : fSumNum = fNum[12];
517 810 : fSumDenom = fDenom[12];
518 10530 : for (nI = 11; nI >= 0; --nI)
519 : {
520 9720 : fSumNum *= fZ;
521 9720 : fSumNum += fNum[nI];
522 9720 : fSumDenom *= fZ;
523 9720 : fSumDenom += fDenom[nI];
524 : }
525 : }
526 : else
527 : // Cancel down with fZ^12; Horner scheme with reverse coefficients
528 : {
529 3465 : double fZInv = 1/fZ;
530 3465 : fSumNum = fNum[0];
531 3465 : fSumDenom = fDenom[0];
532 45045 : for (nI = 1; nI <=12; ++nI)
533 : {
534 41580 : fSumNum *= fZInv;
535 41580 : fSumNum += fNum[nI];
536 41580 : fSumDenom *= fZInv;
537 41580 : fSumDenom += fDenom[nI];
538 : }
539 : }
540 4275 : return fSumNum/fSumDenom;
541 : }
542 :
543 : // The algorithm is based on tgamma in gamma.hpp
544 : // in math library from http://www.boost.org
545 : /** You must ensure fZ>0; fZ>171.624376956302 will overflow. */
546 1590 : static double lcl_GetGammaHelper(double fZ)
547 : {
548 1590 : double fGamma = lcl_getLanczosSum(fZ);
549 1590 : const double fg = 6.024680040776729583740234375;
550 1590 : double fZgHelp = fZ + fg - 0.5;
551 : // avoid intermediate overflow
552 1590 : double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25);
553 1590 : fGamma *= fHalfpower;
554 1590 : fGamma /= exp(fZgHelp);
555 1590 : fGamma *= fHalfpower;
556 1590 : if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ))
557 822 : fGamma = ::rtl::math::round(fGamma);
558 1590 : return fGamma;
559 : }
560 :
561 : // The algorithm is based on tgamma in gamma.hpp
562 : // in math library from http://www.boost.org
563 : /** You must ensure fZ>0 */
564 12 : static double lcl_GetLogGammaHelper(double fZ)
565 : {
566 12 : const double fg = 6.024680040776729583740234375;
567 12 : double fZgHelp = fZ + fg - 0.5;
568 12 : return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;
569 : }
570 :
571 : /** You must ensure non integer arguments for fZ<1 */
572 1053 : double ScInterpreter::GetGamma(double fZ)
573 : {
574 1053 : const double fLogPi = log(F_PI);
575 1053 : const double fLogDblMax = log( ::std::numeric_limits<double>::max());
576 :
577 1053 : if (fZ > fMaxGammaArgument)
578 : {
579 0 : SetError(errIllegalFPOperation);
580 0 : return HUGE_VAL;
581 : }
582 :
583 1053 : if (fZ >= 1.0)
584 1041 : return lcl_GetGammaHelper(fZ);
585 :
586 12 : if (fZ >= 0.5) // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
587 3 : return lcl_GetGammaHelper(fZ+1) / fZ;
588 :
589 9 : if (fZ >= -0.5) // shift to x>=1, might overflow
590 : {
591 9 : double fLogTest = lcl_GetLogGammaHelper(fZ+2) - boost::math::log1p(fZ) - log( fabs(fZ));
592 9 : if (fLogTest >= fLogDblMax)
593 : {
594 0 : SetError( errIllegalFPOperation);
595 0 : return HUGE_VAL;
596 : }
597 9 : return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
598 : }
599 : // fZ<-0.5
600 : // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) )
601 0 : double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( fabs( ::rtl::math::sin( F_PI*fZ)));
602 0 : if (fLogDivisor - fLogPi >= fLogDblMax) // underflow
603 0 : return 0.0;
604 :
605 0 : if (fLogDivisor<0.0)
606 0 : if (fLogPi - fLogDivisor > fLogDblMax) // overflow
607 : {
608 0 : SetError(errIllegalFPOperation);
609 0 : return HUGE_VAL;
610 : }
611 :
612 0 : return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( F_PI*fZ) < 0.0) ? -1.0 : 1.0);
613 : }
614 :
615 : /** You must ensure fZ>0 */
616 540 : double ScInterpreter::GetLogGamma(double fZ)
617 : {
618 540 : if (fZ >= fMaxGammaArgument)
619 0 : return lcl_GetLogGammaHelper(fZ);
620 540 : if (fZ >= 1.0)
621 489 : return log(lcl_GetGammaHelper(fZ));
622 51 : if (fZ >= 0.5)
623 48 : return log( lcl_GetGammaHelper(fZ+1) / fZ);
624 3 : return lcl_GetLogGammaHelper(fZ+2) - boost::math::log1p(fZ) - log(fZ);
625 : }
626 :
627 264 : double ScInterpreter::GetFDist(double x, double fF1, double fF2)
628 : {
629 264 : double arg = fF2/(fF2+fF1*x);
630 264 : double alpha = fF2/2.0;
631 264 : double beta = fF1/2.0;
632 264 : return GetBetaDist(arg, alpha, beta);
633 : }
634 :
635 465 : double ScInterpreter::GetTDist( double T, double fDF, int nType )
636 : {
637 465 : switch ( nType )
638 : {
639 : case 1 : // 1-tailed T-distribution
640 12 : return 0.5 * GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5 );
641 : case 2 : // 2-tailed T-distribution
642 321 : return GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5);
643 : case 3 : // left-tailed T-distribution (probability density function)
644 3 : return pow( 1 + ( T * T / fDF ), -( fDF + 1 ) / 2 ) / ( sqrt( fDF ) * GetBeta( 0.5, fDF / 2.0 ) );
645 : case 4 : // left-tailed T-distribution (cumulative distribution function)
646 129 : double X = fDF / ( T * T + fDF );
647 129 : double R = 0.5 * GetBetaDist( X, 0.5 * fDF, 0.5 );
648 129 : return ( T < 0 ? R : 1 - R );
649 : }
650 0 : SetError( errIllegalArgument );
651 0 : return HUGE_VAL;
652 : }
653 :
654 : // for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom
655 : /** You must ensure fDF>0.0 */
656 204 : double ScInterpreter::GetChiDist(double fX, double fDF)
657 : {
658 204 : if (fX <= 0.0)
659 0 : return 1.0; // see ODFF
660 : else
661 204 : return GetUpRegIGamma( fDF/2.0, fX/2.0);
662 : }
663 :
664 : // ready for ODF 1.2
665 : // for ODF CHISQDIST; cumulative distribution function, fDF=degrees of freedom
666 : // returns left tail
667 : /** You must ensure fDF>0.0 */
668 90 : double ScInterpreter::GetChiSqDistCDF(double fX, double fDF)
669 : {
670 90 : if (fX <= 0.0)
671 6 : return 0.0; // see ODFF
672 : else
673 84 : return GetLowRegIGamma( fDF/2.0, fX/2.0);
674 : }
675 :
676 6 : double ScInterpreter::GetChiSqDistPDF(double fX, double fDF)
677 : {
678 : // you must ensure fDF is positive integer
679 : double fValue;
680 6 : if (fX <= 0.0)
681 0 : return 0.0; // see ODFF
682 6 : if (fDF*fX > 1391000.0)
683 : {
684 : // intermediate invalid values, use log
685 0 : fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0) - GetLogGamma(0.5*fDF));
686 : }
687 : else // fDF is small in most cases, we can iterate
688 : {
689 : double fCount;
690 6 : if (fmod(fDF,2.0)<0.5)
691 : {
692 : // even
693 3 : fValue = 0.5;
694 3 : fCount = 2.0;
695 : }
696 : else
697 : {
698 3 : fValue = 1/sqrt(fX*2*F_PI);
699 3 : fCount = 1.0;
700 : }
701 33 : while ( fCount < fDF)
702 : {
703 21 : fValue *= (fX / fCount);
704 21 : fCount += 2.0;
705 : }
706 6 : if (fX>=1425.0) // underflow in e^(-x/2)
707 0 : fValue = exp(log(fValue)-fX/2);
708 : else
709 6 : fValue *= exp(-fX/2);
710 : }
711 6 : return fValue;
712 : }
713 :
714 6 : void ScInterpreter::ScChiSqDist()
715 : {
716 6 : sal_uInt8 nParamCount = GetByte();
717 6 : if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
718 6 : return;
719 : bool bCumulative;
720 6 : if (nParamCount == 3)
721 3 : bCumulative = GetBool();
722 : else
723 3 : bCumulative = true;
724 6 : double fDF = ::rtl::math::approxFloor(GetDouble());
725 6 : if (fDF < 1.0)
726 0 : PushIllegalArgument();
727 : else
728 : {
729 6 : double fX = GetDouble();
730 6 : if (bCumulative)
731 3 : PushDouble(GetChiSqDistCDF(fX,fDF));
732 : else
733 3 : PushDouble(GetChiSqDistPDF(fX,fDF));
734 : }
735 : }
736 :
737 9 : void ScInterpreter::ScChiSqDist_MS()
738 : {
739 9 : sal_uInt8 nParamCount = GetByte();
740 9 : if ( !MustHaveParamCount( nParamCount, 3, 3 ) )
741 9 : return;
742 9 : bool bCumulative = GetBool();
743 9 : double fDF = ::rtl::math::approxFloor( GetDouble() );
744 9 : if ( fDF < 1.0 || fDF > 1E10 )
745 0 : PushIllegalArgument();
746 : else
747 : {
748 9 : double fX = GetDouble();
749 9 : if ( fX < 0 )
750 0 : PushIllegalArgument();
751 : else
752 : {
753 9 : if ( bCumulative )
754 6 : PushDouble( GetChiSqDistCDF( fX, fDF ) );
755 : else
756 3 : PushDouble( GetChiSqDistPDF( fX, fDF ) );
757 : }
758 : }
759 : }
760 :
761 0 : void ScInterpreter::ScGamma()
762 : {
763 0 : double x = GetDouble();
764 0 : if (x <= 0.0 && x == ::rtl::math::approxFloor(x))
765 0 : PushIllegalArgument();
766 : else
767 : {
768 0 : double fResult = GetGamma(x);
769 0 : if (nGlobalError)
770 : {
771 0 : PushError( nGlobalError);
772 0 : return;
773 : }
774 0 : PushDouble(fResult);
775 : }
776 : }
777 :
778 9 : void ScInterpreter::ScLogGamma()
779 : {
780 9 : double x = GetDouble();
781 9 : if (x > 0.0) // constraint from ODFF
782 9 : PushDouble( GetLogGamma(x));
783 : else
784 0 : PushIllegalArgument();
785 9 : }
786 :
787 351 : double ScInterpreter::GetBeta(double fAlpha, double fBeta)
788 : {
789 : double fA;
790 : double fB;
791 351 : if (fAlpha > fBeta)
792 : {
793 153 : fA = fAlpha; fB = fBeta;
794 : }
795 : else
796 : {
797 198 : fA = fBeta; fB = fAlpha;
798 : }
799 351 : if (fA+fB < fMaxGammaArgument) // simple case
800 351 : return GetGamma(fA)/GetGamma(fA+fB)*GetGamma(fB);
801 : // need logarithm
802 : // GetLogGamma is not accurate enough, back to Lanczos for all three
803 : // GetGamma and arrange factors newly.
804 0 : const double fg = 6.024680040776729583740234375; //see GetGamma
805 0 : double fgm = fg - 0.5;
806 0 : double fLanczos = lcl_getLanczosSum(fA);
807 0 : fLanczos /= lcl_getLanczosSum(fA+fB);
808 0 : fLanczos *= lcl_getLanczosSum(fB);
809 0 : double fABgm = fA+fB+fgm;
810 0 : fLanczos *= sqrt((fABgm/(fA+fgm))/(fB+fgm));
811 0 : double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
812 0 : double fTempB = fA/(fB+fgm);
813 0 : double fResult = exp(-fA * ::rtl::math::log1p(fTempA)
814 0 : -fB * ::rtl::math::log1p(fTempB)-fgm);
815 0 : fResult *= fLanczos;
816 0 : return fResult;
817 : }
818 :
819 : // Same as GetBeta but with logarithm
820 891 : double ScInterpreter::GetLogBeta(double fAlpha, double fBeta)
821 : {
822 : double fA;
823 : double fB;
824 891 : if (fAlpha > fBeta)
825 : {
826 591 : fA = fAlpha; fB = fBeta;
827 : }
828 : else
829 : {
830 300 : fA = fBeta; fB = fAlpha;
831 : }
832 891 : const double fg = 6.024680040776729583740234375; //see GetGamma
833 891 : double fgm = fg - 0.5;
834 891 : double fLanczos = lcl_getLanczosSum(fA);
835 891 : fLanczos /= lcl_getLanczosSum(fA+fB);
836 891 : fLanczos *= lcl_getLanczosSum(fB);
837 891 : double fLogLanczos = log(fLanczos);
838 891 : double fABgm = fA+fB+fgm;
839 891 : fLogLanczos += 0.5*(log(fABgm)-log(fA+fgm)-log(fB+fgm));
840 891 : double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
841 891 : double fTempB = fA/(fB+fgm);
842 891 : double fResult = -fA * ::rtl::math::log1p(fTempA)
843 891 : -fB * ::rtl::math::log1p(fTempB)-fgm;
844 891 : fResult += fLogLanczos;
845 891 : return fResult;
846 : }
847 :
848 : // beta distribution probability density function
849 342 : double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB)
850 : {
851 : // special cases
852 342 : if (fA == 1.0) // result b*(1-x)^(b-1)
853 : {
854 0 : if (fB == 1.0)
855 0 : return 1.0;
856 0 : if (fB == 2.0)
857 0 : return -2.0*fX + 2.0;
858 0 : if (fX == 1.0 && fB < 1.0)
859 : {
860 0 : SetError(errIllegalArgument);
861 0 : return HUGE_VAL;
862 : }
863 0 : if (fX <= 0.01)
864 0 : return fB + fB * ::rtl::math::expm1((fB-1.0) * ::rtl::math::log1p(-fX));
865 : else
866 0 : return fB * pow(0.5-fX+0.5,fB-1.0);
867 : }
868 342 : if (fB == 1.0) // result a*x^(a-1)
869 : {
870 0 : if (fA == 2.0)
871 0 : return fA * fX;
872 0 : if (fX == 0.0 && fA < 1.0)
873 : {
874 0 : SetError(errIllegalArgument);
875 0 : return HUGE_VAL;
876 : }
877 0 : return fA * pow(fX,fA-1);
878 : }
879 342 : if (fX <= 0.0)
880 : {
881 0 : if (fA < 1.0 && fX == 0.0)
882 : {
883 0 : SetError(errIllegalArgument);
884 0 : return HUGE_VAL;
885 : }
886 : else
887 0 : return 0.0;
888 : }
889 342 : if (fX >= 1.0)
890 : {
891 0 : if (fB < 1.0 && fX == 1.0)
892 : {
893 0 : SetError(errIllegalArgument);
894 0 : return HUGE_VAL;
895 : }
896 : else
897 0 : return 0.0;
898 : }
899 :
900 : // normal cases; result x^(a-1)*(1-x)^(b-1)/Beta(a,b)
901 342 : const double fLogDblMax = log( ::std::numeric_limits<double>::max());
902 342 : const double fLogDblMin = log( ::std::numeric_limits<double>::min());
903 342 : double fLogY = (fX < 0.1) ? ::rtl::math::log1p(-fX) : log(0.5-fX+0.5);
904 342 : double fLogX = log(fX);
905 342 : double fAm1LogX = (fA-1.0) * fLogX;
906 342 : double fBm1LogY = (fB-1.0) * fLogY;
907 342 : double fLogBeta = GetLogBeta(fA,fB);
908 : // check whether parts over- or underflow
909 342 : if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin
910 342 : && fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin
911 342 : && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin
912 342 : && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > fLogDblMin)
913 342 : return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);
914 : else // need logarithm;
915 : // might overflow as a whole, but seldom, not worth to pre-detect it
916 0 : return exp( fAm1LogX + fBm1LogY - fLogBeta);
917 : }
918 :
919 : /*
920 : x^a * (1-x)^b
921 : I_x(a,b) = ---------------- * result of ContFrac
922 : a * Beta(a,b)
923 : */
924 885 : static double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)
925 : { // like old version
926 : double a1, b1, a2, b2, fnorm, cfnew, cf;
927 885 : a1 = 1.0; b1 = 1.0;
928 885 : b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;
929 885 : if (b2 == 0.0)
930 : {
931 0 : a2 = 0.0;
932 0 : fnorm = 1.0;
933 0 : cf = 1.0;
934 : }
935 : else
936 : {
937 885 : a2 = 1.0;
938 885 : fnorm = 1.0/b2;
939 885 : cf = a2*fnorm;
940 : }
941 885 : cfnew = 1.0;
942 885 : double rm = 1.0;
943 :
944 885 : const double fMaxIter = 50000.0;
945 : // loop security, normal cases converge in less than 100 iterations.
946 : // FIXME: You will get so much iteratons for fX near mean,
947 : // I do not know a better algorithm.
948 885 : bool bfinished = false;
949 11403 : do
950 : {
951 11403 : const double apl2m = fA + 2.0*rm;
952 11403 : const double d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m);
953 11403 : const double d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0));
954 11403 : a1 = (a2+d2m*a1)*fnorm;
955 11403 : b1 = (b2+d2m*b1)*fnorm;
956 11403 : a2 = a1 + d2m1*a2*fnorm;
957 11403 : b2 = b1 + d2m1*b2*fnorm;
958 11403 : if (b2 != 0.0)
959 : {
960 11403 : fnorm = 1.0/b2;
961 11403 : cfnew = a2*fnorm;
962 11403 : bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);
963 : }
964 11403 : cf = cfnew;
965 11403 : rm += 1.0;
966 : }
967 11403 : while (rm < fMaxIter && !bfinished);
968 885 : return cf;
969 : }
970 :
971 : // cumulative distribution function, normalized
972 936 : double ScInterpreter::GetBetaDist(double fXin, double fAlpha, double fBeta)
973 : {
974 : // special cases
975 936 : if (fXin <= 0.0) // values are valid, see spec
976 9 : return 0.0;
977 927 : if (fXin >= 1.0) // values are valid, see spec
978 39 : return 1.0;
979 888 : if (fBeta == 1.0)
980 3 : return pow(fXin, fAlpha);
981 885 : if (fAlpha == 1.0)
982 : // 1.0 - pow(1.0-fX,fBeta) is not accurate enough
983 0 : return -::rtl::math::expm1(fBeta * ::rtl::math::log1p(-fXin));
984 : //FIXME: need special algorithm for fX near fP for large fA,fB
985 : double fResult;
986 : // I use always continued fraction, power series are neither
987 : // faster nor more accurate.
988 885 : double fY = (0.5-fXin)+0.5;
989 885 : double flnY = ::rtl::math::log1p(-fXin);
990 885 : double fX = fXin;
991 885 : double flnX = log(fXin);
992 885 : double fA = fAlpha;
993 885 : double fB = fBeta;
994 885 : bool bReflect = fXin > fAlpha/(fAlpha+fBeta);
995 885 : if (bReflect)
996 : {
997 156 : fA = fBeta;
998 156 : fB = fAlpha;
999 156 : fX = fY;
1000 156 : fY = fXin;
1001 156 : flnX = flnY;
1002 156 : flnY = log(fXin);
1003 : }
1004 885 : fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);
1005 885 : fResult = fResult/fA;
1006 885 : double fP = fA/(fA+fB);
1007 885 : double fQ = fB/(fA+fB);
1008 : double fTemp;
1009 885 : if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental
1010 336 : fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;
1011 : else
1012 549 : fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));
1013 885 : fResult *= fTemp;
1014 885 : if (bReflect)
1015 156 : fResult = 0.5 - fResult + 0.5;
1016 885 : if (fResult > 1.0) // ensure valid range
1017 0 : fResult = 1.0;
1018 885 : if (fResult < 0.0)
1019 0 : fResult = 0.0;
1020 885 : return fResult;
1021 : }
1022 :
1023 3 : void ScInterpreter::ScBetaDist()
1024 : {
1025 3 : sal_uInt8 nParamCount = GetByte();
1026 3 : if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547#
1027 0 : return;
1028 : double fLowerBound, fUpperBound;
1029 : double alpha, beta, x;
1030 : bool bIsCumulative;
1031 3 : if (nParamCount == 6)
1032 0 : bIsCumulative = GetBool();
1033 : else
1034 3 : bIsCumulative = true;
1035 3 : if (nParamCount >= 5)
1036 0 : fUpperBound = GetDouble();
1037 : else
1038 3 : fUpperBound = 1.0;
1039 3 : if (nParamCount >= 4)
1040 0 : fLowerBound = GetDouble();
1041 : else
1042 3 : fLowerBound = 0.0;
1043 3 : beta = GetDouble();
1044 3 : alpha = GetDouble();
1045 3 : x = GetDouble();
1046 3 : double fScale = fUpperBound - fLowerBound;
1047 3 : if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0)
1048 : {
1049 0 : PushIllegalArgument();
1050 0 : return;
1051 : }
1052 3 : if (bIsCumulative) // cumulative distribution function
1053 : {
1054 : // special cases
1055 3 : if (x < fLowerBound)
1056 : {
1057 0 : PushDouble(0.0); return; //see spec
1058 : }
1059 3 : if (x > fUpperBound)
1060 : {
1061 0 : PushDouble(1.0); return; //see spec
1062 : }
1063 : // normal cases
1064 3 : x = (x-fLowerBound)/fScale; // convert to standard form
1065 3 : PushDouble(GetBetaDist(x, alpha, beta));
1066 3 : return;
1067 : }
1068 : else // probability density function
1069 : {
1070 0 : if (x < fLowerBound || x > fUpperBound)
1071 : {
1072 0 : PushDouble(0.0);
1073 0 : return;
1074 : }
1075 0 : x = (x-fLowerBound)/fScale;
1076 0 : PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
1077 0 : return;
1078 : }
1079 : }
1080 :
1081 : /**
1082 : fdo#71008
1083 : Microsoft version has parameters in different order
1084 : Also, upper and lowerbound are optional and have default values
1085 : otherwise, function is identical with ScInterpreter::ScBetaDist()
1086 : */
1087 12 : void ScInterpreter::ScBetaDist_MS()
1088 : {
1089 12 : sal_uInt8 nParamCount = GetByte();
1090 12 : if ( !MustHaveParamCount( nParamCount, 4, 6 ) )
1091 0 : return;
1092 : double fLowerBound, fUpperBound;
1093 : double alpha, beta, x;
1094 : bool bIsCumulative;
1095 12 : if (nParamCount == 6)
1096 12 : fUpperBound = GetDouble();
1097 : else
1098 0 : fUpperBound = 1.0;
1099 12 : if (nParamCount >= 4)
1100 12 : fLowerBound = GetDouble();
1101 : else
1102 0 : fLowerBound = 0.0;
1103 12 : bIsCumulative = GetBool();
1104 12 : beta = GetDouble();
1105 12 : alpha = GetDouble();
1106 12 : x = GetDouble();
1107 12 : double fScale = fUpperBound - fLowerBound;
1108 12 : if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0)
1109 : {
1110 0 : PushIllegalArgument();
1111 0 : return;
1112 : }
1113 12 : if (bIsCumulative) // cumulative distribution function
1114 : {
1115 : // special cases
1116 6 : if (x < fLowerBound)
1117 : {
1118 0 : PushDouble(0.0); return; //see spec
1119 : }
1120 6 : if (x > fUpperBound)
1121 : {
1122 0 : PushDouble(1.0); return; //see spec
1123 : }
1124 : // normal cases
1125 6 : x = (x-fLowerBound)/fScale; // convert to standard form
1126 6 : PushDouble(GetBetaDist(x, alpha, beta));
1127 6 : return;
1128 : }
1129 : else // probability density function
1130 : {
1131 6 : if (x < fLowerBound || x > fUpperBound)
1132 : {
1133 0 : PushDouble(0.0);
1134 0 : return;
1135 : }
1136 6 : x = (x-fLowerBound)/fScale;
1137 6 : PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
1138 6 : return;
1139 : }
1140 : }
1141 :
1142 9 : void ScInterpreter::ScPhi()
1143 : {
1144 9 : PushDouble(phi(GetDouble()));
1145 9 : }
1146 :
1147 6 : void ScInterpreter::ScGauss()
1148 : {
1149 6 : PushDouble(gauss(GetDouble()));
1150 6 : }
1151 :
1152 3 : void ScInterpreter::ScFisher()
1153 : {
1154 3 : double fVal = GetDouble();
1155 3 : if (fabs(fVal) >= 1.0)
1156 0 : PushIllegalArgument();
1157 : else
1158 3 : PushDouble( ::rtl::math::atanh( fVal));
1159 3 : }
1160 :
1161 3 : void ScInterpreter::ScFisherInv()
1162 : {
1163 3 : PushDouble( tanh( GetDouble()));
1164 3 : }
1165 :
1166 6 : void ScInterpreter::ScFact()
1167 : {
1168 6 : double nVal = GetDouble();
1169 6 : if (nVal < 0.0)
1170 0 : PushIllegalArgument();
1171 : else
1172 6 : PushDouble(Fakultaet(nVal));
1173 6 : }
1174 :
1175 3 : void ScInterpreter::ScCombin()
1176 : {
1177 3 : if ( MustHaveParamCount( GetByte(), 2 ) )
1178 : {
1179 3 : double k = ::rtl::math::approxFloor(GetDouble());
1180 3 : double n = ::rtl::math::approxFloor(GetDouble());
1181 3 : if (k < 0.0 || n < 0.0 || k > n)
1182 0 : PushIllegalArgument();
1183 : else
1184 3 : PushDouble(BinomKoeff(n, k));
1185 : }
1186 3 : }
1187 :
1188 3 : void ScInterpreter::ScCombinA()
1189 : {
1190 3 : if ( MustHaveParamCount( GetByte(), 2 ) )
1191 : {
1192 3 : double k = ::rtl::math::approxFloor(GetDouble());
1193 3 : double n = ::rtl::math::approxFloor(GetDouble());
1194 3 : if (k < 0.0 || n < 0.0 || k > n)
1195 0 : PushIllegalArgument();
1196 : else
1197 3 : PushDouble(BinomKoeff(n + k - 1, k));
1198 : }
1199 3 : }
1200 :
1201 3 : void ScInterpreter::ScPermut()
1202 : {
1203 3 : if ( MustHaveParamCount( GetByte(), 2 ) )
1204 : {
1205 3 : double k = ::rtl::math::approxFloor(GetDouble());
1206 3 : double n = ::rtl::math::approxFloor(GetDouble());
1207 3 : if (n < 0.0 || k < 0.0 || k > n)
1208 0 : PushIllegalArgument();
1209 3 : else if (k == 0.0)
1210 0 : PushInt(1); // (n! / (n - 0)!) == 1
1211 : else
1212 : {
1213 3 : double nVal = n;
1214 9 : for (sal_uLong i = (sal_uLong)k-1; i >= 1; i--)
1215 6 : nVal *= n-(double)i;
1216 3 : PushDouble(nVal);
1217 : }
1218 : }
1219 3 : }
1220 :
1221 6 : void ScInterpreter::ScPermutationA()
1222 : {
1223 6 : if ( MustHaveParamCount( GetByte(), 2 ) )
1224 : {
1225 6 : double k = ::rtl::math::approxFloor(GetDouble());
1226 6 : double n = ::rtl::math::approxFloor(GetDouble());
1227 6 : if (n < 0.0 || k < 0.0 || k > n)
1228 0 : PushIllegalArgument();
1229 : else
1230 6 : PushDouble(pow(n,k));
1231 : }
1232 6 : }
1233 :
1234 12 : double ScInterpreter::GetBinomDistPMF(double x, double n, double p)
1235 : // used in ScB and ScBinomDist
1236 : // preconditions: 0.0 <= x <= n, 0.0 < p < 1.0; x,n integral although double
1237 : {
1238 12 : double q = (0.5 - p) + 0.5;
1239 12 : double fFactor = pow(q, n);
1240 12 : if (fFactor <=::std::numeric_limits<double>::min())
1241 : {
1242 0 : fFactor = pow(p, n);
1243 0 : if (fFactor <= ::std::numeric_limits<double>::min())
1244 0 : return GetBetaDistPDF(p, x+1.0, n-x+1.0)/(n+1.0);
1245 : else
1246 : {
1247 0 : sal_uInt32 max = static_cast<sal_uInt32>(n - x);
1248 0 : for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
1249 0 : fFactor *= (n-i)/(i+1)*q/p;
1250 0 : return fFactor;
1251 : }
1252 : }
1253 : else
1254 : {
1255 12 : sal_uInt32 max = static_cast<sal_uInt32>(x);
1256 45 : for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
1257 33 : fFactor *= (n-i)/(i+1)*p/q;
1258 12 : return fFactor;
1259 : }
1260 : }
1261 :
1262 9 : double lcl_GetBinomDistRange(double n, double xs,double xe,
1263 : double fFactor /* q^n */, double p, double q)
1264 : //preconditions: 0.0 <= xs < xe <= n; xs,xe,n integral although double
1265 : {
1266 : sal_uInt32 i;
1267 : double fSum;
1268 : // skip summands index 0 to xs-1, start sum with index xs
1269 9 : sal_uInt32 nXs = static_cast<sal_uInt32>( xs );
1270 9 : for (i = 1; i <= nXs && fFactor > 0.0; i++)
1271 0 : fFactor *= (n-i+1)/i * p/q;
1272 9 : fSum = fFactor; // Summand xs
1273 9 : sal_uInt32 nXe = static_cast<sal_uInt32>(xe);
1274 36 : for (i = nXs+1; i <= nXe && fFactor > 0.0; i++)
1275 : {
1276 27 : fFactor *= (n-i+1)/i * p/q;
1277 27 : fSum += fFactor;
1278 : }
1279 9 : return (fSum>1.0) ? 1.0 : fSum;
1280 : }
1281 :
1282 3 : void ScInterpreter::ScB()
1283 : {
1284 3 : sal_uInt8 nParamCount = GetByte();
1285 3 : if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1286 3 : return ;
1287 3 : if (nParamCount == 3) // mass function
1288 : {
1289 3 : double x = ::rtl::math::approxFloor(GetDouble());
1290 3 : double p = GetDouble();
1291 3 : double n = ::rtl::math::approxFloor(GetDouble());
1292 3 : if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1293 0 : PushIllegalArgument();
1294 3 : else if (p == 0.0)
1295 0 : PushDouble( (x == 0.0) ? 1.0 : 0.0 );
1296 3 : else if ( p == 1.0)
1297 0 : PushDouble( (x == n) ? 1.0 : 0.0);
1298 : else
1299 3 : PushDouble(GetBinomDistPMF(x,n,p));
1300 : }
1301 : else
1302 : { // nParamCount == 4
1303 0 : double xe = ::rtl::math::approxFloor(GetDouble());
1304 0 : double xs = ::rtl::math::approxFloor(GetDouble());
1305 0 : double p = GetDouble();
1306 0 : double n = ::rtl::math::approxFloor(GetDouble());
1307 0 : double q = (0.5 - p) + 0.5;
1308 0 : bool bIsValidX = ( 0.0 <= xs && xs <= xe && xe <= n);
1309 0 : if ( bIsValidX && 0.0 < p && p < 1.0)
1310 : {
1311 0 : if (xs == xe) // mass function
1312 0 : PushDouble(GetBinomDistPMF(xs,n,p));
1313 : else
1314 : {
1315 0 : double fFactor = pow(q, n);
1316 0 : if (fFactor > ::std::numeric_limits<double>::min())
1317 0 : PushDouble(lcl_GetBinomDistRange(n,xs,xe,fFactor,p,q));
1318 : else
1319 : {
1320 0 : fFactor = pow(p, n);
1321 0 : if (fFactor > ::std::numeric_limits<double>::min())
1322 : {
1323 : // sum from j=xs to xe {(n choose j) * p^j * q^(n-j)}
1324 : // = sum from i = n-xe to n-xs { (n choose i) * q^i * p^(n-i)}
1325 0 : PushDouble(lcl_GetBinomDistRange(n,n-xe,n-xs,fFactor,q,p));
1326 : }
1327 : else
1328 0 : PushDouble(GetBetaDist(q,n-xe,xe+1.0)-GetBetaDist(q,n-xs+1,xs) );
1329 : }
1330 0 : }
1331 : }
1332 : else
1333 : {
1334 0 : if ( bIsValidX ) // not(0<p<1)
1335 : {
1336 0 : if ( p == 0.0 )
1337 0 : PushDouble( (xs == 0.0) ? 1.0 : 0.0 );
1338 0 : else if ( p == 1.0 )
1339 0 : PushDouble( (xe == n) ? 1.0 : 0.0 );
1340 : else
1341 0 : PushIllegalArgument();
1342 : }
1343 : else
1344 0 : PushIllegalArgument();
1345 : }
1346 : }
1347 : }
1348 :
1349 18 : void ScInterpreter::ScBinomDist()
1350 : {
1351 18 : if ( MustHaveParamCount( GetByte(), 4 ) )
1352 : {
1353 18 : bool bIsCum = GetBool(); // false=mass function; true=cumulative
1354 18 : double p = GetDouble();
1355 18 : double n = ::rtl::math::approxFloor(GetDouble());
1356 18 : double x = ::rtl::math::approxFloor(GetDouble());
1357 18 : double q = (0.5 - p) + 0.5; // get one bit more for p near 1.0
1358 18 : if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1359 : {
1360 0 : PushIllegalArgument();
1361 0 : return;
1362 : }
1363 18 : if ( p == 0.0)
1364 : {
1365 0 : PushDouble( (x==0.0 || bIsCum) ? 1.0 : 0.0 );
1366 0 : return;
1367 : }
1368 18 : if ( p == 1.0)
1369 : {
1370 0 : PushDouble( (x==n) ? 1.0 : 0.0);
1371 0 : return;
1372 : }
1373 18 : if (!bIsCum)
1374 9 : PushDouble( GetBinomDistPMF(x,n,p));
1375 : else
1376 : {
1377 9 : if (x == n)
1378 0 : PushDouble(1.0);
1379 : else
1380 : {
1381 9 : double fFactor = pow(q, n);
1382 9 : if (x == 0.0)
1383 0 : PushDouble(fFactor);
1384 9 : else if (fFactor <= ::std::numeric_limits<double>::min())
1385 : {
1386 0 : fFactor = pow(p, n);
1387 0 : if (fFactor <= ::std::numeric_limits<double>::min())
1388 0 : PushDouble(GetBetaDist(q,n-x,x+1.0));
1389 : else
1390 : {
1391 0 : if (fFactor > fMachEps)
1392 : {
1393 0 : double fSum = 1.0 - fFactor;
1394 0 : sal_uInt32 max = static_cast<sal_uInt32> (n - x) - 1;
1395 0 : for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
1396 : {
1397 0 : fFactor *= (n-i)/(i+1)*q/p;
1398 0 : fSum -= fFactor;
1399 : }
1400 0 : PushDouble( (fSum < 0.0) ? 0.0 : fSum );
1401 : }
1402 : else
1403 0 : PushDouble(lcl_GetBinomDistRange(n,n-x,n,fFactor,q,p));
1404 : }
1405 : }
1406 : else
1407 9 : PushDouble( lcl_GetBinomDistRange(n,0.0,x,fFactor,p,q)) ;
1408 : }
1409 : }
1410 : }
1411 : }
1412 :
1413 6 : void ScInterpreter::ScCritBinom()
1414 : {
1415 6 : if ( MustHaveParamCount( GetByte(), 3 ) )
1416 : {
1417 6 : double alpha = GetDouble();
1418 6 : double p = GetDouble();
1419 6 : double n = ::rtl::math::approxFloor(GetDouble());
1420 6 : if (n < 0.0 || alpha <= 0.0 || alpha >= 1.0 || p < 0.0 || p > 1.0)
1421 0 : PushIllegalArgument();
1422 : else
1423 : {
1424 : double fFactor;
1425 6 : double q = (0.5 - p) + 0.5; // get one bit more for p near 1.0
1426 6 : if ( q > p ) // work from the side where the cumulative curve is
1427 : {
1428 : // work from 0 upwards
1429 3 : fFactor = pow(q,n);
1430 3 : if (fFactor > ::std::numeric_limits<double>::min())
1431 : {
1432 3 : double fSum = fFactor;
1433 3 : sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1434 12 : for (i = 0; i < max && fSum < alpha; i++)
1435 : {
1436 9 : fFactor *= (n-i)/(i+1)*p/q;
1437 9 : fSum += fFactor;
1438 : }
1439 3 : PushDouble(i);
1440 : }
1441 : else
1442 : {
1443 : // accumulate BinomDist until accumulated BinomDist reaches alpha
1444 0 : double fSum = 0.0;
1445 0 : sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1446 0 : for (i = 0; i < max && fSum < alpha; i++)
1447 : {
1448 0 : const double x = GetBetaDistPDF( p, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
1449 0 : if ( !nGlobalError )
1450 : {
1451 0 : fSum += x;
1452 : }
1453 : else
1454 : {
1455 0 : PushNoValue();
1456 0 : return;
1457 : }
1458 : }
1459 0 : PushDouble( i - 1 );
1460 : }
1461 : }
1462 : else
1463 : {
1464 : // work from n backwards
1465 3 : fFactor = pow(p, n);
1466 3 : if (fFactor > ::std::numeric_limits<double>::min())
1467 : {
1468 3 : double fSum = 1.0 - fFactor;
1469 3 : sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1470 6 : for (i = 0; i < max && fSum >= alpha; i++)
1471 : {
1472 3 : fFactor *= (n-i)/(i+1)*q/p;
1473 3 : fSum -= fFactor;
1474 : }
1475 3 : PushDouble(n-i);
1476 : }
1477 : else
1478 : {
1479 : // accumulate BinomDist until accumulated BinomDist reaches alpha
1480 0 : double fSum = 0.0;
1481 0 : sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1482 0 : alpha = 1 - alpha;
1483 0 : for (i = 0; i < max && fSum < alpha; i++)
1484 : {
1485 0 : const double x = GetBetaDistPDF( q, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
1486 0 : if ( !nGlobalError )
1487 : {
1488 0 : fSum += x;
1489 : }
1490 : else
1491 : {
1492 0 : PushNoValue();
1493 0 : return;
1494 : }
1495 : }
1496 0 : PushDouble( n - i + 1 );
1497 : }
1498 : }
1499 : }
1500 : }
1501 : }
1502 :
1503 3 : void ScInterpreter::ScNegBinomDist()
1504 : {
1505 3 : if ( MustHaveParamCount( GetByte(), 3 ) )
1506 : {
1507 3 : double p = GetDouble(); // p
1508 3 : double r = GetDouble(); // r
1509 3 : double x = GetDouble(); // x
1510 3 : if (r < 0.0 || x < 0.0 || p < 0.0 || p > 1.0)
1511 0 : PushIllegalArgument();
1512 : else
1513 : {
1514 3 : double q = 1.0 - p;
1515 3 : double fFactor = pow(p,r);
1516 6 : for (double i = 0.0; i < x; i++)
1517 3 : fFactor *= (i+r)/(i+1.0)*q;
1518 3 : PushDouble(fFactor);
1519 : }
1520 : }
1521 3 : }
1522 :
1523 12 : void ScInterpreter::ScNegBinomDist_MS()
1524 : {
1525 12 : if ( MustHaveParamCount( GetByte(), 4 ) )
1526 : {
1527 12 : bool bCumulative = GetBool();
1528 12 : double p = GetDouble(); // p
1529 12 : double r = GetDouble(); // r
1530 12 : double x = GetDouble(); // x
1531 12 : if ( r < 0.0 || x < 0.0 || p < 0.0 || p > 1.0 )
1532 0 : PushIllegalArgument();
1533 : else
1534 : {
1535 12 : double q = 1.0 - p;
1536 12 : if ( bCumulative )
1537 6 : PushDouble( 1.0 - GetBetaDist( q, x + 1, r ) );
1538 : else
1539 : {
1540 6 : double fFactor = pow( p, r );
1541 15 : for ( double i = 0.0; i < x; i++ )
1542 9 : fFactor *= ( i + r ) / ( i + 1.0 ) * q;
1543 6 : PushDouble( fFactor );
1544 : }
1545 : }
1546 : }
1547 12 : }
1548 :
1549 18 : void ScInterpreter::ScNormDist( int nMinParamCount )
1550 : {
1551 18 : sal_uInt8 nParamCount = GetByte();
1552 18 : if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
1553 0 : return;
1554 18 : bool bCumulative = nParamCount != 4 || GetBool();
1555 18 : double sigma = GetDouble(); // standard deviation
1556 18 : double mue = GetDouble(); // mean
1557 18 : double x = GetDouble(); // x
1558 18 : if (sigma <= 0.0)
1559 : {
1560 0 : PushIllegalArgument();
1561 0 : return;
1562 : }
1563 18 : if (bCumulative)
1564 9 : PushDouble(integralPhi((x-mue)/sigma));
1565 : else
1566 9 : PushDouble(phi((x-mue)/sigma)/sigma);
1567 : }
1568 :
1569 12 : void ScInterpreter::ScLogNormDist( int nMinParamCount ) //expanded, see #i100119# and fdo72158
1570 : {
1571 12 : sal_uInt8 nParamCount = GetByte();
1572 12 : if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
1573 0 : return;
1574 12 : bool bCumulative = nParamCount != 4 || GetBool(); // cumulative
1575 12 : double sigma = nParamCount >= 3 ? GetDouble() : 1.0; // standard deviation
1576 12 : double mue = nParamCount >= 2 ? GetDouble() : 0.0; // mean
1577 12 : double x = GetDouble(); // x
1578 12 : if (sigma <= 0.0)
1579 : {
1580 0 : PushIllegalArgument();
1581 0 : return;
1582 : }
1583 12 : if (bCumulative)
1584 : { // cumulative
1585 9 : if (x <= 0.0)
1586 0 : PushDouble(0.0);
1587 : else
1588 9 : PushDouble(integralPhi((log(x)-mue)/sigma));
1589 : }
1590 : else
1591 : { // density
1592 3 : if (x <= 0.0)
1593 0 : PushIllegalArgument();
1594 : else
1595 3 : PushDouble(phi((log(x)-mue)/sigma)/sigma/x);
1596 : }
1597 : }
1598 :
1599 3 : void ScInterpreter::ScStdNormDist()
1600 : {
1601 3 : PushDouble(integralPhi(GetDouble()));
1602 3 : }
1603 :
1604 12 : void ScInterpreter::ScStdNormDist_MS()
1605 : {
1606 12 : sal_uInt8 nParamCount = GetByte();
1607 12 : if ( !MustHaveParamCount( nParamCount, 2 ) )
1608 12 : return;
1609 12 : bool bCumulative = GetBool(); // cumulative
1610 12 : double x = GetDouble(); // x
1611 :
1612 12 : if ( bCumulative )
1613 6 : PushDouble( integralPhi( x ) );
1614 : else
1615 6 : PushDouble( exp( - pow( x, 2 ) / 2 ) / sqrt( 2 * F_PI ) );
1616 : }
1617 :
1618 12 : void ScInterpreter::ScExpDist()
1619 : {
1620 12 : if ( MustHaveParamCount( GetByte(), 3 ) )
1621 : {
1622 12 : double kum = GetDouble(); // 0 oder 1
1623 12 : double lambda = GetDouble(); // lambda
1624 12 : double x = GetDouble(); // x
1625 12 : if (lambda <= 0.0)
1626 0 : PushIllegalArgument();
1627 12 : else if (kum == 0.0) // Dichte
1628 : {
1629 3 : if (x >= 0.0)
1630 3 : PushDouble(lambda * exp(-lambda*x));
1631 : else
1632 0 : PushInt(0);
1633 : }
1634 : else // Verteilung
1635 : {
1636 9 : if (x > 0.0)
1637 9 : PushDouble(1.0 - exp(-lambda*x));
1638 : else
1639 0 : PushInt(0);
1640 : }
1641 : }
1642 12 : }
1643 :
1644 3 : void ScInterpreter::ScTDist()
1645 : {
1646 3 : if ( !MustHaveParamCount( GetByte(), 3 ) )
1647 0 : return;
1648 3 : double fFlag = ::rtl::math::approxFloor(GetDouble());
1649 3 : double fDF = ::rtl::math::approxFloor(GetDouble());
1650 3 : double T = GetDouble();
1651 3 : if (fDF < 1.0 || T < 0.0 || (fFlag != 1.0 && fFlag != 2.0) )
1652 : {
1653 0 : PushIllegalArgument();
1654 0 : return;
1655 : }
1656 3 : PushDouble( GetTDist( T, fDF, ( int )fFlag ) );
1657 : }
1658 :
1659 12 : void ScInterpreter::ScTDist_T( int nTails )
1660 : {
1661 12 : if ( !MustHaveParamCount( GetByte(), 2 ) )
1662 0 : return;
1663 12 : double fDF = ::rtl::math::approxFloor( GetDouble() );
1664 12 : double T = GetDouble();
1665 12 : if ( fDF < 1.0 || T < 0.0 )
1666 : {
1667 0 : PushIllegalArgument();
1668 0 : return;
1669 : }
1670 12 : PushDouble( GetTDist( T, fDF, nTails ) );
1671 : }
1672 :
1673 9 : void ScInterpreter::ScTDist_MS()
1674 : {
1675 9 : if ( !MustHaveParamCount( GetByte(), 3 ) )
1676 0 : return;
1677 9 : bool bCumulative = GetBool();
1678 9 : double fDF = ::rtl::math::approxFloor( GetDouble() );
1679 9 : double T = GetDouble();
1680 9 : if ( fDF < 1.0 )
1681 : {
1682 0 : PushIllegalArgument();
1683 0 : return;
1684 : }
1685 9 : PushDouble( GetTDist( T, fDF, ( bCumulative ? 4 : 3 ) ) );
1686 : }
1687 :
1688 9 : void ScInterpreter::ScFDist()
1689 : {
1690 9 : if ( !MustHaveParamCount( GetByte(), 3 ) )
1691 0 : return;
1692 9 : double fF2 = ::rtl::math::approxFloor(GetDouble());
1693 9 : double fF1 = ::rtl::math::approxFloor(GetDouble());
1694 9 : double fF = GetDouble();
1695 9 : if (fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10)
1696 : {
1697 0 : PushIllegalArgument();
1698 0 : return;
1699 : }
1700 9 : PushDouble(GetFDist(fF, fF1, fF2));
1701 : }
1702 :
1703 12 : void ScInterpreter::ScFDist_LT()
1704 : {
1705 12 : int nParamCount = GetByte();
1706 12 : if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1707 0 : return;
1708 : bool bCum;
1709 12 : if ( nParamCount == 3 )
1710 0 : bCum = true;
1711 12 : else if ( IsMissing() )
1712 : {
1713 0 : bCum = true;
1714 0 : Pop();
1715 : }
1716 : else
1717 12 : bCum = GetBool();
1718 12 : double fF2 = ::rtl::math::approxFloor( GetDouble() );
1719 12 : double fF1 = ::rtl::math::approxFloor( GetDouble() );
1720 12 : double fF = GetDouble();
1721 12 : if ( fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 )
1722 : {
1723 0 : PushIllegalArgument();
1724 0 : return;
1725 : }
1726 12 : if ( bCum )
1727 : {
1728 : // left tail cumulative distribution
1729 6 : PushDouble( 1.0 - GetFDist( fF, fF1, fF2 ) );
1730 : }
1731 : else
1732 : {
1733 : // probability density function
1734 6 : PushDouble( pow( fF1 / fF2, fF1 / 2 ) * pow( fF, ( fF1 / 2 ) - 1 ) /
1735 12 : ( pow( ( 1 + ( fF * fF1 / fF2 ) ), ( fF1 + fF2 ) / 2 ) *
1736 12 : GetBeta( fF1 / 2, fF2 / 2 ) ) );
1737 : }
1738 : }
1739 :
1740 9 : void ScInterpreter::ScChiDist()
1741 : {
1742 : double fResult;
1743 9 : if ( !MustHaveParamCount( GetByte(), 2 ) )
1744 0 : return;
1745 9 : double fDF = ::rtl::math::approxFloor(GetDouble());
1746 9 : double fChi = GetDouble();
1747 9 : if (fDF < 1.0) // x<=0 returns 1, see ODFF 6.17.10
1748 : {
1749 0 : PushIllegalArgument();
1750 0 : return;
1751 : }
1752 9 : fResult = GetChiDist( fChi, fDF);
1753 9 : if (nGlobalError)
1754 : {
1755 0 : PushError( nGlobalError);
1756 0 : return;
1757 : }
1758 9 : PushDouble(fResult);
1759 : }
1760 :
1761 12 : void ScInterpreter::ScWeibull()
1762 : {
1763 12 : if ( MustHaveParamCount( GetByte(), 4 ) )
1764 : {
1765 12 : double kum = GetDouble(); // 0 oder 1
1766 12 : double beta = GetDouble(); // beta
1767 12 : double alpha = GetDouble(); // alpha
1768 12 : double x = GetDouble(); // x
1769 12 : if (alpha <= 0.0 || beta <= 0.0 || x < 0.0)
1770 0 : PushIllegalArgument();
1771 12 : else if (kum == 0.0) // Dichte
1772 3 : PushDouble(alpha/pow(beta,alpha)*pow(x,alpha-1.0)*
1773 3 : exp(-pow(x/beta,alpha)));
1774 : else // Verteilung
1775 9 : PushDouble(1.0 - exp(-pow(x/beta,alpha)));
1776 : }
1777 12 : }
1778 :
1779 12 : void ScInterpreter::ScPoissonDist()
1780 : {
1781 12 : sal_uInt8 nParamCount = GetByte();
1782 12 : if ( MustHaveParamCount( nParamCount, 2, 3 ) )
1783 : {
1784 12 : bool bCumulative = nParamCount != 3 || GetBool(); // default cumulative
1785 12 : double lambda = GetDouble(); // Mean
1786 12 : double x = ::rtl::math::approxFloor(GetDouble()); // discrete distribution
1787 12 : if (lambda < 0.0 || x < 0.0)
1788 0 : PushIllegalArgument();
1789 12 : else if (!bCumulative) // Probability mass function
1790 : {
1791 3 : if (lambda == 0.0)
1792 0 : PushInt(0);
1793 : else
1794 : {
1795 3 : if (lambda >712) // underflow in exp(-lambda)
1796 : { // accuracy 11 Digits
1797 0 : PushDouble( exp(x*log(lambda)-lambda-GetLogGamma(x+1.0)));
1798 : }
1799 : else
1800 : {
1801 3 : double fPoissonVar = 1.0;
1802 6 : for ( double f = 0.0; f < x; ++f )
1803 3 : fPoissonVar *= lambda / ( f + 1.0 );
1804 3 : PushDouble( fPoissonVar * exp( -lambda ) );
1805 : }
1806 : }
1807 : }
1808 : else // Cumulative distribution function
1809 : {
1810 9 : if (lambda == 0.0)
1811 0 : PushInt(1);
1812 : else
1813 : {
1814 9 : if (lambda > 712 ) // underflow in exp(-lambda)
1815 : { // accuracy 12 Digits
1816 0 : PushDouble(GetUpRegIGamma(x+1.0,lambda));
1817 : }
1818 : else
1819 : {
1820 9 : if (x >= 936.0) // result is always undistinghable from 1
1821 0 : PushDouble (1.0);
1822 : else
1823 : {
1824 9 : double fSummand = exp(-lambda);
1825 9 : double fSum = fSummand;
1826 9 : int nEnd = sal::static_int_cast<int>( x );
1827 372 : for (int i = 1; i <= nEnd; i++)
1828 : {
1829 363 : fSummand = (fSummand * lambda)/(double)i;
1830 363 : fSum += fSummand;
1831 : }
1832 9 : PushDouble(fSum);
1833 : }
1834 : }
1835 : }
1836 : }
1837 : }
1838 12 : }
1839 :
1840 : /** Local function used in the calculation of the hypergeometric distribution.
1841 : */
1842 135 : static void lcl_PutFactorialElements( ::std::vector< double >& cn, double fLower, double fUpper, double fBase )
1843 : {
1844 375 : for ( double i = fLower; i <= fUpper; ++i )
1845 : {
1846 240 : double fVal = fBase - i;
1847 240 : if ( fVal > 1.0 )
1848 228 : cn.push_back( fVal );
1849 : }
1850 135 : }
1851 :
1852 : /** Calculates a value of the hypergeometric distribution.
1853 :
1854 : @author Kohei Yoshida <kohei@openoffice.org>
1855 :
1856 : @see #i47296#
1857 :
1858 : */
1859 3 : void ScInterpreter::ScHypGeomDist()
1860 : {
1861 3 : if ( !MustHaveParamCount( GetByte(), 4 ) )
1862 0 : return;
1863 :
1864 3 : double N = ::rtl::math::approxFloor(GetDouble());
1865 3 : double M = ::rtl::math::approxFloor(GetDouble());
1866 3 : double n = ::rtl::math::approxFloor(GetDouble());
1867 3 : double x = ::rtl::math::approxFloor(GetDouble());
1868 :
1869 3 : if( (x < 0.0) || (n < x) || (M < x) || (N < n) || (N < M) || (x < n - N + M) )
1870 : {
1871 0 : PushIllegalArgument();
1872 0 : return;
1873 : }
1874 :
1875 3 : PushDouble( GetHypGeomDist( x, n, M, N ) );
1876 : }
1877 :
1878 : /** Calculates a value of the hypergeometric distribution (Excel 2010 function).
1879 :
1880 : This function has an extra argument bCumulative as compared to ScHypGeomDist(),
1881 : which only calculates the non-cumulative distribution.
1882 :
1883 : @see fdo#71722
1884 : */
1885 12 : void ScInterpreter::ScHypGeomDist_MS()
1886 : {
1887 12 : if ( !MustHaveParamCount( GetByte(), 5 ) )
1888 0 : return;
1889 :
1890 12 : bool bCumulative = GetBool();
1891 12 : double N = ::rtl::math::approxFloor(GetDouble());
1892 12 : double M = ::rtl::math::approxFloor(GetDouble());
1893 12 : double n = ::rtl::math::approxFloor(GetDouble());
1894 12 : double x = ::rtl::math::approxFloor(GetDouble());
1895 :
1896 12 : if( (x < 0.0) || (n < x) || (M < x) || (N < n) || (N < M) || (x < n - N + M) )
1897 : {
1898 0 : PushIllegalArgument();
1899 0 : return;
1900 : }
1901 :
1902 12 : if ( bCumulative )
1903 : {
1904 6 : double fVal = 0.0;
1905 :
1906 24 : for ( int i = 0; i <= x && !nGlobalError; i++ )
1907 18 : fVal += GetHypGeomDist( i, n, M, N );
1908 :
1909 6 : PushDouble( fVal );
1910 : }
1911 : else
1912 6 : PushDouble( GetHypGeomDist( x, n, M, N ) );
1913 : }
1914 :
1915 : /** Calculates a value of the hypergeometric distribution.
1916 :
1917 : The algorithm is designed to avoid unnecessary multiplications and division
1918 : by expanding all factorial elements (9 of them). It is done by excluding
1919 : those ranges that overlap in the numerator and the denominator. This allows
1920 : for a fast calculation for large values which would otherwise cause an overflow
1921 : in the intermediate values.
1922 :
1923 : @author Kohei Yoshida <kohei@openoffice.org>
1924 :
1925 : @see #i47296#
1926 :
1927 : */
1928 27 : double ScInterpreter::GetHypGeomDist( double x, double n, double M, double N )
1929 : {
1930 27 : const size_t nMaxArraySize = 500000; // arbitrary max array size
1931 :
1932 : typedef ::std::vector< double > HypContainer;
1933 54 : HypContainer cnNumer, cnDenom;
1934 :
1935 27 : size_t nEstContainerSize = static_cast<size_t>( x + ::std::min( n, M ) );
1936 27 : size_t nMaxSize = ::std::min( cnNumer.max_size(), nMaxArraySize );
1937 27 : if ( nEstContainerSize > nMaxSize )
1938 : {
1939 0 : PushNoValue();
1940 0 : return 0;
1941 : }
1942 27 : cnNumer.reserve( nEstContainerSize + 10 );
1943 27 : cnDenom.reserve( nEstContainerSize + 10 );
1944 :
1945 : // Trim coefficient C first
1946 27 : double fCNumVarUpper = N - n - M + x - 1.0;
1947 27 : double fCDenomVarLower = 1.0;
1948 27 : if ( N - n - M + x >= M - x + 1.0 )
1949 : {
1950 12 : fCNumVarUpper = M - x - 1.0;
1951 12 : fCDenomVarLower = N - n - 2.0*(M - x) + 1.0;
1952 : }
1953 :
1954 : #if OSL_DEBUG_LEVEL > 0
1955 : double fCNumLower = N - n - fCNumVarUpper;
1956 : #endif
1957 27 : double fCDenomUpper = N - n - M + x + 1.0 - fCDenomVarLower;
1958 :
1959 27 : double fDNumVarLower = n - M;
1960 :
1961 27 : if ( n >= M + 1.0 )
1962 : {
1963 0 : if ( N - M < n + 1.0 )
1964 : {
1965 : // Case 1
1966 :
1967 0 : if ( N - n < n + 1.0 )
1968 : {
1969 : // no overlap
1970 0 : lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1971 0 : lcl_PutFactorialElements( cnDenom, 0.0, N - n - 1.0, N );
1972 : }
1973 : else
1974 : {
1975 : // overlap
1976 : OSL_ENSURE( fCNumLower < n + 1.0, "ScHypGeomDist: wrong assertion" );
1977 0 : lcl_PutFactorialElements( cnNumer, N - 2.0*n, fCNumVarUpper, N - n );
1978 0 : lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1979 : }
1980 :
1981 : OSL_ENSURE( fCDenomUpper <= N - M, "ScHypGeomDist: wrong assertion" );
1982 :
1983 0 : if ( fCDenomUpper < n - x + 1.0 )
1984 : // no overlap
1985 0 : lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
1986 : else
1987 : {
1988 : // overlap
1989 0 : lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1990 :
1991 0 : fCDenomUpper = n - x;
1992 0 : fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1993 : }
1994 : }
1995 : else
1996 : {
1997 : // Case 2
1998 :
1999 0 : if ( n > M - 1.0 )
2000 : {
2001 : // no overlap
2002 0 : lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
2003 0 : lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
2004 : }
2005 : else
2006 : {
2007 0 : lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
2008 0 : lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
2009 : }
2010 :
2011 : OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
2012 :
2013 0 : if ( fCDenomUpper < n - x + 1.0 )
2014 : // no overlap
2015 0 : lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - n + x, N - M + 1.0 );
2016 : else
2017 : {
2018 0 : lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
2019 0 : fCDenomUpper = n - x;
2020 0 : fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
2021 : }
2022 : }
2023 :
2024 : OSL_ENSURE( fCDenomUpper <= M, "ScHypGeomDist: wrong assertion" );
2025 : }
2026 : else
2027 : {
2028 27 : if ( N - M < M + 1.0 )
2029 : {
2030 : // Case 3
2031 :
2032 15 : if ( N - n < M + 1.0 )
2033 : {
2034 : // No overlap
2035 0 : lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
2036 0 : lcl_PutFactorialElements( cnDenom, 0.0, N - M - 1.0, N );
2037 : }
2038 : else
2039 : {
2040 15 : lcl_PutFactorialElements( cnNumer, N - n - M, fCNumVarUpper, N - n );
2041 15 : lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
2042 : }
2043 :
2044 15 : if ( n - x + 1.0 > fCDenomUpper )
2045 : // No overlap
2046 0 : lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
2047 : else
2048 : {
2049 : // Overlap
2050 15 : lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
2051 :
2052 15 : fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
2053 15 : fCDenomUpper = n - x;
2054 : }
2055 : }
2056 : else
2057 : {
2058 : // Case 4
2059 :
2060 : OSL_ENSURE( M >= n - x, "ScHypGeomDist: wrong assertion" );
2061 : OSL_ENSURE( M - x <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
2062 :
2063 12 : if ( N - n < N - M + 1.0 )
2064 : {
2065 : // No overlap
2066 0 : lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
2067 0 : lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
2068 : }
2069 : else
2070 : {
2071 : // Overlap
2072 : OSL_ENSURE( fCNumLower <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
2073 12 : lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
2074 12 : lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
2075 : }
2076 :
2077 12 : if ( n - x + 1.0 > fCDenomUpper )
2078 : // No overlap
2079 0 : lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - n + x, N - M + 1.0 );
2080 12 : else if ( M >= fCDenomUpper )
2081 : {
2082 12 : lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
2083 :
2084 12 : fCDenomUpper = n - x;
2085 12 : fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
2086 : }
2087 : else
2088 : {
2089 : OSL_ENSURE( M <= fCDenomUpper, "ScHypGeomDist: wrong assertion" );
2090 0 : lcl_PutFactorialElements( cnDenom, fCDenomVarLower, N - n - 2.0*M + x,
2091 0 : N - n - M + x + 1.0 );
2092 :
2093 0 : fCDenomUpper = n - x;
2094 0 : fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
2095 : }
2096 : }
2097 :
2098 : OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
2099 :
2100 27 : fDNumVarLower = 0.0;
2101 : }
2102 :
2103 27 : double nDNumVarUpper = fCDenomUpper < x + 1.0 ? n - x - 1.0 : n - fCDenomUpper - 1.0;
2104 27 : double nDDenomVarLower = fCDenomUpper < x + 1.0 ? fCDenomVarLower : N - n - M + 1.0;
2105 27 : lcl_PutFactorialElements( cnNumer, fDNumVarLower, nDNumVarUpper, n );
2106 27 : lcl_PutFactorialElements( cnDenom, nDDenomVarLower, N - n - M + x, N - n - M + x + 1.0 );
2107 :
2108 27 : ::std::sort( cnNumer.begin(), cnNumer.end() );
2109 27 : ::std::sort( cnDenom.begin(), cnDenom.end() );
2110 27 : HypContainer::reverse_iterator it1 = cnNumer.rbegin(), it1End = cnNumer.rend();
2111 27 : HypContainer::reverse_iterator it2 = cnDenom.rbegin(), it2End = cnDenom.rend();
2112 :
2113 27 : double fFactor = 1.0;
2114 174 : for ( ; it1 != it1End || it2 != it2End; )
2115 : {
2116 120 : double fEnum = 1.0, fDenom = 1.0;
2117 120 : if ( it1 != it1End )
2118 120 : fEnum = *it1++;
2119 120 : if ( it2 != it2End )
2120 108 : fDenom = *it2++;
2121 120 : fFactor *= fEnum / fDenom;
2122 : }
2123 :
2124 54 : return fFactor;
2125 : }
2126 :
2127 9 : void ScInterpreter::ScGammaDist( int nMinParamCount )
2128 : {
2129 9 : sal_uInt8 nParamCount = GetByte();
2130 9 : if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
2131 9 : return;
2132 : bool bCumulative;
2133 9 : if (nParamCount == 4)
2134 9 : bCumulative = GetBool();
2135 : else
2136 0 : bCumulative = true;
2137 9 : double fBeta = GetDouble(); // scale
2138 9 : double fAlpha = GetDouble(); // shape
2139 9 : double fX = GetDouble(); // x
2140 9 : if (fAlpha <= 0.0 || fBeta <= 0.0)
2141 0 : PushIllegalArgument();
2142 : else
2143 : {
2144 9 : if (bCumulative) // distribution
2145 9 : PushDouble( GetGammaDist( fX, fAlpha, fBeta));
2146 : else // density
2147 0 : PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta));
2148 : }
2149 : }
2150 :
2151 9 : void ScInterpreter::ScNormInv()
2152 : {
2153 9 : if ( MustHaveParamCount( GetByte(), 3 ) )
2154 : {
2155 9 : double sigma = GetDouble();
2156 9 : double mue = GetDouble();
2157 9 : double x = GetDouble();
2158 9 : if (sigma <= 0.0 || x < 0.0 || x > 1.0)
2159 0 : PushIllegalArgument();
2160 9 : else if (x == 0.0 || x == 1.0)
2161 0 : PushNoValue();
2162 : else
2163 9 : PushDouble(gaussinv(x)*sigma + mue);
2164 : }
2165 9 : }
2166 :
2167 9 : void ScInterpreter::ScSNormInv()
2168 : {
2169 9 : double x = GetDouble();
2170 9 : if (x < 0.0 || x > 1.0)
2171 0 : PushIllegalArgument();
2172 9 : else if (x == 0.0 || x == 1.0)
2173 0 : PushNoValue();
2174 : else
2175 9 : PushDouble(gaussinv(x));
2176 9 : }
2177 :
2178 9 : void ScInterpreter::ScLogNormInv()
2179 : {
2180 9 : if ( MustHaveParamCount( GetByte(), 3 ) )
2181 : {
2182 9 : double sigma = GetDouble(); // Stdabw
2183 9 : double mue = GetDouble(); // Mittelwert
2184 9 : double y = GetDouble(); // y
2185 9 : if (sigma <= 0.0 || y <= 0.0 || y >= 1.0)
2186 0 : PushIllegalArgument();
2187 : else
2188 9 : PushDouble(exp(mue+sigma*gaussinv(y)));
2189 : }
2190 9 : }
2191 :
2192 : class ScGammaDistFunction : public ScDistFunc
2193 : {
2194 : ScInterpreter& rInt;
2195 : double fp, fAlpha, fBeta;
2196 :
2197 : public:
2198 9 : ScGammaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2199 9 : rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2200 :
2201 9 : virtual ~ScGammaDistFunction() {}
2202 :
2203 237 : double GetValue( double x ) const SAL_OVERRIDE { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); }
2204 : };
2205 :
2206 9 : void ScInterpreter::ScGammaInv()
2207 : {
2208 9 : if ( !MustHaveParamCount( GetByte(), 3 ) )
2209 0 : return;
2210 9 : double fBeta = GetDouble();
2211 9 : double fAlpha = GetDouble();
2212 9 : double fP = GetDouble();
2213 9 : if (fAlpha <= 0.0 || fBeta <= 0.0 || fP < 0.0 || fP >= 1.0 )
2214 : {
2215 0 : PushIllegalArgument();
2216 0 : return;
2217 : }
2218 9 : if (fP == 0.0)
2219 0 : PushInt(0);
2220 : else
2221 : {
2222 : bool bConvError;
2223 9 : ScGammaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2224 9 : double fStart = fAlpha * fBeta;
2225 9 : double fVal = lcl_IterateInverse( aFunc, fStart*0.5, fStart, bConvError );
2226 9 : if (bConvError)
2227 0 : SetError(errNoConvergence);
2228 9 : PushDouble(fVal);
2229 : }
2230 : }
2231 :
2232 : class ScBetaDistFunction : public ScDistFunc
2233 : {
2234 : ScInterpreter& rInt;
2235 : double fp, fAlpha, fBeta;
2236 :
2237 : public:
2238 9 : ScBetaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2239 9 : rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2240 :
2241 9 : virtual ~ScBetaDistFunction() {}
2242 :
2243 195 : double GetValue( double x ) const SAL_OVERRIDE { return fp - rInt.GetBetaDist(x, fAlpha, fBeta); }
2244 : };
2245 :
2246 9 : void ScInterpreter::ScBetaInv()
2247 : {
2248 9 : sal_uInt8 nParamCount = GetByte();
2249 9 : if ( !MustHaveParamCount( nParamCount, 3, 5 ) )
2250 0 : return;
2251 : double fP, fA, fB, fAlpha, fBeta;
2252 9 : if (nParamCount == 5)
2253 3 : fB = GetDouble();
2254 : else
2255 6 : fB = 1.0;
2256 9 : if (nParamCount >= 4)
2257 3 : fA = GetDouble();
2258 : else
2259 6 : fA = 0.0;
2260 9 : fBeta = GetDouble();
2261 9 : fAlpha = GetDouble();
2262 9 : fP = GetDouble();
2263 9 : if (fP < 0.0 || fP >= 1.0 || fA == fB || fAlpha <= 0.0 || fBeta <= 0.0)
2264 : {
2265 0 : PushIllegalArgument();
2266 0 : return;
2267 : }
2268 9 : if (fP == 0.0)
2269 0 : PushInt(0);
2270 : else
2271 : {
2272 : bool bConvError;
2273 9 : ScBetaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2274 : // 0..1 as range for iteration so it isn't extended beyond the valid range
2275 9 : double fVal = lcl_IterateInverse( aFunc, 0.0, 1.0, bConvError );
2276 9 : if (bConvError)
2277 0 : PushError( errNoConvergence);
2278 : else
2279 9 : PushDouble(fA + fVal*(fB-fA)); // scale to (A,B)
2280 : }
2281 : }
2282 :
2283 : // Achtung: T, F und Chi
2284 : // sind monoton fallend,
2285 : // deshalb 1-Dist als Funktion
2286 :
2287 : class ScTDistFunction : public ScDistFunc
2288 : {
2289 : ScInterpreter& rInt;
2290 : double fp, fDF;
2291 : int nT;
2292 :
2293 : public:
2294 21 : ScTDistFunction( ScInterpreter& rI, double fpVal, double fDFVal, int nType ) :
2295 21 : rInt( rI ), fp( fpVal ), fDF( fDFVal ), nT( nType ) {}
2296 :
2297 21 : virtual ~ScTDistFunction() {}
2298 :
2299 432 : double GetValue( double x ) const SAL_OVERRIDE { return fp - rInt.GetTDist( x, fDF, nT ); }
2300 : };
2301 :
2302 15 : void ScInterpreter::ScTInv( int nType )
2303 : {
2304 15 : if ( !MustHaveParamCount( GetByte(), 2 ) )
2305 0 : return;
2306 15 : double fDF = ::rtl::math::approxFloor(GetDouble());
2307 15 : double fP = GetDouble();
2308 15 : if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
2309 : {
2310 0 : PushIllegalArgument();
2311 0 : return;
2312 : }
2313 15 : if ( nType == 4 ) // left-tailed cumulative t-distribution
2314 : {
2315 6 : if ( fP < 0.5 )
2316 6 : PushDouble( -GetTInv( 1 - fP, fDF, nType ) );
2317 : else
2318 0 : PushDouble( GetTInv( fP, fDF, nType ) );
2319 : }
2320 : else
2321 9 : PushDouble( GetTInv( fP, fDF, nType ) );
2322 : };
2323 :
2324 21 : double ScInterpreter::GetTInv( double fAlpha, double fSize, int nType )
2325 : {
2326 : bool bConvError;
2327 21 : ScTDistFunction aFunc( *this, fAlpha, fSize, nType );
2328 21 : double fVal = lcl_IterateInverse( aFunc, fSize * 0.5, fSize, bConvError );
2329 21 : if (bConvError)
2330 0 : SetError(errNoConvergence);
2331 21 : return fVal;
2332 : }
2333 :
2334 : class ScFDistFunction : public ScDistFunc
2335 : {
2336 : ScInterpreter& rInt;
2337 : double fp, fF1, fF2;
2338 :
2339 : public:
2340 12 : ScFDistFunction( ScInterpreter& rI, double fpVal, double fF1Val, double fF2Val ) :
2341 12 : rInt(rI), fp(fpVal), fF1(fF1Val), fF2(fF2Val) {}
2342 :
2343 12 : virtual ~ScFDistFunction() {}
2344 :
2345 240 : double GetValue( double x ) const SAL_OVERRIDE { return fp - rInt.GetFDist(x, fF1, fF2); }
2346 : };
2347 :
2348 9 : void ScInterpreter::ScFInv()
2349 : {
2350 9 : if ( !MustHaveParamCount( GetByte(), 3 ) )
2351 0 : return;
2352 9 : double fF2 = ::rtl::math::approxFloor(GetDouble());
2353 9 : double fF1 = ::rtl::math::approxFloor(GetDouble());
2354 9 : double fP = GetDouble();
2355 9 : if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
2356 : {
2357 0 : PushIllegalArgument();
2358 0 : return;
2359 : }
2360 :
2361 : bool bConvError;
2362 9 : ScFDistFunction aFunc( *this, fP, fF1, fF2 );
2363 9 : double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
2364 9 : if (bConvError)
2365 0 : SetError(errNoConvergence);
2366 9 : PushDouble(fVal);
2367 : }
2368 :
2369 3 : void ScInterpreter::ScFInv_LT()
2370 : {
2371 3 : if ( !MustHaveParamCount( GetByte(), 3 ) )
2372 0 : return;
2373 3 : double fF2 = ::rtl::math::approxFloor(GetDouble());
2374 3 : double fF1 = ::rtl::math::approxFloor(GetDouble());
2375 3 : double fP = GetDouble();
2376 3 : if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
2377 : {
2378 0 : PushIllegalArgument();
2379 0 : return;
2380 : }
2381 :
2382 : bool bConvError;
2383 3 : ScFDistFunction aFunc( *this, ( 1.0 - fP ), fF1, fF2 );
2384 3 : double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
2385 3 : if (bConvError)
2386 0 : SetError(errNoConvergence);
2387 3 : PushDouble(fVal);
2388 : }
2389 :
2390 : class ScChiDistFunction : public ScDistFunc
2391 : {
2392 : ScInterpreter& rInt;
2393 : double fp, fDF;
2394 :
2395 : public:
2396 12 : ScChiDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2397 12 : rInt(rI), fp(fpVal), fDF(fDFVal) {}
2398 :
2399 12 : virtual ~ScChiDistFunction() {}
2400 :
2401 183 : double GetValue( double x ) const SAL_OVERRIDE { return fp - rInt.GetChiDist(x, fDF); }
2402 : };
2403 :
2404 12 : void ScInterpreter::ScChiInv()
2405 : {
2406 12 : if ( !MustHaveParamCount( GetByte(), 2 ) )
2407 0 : return;
2408 12 : double fDF = ::rtl::math::approxFloor(GetDouble());
2409 12 : double fP = GetDouble();
2410 12 : if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
2411 : {
2412 0 : PushIllegalArgument();
2413 0 : return;
2414 : }
2415 :
2416 : bool bConvError;
2417 12 : ScChiDistFunction aFunc( *this, fP, fDF );
2418 12 : double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2419 12 : if (bConvError)
2420 0 : SetError(errNoConvergence);
2421 12 : PushDouble(fVal);
2422 : }
2423 :
2424 : /***********************************************/
2425 : class ScChiSqDistFunction : public ScDistFunc
2426 : {
2427 : ScInterpreter& rInt;
2428 : double fp, fDF;
2429 :
2430 : public:
2431 7 : ScChiSqDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2432 7 : rInt(rI), fp(fpVal), fDF(fDFVal) {}
2433 :
2434 7 : virtual ~ScChiSqDistFunction() {}
2435 :
2436 81 : double GetValue( double x ) const SAL_OVERRIDE { return fp - rInt.GetChiSqDistCDF(x, fDF); }
2437 : };
2438 :
2439 7 : void ScInterpreter::ScChiSqInv()
2440 : {
2441 7 : if ( !MustHaveParamCount( GetByte(), 2 ) )
2442 0 : return;
2443 7 : double fDF = ::rtl::math::approxFloor(GetDouble());
2444 7 : double fP = GetDouble();
2445 7 : if (fDF < 1.0 || fP < 0.0 || fP >= 1.0 )
2446 : {
2447 0 : PushIllegalArgument();
2448 0 : return;
2449 : }
2450 :
2451 : bool bConvError;
2452 7 : ScChiSqDistFunction aFunc( *this, fP, fDF );
2453 7 : double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2454 7 : if (bConvError)
2455 0 : SetError(errNoConvergence);
2456 7 : PushDouble(fVal);
2457 : }
2458 :
2459 9 : void ScInterpreter::ScConfidence()
2460 : {
2461 9 : if ( MustHaveParamCount( GetByte(), 3 ) )
2462 : {
2463 9 : double n = ::rtl::math::approxFloor(GetDouble());
2464 9 : double sigma = GetDouble();
2465 9 : double alpha = GetDouble();
2466 9 : if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
2467 0 : PushIllegalArgument();
2468 : else
2469 9 : PushDouble( gaussinv(1.0-alpha/2.0) * sigma/sqrt(n) );
2470 : }
2471 9 : }
2472 :
2473 6 : void ScInterpreter::ScConfidenceT()
2474 : {
2475 6 : if ( MustHaveParamCount( GetByte(), 3 ) )
2476 : {
2477 6 : double n = ::rtl::math::approxFloor(GetDouble());
2478 6 : double sigma = GetDouble();
2479 6 : double alpha = GetDouble();
2480 6 : if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
2481 0 : PushIllegalArgument();
2482 : else
2483 6 : PushDouble( sigma * GetTInv( alpha, n - 1, 2 ) / sqrt( n ) );
2484 : }
2485 6 : }
2486 :
2487 3 : void ScInterpreter::ScZTest()
2488 : {
2489 3 : sal_uInt8 nParamCount = GetByte();
2490 3 : if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
2491 0 : return;
2492 3 : double sigma = 0.0, x;
2493 3 : if (nParamCount == 3)
2494 : {
2495 3 : sigma = GetDouble();
2496 3 : if (sigma <= 0.0)
2497 : {
2498 0 : PushIllegalArgument();
2499 0 : return;
2500 : }
2501 : }
2502 3 : x = GetDouble();
2503 :
2504 3 : double fSum = 0.0;
2505 3 : double fSumSqr = 0.0;
2506 : double fVal;
2507 3 : double rValCount = 0.0;
2508 3 : switch (GetStackType())
2509 : {
2510 : case formula::svDouble :
2511 : {
2512 0 : fVal = GetDouble();
2513 0 : fSum += fVal;
2514 0 : fSumSqr += fVal*fVal;
2515 0 : rValCount++;
2516 : }
2517 0 : break;
2518 : case svSingleRef :
2519 : {
2520 0 : ScAddress aAdr;
2521 0 : PopSingleRef( aAdr );
2522 0 : ScRefCellValue aCell;
2523 0 : aCell.assign(*pDok, aAdr);
2524 0 : if (aCell.hasNumeric())
2525 : {
2526 0 : fVal = GetCellValue(aAdr, aCell);
2527 0 : fSum += fVal;
2528 0 : fSumSqr += fVal*fVal;
2529 0 : rValCount++;
2530 0 : }
2531 : }
2532 0 : break;
2533 : case svRefList :
2534 : case formula::svDoubleRef :
2535 : {
2536 3 : short nParam = 1;
2537 3 : size_t nRefInList = 0;
2538 9 : while (nParam-- > 0)
2539 : {
2540 3 : ScRange aRange;
2541 3 : sal_uInt16 nErr = 0;
2542 3 : PopDoubleRef( aRange, nParam, nRefInList);
2543 3 : ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
2544 3 : if (aValIter.GetFirst(fVal, nErr))
2545 : {
2546 3 : fSum += fVal;
2547 3 : fSumSqr += fVal*fVal;
2548 3 : rValCount++;
2549 30 : while ((nErr == 0) && aValIter.GetNext(fVal, nErr))
2550 : {
2551 24 : fSum += fVal;
2552 24 : fSumSqr += fVal*fVal;
2553 24 : rValCount++;
2554 : }
2555 3 : SetError(nErr);
2556 : }
2557 : }
2558 : }
2559 3 : break;
2560 : case svMatrix :
2561 : case svExternalSingleRef:
2562 : case svExternalDoubleRef:
2563 : {
2564 0 : ScMatrixRef pMat = GetMatrix();
2565 0 : if (pMat)
2566 : {
2567 0 : SCSIZE nCount = pMat->GetElementCount();
2568 0 : if (pMat->IsNumeric())
2569 : {
2570 0 : for ( SCSIZE i = 0; i < nCount; i++ )
2571 : {
2572 0 : fVal= pMat->GetDouble(i);
2573 0 : fSum += fVal;
2574 0 : fSumSqr += fVal * fVal;
2575 0 : rValCount++;
2576 : }
2577 : }
2578 : else
2579 : {
2580 0 : for (SCSIZE i = 0; i < nCount; i++)
2581 0 : if (!pMat->IsString(i))
2582 : {
2583 0 : fVal= pMat->GetDouble(i);
2584 0 : fSum += fVal;
2585 0 : fSumSqr += fVal * fVal;
2586 0 : rValCount++;
2587 : }
2588 : }
2589 0 : }
2590 : }
2591 0 : break;
2592 0 : default : SetError(errIllegalParameter); break;
2593 : }
2594 3 : if (rValCount <= 1.0)
2595 0 : PushError( errDivisionByZero);
2596 : else
2597 : {
2598 3 : double mue = fSum/rValCount;
2599 :
2600 3 : if (nParamCount != 3)
2601 : {
2602 0 : sigma = (fSumSqr - fSum*fSum/rValCount)/(rValCount-1.0);
2603 0 : PushDouble(0.5 - gauss((mue-x)/sqrt(sigma/rValCount)));
2604 : }
2605 : else
2606 3 : PushDouble(0.5 - gauss((mue-x)*sqrt(rValCount)/sigma));
2607 : }
2608 : }
2609 :
2610 6 : bool ScInterpreter::CalculateTest(bool _bTemplin
2611 : ,const SCSIZE nC1, const SCSIZE nC2,const SCSIZE nR1,const SCSIZE nR2
2612 : ,const ScMatrixRef& pMat1,const ScMatrixRef& pMat2
2613 : ,double& fT,double& fF)
2614 : {
2615 6 : double fCount1 = 0.0;
2616 6 : double fCount2 = 0.0;
2617 6 : double fSum1 = 0.0;
2618 6 : double fSumSqr1 = 0.0;
2619 6 : double fSum2 = 0.0;
2620 6 : double fSumSqr2 = 0.0;
2621 : double fVal;
2622 : SCSIZE i,j;
2623 12 : for (i = 0; i < nC1; i++)
2624 48 : for (j = 0; j < nR1; j++)
2625 : {
2626 42 : if (!pMat1->IsString(i,j))
2627 : {
2628 42 : fVal = pMat1->GetDouble(i,j);
2629 42 : fSum1 += fVal;
2630 42 : fSumSqr1 += fVal * fVal;
2631 42 : fCount1++;
2632 : }
2633 : }
2634 12 : for (i = 0; i < nC2; i++)
2635 48 : for (j = 0; j < nR2; j++)
2636 : {
2637 42 : if (!pMat2->IsString(i,j))
2638 : {
2639 42 : fVal = pMat2->GetDouble(i,j);
2640 42 : fSum2 += fVal;
2641 42 : fSumSqr2 += fVal * fVal;
2642 42 : fCount2++;
2643 : }
2644 : }
2645 6 : if (fCount1 < 2.0 || fCount2 < 2.0)
2646 : {
2647 0 : PushNoValue();
2648 0 : return false;
2649 : } // if (fCount1 < 2.0 || fCount2 < 2.0)
2650 6 : if ( _bTemplin )
2651 : {
2652 0 : double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0)/fCount1;
2653 0 : double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0)/fCount2;
2654 0 : if (fS1 + fS2 == 0.0)
2655 : {
2656 0 : PushNoValue();
2657 0 : return false;
2658 : }
2659 0 : fT = fabs(fSum1/fCount1 - fSum2/fCount2)/sqrt(fS1+fS2);
2660 0 : double c = fS1/(fS1+fS2);
2661 : // GetTDist wird mit GetBetaDist berechnet und kommt auch mit nicht ganzzahligen
2662 : // Freiheitsgraden klar. Dann stimmt das Ergebnis auch mit Excel ueberein (#52406#):
2663 0 : fF = 1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0));
2664 : }
2665 : else
2666 : {
2667 : // laut Bronstein-Semendjajew
2668 6 : double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1) / (fCount1 - 1.0); // Varianz
2669 6 : double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2) / (fCount2 - 1.0);
2670 12 : fT = fabs( fSum1/fCount1 - fSum2/fCount2 ) /
2671 12 : sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
2672 12 : sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
2673 6 : fF = fCount1 + fCount2 - 2;
2674 : }
2675 6 : return true;
2676 : }
2677 9 : void ScInterpreter::ScTTest()
2678 : {
2679 9 : if ( !MustHaveParamCount( GetByte(), 4 ) )
2680 0 : return;
2681 9 : double fTyp = ::rtl::math::approxFloor(GetDouble());
2682 9 : double fTails = ::rtl::math::approxFloor(GetDouble());
2683 9 : if (fTails != 1.0 && fTails != 2.0)
2684 : {
2685 0 : PushIllegalArgument();
2686 0 : return;
2687 : }
2688 :
2689 9 : ScMatrixRef pMat2 = GetMatrix();
2690 18 : ScMatrixRef pMat1 = GetMatrix();
2691 9 : if (!pMat1 || !pMat2)
2692 : {
2693 0 : PushIllegalParameter();
2694 0 : return;
2695 : }
2696 : double fT, fF;
2697 : SCSIZE nC1, nC2;
2698 : SCSIZE nR1, nR2;
2699 : SCSIZE i, j;
2700 9 : pMat1->GetDimensions(nC1, nR1);
2701 9 : pMat2->GetDimensions(nC2, nR2);
2702 9 : if (fTyp == 1.0)
2703 : {
2704 3 : if (nC1 != nC2 || nR1 != nR2)
2705 : {
2706 0 : PushIllegalArgument();
2707 0 : return;
2708 : }
2709 3 : double fCount = 0.0;
2710 3 : double fSum1 = 0.0;
2711 3 : double fSum2 = 0.0;
2712 3 : double fSumSqrD = 0.0;
2713 : double fVal1, fVal2;
2714 6 : for (i = 0; i < nC1; i++)
2715 30 : for (j = 0; j < nR1; j++)
2716 : {
2717 27 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
2718 : {
2719 27 : fVal1 = pMat1->GetDouble(i,j);
2720 27 : fVal2 = pMat2->GetDouble(i,j);
2721 27 : fSum1 += fVal1;
2722 27 : fSum2 += fVal2;
2723 27 : fSumSqrD += (fVal1 - fVal2)*(fVal1 - fVal2);
2724 27 : fCount++;
2725 : }
2726 : }
2727 3 : if (fCount < 1.0)
2728 : {
2729 0 : PushNoValue();
2730 0 : return;
2731 : }
2732 6 : fT = sqrt(fCount-1.0) * fabs(fSum1 - fSum2) /
2733 6 : sqrt(fCount * fSumSqrD - (fSum1-fSum2)*(fSum1-fSum2));
2734 3 : fF = fCount - 1.0;
2735 : }
2736 6 : else if (fTyp == 2.0)
2737 : {
2738 6 : CalculateTest(false,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF);
2739 : }
2740 0 : else if (fTyp == 3.0)
2741 : {
2742 0 : CalculateTest(true,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF);
2743 : }
2744 :
2745 : else
2746 : {
2747 0 : PushIllegalArgument();
2748 0 : return;
2749 : }
2750 18 : PushDouble( GetTDist( fT, fF, ( int )fTails ) );
2751 : }
2752 :
2753 9 : void ScInterpreter::ScFTest()
2754 : {
2755 9 : if ( !MustHaveParamCount( GetByte(), 2 ) )
2756 0 : return;
2757 9 : ScMatrixRef pMat2 = GetMatrix();
2758 18 : ScMatrixRef pMat1 = GetMatrix();
2759 9 : if (!pMat1 || !pMat2)
2760 : {
2761 0 : PushIllegalParameter();
2762 0 : return;
2763 : }
2764 : SCSIZE nC1, nC2;
2765 : SCSIZE nR1, nR2;
2766 : SCSIZE i, j;
2767 9 : pMat1->GetDimensions(nC1, nR1);
2768 9 : pMat2->GetDimensions(nC2, nR2);
2769 9 : double fCount1 = 0.0;
2770 9 : double fCount2 = 0.0;
2771 9 : double fSum1 = 0.0;
2772 9 : double fSumSqr1 = 0.0;
2773 9 : double fSum2 = 0.0;
2774 9 : double fSumSqr2 = 0.0;
2775 : double fVal;
2776 18 : for (i = 0; i < nC1; i++)
2777 63 : for (j = 0; j < nR1; j++)
2778 : {
2779 54 : if (!pMat1->IsString(i,j))
2780 : {
2781 54 : fVal = pMat1->GetDouble(i,j);
2782 54 : fSum1 += fVal;
2783 54 : fSumSqr1 += fVal * fVal;
2784 54 : fCount1++;
2785 : }
2786 : }
2787 18 : for (i = 0; i < nC2; i++)
2788 66 : for (j = 0; j < nR2; j++)
2789 : {
2790 57 : if (!pMat2->IsString(i,j))
2791 : {
2792 57 : fVal = pMat2->GetDouble(i,j);
2793 57 : fSum2 += fVal;
2794 57 : fSumSqr2 += fVal * fVal;
2795 57 : fCount2++;
2796 : }
2797 : }
2798 9 : if (fCount1 < 2.0 || fCount2 < 2.0)
2799 : {
2800 0 : PushNoValue();
2801 0 : return;
2802 : }
2803 9 : double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0);
2804 9 : double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0);
2805 9 : if (fS1 == 0.0 || fS2 == 0.0)
2806 : {
2807 0 : PushNoValue();
2808 0 : return;
2809 : }
2810 : double fF, fF1, fF2;
2811 9 : if (fS1 > fS2)
2812 : {
2813 3 : fF = fS1/fS2;
2814 3 : fF1 = fCount1-1.0;
2815 3 : fF2 = fCount2-1.0;
2816 : }
2817 : else
2818 : {
2819 6 : fF = fS2/fS1;
2820 6 : fF1 = fCount2-1.0;
2821 6 : fF2 = fCount1-1.0;
2822 : }
2823 18 : PushDouble(2.0*GetFDist(fF, fF1, fF2));
2824 : }
2825 :
2826 12 : void ScInterpreter::ScChiTest()
2827 : {
2828 12 : if ( !MustHaveParamCount( GetByte(), 2 ) )
2829 0 : return;
2830 12 : ScMatrixRef pMat2 = GetMatrix();
2831 24 : ScMatrixRef pMat1 = GetMatrix();
2832 12 : if (!pMat1 || !pMat2)
2833 : {
2834 0 : PushIllegalParameter();
2835 0 : return;
2836 : }
2837 : SCSIZE nC1, nC2;
2838 : SCSIZE nR1, nR2;
2839 12 : pMat1->GetDimensions(nC1, nR1);
2840 12 : pMat2->GetDimensions(nC2, nR2);
2841 12 : if (nR1 != nR2 || nC1 != nC2)
2842 : {
2843 0 : PushIllegalArgument();
2844 0 : return;
2845 : }
2846 12 : double fChi = 0.0;
2847 24 : for (SCSIZE i = 0; i < nC1; i++)
2848 : {
2849 93 : for (SCSIZE j = 0; j < nR1; j++)
2850 : {
2851 81 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
2852 : {
2853 81 : double fValX = pMat1->GetDouble(i,j);
2854 81 : double fValE = pMat2->GetDouble(i,j);
2855 81 : fChi += (fValX - fValE) * (fValX - fValE) / fValE;
2856 : }
2857 : else
2858 : {
2859 0 : PushIllegalArgument();
2860 0 : return;
2861 : }
2862 : }
2863 : }
2864 : double fDF;
2865 12 : if (nC1 == 1 || nR1 == 1)
2866 : {
2867 12 : fDF = (double)(nC1*nR1 - 1);
2868 24 : if (fDF == 0.0)
2869 : {
2870 0 : PushNoValue();
2871 0 : return;
2872 : }
2873 : }
2874 : else
2875 0 : fDF = (double)(nC1-1)*(double)(nR1-1);
2876 24 : PushDouble(GetChiDist(fChi, fDF));
2877 : }
2878 :
2879 3 : void ScInterpreter::ScKurt()
2880 : {
2881 : double fSum,fCount,vSum;
2882 3 : std::vector<double> values;
2883 3 : if ( !CalculateSkew(fSum,fCount,vSum,values) )
2884 0 : return;
2885 :
2886 3 : if (fCount == 0.0)
2887 : {
2888 0 : PushError( errDivisionByZero);
2889 0 : return;
2890 : }
2891 :
2892 3 : double fMean = fSum / fCount;
2893 :
2894 21 : for (size_t i = 0; i < values.size(); i++)
2895 18 : vSum += (values[i] - fMean) * (values[i] - fMean);
2896 :
2897 3 : double fStdDev = sqrt(vSum / (fCount - 1.0));
2898 3 : double xpower4 = 0.0;
2899 :
2900 3 : if (fStdDev == 0.0)
2901 : {
2902 0 : PushError( errDivisionByZero);
2903 0 : return;
2904 : }
2905 :
2906 21 : for (size_t i = 0; i < values.size(); i++)
2907 : {
2908 18 : double dx = (values[i] - fMean) / fStdDev;
2909 18 : xpower4 = xpower4 + (dx * dx * dx * dx);
2910 : }
2911 :
2912 3 : double k_d = (fCount - 2.0) * (fCount - 3.0);
2913 3 : double k_l = fCount * (fCount + 1.0) / ((fCount - 1.0) * k_d);
2914 3 : double k_t = 3.0 * (fCount - 1.0) * (fCount - 1.0) / k_d;
2915 :
2916 3 : PushDouble(xpower4 * k_l - k_t);
2917 : }
2918 :
2919 3 : void ScInterpreter::ScHarMean()
2920 : {
2921 3 : short nParamCount = GetByte();
2922 3 : double nVal = 0.0;
2923 3 : double nValCount = 0.0;
2924 3 : ScAddress aAdr;
2925 3 : ScRange aRange;
2926 3 : size_t nRefInList = 0;
2927 15 : while ((nGlobalError == 0) && (nParamCount-- > 0))
2928 : {
2929 9 : switch (GetStackType())
2930 : {
2931 : case formula::svDouble :
2932 : {
2933 9 : double x = GetDouble();
2934 9 : if (x > 0.0)
2935 : {
2936 9 : nVal += 1.0/x;
2937 9 : nValCount++;
2938 : }
2939 : else
2940 0 : SetError( errIllegalArgument);
2941 9 : break;
2942 : }
2943 : case svSingleRef :
2944 : {
2945 0 : PopSingleRef( aAdr );
2946 0 : ScRefCellValue aCell;
2947 0 : aCell.assign(*pDok, aAdr);
2948 0 : if (aCell.hasNumeric())
2949 : {
2950 0 : double x = GetCellValue(aAdr, aCell);
2951 0 : if (x > 0.0)
2952 : {
2953 0 : nVal += 1.0/x;
2954 0 : nValCount++;
2955 : }
2956 : else
2957 0 : SetError( errIllegalArgument);
2958 : }
2959 0 : break;
2960 : }
2961 : case formula::svDoubleRef :
2962 : case svRefList :
2963 : {
2964 0 : sal_uInt16 nErr = 0;
2965 0 : PopDoubleRef( aRange, nParamCount, nRefInList);
2966 : double nCellVal;
2967 0 : ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
2968 0 : if (aValIter.GetFirst(nCellVal, nErr))
2969 : {
2970 0 : if (nCellVal > 0.0)
2971 : {
2972 0 : nVal += 1.0/nCellVal;
2973 0 : nValCount++;
2974 : }
2975 : else
2976 0 : SetError( errIllegalArgument);
2977 0 : SetError(nErr);
2978 0 : while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
2979 : {
2980 0 : if (nCellVal > 0.0)
2981 : {
2982 0 : nVal += 1.0/nCellVal;
2983 0 : nValCount++;
2984 : }
2985 : else
2986 0 : SetError( errIllegalArgument);
2987 : }
2988 0 : SetError(nErr);
2989 : }
2990 : }
2991 0 : break;
2992 : case svMatrix :
2993 : case svExternalSingleRef:
2994 : case svExternalDoubleRef:
2995 : {
2996 0 : ScMatrixRef pMat = GetMatrix();
2997 0 : if (pMat)
2998 : {
2999 0 : SCSIZE nCount = pMat->GetElementCount();
3000 0 : if (pMat->IsNumeric())
3001 : {
3002 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3003 : {
3004 0 : double x = pMat->GetDouble(nElem);
3005 0 : if (x > 0.0)
3006 : {
3007 0 : nVal += 1.0/x;
3008 0 : nValCount++;
3009 : }
3010 : else
3011 0 : SetError( errIllegalArgument);
3012 : }
3013 : }
3014 : else
3015 : {
3016 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3017 0 : if (!pMat->IsString(nElem))
3018 : {
3019 0 : double x = pMat->GetDouble(nElem);
3020 0 : if (x > 0.0)
3021 : {
3022 0 : nVal += 1.0/x;
3023 0 : nValCount++;
3024 : }
3025 : else
3026 0 : SetError( errIllegalArgument);
3027 : }
3028 : }
3029 0 : }
3030 : }
3031 0 : break;
3032 0 : default : SetError(errIllegalParameter); break;
3033 : }
3034 : }
3035 3 : if (nGlobalError == 0)
3036 3 : PushDouble((double)nValCount/nVal);
3037 : else
3038 0 : PushError( nGlobalError);
3039 3 : }
3040 :
3041 3 : void ScInterpreter::ScGeoMean()
3042 : {
3043 3 : short nParamCount = GetByte();
3044 3 : double nVal = 0.0;
3045 3 : double nValCount = 0.0;
3046 3 : ScAddress aAdr;
3047 3 : ScRange aRange;
3048 :
3049 3 : size_t nRefInList = 0;
3050 15 : while ((nGlobalError == 0) && (nParamCount-- > 0))
3051 : {
3052 9 : switch (GetStackType())
3053 : {
3054 : case formula::svDouble :
3055 : {
3056 9 : double x = GetDouble();
3057 9 : if (x > 0.0)
3058 : {
3059 9 : nVal += log(x);
3060 9 : nValCount++;
3061 : }
3062 : else
3063 0 : SetError( errIllegalArgument);
3064 9 : break;
3065 : }
3066 : case svSingleRef :
3067 : {
3068 0 : PopSingleRef( aAdr );
3069 0 : ScRefCellValue aCell;
3070 0 : aCell.assign(*pDok, aAdr);
3071 0 : if (aCell.hasNumeric())
3072 : {
3073 0 : double x = GetCellValue(aAdr, aCell);
3074 0 : if (x > 0.0)
3075 : {
3076 0 : nVal += log(x);
3077 0 : nValCount++;
3078 : }
3079 : else
3080 0 : SetError( errIllegalArgument);
3081 : }
3082 0 : break;
3083 : }
3084 : case formula::svDoubleRef :
3085 : case svRefList :
3086 : {
3087 0 : sal_uInt16 nErr = 0;
3088 0 : PopDoubleRef( aRange, nParamCount, nRefInList);
3089 : double nCellVal;
3090 0 : ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
3091 0 : if (aValIter.GetFirst(nCellVal, nErr))
3092 : {
3093 0 : if (nCellVal > 0.0)
3094 : {
3095 0 : nVal += log(nCellVal);
3096 0 : nValCount++;
3097 : }
3098 : else
3099 0 : SetError( errIllegalArgument);
3100 0 : SetError(nErr);
3101 0 : while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
3102 : {
3103 0 : if (nCellVal > 0.0)
3104 : {
3105 0 : nVal += log(nCellVal);
3106 0 : nValCount++;
3107 : }
3108 : else
3109 0 : SetError( errIllegalArgument);
3110 : }
3111 0 : SetError(nErr);
3112 : }
3113 : }
3114 0 : break;
3115 : case svMatrix :
3116 : case svExternalSingleRef:
3117 : case svExternalDoubleRef:
3118 : {
3119 0 : ScMatrixRef pMat = GetMatrix();
3120 0 : if (pMat)
3121 : {
3122 0 : SCSIZE nCount = pMat->GetElementCount();
3123 0 : if (pMat->IsNumeric())
3124 : {
3125 0 : for (SCSIZE ui = 0; ui < nCount; ui++)
3126 : {
3127 0 : double x = pMat->GetDouble(ui);
3128 0 : if (x > 0.0)
3129 : {
3130 0 : nVal += log(x);
3131 0 : nValCount++;
3132 : }
3133 : else
3134 0 : SetError( errIllegalArgument);
3135 : }
3136 : }
3137 : else
3138 : {
3139 0 : for (SCSIZE ui = 0; ui < nCount; ui++)
3140 0 : if (!pMat->IsString(ui))
3141 : {
3142 0 : double x = pMat->GetDouble(ui);
3143 0 : if (x > 0.0)
3144 : {
3145 0 : nVal += log(x);
3146 0 : nValCount++;
3147 : }
3148 : else
3149 0 : SetError( errIllegalArgument);
3150 : }
3151 : }
3152 0 : }
3153 : }
3154 0 : break;
3155 0 : default : SetError(errIllegalParameter); break;
3156 : }
3157 : }
3158 3 : if (nGlobalError == 0)
3159 3 : PushDouble(exp(nVal / nValCount));
3160 : else
3161 0 : PushError( nGlobalError);
3162 3 : }
3163 :
3164 3 : void ScInterpreter::ScStandard()
3165 : {
3166 3 : if ( MustHaveParamCount( GetByte(), 3 ) )
3167 : {
3168 3 : double sigma = GetDouble();
3169 3 : double mue = GetDouble();
3170 3 : double x = GetDouble();
3171 3 : if (sigma < 0.0)
3172 0 : PushError( errIllegalArgument);
3173 3 : else if (sigma == 0.0)
3174 0 : PushError( errDivisionByZero);
3175 : else
3176 3 : PushDouble((x-mue)/sigma);
3177 : }
3178 3 : }
3179 4 : bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values)
3180 : {
3181 4 : short nParamCount = GetByte();
3182 4 : if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3183 0 : return false;
3184 :
3185 4 : fSum = 0.0;
3186 4 : fCount = 0.0;
3187 4 : vSum = 0.0;
3188 4 : double fVal = 0.0;
3189 4 : ScAddress aAdr;
3190 4 : ScRange aRange;
3191 4 : size_t nRefInList = 0;
3192 27 : while (nParamCount-- > 0)
3193 : {
3194 19 : switch (GetStackType())
3195 : {
3196 : case formula::svDouble :
3197 : {
3198 0 : fVal = GetDouble();
3199 0 : fSum += fVal;
3200 0 : values.push_back(fVal);
3201 0 : fCount++;
3202 : }
3203 0 : break;
3204 : case svSingleRef :
3205 : {
3206 18 : PopSingleRef( aAdr );
3207 18 : ScRefCellValue aCell;
3208 18 : aCell.assign(*pDok, aAdr);
3209 18 : if (aCell.hasNumeric())
3210 : {
3211 18 : fVal = GetCellValue(aAdr, aCell);
3212 18 : fSum += fVal;
3213 18 : values.push_back(fVal);
3214 18 : fCount++;
3215 18 : }
3216 : }
3217 18 : break;
3218 : case formula::svDoubleRef :
3219 : case svRefList :
3220 : {
3221 1 : PopDoubleRef( aRange, nParamCount, nRefInList);
3222 1 : sal_uInt16 nErr = 0;
3223 1 : ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
3224 1 : if (aValIter.GetFirst(fVal, nErr))
3225 : {
3226 1 : fSum += fVal;
3227 1 : values.push_back(fVal);
3228 1 : fCount++;
3229 1 : SetError(nErr);
3230 11 : while ((nErr == 0) && aValIter.GetNext(fVal, nErr))
3231 : {
3232 9 : fSum += fVal;
3233 9 : values.push_back(fVal);
3234 9 : fCount++;
3235 : }
3236 1 : SetError(nErr);
3237 : }
3238 : }
3239 1 : break;
3240 : case svMatrix :
3241 : case svExternalSingleRef:
3242 : case svExternalDoubleRef:
3243 : {
3244 0 : ScMatrixRef pMat = GetMatrix();
3245 0 : if (pMat)
3246 : {
3247 0 : SCSIZE nCount = pMat->GetElementCount();
3248 0 : if (pMat->IsNumeric())
3249 : {
3250 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3251 : {
3252 0 : fVal = pMat->GetDouble(nElem);
3253 0 : fSum += fVal;
3254 0 : values.push_back(fVal);
3255 0 : fCount++;
3256 : }
3257 : }
3258 : else
3259 : {
3260 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3261 0 : if (!pMat->IsString(nElem))
3262 : {
3263 0 : fVal = pMat->GetDouble(nElem);
3264 0 : fSum += fVal;
3265 0 : values.push_back(fVal);
3266 0 : fCount++;
3267 : }
3268 : }
3269 0 : }
3270 : }
3271 0 : break;
3272 : default :
3273 0 : SetError(errIllegalParameter);
3274 0 : break;
3275 : }
3276 : }
3277 :
3278 4 : if (nGlobalError)
3279 : {
3280 0 : PushError( nGlobalError);
3281 0 : return false;
3282 : } // if (nGlobalError)
3283 4 : return true;
3284 : }
3285 :
3286 1 : void ScInterpreter::CalculateSkewOrSkewp( bool bSkewp )
3287 : {
3288 : double fSum, fCount, vSum;
3289 1 : std::vector<double> values;
3290 1 : if (!CalculateSkew( fSum, fCount, vSum, values))
3291 0 : return;
3292 :
3293 1 : double fMean = fSum / fCount;
3294 :
3295 11 : for (size_t i = 0; i < values.size(); ++i)
3296 10 : vSum += (values[i] - fMean) * (values[i] - fMean);
3297 :
3298 1 : double fStdDev = sqrt( vSum / (bSkewp ? fCount : (fCount - 1.0)));
3299 1 : double xcube = 0.0;
3300 :
3301 1 : if (fStdDev == 0)
3302 : {
3303 0 : PushIllegalArgument();
3304 0 : return;
3305 : }
3306 :
3307 11 : for (size_t i = 0; i < values.size(); ++i)
3308 : {
3309 10 : double dx = (values[i] - fMean) / fStdDev;
3310 10 : xcube = xcube + (dx * dx * dx);
3311 : }
3312 :
3313 1 : if (bSkewp)
3314 0 : PushDouble( xcube / fCount );
3315 : else
3316 1 : PushDouble( ((xcube * fCount) / (fCount - 1.0)) / (fCount - 2.0) );
3317 : }
3318 :
3319 1 : void ScInterpreter::ScSkew()
3320 : {
3321 1 : CalculateSkewOrSkewp( false );
3322 1 : }
3323 :
3324 0 : void ScInterpreter::ScSkewp()
3325 : {
3326 0 : CalculateSkewOrSkewp( true );
3327 0 : }
3328 :
3329 103 : double ScInterpreter::GetMedian( vector<double> & rArray )
3330 : {
3331 103 : size_t nSize = rArray.size();
3332 103 : if (rArray.empty() || nSize == 0 || nGlobalError)
3333 : {
3334 0 : SetError( errNoValue);
3335 0 : return 0.0;
3336 : }
3337 :
3338 : // Upper median.
3339 103 : size_t nMid = nSize / 2;
3340 103 : vector<double>::iterator iMid = rArray.begin() + nMid;
3341 103 : ::std::nth_element( rArray.begin(), iMid, rArray.end());
3342 103 : if (nSize & 1)
3343 40 : return *iMid; // Lower and upper median are equal.
3344 : else
3345 : {
3346 63 : double fUp = *iMid;
3347 : // Lower median.
3348 63 : iMid = rArray.begin() + nMid - 1;
3349 63 : ::std::nth_element( rArray.begin(), iMid, rArray.end());
3350 63 : return (fUp + *iMid) / 2;
3351 : }
3352 : }
3353 :
3354 60 : void ScInterpreter::ScMedian()
3355 : {
3356 60 : sal_uInt8 nParamCount = GetByte();
3357 60 : if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3358 60 : return;
3359 60 : vector<double> aArray;
3360 60 : GetNumberSequenceArray( nParamCount, aArray, false );
3361 60 : PushDouble( GetMedian( aArray));
3362 : }
3363 :
3364 46 : double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile )
3365 : {
3366 46 : size_t nSize = rArray.size();
3367 46 : if (rArray.empty() || nSize == 0 || nGlobalError)
3368 : {
3369 0 : SetError( errNoValue);
3370 0 : return 0.0;
3371 : }
3372 :
3373 46 : if (nSize == 1)
3374 0 : return rArray[0];
3375 : else
3376 : {
3377 46 : size_t nIndex = (size_t)::rtl::math::approxFloor( fPercentile * (nSize-1));
3378 46 : double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1));
3379 : OSL_ENSURE(nIndex < nSize, "GetPercentile: wrong index(1)");
3380 46 : vector<double>::iterator iter = rArray.begin() + nIndex;
3381 46 : ::std::nth_element( rArray.begin(), iter, rArray.end());
3382 46 : if (fDiff == 0.0)
3383 21 : return *iter;
3384 : else
3385 : {
3386 : OSL_ENSURE(nIndex < nSize-1, "GetPercentile: wrong index(2)");
3387 25 : double fVal = *iter;
3388 25 : iter = rArray.begin() + nIndex+1;
3389 25 : ::std::nth_element( rArray.begin(), iter, rArray.end());
3390 25 : return fVal + fDiff * (*iter - fVal);
3391 : }
3392 : }
3393 : }
3394 :
3395 48 : double ScInterpreter::GetPercentileExclusive( vector<double> & rArray, double fPercentile )
3396 : {
3397 48 : size_t nSize1 = rArray.size() + 1;
3398 48 : if ( rArray.empty() || nSize1 == 1 || nGlobalError )
3399 : {
3400 0 : SetError( errNoValue );
3401 0 : return 0.0;
3402 : }
3403 48 : if ( fPercentile * nSize1 < 1.0 || fPercentile * nSize1 > (double) ( nSize1 - 1 ) )
3404 : {
3405 0 : SetError( errIllegalParameter );
3406 0 : return 0.0;
3407 : }
3408 :
3409 48 : size_t nIndex = (size_t)::rtl::math::approxFloor( fPercentile * nSize1 - 1 );
3410 48 : double fDiff = fPercentile * nSize1 - 1 - ::rtl::math::approxFloor( fPercentile * nSize1 - 1 );
3411 : OSL_ENSURE(nIndex < ( nSize1 - 1 ), "GetPercentile: wrong index(1)");
3412 48 : vector<double>::iterator iter = rArray.begin() + nIndex;
3413 48 : ::std::nth_element( rArray.begin(), iter, rArray.end());
3414 48 : if (fDiff == 0.0)
3415 12 : return *iter;
3416 : else
3417 : {
3418 : OSL_ENSURE(nIndex < nSize1, "GetPercentile: wrong index(2)");
3419 36 : double fVal = *iter;
3420 36 : iter = rArray.begin() + nIndex + 1;
3421 36 : ::std::nth_element( rArray.begin(), iter, rArray.end());
3422 36 : return fVal + fDiff * (*iter - fVal);
3423 : }
3424 : }
3425 :
3426 70 : void ScInterpreter::ScPercentile( bool bInclusive )
3427 : {
3428 70 : if ( !MustHaveParamCount( GetByte(), 2 ) )
3429 0 : return;
3430 70 : double alpha = GetDouble();
3431 70 : if ( bInclusive ? ( alpha < 0.0 || alpha > 1.0 ) : ( alpha <= 0.0 || alpha >= 1.0 ) )
3432 : {
3433 0 : PushIllegalArgument();
3434 0 : return;
3435 : }
3436 70 : vector<double> aArray;
3437 70 : GetNumberSequenceArray( 1, aArray, false );
3438 70 : if ( bInclusive )
3439 34 : PushDouble( GetPercentile( aArray, alpha ));
3440 : else
3441 36 : PushDouble( GetPercentileExclusive( aArray, alpha ));
3442 : }
3443 :
3444 67 : void ScInterpreter::ScQuartile( bool bInclusive )
3445 : {
3446 67 : if ( !MustHaveParamCount( GetByte(), 2 ) )
3447 0 : return;
3448 67 : double fFlag = ::rtl::math::approxFloor(GetDouble());
3449 67 : if ( bInclusive ? ( fFlag < 0.0 || fFlag > 4.0 ) : ( fFlag <= 0.0 || fFlag >= 4.0 ) )
3450 : {
3451 0 : PushIllegalArgument();
3452 0 : return;
3453 : }
3454 67 : vector<double> aArray;
3455 67 : GetNumberSequenceArray( 1, aArray, false );
3456 67 : if ( bInclusive )
3457 34 : PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentile( aArray, 0.25 * fFlag ) );
3458 : else
3459 33 : PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentileExclusive( aArray, 0.25 * fFlag ) );
3460 : }
3461 :
3462 61 : void ScInterpreter::ScModalValue()
3463 : {
3464 61 : sal_uInt8 nParamCount = GetByte();
3465 61 : if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3466 61 : return;
3467 61 : vector<double> aSortArray;
3468 61 : GetSortArray( nParamCount, aSortArray, NULL, false );
3469 61 : SCSIZE nSize = aSortArray.size();
3470 61 : if (aSortArray.empty() || nSize == 0 || nGlobalError)
3471 0 : PushNoValue();
3472 : else
3473 : {
3474 61 : SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
3475 61 : double nOldVal = aSortArray[0];
3476 : SCSIZE i;
3477 :
3478 604 : for ( i = 1; i < nSize; i++)
3479 : {
3480 543 : if (aSortArray[i] == nOldVal)
3481 109 : nCount++;
3482 : else
3483 : {
3484 434 : nOldVal = aSortArray[i];
3485 434 : if (nCount > nMax)
3486 : {
3487 52 : nMax = nCount;
3488 52 : nMaxIndex = i-1;
3489 : }
3490 434 : nCount = 1;
3491 : }
3492 : }
3493 61 : if (nCount > nMax)
3494 : {
3495 12 : nMax = nCount;
3496 12 : nMaxIndex = i-1;
3497 : }
3498 61 : if (nMax == 1 && nCount == 1)
3499 0 : PushNoValue();
3500 61 : else if (nMax == 1)
3501 0 : PushDouble(nOldVal);
3502 : else
3503 61 : PushDouble(aSortArray[nMaxIndex]);
3504 61 : }
3505 : }
3506 :
3507 78 : void ScInterpreter::CalculateSmallLarge(bool bSmall)
3508 : {
3509 78 : if ( !MustHaveParamCount( GetByte(), 2 ) )
3510 0 : return;
3511 78 : double f = ::rtl::math::approxFloor(GetDouble());
3512 78 : if (f < 1.0)
3513 : {
3514 0 : PushIllegalArgument();
3515 0 : return;
3516 : }
3517 78 : SCSIZE k = static_cast<SCSIZE>(f);
3518 78 : vector<double> aSortArray;
3519 : /* TODO: using nth_element() is best for one single value, but LARGE/SMALL
3520 : * actually are defined to return an array of values if an array of
3521 : * positions was passed, in which case, depending on the number of values,
3522 : * we may or will need a real sorted array again, see #i32345. */
3523 78 : GetNumberSequenceArray(1, aSortArray, false );
3524 78 : SCSIZE nSize = aSortArray.size();
3525 78 : if (aSortArray.empty() || nSize == 0 || nGlobalError || nSize < k)
3526 0 : PushNoValue();
3527 : else
3528 : {
3529 : // TODO: the sorted case for array: PushDouble( aSortArray[ bSmall ? k-1 : nSize-k ] );
3530 78 : vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
3531 78 : ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
3532 78 : PushDouble( *iPos);
3533 78 : }
3534 : }
3535 :
3536 39 : void ScInterpreter::ScLarge()
3537 : {
3538 39 : CalculateSmallLarge(false);
3539 39 : }
3540 :
3541 39 : void ScInterpreter::ScSmall()
3542 : {
3543 39 : CalculateSmallLarge(true);
3544 39 : }
3545 :
3546 8 : void ScInterpreter::ScPercentrank( bool bInclusive )
3547 : {
3548 8 : sal_uInt8 nParamCount = GetByte();
3549 8 : if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
3550 8 : return;
3551 8 : double fSignificance = ( nParamCount == 3 ? ::rtl::math::approxFloor( GetDouble() ) : 3.0 );
3552 8 : double fNum = GetDouble();
3553 8 : vector<double> aSortArray;
3554 8 : GetSortArray( 1, aSortArray, NULL, false );
3555 8 : SCSIZE nSize = aSortArray.size();
3556 8 : if ( aSortArray.empty() || nSize == 0 || nGlobalError )
3557 0 : PushNoValue();
3558 : else
3559 : {
3560 8 : if ( fNum < aSortArray[ 0 ] || fNum > aSortArray[ nSize - 1 ] )
3561 0 : PushNoValue();
3562 : else
3563 : {
3564 : double fRes;
3565 8 : if ( nSize == 1 )
3566 0 : fRes = 1.0; // fNum == aSortArray[ 0 ], see test above
3567 : else
3568 8 : fRes = GetPercentrank( aSortArray, fNum, bInclusive );
3569 8 : if ( fRes != 0.0 )
3570 : {
3571 5 : double fExp = ::rtl::math::approxFloor( log( fRes ) );
3572 5 : fRes = ::rtl::math::round( fRes * pow( 10, -fExp + fSignificance - 1 ) ) / pow( 10, -fExp + fSignificance - 1 );
3573 : }
3574 8 : PushDouble( fRes );
3575 : }
3576 8 : }
3577 : }
3578 :
3579 8 : double ScInterpreter::GetPercentrank( ::std::vector<double> & rArray, double fVal, bool bInclusive )
3580 : {
3581 8 : SCSIZE nSize = rArray.size();
3582 : double fRes;
3583 8 : if ( fVal == rArray[ 0 ] )
3584 : {
3585 3 : if ( bInclusive )
3586 3 : fRes = 0.0;
3587 : else
3588 0 : fRes = 1.0 / ( double )( nSize + 1 );
3589 : }
3590 : else
3591 : {
3592 5 : SCSIZE nOldCount = 0;
3593 5 : double fOldVal = rArray[ 0 ];
3594 : SCSIZE i;
3595 20 : for ( i = 1; i < nSize && rArray[ i ] < fVal; i++ )
3596 : {
3597 15 : if ( rArray[ i ] != fOldVal )
3598 : {
3599 15 : nOldCount = i;
3600 15 : fOldVal = rArray[ i ];
3601 : }
3602 : }
3603 5 : if ( rArray[ i ] != fOldVal )
3604 5 : nOldCount = i;
3605 5 : if ( fVal == rArray[ i ] )
3606 : {
3607 5 : if ( bInclusive )
3608 4 : fRes = div( nOldCount, nSize - 1 );
3609 : else
3610 1 : fRes = ( double )( i + 1 ) / ( double )( nSize + 1 );
3611 : }
3612 : else
3613 : {
3614 : // nOldCount is the count of smaller entries
3615 : // fVal is between rArray[ nOldCount - 1 ] and rArray[ nOldCount ]
3616 : // use linear interpolation to find a position between the entries
3617 0 : if ( nOldCount == 0 )
3618 : {
3619 : OSL_FAIL( "should not happen" );
3620 0 : fRes = 0.0;
3621 : }
3622 : else
3623 : {
3624 0 : double fFract = ( fVal - rArray[ nOldCount - 1 ] ) /
3625 0 : ( rArray[ nOldCount ] - rArray[ nOldCount - 1 ] );
3626 0 : if ( bInclusive )
3627 0 : fRes = div( ( double )( nOldCount - 1 ) + fFract, nSize - 1 );
3628 : else
3629 0 : fRes = ( ( double )nOldCount + fFract ) / ( double )( nSize + 1 );
3630 : }
3631 : }
3632 : }
3633 8 : return fRes;
3634 : }
3635 :
3636 3 : void ScInterpreter::ScTrimMean()
3637 : {
3638 3 : if ( !MustHaveParamCount( GetByte(), 2 ) )
3639 0 : return;
3640 3 : double alpha = GetDouble();
3641 3 : if (alpha < 0.0 || alpha >= 1.0)
3642 : {
3643 0 : PushIllegalArgument();
3644 0 : return;
3645 : }
3646 3 : vector<double> aSortArray;
3647 3 : GetSortArray( 1, aSortArray, NULL, false );
3648 3 : SCSIZE nSize = aSortArray.size();
3649 3 : if (aSortArray.empty() || nSize == 0 || nGlobalError)
3650 0 : PushNoValue();
3651 : else
3652 : {
3653 3 : sal_uLong nIndex = (sal_uLong) ::rtl::math::approxFloor(alpha*(double)nSize);
3654 3 : if (nIndex % 2 != 0)
3655 3 : nIndex--;
3656 3 : nIndex /= 2;
3657 : OSL_ENSURE(nIndex < nSize, "ScTrimMean: falscher Index");
3658 3 : double fSum = 0.0;
3659 33 : for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
3660 30 : fSum += aSortArray[i];
3661 3 : PushDouble(fSum/(double)(nSize-2*nIndex));
3662 3 : }
3663 : }
3664 :
3665 390 : void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray, bool bConvertTextInArray )
3666 : {
3667 390 : ScAddress aAdr;
3668 390 : ScRange aRange;
3669 390 : short nParam = nParamCount;
3670 390 : size_t nRefInList = 0;
3671 1191 : while (nParam-- > 0)
3672 : {
3673 411 : const StackVar eStackType = GetStackType();
3674 411 : switch (eStackType)
3675 : {
3676 : case formula::svDouble :
3677 27 : rArray.push_back( PopDouble());
3678 27 : break;
3679 : case svSingleRef :
3680 : {
3681 9 : PopSingleRef( aAdr );
3682 9 : ScRefCellValue aCell;
3683 9 : aCell.assign(*pDok, aAdr);
3684 9 : if (aCell.hasNumeric())
3685 9 : rArray.push_back(GetCellValue(aAdr, aCell));
3686 : }
3687 9 : break;
3688 : case formula::svDoubleRef :
3689 : case svRefList :
3690 : {
3691 362 : PopDoubleRef( aRange, nParam, nRefInList);
3692 362 : if (nGlobalError)
3693 0 : break;
3694 :
3695 362 : aRange.Justify();
3696 362 : SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3697 362 : nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3698 362 : rArray.reserve( rArray.size() + nCellCount);
3699 :
3700 362 : sal_uInt16 nErr = 0;
3701 : double fCellVal;
3702 362 : ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
3703 362 : if (aValIter.GetFirst( fCellVal, nErr))
3704 : {
3705 362 : rArray.push_back( fCellVal);
3706 362 : SetError(nErr);
3707 5174 : while ((nErr == 0) && aValIter.GetNext( fCellVal, nErr))
3708 4450 : rArray.push_back( fCellVal);
3709 362 : SetError(nErr);
3710 : }
3711 : }
3712 362 : break;
3713 : case svMatrix :
3714 : case svExternalSingleRef:
3715 : case svExternalDoubleRef:
3716 : {
3717 13 : ScMatrixRef pMat = GetMatrix();
3718 13 : if (!pMat)
3719 0 : break;
3720 :
3721 13 : SCSIZE nCount = pMat->GetElementCount();
3722 13 : rArray.reserve( rArray.size() + nCount);
3723 13 : if (pMat->IsNumeric())
3724 : {
3725 164 : for (SCSIZE i = 0; i < nCount; ++i)
3726 151 : rArray.push_back( pMat->GetDouble(i));
3727 : }
3728 0 : else if (bConvertTextInArray && eStackType == svMatrix)
3729 : {
3730 0 : for (SCSIZE i = 0; i < nCount; ++i)
3731 : {
3732 0 : if ( pMat->IsValue( i ) )
3733 0 : rArray.push_back( pMat->GetDouble(i));
3734 : else
3735 : {
3736 : // tdf#88547 try to convert string to (date)value
3737 0 : OUString aStr = pMat->GetString( i ).getString();
3738 0 : if ( aStr.getLength() > 0 )
3739 : {
3740 0 : sal_uInt16 nErr = nGlobalError;
3741 0 : nGlobalError = 0;
3742 0 : double fVal = ConvertStringToValue( aStr );
3743 0 : if ( !nGlobalError )
3744 : {
3745 0 : rArray.push_back( fVal );
3746 0 : nGlobalError = nErr;
3747 : }
3748 : else
3749 : {
3750 0 : rArray.push_back( CreateDoubleError( errNoValue));
3751 : // Propagate previous error if any, else
3752 : // the current #VALUE! error.
3753 0 : if (nErr)
3754 0 : nGlobalError = nErr;
3755 : else
3756 0 : nGlobalError = errNoValue;
3757 : }
3758 0 : }
3759 : }
3760 0 : }
3761 : }
3762 : else
3763 : {
3764 0 : for (SCSIZE i = 0; i < nCount; ++i)
3765 : {
3766 0 : if ( pMat->IsValue( i ) )
3767 0 : rArray.push_back( pMat->GetDouble(i));
3768 : }
3769 13 : }
3770 : }
3771 13 : break;
3772 : default :
3773 0 : PopError();
3774 0 : SetError( errIllegalParameter);
3775 0 : break;
3776 : }
3777 411 : if (nGlobalError)
3778 0 : break; // while
3779 : }
3780 : // nParam > 0 in case of error, clean stack environment and obtain earlier
3781 : // error if there was one.
3782 780 : while (nParam-- > 0)
3783 0 : PopError();
3784 390 : }
3785 :
3786 115 : void ScInterpreter::GetSortArray( sal_uInt8 nParamCount, vector<double>& rSortArray, vector<long>* pIndexOrder, bool bConvertTextInArray )
3787 : {
3788 115 : GetNumberSequenceArray( nParamCount, rSortArray, bConvertTextInArray );
3789 115 : if (rSortArray.size() > MAX_ANZ_DOUBLE_FOR_SORT)
3790 0 : SetError( errStackOverflow);
3791 115 : else if (rSortArray.empty())
3792 0 : SetError( errNoValue);
3793 :
3794 115 : if (nGlobalError == 0)
3795 115 : QuickSort( rSortArray, pIndexOrder);
3796 115 : }
3797 :
3798 2738 : static void lcl_QuickSort( long nLo, long nHi, vector<double>& rSortArray, vector<long>* pIndexOrder )
3799 : {
3800 : // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().
3801 :
3802 : using ::std::swap;
3803 :
3804 2738 : if (nHi - nLo == 1)
3805 : {
3806 1039 : if (rSortArray[nLo] > rSortArray[nHi])
3807 : {
3808 281 : swap(rSortArray[nLo], rSortArray[nHi]);
3809 281 : if (pIndexOrder)
3810 0 : swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
3811 : }
3812 3777 : return;
3813 : }
3814 :
3815 1699 : long ni = nLo;
3816 1699 : long nj = nHi;
3817 4511 : do
3818 : {
3819 4511 : double fLo = rSortArray[nLo];
3820 4511 : while (ni <= nHi && rSortArray[ni] < fLo) ni++;
3821 4511 : while (nj >= nLo && fLo < rSortArray[nj]) nj--;
3822 4511 : if (ni <= nj)
3823 : {
3824 3923 : if (ni != nj)
3825 : {
3826 3575 : swap(rSortArray[ni], rSortArray[nj]);
3827 3575 : if (pIndexOrder)
3828 6 : swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
3829 : }
3830 :
3831 3923 : ++ni;
3832 3923 : --nj;
3833 : }
3834 : }
3835 : while (ni < nj);
3836 :
3837 1699 : if ((nj - nLo) < (nHi - ni))
3838 : {
3839 1144 : if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
3840 1144 : if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
3841 : }
3842 : else
3843 : {
3844 555 : if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
3845 555 : if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
3846 : }
3847 : }
3848 :
3849 115 : void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<long>* pIndexOrder )
3850 : {
3851 115 : long n = static_cast<long>(rSortArray.size());
3852 :
3853 115 : if (pIndexOrder)
3854 : {
3855 3 : pIndexOrder->clear();
3856 3 : pIndexOrder->reserve(n);
3857 33 : for (long i = 0; i < n; ++i)
3858 30 : pIndexOrder->push_back(i);
3859 : }
3860 :
3861 115 : if (n < 2)
3862 124 : return;
3863 :
3864 106 : size_t nValCount = rSortArray.size();
3865 647 : for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
3866 : {
3867 541 : size_t nInd = comphelper::rng::uniform_size_distribution(0, nValCount-2);
3868 541 : ::std::swap( rSortArray[i], rSortArray[nInd]);
3869 541 : if (pIndexOrder)
3870 6 : ::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
3871 : }
3872 :
3873 106 : lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
3874 : }
3875 :
3876 13 : void ScInterpreter::ScRank( bool bAverage )
3877 : {
3878 13 : sal_uInt8 nParamCount = GetByte();
3879 13 : if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
3880 13 : return;
3881 : bool bAscending;
3882 13 : if ( nParamCount == 3 )
3883 6 : bAscending = GetBool();
3884 : else
3885 7 : bAscending = false;
3886 :
3887 13 : vector<double> aSortArray;
3888 13 : GetSortArray( 1, aSortArray, NULL, false );
3889 13 : double fVal = GetDouble();
3890 13 : SCSIZE nSize = aSortArray.size();
3891 13 : if ( aSortArray.empty() || nSize == 0 || nGlobalError )
3892 0 : PushNoValue();
3893 : else
3894 : {
3895 13 : if ( fVal < aSortArray[ 0 ] || fVal > aSortArray[ nSize - 1 ] )
3896 0 : PushNoValue();
3897 : else
3898 : {
3899 13 : double fLastPos = 0;
3900 13 : double fFirstPos = -1.0;
3901 13 : bool bFinished = false;
3902 : SCSIZE i;
3903 126 : for ( i = 0; i < nSize && !bFinished && !nGlobalError; i++ )
3904 : {
3905 113 : if ( aSortArray[ i ] == fVal )
3906 : {
3907 6 : if ( fFirstPos < 0 )
3908 6 : fFirstPos = i + 1.0;
3909 : }
3910 : else
3911 : {
3912 107 : if ( aSortArray[ i ] > fVal )
3913 : {
3914 10 : fLastPos = i;
3915 10 : bFinished = true;
3916 : }
3917 : }
3918 : }
3919 13 : if ( !bFinished )
3920 3 : fLastPos = i;
3921 13 : if ( !bAverage )
3922 : {
3923 7 : if ( bAscending )
3924 0 : PushDouble( fFirstPos );
3925 : else
3926 7 : PushDouble( nSize + 1.0 - fLastPos );
3927 : }
3928 : else
3929 : {
3930 6 : if ( bAscending )
3931 0 : PushDouble( ( fFirstPos + fLastPos ) / 2.0 );
3932 : else
3933 6 : PushDouble( nSize + 1.0 - ( fFirstPos + fLastPos ) / 2.0 );
3934 : }
3935 : }
3936 13 : }
3937 : }
3938 :
3939 3 : void ScInterpreter::ScAveDev()
3940 : {
3941 3 : sal_uInt8 nParamCount = GetByte();
3942 3 : if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3943 0 : return;
3944 3 : sal_uInt16 SaveSP = sp;
3945 3 : double nMiddle = 0.0;
3946 3 : double rVal = 0.0;
3947 3 : double rValCount = 0.0;
3948 3 : ScAddress aAdr;
3949 3 : ScRange aRange;
3950 3 : short nParam = nParamCount;
3951 3 : size_t nRefInList = 0;
3952 9 : while (nParam-- > 0)
3953 : {
3954 3 : switch (GetStackType())
3955 : {
3956 : case formula::svDouble :
3957 0 : rVal += GetDouble();
3958 0 : rValCount++;
3959 0 : break;
3960 : case svSingleRef :
3961 : {
3962 0 : PopSingleRef( aAdr );
3963 0 : ScRefCellValue aCell;
3964 0 : aCell.assign(*pDok, aAdr);
3965 0 : if (aCell.hasNumeric())
3966 : {
3967 0 : rVal += GetCellValue(aAdr, aCell);
3968 0 : rValCount++;
3969 0 : }
3970 : }
3971 0 : break;
3972 : case formula::svDoubleRef :
3973 : case svRefList :
3974 : {
3975 3 : sal_uInt16 nErr = 0;
3976 : double nCellVal;
3977 3 : PopDoubleRef( aRange, nParam, nRefInList);
3978 3 : ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
3979 3 : if (aValIter.GetFirst(nCellVal, nErr))
3980 : {
3981 3 : rVal += nCellVal;
3982 3 : rValCount++;
3983 3 : SetError(nErr);
3984 18 : while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
3985 : {
3986 12 : rVal += nCellVal;
3987 12 : rValCount++;
3988 : }
3989 3 : SetError(nErr);
3990 : }
3991 : }
3992 3 : break;
3993 : case svMatrix :
3994 : case svExternalSingleRef:
3995 : case svExternalDoubleRef:
3996 : {
3997 0 : ScMatrixRef pMat = GetMatrix();
3998 0 : if (pMat)
3999 : {
4000 0 : SCSIZE nCount = pMat->GetElementCount();
4001 0 : if (pMat->IsNumeric())
4002 : {
4003 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4004 : {
4005 0 : rVal += pMat->GetDouble(nElem);
4006 0 : rValCount++;
4007 : }
4008 : }
4009 : else
4010 : {
4011 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4012 0 : if (!pMat->IsString(nElem))
4013 : {
4014 0 : rVal += pMat->GetDouble(nElem);
4015 0 : rValCount++;
4016 : }
4017 : }
4018 0 : }
4019 : }
4020 0 : break;
4021 : default :
4022 0 : SetError(errIllegalParameter);
4023 0 : break;
4024 : }
4025 : }
4026 3 : if (nGlobalError)
4027 : {
4028 0 : PushError( nGlobalError);
4029 0 : return;
4030 : }
4031 3 : nMiddle = rVal / rValCount;
4032 3 : sp = SaveSP;
4033 3 : rVal = 0.0;
4034 3 : nParam = nParamCount;
4035 3 : nRefInList = 0;
4036 9 : while (nParam-- > 0)
4037 : {
4038 3 : switch (GetStackType())
4039 : {
4040 : case formula::svDouble :
4041 0 : rVal += fabs(GetDouble() - nMiddle);
4042 0 : break;
4043 : case svSingleRef :
4044 : {
4045 0 : PopSingleRef( aAdr );
4046 0 : ScRefCellValue aCell;
4047 0 : aCell.assign(*pDok, aAdr);
4048 0 : if (aCell.hasNumeric())
4049 0 : rVal += fabs(GetCellValue(aAdr, aCell) - nMiddle);
4050 : }
4051 0 : break;
4052 : case formula::svDoubleRef :
4053 : case svRefList :
4054 : {
4055 3 : sal_uInt16 nErr = 0;
4056 : double nCellVal;
4057 3 : PopDoubleRef( aRange, nParam, nRefInList);
4058 3 : ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
4059 3 : if (aValIter.GetFirst(nCellVal, nErr))
4060 : {
4061 3 : rVal += (fabs(nCellVal - nMiddle));
4062 18 : while (aValIter.GetNext(nCellVal, nErr))
4063 12 : rVal += fabs(nCellVal - nMiddle);
4064 : }
4065 : }
4066 3 : break;
4067 : case svMatrix :
4068 : case svExternalSingleRef:
4069 : case svExternalDoubleRef:
4070 : {
4071 0 : ScMatrixRef pMat = GetMatrix();
4072 0 : if (pMat)
4073 : {
4074 0 : SCSIZE nCount = pMat->GetElementCount();
4075 0 : if (pMat->IsNumeric())
4076 : {
4077 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4078 : {
4079 0 : rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
4080 : }
4081 : }
4082 : else
4083 : {
4084 0 : for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4085 : {
4086 0 : if (!pMat->IsString(nElem))
4087 0 : rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
4088 : }
4089 : }
4090 0 : }
4091 : }
4092 0 : break;
4093 0 : default : SetError(errIllegalParameter); break;
4094 : }
4095 : }
4096 3 : PushDouble(rVal / rValCount);
4097 : }
4098 :
4099 3 : void ScInterpreter::ScDevSq()
4100 : {
4101 : double nVal;
4102 : double nValCount;
4103 3 : GetStVarParams(nVal, nValCount);
4104 3 : PushDouble(nVal);
4105 3 : }
4106 :
4107 2 : void ScInterpreter::ScProbability()
4108 : {
4109 2 : sal_uInt8 nParamCount = GetByte();
4110 2 : if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
4111 2 : return;
4112 : double fUp, fLo;
4113 2 : fUp = GetDouble();
4114 2 : if (nParamCount == 4)
4115 1 : fLo = GetDouble();
4116 : else
4117 1 : fLo = fUp;
4118 2 : if (fLo > fUp)
4119 : {
4120 0 : double fTemp = fLo;
4121 0 : fLo = fUp;
4122 0 : fUp = fTemp;
4123 : }
4124 2 : ScMatrixRef pMatP = GetMatrix();
4125 4 : ScMatrixRef pMatW = GetMatrix();
4126 2 : if (!pMatP || !pMatW)
4127 0 : PushIllegalParameter();
4128 : else
4129 : {
4130 : SCSIZE nC1, nC2;
4131 : SCSIZE nR1, nR2;
4132 2 : pMatP->GetDimensions(nC1, nR1);
4133 2 : pMatW->GetDimensions(nC2, nR2);
4134 4 : if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
4135 4 : nC2 == 0 || nR2 == 0)
4136 0 : PushNA();
4137 : else
4138 : {
4139 2 : double fSum = 0.0;
4140 2 : double fRes = 0.0;
4141 2 : bool bStop = false;
4142 : double fP, fW;
4143 4 : for ( SCSIZE i = 0; i < nC1 && !bStop; i++ )
4144 : {
4145 10 : for (SCSIZE j = 0; j < nR1 && !bStop; ++j )
4146 : {
4147 8 : if (pMatP->IsValue(i,j) && pMatW->IsValue(i,j))
4148 : {
4149 8 : fP = pMatP->GetDouble(i,j);
4150 8 : fW = pMatW->GetDouble(i,j);
4151 8 : if (fP < 0.0 || fP > 1.0)
4152 0 : bStop = true;
4153 : else
4154 : {
4155 8 : fSum += fP;
4156 8 : if (fW >= fLo && fW <= fUp)
4157 3 : fRes += fP;
4158 : }
4159 : }
4160 : else
4161 0 : SetError( errIllegalArgument);
4162 : }
4163 : }
4164 2 : if (bStop || fabs(fSum -1.0) > 1.0E-7)
4165 0 : PushNoValue();
4166 : else
4167 2 : PushDouble(fRes);
4168 : }
4169 2 : }
4170 : }
4171 :
4172 3 : void ScInterpreter::ScCorrel()
4173 : {
4174 : // This is identical to ScPearson()
4175 3 : ScPearson();
4176 3 : }
4177 :
4178 9 : void ScInterpreter::ScCovarianceP()
4179 : {
4180 9 : CalculatePearsonCovar( false, false, false );
4181 9 : }
4182 :
4183 6 : void ScInterpreter::ScCovarianceS()
4184 : {
4185 6 : CalculatePearsonCovar( false, false, true );
4186 6 : }
4187 :
4188 12 : void ScInterpreter::ScPearson()
4189 : {
4190 12 : CalculatePearsonCovar( true, false, false );
4191 12 : }
4192 :
4193 30 : void ScInterpreter::CalculatePearsonCovar( bool _bPearson, bool _bStexy, bool _bSample )
4194 : {
4195 30 : if ( !MustHaveParamCount( GetByte(), 2 ) )
4196 0 : return;
4197 30 : ScMatrixRef pMat1 = GetMatrix();
4198 60 : ScMatrixRef pMat2 = GetMatrix();
4199 30 : if (!pMat1 || !pMat2)
4200 : {
4201 0 : PushIllegalParameter();
4202 0 : return;
4203 : }
4204 : SCSIZE nC1, nC2;
4205 : SCSIZE nR1, nR2;
4206 30 : pMat1->GetDimensions(nC1, nR1);
4207 30 : pMat2->GetDimensions(nC2, nR2);
4208 30 : if (nR1 != nR2 || nC1 != nC2)
4209 : {
4210 0 : PushIllegalArgument();
4211 0 : return;
4212 : }
4213 : /* #i78250#
4214 : * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
4215 : * but the latter produces wrong results if the absolute values are high,
4216 : * for example above 10^8
4217 : */
4218 30 : double fCount = 0.0;
4219 30 : double fSumX = 0.0;
4220 30 : double fSumY = 0.0;
4221 :
4222 60 : for (SCSIZE i = 0; i < nC1; i++)
4223 : {
4224 258 : for (SCSIZE j = 0; j < nR1; j++)
4225 : {
4226 228 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4227 : {
4228 219 : double fValX = pMat1->GetDouble(i,j);
4229 219 : double fValY = pMat2->GetDouble(i,j);
4230 219 : fSumX += fValX;
4231 219 : fSumY += fValY;
4232 219 : fCount++;
4233 : }
4234 : }
4235 : }
4236 30 : if (fCount < (_bStexy ? 3.0 : (_bSample ? 2.0 : 1.0)))
4237 0 : PushNoValue();
4238 : else
4239 : {
4240 30 : double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4241 30 : double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4242 30 : double fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2
4243 30 : const double fMeanX = fSumX / fCount;
4244 30 : const double fMeanY = fSumY / fCount;
4245 60 : for (SCSIZE i = 0; i < nC1; i++)
4246 : {
4247 258 : for (SCSIZE j = 0; j < nR1; j++)
4248 : {
4249 228 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4250 : {
4251 219 : const double fValX = pMat1->GetDouble(i,j);
4252 219 : const double fValY = pMat2->GetDouble(i,j);
4253 219 : fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4254 219 : if ( _bPearson )
4255 : {
4256 105 : fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4257 105 : fSumSqrDeltaY += (fValY - fMeanY) * (fValY - fMeanY);
4258 : }
4259 : }
4260 : }
4261 : }
4262 30 : if ( _bPearson )
4263 : {
4264 15 : if (fSumSqrDeltaX == 0.0 || ( !_bStexy && fSumSqrDeltaY == 0.0) )
4265 0 : PushError( errDivisionByZero);
4266 15 : else if ( _bStexy )
4267 6 : PushDouble( sqrt( (fSumSqrDeltaY - fSumDeltaXDeltaY *
4268 6 : fSumDeltaXDeltaY / fSumSqrDeltaX) / (fCount-2)));
4269 : else
4270 12 : PushDouble( fSumDeltaXDeltaY / sqrt( fSumSqrDeltaX * fSumSqrDeltaY));
4271 : }
4272 : else
4273 : {
4274 15 : if ( _bSample )
4275 6 : PushDouble( fSumDeltaXDeltaY / ( fCount - 1 ) );
4276 : else
4277 9 : PushDouble( fSumDeltaXDeltaY / fCount);
4278 : }
4279 30 : }
4280 : }
4281 :
4282 6 : void ScInterpreter::ScRSQ()
4283 : {
4284 : // Same as ScPearson()*ScPearson()
4285 6 : ScPearson();
4286 6 : if (!nGlobalError)
4287 : {
4288 6 : switch (GetStackType())
4289 : {
4290 : case formula::svDouble:
4291 : {
4292 6 : double fVal = PopDouble();
4293 6 : PushDouble( fVal * fVal);
4294 : }
4295 6 : break;
4296 : default:
4297 0 : PopError();
4298 0 : PushNoValue();
4299 : }
4300 : }
4301 6 : }
4302 :
4303 3 : void ScInterpreter::ScSTEYX()
4304 : {
4305 3 : CalculatePearsonCovar( true, true, false );
4306 3 : }
4307 9 : void ScInterpreter::CalculateSlopeIntercept(bool bSlope)
4308 : {
4309 9 : if ( !MustHaveParamCount( GetByte(), 2 ) )
4310 0 : return;
4311 9 : ScMatrixRef pMat1 = GetMatrix();
4312 18 : ScMatrixRef pMat2 = GetMatrix();
4313 9 : if (!pMat1 || !pMat2)
4314 : {
4315 0 : PushIllegalParameter();
4316 0 : return;
4317 : }
4318 : SCSIZE nC1, nC2;
4319 : SCSIZE nR1, nR2;
4320 9 : pMat1->GetDimensions(nC1, nR1);
4321 9 : pMat2->GetDimensions(nC2, nR2);
4322 9 : if (nR1 != nR2 || nC1 != nC2)
4323 : {
4324 0 : PushIllegalArgument();
4325 0 : return;
4326 : }
4327 : // #i78250# numerical stability improved
4328 9 : double fCount = 0.0;
4329 9 : double fSumX = 0.0;
4330 9 : double fSumY = 0.0;
4331 :
4332 18 : for (SCSIZE i = 0; i < nC1; i++)
4333 : {
4334 81 : for (SCSIZE j = 0; j < nR1; j++)
4335 : {
4336 72 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4337 : {
4338 63 : double fValX = pMat1->GetDouble(i,j);
4339 63 : double fValY = pMat2->GetDouble(i,j);
4340 63 : fSumX += fValX;
4341 63 : fSumY += fValY;
4342 63 : fCount++;
4343 : }
4344 : }
4345 : }
4346 9 : if (fCount < 1.0)
4347 0 : PushNoValue();
4348 : else
4349 : {
4350 9 : double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4351 9 : double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4352 9 : double fMeanX = fSumX / fCount;
4353 9 : double fMeanY = fSumY / fCount;
4354 18 : for (SCSIZE i = 0; i < nC1; i++)
4355 : {
4356 81 : for (SCSIZE j = 0; j < nR1; j++)
4357 : {
4358 72 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4359 : {
4360 63 : double fValX = pMat1->GetDouble(i,j);
4361 63 : double fValY = pMat2->GetDouble(i,j);
4362 63 : fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4363 63 : fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4364 : }
4365 : }
4366 : }
4367 9 : if (fSumSqrDeltaX == 0.0)
4368 0 : PushError( errDivisionByZero);
4369 : else
4370 : {
4371 9 : if ( bSlope )
4372 3 : PushDouble( fSumDeltaXDeltaY / fSumSqrDeltaX);
4373 : else
4374 6 : PushDouble( fMeanY - fSumDeltaXDeltaY / fSumSqrDeltaX * fMeanX);
4375 : }
4376 9 : }
4377 : }
4378 :
4379 3 : void ScInterpreter::ScSlope()
4380 : {
4381 3 : CalculateSlopeIntercept(true);
4382 3 : }
4383 :
4384 6 : void ScInterpreter::ScIntercept()
4385 : {
4386 6 : CalculateSlopeIntercept(false);
4387 6 : }
4388 :
4389 3 : void ScInterpreter::ScForecast()
4390 : {
4391 3 : if ( !MustHaveParamCount( GetByte(), 3 ) )
4392 0 : return;
4393 3 : ScMatrixRef pMat1 = GetMatrix();
4394 6 : ScMatrixRef pMat2 = GetMatrix();
4395 3 : if (!pMat1 || !pMat2)
4396 : {
4397 0 : PushIllegalParameter();
4398 0 : return;
4399 : }
4400 : SCSIZE nC1, nC2;
4401 : SCSIZE nR1, nR2;
4402 3 : pMat1->GetDimensions(nC1, nR1);
4403 3 : pMat2->GetDimensions(nC2, nR2);
4404 3 : if (nR1 != nR2 || nC1 != nC2)
4405 : {
4406 0 : PushIllegalArgument();
4407 0 : return;
4408 : }
4409 3 : double fVal = GetDouble();
4410 : // #i78250# numerical stability improved
4411 3 : double fCount = 0.0;
4412 3 : double fSumX = 0.0;
4413 3 : double fSumY = 0.0;
4414 :
4415 6 : for (SCSIZE i = 0; i < nC1; i++)
4416 : {
4417 24 : for (SCSIZE j = 0; j < nR1; j++)
4418 : {
4419 21 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4420 : {
4421 21 : double fValX = pMat1->GetDouble(i,j);
4422 21 : double fValY = pMat2->GetDouble(i,j);
4423 21 : fSumX += fValX;
4424 21 : fSumY += fValY;
4425 21 : fCount++;
4426 : }
4427 : }
4428 : }
4429 3 : if (fCount < 1.0)
4430 0 : PushNoValue();
4431 : else
4432 : {
4433 3 : double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4434 3 : double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4435 3 : double fMeanX = fSumX / fCount;
4436 3 : double fMeanY = fSumY / fCount;
4437 6 : for (SCSIZE i = 0; i < nC1; i++)
4438 : {
4439 24 : for (SCSIZE j = 0; j < nR1; j++)
4440 : {
4441 21 : if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4442 : {
4443 21 : double fValX = pMat1->GetDouble(i,j);
4444 21 : double fValY = pMat2->GetDouble(i,j);
4445 21 : fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4446 21 : fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4447 : }
4448 : }
4449 : }
4450 3 : if (fSumSqrDeltaX == 0.0)
4451 0 : PushError( errDivisionByZero);
4452 : else
4453 3 : PushDouble( fMeanY + fSumDeltaXDeltaY / fSumSqrDeltaX * (fVal - fMeanX));
4454 3 : }
4455 156 : }
4456 :
4457 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
