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