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