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