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