LibreOffice Module sc (master) 1
interpr3.cxx
Go to the documentation of this file.
1/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2/*
3 * This file is part of the LibreOffice project.
4 *
5 * This Source Code Form is subject to the terms of the Mozilla Public
6 * License, v. 2.0. If a copy of the MPL was not distributed with this
7 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 *
9 * This file incorporates work covered by the following license notice:
10 *
11 * Licensed to the Apache Software Foundation (ASF) under one or more
12 * contributor license agreements. See the NOTICE file distributed
13 * with this work for additional information regarding copyright
14 * ownership. The ASF licenses this file to you under the Apache
15 * License, Version 2.0 (the "License"); you may not use this file
16 * except in compliance with the License. You may obtain a copy of
17 * the License at http://www.apache.org/licenses/LICENSE-2.0 .
18 */
19
20#include <tools/solar.h>
21#include <stdlib.h>
22
23#include <interpre.hxx>
24#include <global.hxx>
25#include <document.hxx>
26#include <dociter.hxx>
27#include <matrixoperators.hxx>
28#include <scmatrix.hxx>
29
30#include <cassert>
31#include <cmath>
32#include <memory>
33#include <set>
34#include <vector>
35#include <algorithm>
36#include <comphelper/random.hxx>
38#include <osl/diagnose.h>
39
40using ::std::vector;
41using namespace formula;
42
44// This is an arbitrary limit.
45static size_t MAX_COUNT_DOUBLE_FOR_SORT(const ScSheetLimits& rSheetLimits)
46{
47 return rSheetLimits.GetMaxRowCount() * 2;
48}
49
50const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental
51const double fMachEps = ::std::numeric_limits<double>::epsilon();
52
53namespace {
54
55class ScDistFunc
56{
57public:
58 virtual double GetValue(double x) const = 0;
59
60protected:
61 ~ScDistFunc() {}
62};
63
64}
65
66// iteration for inverse distributions
67
68//template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, bool& rConvError )
69
71static bool lcl_HasChangeOfSign( double u, double w )
72{
73 return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);
74}
75
76static double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError )
77{
78 rConvError = false;
79 const double fYEps = 1.0E-307;
80 const double fXEps = ::std::numeric_limits<double>::epsilon();
81
82 OSL_ENSURE(fAx<fBx, "IterateInverse: wrong interval");
83
84 // find enclosing interval
85
86 KahanSum fkAx = fAx;
87 KahanSum fkBx = fBx;
88 double fAy = rFunction.GetValue(fAx);
89 double fBy = rFunction.GetValue(fBx);
90 KahanSum fTemp;
91 unsigned short nCount;
92 for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++)
93 {
94 if (std::abs(fAy) <= std::abs(fBy))
95 {
96 fTemp = fkAx;
97 fkAx += (fkAx - fkBx) * 2.0;
98 if (fkAx < 0.0)
99 fkAx = 0.0;
100 fkBx = fTemp;
101 fBy = fAy;
102 fAy = rFunction.GetValue(fkAx.get());
103 }
104 else
105 {
106 fTemp = fkBx;
107 fkBx += (fkBx - fkAx) * 2.0;
108 fkAx = fTemp;
109 fAy = fBy;
110 fBy = rFunction.GetValue(fkBx.get());
111 }
112 }
113
114 fAx = fkAx.get();
115 fBx = fkBx.get();
116 if (fAy == 0.0)
117 return fAx;
118 if (fBy == 0.0)
119 return fBx;
120 if (!lcl_HasChangeOfSign( fAy, fBy))
121 {
122 rConvError = true;
123 return 0.0;
124 }
125 // inverse quadric interpolation with additional brackets
126 // set three points
127 double fPx = fAx;
128 double fPy = fAy;
129 double fQx = fBx;
130 double fQy = fBy;
131 double fRx = fAx;
132 double fRy = fAy;
133 double fSx = 0.5 * (fAx + fBx); // potential next point
134 bool bHasToInterpolate = true;
135 nCount = 0;
136 while ( nCount < 500 && std::abs(fRy) > fYEps &&
137 (fBx-fAx) > ::std::max( std::abs(fAx), std::abs(fBx)) * fXEps )
138 {
139 if (bHasToInterpolate)
140 {
141 if (fPy!=fQy && fQy!=fRy && fRy!=fPy)
142 {
143 fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)
144 + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)
145 + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);
146 bHasToInterpolate = (fAx < fSx) && (fSx < fBx); // inside the brackets?
147 }
148 else
149 bHasToInterpolate = false;
150 }
151 if(!bHasToInterpolate)
152 {
153 fSx = 0.5 * (fAx + fBx);
154 // reset points
155 fQx = fBx; fQy = fBy;
156 bHasToInterpolate = true;
157 }
158 // shift points for next interpolation
159 fPx = fQx; fQx = fRx; fRx = fSx;
160 fPy = fQy; fQy = fRy; fRy = rFunction.GetValue(fSx);
161 // update brackets
162 if (lcl_HasChangeOfSign( fAy, fRy))
163 {
164 fBx = fRx; fBy = fRy;
165 }
166 else
167 {
168 fAx = fRx; fAy = fRy;
169 }
170 // if last iteration brought too small advance, then do bisection next
171 // time, for safety
172 bHasToInterpolate = bHasToInterpolate && (std::abs(fRy) * 2.0 <= std::abs(fQy));
173 ++nCount;
174 }
175 return fRx;
176}
177
178// General functions
179
181{
182 PushError(FormulaError::NoName);
183}
184
186{
187 short nParamCount = GetByte();
188 while (nParamCount-- > 0)
189 {
190 PopError();
191 }
192 PushError( FormulaError::NoName);
193}
194
195double ScInterpreter::phi(double x)
196{
197 return 0.39894228040143268 * exp(-(x * x) / 2.0);
198}
199
201{ // Using gauss(x)+0.5 has severe cancellation errors for x<-4
202 return 0.5 * std::erfc(-x * 0.7071067811865475); // * 1/sqrt(2)
203}
204
205double ScInterpreter::taylor(const double* pPolynom, sal_uInt16 nMax, double x)
206{
207 KahanSum nVal = pPolynom[nMax];
208 for (short i = nMax-1; i >= 0; i--)
209 {
210 nVal = (nVal * x) + pPolynom[i];
211 }
212 return nVal.get();
213}
214
215double ScInterpreter::gauss(double x)
216{
217
218 double xAbs = std::abs(x);
219 sal_uInt16 xShort = static_cast<sal_uInt16>(::rtl::math::approxFloor(xAbs));
220 double nVal = 0.0;
221 if (xShort == 0)
222 {
223 static const double t0[] =
224 { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,
225 -0.00118732821548045, 0.00011543468761616, -0.00000944465625950,
226 0.00000066596935163, -0.00000004122667415, 0.00000000227352982,
227 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };
228 nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;
229 }
230 else if (xShort <= 2)
231 {
232 static const double t2[] =
233 { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,
234 0.02699548325659403, -0.00449924720943234, -0.00224962360471617,
235 0.00134977416282970, -0.00011783742691370, -0.00011515930357476,
236 0.00003704737285544, 0.00000282690796889, -0.00000354513195524,
237 0.00000037669563126, 0.00000019202407921, -0.00000005226908590,
238 -0.00000000491799345, 0.00000000366377919, -0.00000000015981997,
239 -0.00000000017381238, 0.00000000002624031, 0.00000000000560919,
240 -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };
241 nVal = taylor(t2, 23, (xAbs - 2.0));
242 }
243 else if (xShort <= 4)
244 {
245 static const double t4[] =
246 { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,
247 0.00033457556441221, -0.00028996548915725, 0.00018178605666397,
248 -0.00008252863922168, 0.00002551802519049, -0.00000391665839292,
249 -0.00000074018205222, 0.00000064422023359, -0.00000017370155340,
250 0.00000000909595465, 0.00000000944943118, -0.00000000329957075,
251 0.00000000029492075, 0.00000000011874477, -0.00000000004420396,
252 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };
253 nVal = taylor(t4, 20, (xAbs - 4.0));
254 }
255 else
256 {
257 static const double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };
258 nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs;
259 }
260 if (x < 0.0)
261 return -nVal;
262 else
263 return nVal;
264}
265
266// #i26836# new gaussinv implementation by Martin Eitzenberger <m.eitzenberger@unix.net>
267
269{
270 double q,t,z;
271
272 q=x-0.5;
273
274 if(std::abs(q)<=.425)
275 {
276 t=0.180625-q*q;
277
278 z=
279 q*
280 (
281 (
282 (
283 (
284 (
285 (
286 (
287 t*2509.0809287301226727+33430.575583588128105
288 )
289 *t+67265.770927008700853
290 )
291 *t+45921.953931549871457
292 )
293 *t+13731.693765509461125
294 )
295 *t+1971.5909503065514427
296 )
297 *t+133.14166789178437745
298 )
299 *t+3.387132872796366608
300 )
301 /
302 (
303 (
304 (
305 (
306 (
307 (
308 (
309 t*5226.495278852854561+28729.085735721942674
310 )
311 *t+39307.89580009271061
312 )
313 *t+21213.794301586595867
314 )
315 *t+5394.1960214247511077
316 )
317 *t+687.1870074920579083
318 )
319 *t+42.313330701600911252
320 )
321 *t+1.0
322 );
323
324 }
325 else
326 {
327 if(q>0) t=1-x;
328 else t=x;
329
330 t=sqrt(-log(t));
331
332 if(t<=5.0)
333 {
334 t+=-1.6;
335
336 z=
337 (
338 (
339 (
340 (
341 (
342 (
343 (
344 t*7.7454501427834140764e-4+0.0227238449892691845833
345 )
346 *t+0.24178072517745061177
347 )
348 *t+1.27045825245236838258
349 )
350 *t+3.64784832476320460504
351 )
352 *t+5.7694972214606914055
353 )
354 *t+4.6303378461565452959
355 )
356 *t+1.42343711074968357734
357 )
358 /
359 (
360 (
361 (
362 (
363 (
364 (
365 (
366 t*1.05075007164441684324e-9+5.475938084995344946e-4
367 )
368 *t+0.0151986665636164571966
369 )
370 *t+0.14810397642748007459
371 )
372 *t+0.68976733498510000455
373 )
374 *t+1.6763848301838038494
375 )
376 *t+2.05319162663775882187
377 )
378 *t+1.0
379 );
380
381 }
382 else
383 {
384 t+=-5.0;
385
386 z=
387 (
388 (
389 (
390 (
391 (
392 (
393 (
394 t*2.01033439929228813265e-7+2.71155556874348757815e-5
395 )
396 *t+0.0012426609473880784386
397 )
398 *t+0.026532189526576123093
399 )
400 *t+0.29656057182850489123
401 )
402 *t+1.7848265399172913358
403 )
404 *t+5.4637849111641143699
405 )
406 *t+6.6579046435011037772
407 )
408 /
409 (
410 (
411 (
412 (
413 (
414 (
415 (
416 t*2.04426310338993978564e-15+1.4215117583164458887e-7
417 )
418 *t+1.8463183175100546818e-5
419 )
420 *t+7.868691311456132591e-4
421 )
422 *t+0.0148753612908506148525
423 )
424 *t+0.13692988092273580531
425 )
426 *t+0.59983220655588793769
427 )
428 *t+1.0
429 );
430
431 }
432
433 if(q<0.0) z=-z;
434 }
435
436 return z;
437}
438
440{
441 x = ::rtl::math::approxFloor(x);
442 if (x < 0.0)
443 return 0.0;
444 else if (x == 0.0)
445 return 1.0;
446 else if (x <= 170.0)
447 {
448 double fTemp = x;
449 while (fTemp > 2.0)
450 {
451 fTemp--;
452 x *= fTemp;
453 }
454 }
455 else
456 SetError(FormulaError::NoValue);
457 return x;
458}
459
460double ScInterpreter::BinomKoeff(double n, double k)
461{
462 // this method has been duplicated as BinomialCoefficient()
463 // in scaddins/source/analysis/analysishelper.cxx
464
465 double nVal = 0.0;
466 k = ::rtl::math::approxFloor(k);
467 if (n < k)
468 nVal = 0.0;
469 else if (k == 0.0)
470 nVal = 1.0;
471 else
472 {
473 nVal = n/k;
474 n--;
475 k--;
476 while (k > 0.0)
477 {
478 nVal *= n/k;
479 k--;
480 n--;
481 }
482
483 }
484 return nVal;
485}
486
487// The algorithm is based on lanczos13m53 in lanczos.hpp
488// in math library from http://www.boost.org
491static double lcl_getLanczosSum(double fZ)
492{
493 static const double fNum[13] ={
494 23531376880.41075968857200767445163675473,
495 42919803642.64909876895789904700198885093,
496 35711959237.35566804944018545154716670596,
497 17921034426.03720969991975575445893111267,
498 6039542586.35202800506429164430729792107,
499 1439720407.311721673663223072794912393972,
500 248874557.8620541565114603864132294232163,
501 31426415.58540019438061423162831820536287,
502 2876370.628935372441225409051620849613599,
503 186056.2653952234950402949897160456992822,
504 8071.672002365816210638002902272250613822,
505 210.8242777515793458725097339207133627117,
506 2.506628274631000270164908177133837338626
507 };
508 static const double fDenom[13] = {
509 0,
510 39916800,
511 120543840,
512 150917976,
513 105258076,
514 45995730,
515 13339535,
516 2637558,
517 357423,
518 32670,
519 1925,
520 66,
521 1
522 };
523 // Horner scheme
524 double fSumNum;
525 double fSumDenom;
526 int nI;
527 if (fZ<=1.0)
528 {
529 fSumNum = fNum[12];
530 fSumDenom = fDenom[12];
531 for (nI = 11; nI >= 0; --nI)
532 {
533 fSumNum *= fZ;
534 fSumNum += fNum[nI];
535 fSumDenom *= fZ;
536 fSumDenom += fDenom[nI];
537 }
538 }
539 else
540 // Cancel down with fZ^12; Horner scheme with reverse coefficients
541 {
542 double fZInv = 1/fZ;
543 fSumNum = fNum[0];
544 fSumDenom = fDenom[0];
545 for (nI = 1; nI <=12; ++nI)
546 {
547 fSumNum *= fZInv;
548 fSumNum += fNum[nI];
549 fSumDenom *= fZInv;
550 fSumDenom += fDenom[nI];
551 }
552 }
553 return fSumNum/fSumDenom;
554}
555
556// The algorithm is based on tgamma in gamma.hpp
557// in math library from http://www.boost.org
559static double lcl_GetGammaHelper(double fZ)
560{
561 double fGamma = lcl_getLanczosSum(fZ);
562 const double fg = 6.024680040776729583740234375;
563 double fZgHelp = fZ + fg - 0.5;
564 // avoid intermediate overflow
565 double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25);
566 fGamma *= fHalfpower;
567 fGamma /= exp(fZgHelp);
568 fGamma *= fHalfpower;
569 if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ))
570 fGamma = ::rtl::math::round(fGamma);
571 return fGamma;
572}
573
574// The algorithm is based on tgamma in gamma.hpp
575// in math library from http://www.boost.org
577static double lcl_GetLogGammaHelper(double fZ)
578{
579 const double fg = 6.024680040776729583740234375;
580 double fZgHelp = fZ + fg - 0.5;
581 return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;
582}
583
585double ScInterpreter::GetGamma(double fZ)
586{
587 const double fLogPi = log(M_PI);
588 const double fLogDblMax = log( ::std::numeric_limits<double>::max());
589
590 if (fZ > fMaxGammaArgument)
591 {
592 SetError(FormulaError::IllegalFPOperation);
593 return HUGE_VAL;
594 }
595
596 if (fZ >= 1.0)
597 return lcl_GetGammaHelper(fZ);
598
599 if (fZ >= 0.5) // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
600 return lcl_GetGammaHelper(fZ+1) / fZ;
601
602 if (fZ >= -0.5) // shift to x>=1, might overflow
603 {
604 double fLogTest = lcl_GetLogGammaHelper(fZ+2) - std::log1p(fZ) - log( std::abs(fZ));
605 if (fLogTest >= fLogDblMax)
606 {
607 SetError( FormulaError::IllegalFPOperation);
608 return HUGE_VAL;
609 }
610 return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
611 }
612 // fZ<-0.5
613 // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) )
614 double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( std::abs( ::rtl::math::sin( M_PI*fZ)));
615 if (fLogDivisor - fLogPi >= fLogDblMax) // underflow
616 return 0.0;
617
618 if (fLogDivisor<0.0)
619 if (fLogPi - fLogDivisor > fLogDblMax) // overflow
620 {
621 SetError(FormulaError::IllegalFPOperation);
622 return HUGE_VAL;
623 }
624
625 return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( M_PI*fZ) < 0.0) ? -1.0 : 1.0);
626}
627
630{
631 if (fZ >= fMaxGammaArgument)
632 return lcl_GetLogGammaHelper(fZ);
633 if (fZ >= 1.0)
634 return log(lcl_GetGammaHelper(fZ));
635 if (fZ >= 0.5)
636 return log( lcl_GetGammaHelper(fZ+1) / fZ);
637 return lcl_GetLogGammaHelper(fZ+2) - std::log1p(fZ) - log(fZ);
638}
639
640double ScInterpreter::GetFDist(double x, double fF1, double fF2)
641{
642 double arg = fF2/(fF2+fF1*x);
643 double alpha = fF2/2.0;
644 double beta = fF1/2.0;
645 return GetBetaDist(arg, alpha, beta);
646}
647
648double ScInterpreter::GetTDist( double T, double fDF, int nType )
649{
650 switch ( nType )
651 {
652 case 1 : // 1-tailed T-distribution
653 return 0.5 * GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5 );
654 case 2 : // 2-tailed T-distribution
655 return GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5);
656 case 3 : // left-tailed T-distribution (probability density function)
657 return pow( 1 + ( T * T / fDF ), -( fDF + 1 ) / 2 ) / ( sqrt( fDF ) * GetBeta( 0.5, fDF / 2.0 ) );
658 case 4 : // left-tailed T-distribution (cumulative distribution function)
659 double X = fDF / ( T * T + fDF );
660 double R = 0.5 * GetBetaDist( X, 0.5 * fDF, 0.5 );
661 return ( T < 0 ? R : 1 - R );
662 }
663 SetError( FormulaError::IllegalArgument );
664 return HUGE_VAL;
665}
666
667// for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom
669double ScInterpreter::GetChiDist(double fX, double fDF)
670{
671 if (fX <= 0.0)
672 return 1.0; // see ODFF
673 else
674 return GetUpRegIGamma( fDF/2.0, fX/2.0);
675}
676
677// ready for ODF 1.2
678// for ODF CHISQDIST; cumulative distribution function, fDF=degrees of freedom
679// returns left tail
681double ScInterpreter::GetChiSqDistCDF(double fX, double fDF)
682{
683 if (fX <= 0.0)
684 return 0.0; // see ODFF
685 else
686 return GetLowRegIGamma( fDF/2.0, fX/2.0);
687}
688
689double ScInterpreter::GetChiSqDistPDF(double fX, double fDF)
690{
691 // you must ensure fDF is positive integer
692 double fValue;
693 if (fX <= 0.0)
694 return 0.0; // see ODFF
695 if (fDF*fX > 1391000.0)
696 {
697 // intermediate invalid values, use log
698 fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0) - GetLogGamma(0.5*fDF));
699 }
700 else // fDF is small in most cases, we can iterate
701 {
702 double fCount;
703 if (fmod(fDF,2.0)<0.5)
704 {
705 // even
706 fValue = 0.5;
707 fCount = 2.0;
708 }
709 else
710 {
711 fValue = 1/sqrt(fX*2*M_PI);
712 fCount = 1.0;
713 }
714 while ( fCount < fDF)
715 {
716 fValue *= (fX / fCount);
717 fCount += 2.0;
718 }
719 if (fX>=1425.0) // underflow in e^(-x/2)
720 fValue = exp(log(fValue)-fX/2);
721 else
722 fValue *= exp(-fX/2);
723 }
724 return fValue;
725}
726
728{
729 sal_uInt8 nParamCount = GetByte();
730 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
731 return;
732 bool bCumulative;
733 if (nParamCount == 3)
734 bCumulative = GetBool();
735 else
736 bCumulative = true;
737 double fDF = ::rtl::math::approxFloor(GetDouble());
738 if (fDF < 1.0)
740 else
741 {
742 double fX = GetDouble();
743 if (bCumulative)
745 else
747 }
748}
749
751{
752 sal_uInt8 nParamCount = GetByte();
753 if ( !MustHaveParamCount( nParamCount, 3, 3 ) )
754 return;
755 bool bCumulative = GetBool();
756 double fDF = ::rtl::math::approxFloor( GetDouble() );
757 if ( fDF < 1.0 || fDF > 1E10 )
759 else
760 {
761 double fX = GetDouble();
762 if ( fX < 0 )
764 else
765 {
766 if ( bCumulative )
767 PushDouble( GetChiSqDistCDF( fX, fDF ) );
768 else
769 PushDouble( GetChiSqDistPDF( fX, fDF ) );
770 }
771 }
772}
773
775{
776 double x = GetDouble();
777 if (x <= 0.0 && x == ::rtl::math::approxFloor(x))
779 else
780 {
781 double fResult = GetGamma(x);
782 if (nGlobalError != FormulaError::NONE)
783 {
785 return;
786 }
787 PushDouble(fResult);
788 }
789}
790
792{
793 double x = GetDouble();
794 if (x > 0.0) // constraint from ODFF
796 else
798}
799
800double ScInterpreter::GetBeta(double fAlpha, double fBeta)
801{
802 double fA;
803 double fB;
804 if (fAlpha > fBeta)
805 {
806 fA = fAlpha; fB = fBeta;
807 }
808 else
809 {
810 fA = fBeta; fB = fAlpha;
811 }
812 if (fA+fB < fMaxGammaArgument) // simple case
813 return GetGamma(fA)/GetGamma(fA+fB)*GetGamma(fB);
814 // need logarithm
815 // GetLogGamma is not accurate enough, back to Lanczos for all three
816 // GetGamma and arrange factors newly.
817 const double fg = 6.024680040776729583740234375; //see GetGamma
818 double fgm = fg - 0.5;
819 double fLanczos = lcl_getLanczosSum(fA);
820 fLanczos /= lcl_getLanczosSum(fA+fB);
821 fLanczos *= lcl_getLanczosSum(fB);
822 double fABgm = fA+fB+fgm;
823 fLanczos *= sqrt((fABgm/(fA+fgm))/(fB+fgm));
824 double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
825 double fTempB = fA/(fB+fgm);
826 double fResult = exp(-fA * std::log1p(fTempA)
827 -fB * std::log1p(fTempB)-fgm);
828 fResult *= fLanczos;
829 return fResult;
830}
831
832// Same as GetBeta but with logarithm
833double ScInterpreter::GetLogBeta(double fAlpha, double fBeta)
834{
835 double fA;
836 double fB;
837 if (fAlpha > fBeta)
838 {
839 fA = fAlpha; fB = fBeta;
840 }
841 else
842 {
843 fA = fBeta; fB = fAlpha;
844 }
845 const double fg = 6.024680040776729583740234375; //see GetGamma
846 double fgm = fg - 0.5;
847 double fLanczos = lcl_getLanczosSum(fA);
848 fLanczos /= lcl_getLanczosSum(fA+fB);
849 fLanczos *= lcl_getLanczosSum(fB);
850 double fLogLanczos = log(fLanczos);
851 double fABgm = fA+fB+fgm;
852 fLogLanczos += 0.5*(log(fABgm)-log(fA+fgm)-log(fB+fgm));
853 double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
854 double fTempB = fA/(fB+fgm);
855 double fResult = -fA * std::log1p(fTempA)
856 -fB * std::log1p(fTempB)-fgm;
857 fResult += fLogLanczos;
858 return fResult;
859}
860
861// beta distribution probability density function
862double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB)
863{
864 // special cases
865 if (fA == 1.0) // result b*(1-x)^(b-1)
866 {
867 if (fB == 1.0)
868 return 1.0;
869 if (fB == 2.0)
870 return -2.0*fX + 2.0;
871 if (fX == 1.0 && fB < 1.0)
872 {
873 SetError(FormulaError::IllegalArgument);
874 return HUGE_VAL;
875 }
876 if (fX <= 0.01)
877 return fB + fB * std::expm1((fB-1.0) * std::log1p(-fX));
878 else
879 return fB * pow(0.5-fX+0.5,fB-1.0);
880 }
881 if (fB == 1.0) // result a*x^(a-1)
882 {
883 if (fA == 2.0)
884 return fA * fX;
885 if (fX == 0.0 && fA < 1.0)
886 {
887 SetError(FormulaError::IllegalArgument);
888 return HUGE_VAL;
889 }
890 return fA * pow(fX,fA-1);
891 }
892 if (fX <= 0.0)
893 {
894 if (fA < 1.0 && fX == 0.0)
895 {
896 SetError(FormulaError::IllegalArgument);
897 return HUGE_VAL;
898 }
899 else
900 return 0.0;
901 }
902 if (fX >= 1.0)
903 {
904 if (fB < 1.0 && fX == 1.0)
905 {
906 SetError(FormulaError::IllegalArgument);
907 return HUGE_VAL;
908 }
909 else
910 return 0.0;
911 }
912
913 // normal cases; result x^(a-1)*(1-x)^(b-1)/Beta(a,b)
914 const double fLogDblMax = log( ::std::numeric_limits<double>::max());
915 const double fLogDblMin = log( ::std::numeric_limits<double>::min());
916 double fLogY = (fX < 0.1) ? std::log1p(-fX) : log(0.5-fX+0.5);
917 double fLogX = log(fX);
918 double fAm1LogX = (fA-1.0) * fLogX;
919 double fBm1LogY = (fB-1.0) * fLogY;
920 double fLogBeta = GetLogBeta(fA,fB);
921 // check whether parts over- or underflow
922 if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin
923 && fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin
924 && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin
925 && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > fLogDblMin)
926 return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);
927 else // need logarithm;
928 // might overflow as a whole, but seldom, not worth to pre-detect it
929 return exp( fAm1LogX + fBm1LogY - fLogBeta);
930}
931
932/*
933 x^a * (1-x)^b
934 I_x(a,b) = ---------------- * result of ContFrac
935 a * Beta(a,b)
936*/
937static double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)
938{ // like old version
939 double a1, b1, a2, b2, fnorm, cfnew, cf;
940 a1 = 1.0; b1 = 1.0;
941 b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;
942 if (b2 == 0.0)
943 {
944 a2 = 0.0;
945 fnorm = 1.0;
946 cf = 1.0;
947 }
948 else
949 {
950 a2 = 1.0;
951 fnorm = 1.0/b2;
952 cf = a2*fnorm;
953 }
954 cfnew = 1.0;
955 double rm = 1.0;
956
957 const double fMaxIter = 50000.0;
958 // loop security, normal cases converge in less than 100 iterations.
959 // FIXME: You will get so much iterations for fX near mean,
960 // I do not know a better algorithm.
961 bool bfinished = false;
962 do
963 {
964 const double apl2m = fA + 2.0*rm;
965 const double d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m);
966 const double d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0));
967 a1 = (a2+d2m*a1)*fnorm;
968 b1 = (b2+d2m*b1)*fnorm;
969 a2 = a1 + d2m1*a2*fnorm;
970 b2 = b1 + d2m1*b2*fnorm;
971 if (b2 != 0.0)
972 {
973 fnorm = 1.0/b2;
974 cfnew = a2*fnorm;
975 bfinished = (std::abs(cf-cfnew) < std::abs(cf)*fMachEps);
976 }
977 cf = cfnew;
978 rm += 1.0;
979 }
980 while (rm < fMaxIter && !bfinished);
981 return cf;
982}
983
984// cumulative distribution function, normalized
985double ScInterpreter::GetBetaDist(double fXin, double fAlpha, double fBeta)
986{
987// special cases
988 if (fXin <= 0.0) // values are valid, see spec
989 return 0.0;
990 if (fXin >= 1.0) // values are valid, see spec
991 return 1.0;
992 if (fBeta == 1.0)
993 return pow(fXin, fAlpha);
994 if (fAlpha == 1.0)
995 // 1.0 - pow(1.0-fX,fBeta) is not accurate enough
996 return -std::expm1(fBeta * std::log1p(-fXin));
997 //FIXME: need special algorithm for fX near fP for large fA,fB
998 double fResult;
999 // I use always continued fraction, power series are neither
1000 // faster nor more accurate.
1001 double fY = (0.5-fXin)+0.5;
1002 double flnY = std::log1p(-fXin);
1003 double fX = fXin;
1004 double flnX = log(fXin);
1005 double fA = fAlpha;
1006 double fB = fBeta;
1007 bool bReflect = fXin > fAlpha/(fAlpha+fBeta);
1008 if (bReflect)
1009 {
1010 fA = fBeta;
1011 fB = fAlpha;
1012 fX = fY;
1013 fY = fXin;
1014 flnX = flnY;
1015 flnY = log(fXin);
1016 }
1017 fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);
1018 fResult = fResult/fA;
1019 double fP = fA/(fA+fB);
1020 double fQ = fB/(fA+fB);
1021 double fTemp;
1022 if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental
1023 fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;
1024 else
1025 fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));
1026 fResult *= fTemp;
1027 if (bReflect)
1028 fResult = 0.5 - fResult + 0.5;
1029 if (fResult > 1.0) // ensure valid range
1030 fResult = 1.0;
1031 if (fResult < 0.0)
1032 fResult = 0.0;
1033 return fResult;
1034}
1035
1037{
1038 sal_uInt8 nParamCount = GetByte();
1039 if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547#
1040 return;
1041 double fLowerBound, fUpperBound;
1042 double alpha, beta, x;
1043 bool bIsCumulative;
1044 if (nParamCount == 6)
1045 bIsCumulative = GetBool();
1046 else
1047 bIsCumulative = true;
1048 if (nParamCount >= 5)
1049 fUpperBound = GetDouble();
1050 else
1051 fUpperBound = 1.0;
1052 if (nParamCount >= 4)
1053 fLowerBound = GetDouble();
1054 else
1055 fLowerBound = 0.0;
1056 beta = GetDouble();
1057 alpha = GetDouble();
1058 x = GetDouble();
1059 double fScale = fUpperBound - fLowerBound;
1060 if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0)
1061 {
1063 return;
1064 }
1065 if (bIsCumulative) // cumulative distribution function
1066 {
1067 // special cases
1068 if (x < fLowerBound)
1069 {
1070 PushDouble(0.0); return; //see spec
1071 }
1072 if (x > fUpperBound)
1073 {
1074 PushDouble(1.0); return; //see spec
1075 }
1076 // normal cases
1077 x = (x-fLowerBound)/fScale; // convert to standard form
1078 PushDouble(GetBetaDist(x, alpha, beta));
1079 return;
1080 }
1081 else // probability density function
1082 {
1083 if (x < fLowerBound || x > fUpperBound)
1084 {
1085 PushDouble(0.0);
1086 return;
1087 }
1088 x = (x-fLowerBound)/fScale;
1089 PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
1090 return;
1091 }
1092}
1093
1101{
1102 sal_uInt8 nParamCount = GetByte();
1103 if ( !MustHaveParamCount( nParamCount, 4, 6 ) )
1104 return;
1105 double fLowerBound, fUpperBound;
1106 double alpha, beta, x;
1107 bool bIsCumulative;
1108 if (nParamCount == 6)
1109 fUpperBound = GetDouble();
1110 else
1111 fUpperBound = 1.0;
1112 if (nParamCount >= 5)
1113 fLowerBound = GetDouble();
1114 else
1115 fLowerBound = 0.0;
1116 bIsCumulative = GetBool();
1117 beta = GetDouble();
1118 alpha = GetDouble();
1119 x = GetDouble();
1120 if (alpha <= 0.0 || beta <= 0.0 || x < fLowerBound || x > fUpperBound)
1121 {
1123 return;
1124 }
1125 double fScale = fUpperBound - fLowerBound;
1126 if (bIsCumulative) // cumulative distribution function
1127 {
1128 x = (x-fLowerBound)/fScale; // convert to standard form
1129 PushDouble(GetBetaDist(x, alpha, beta));
1130 return;
1131 }
1132 else // probability density function
1133 {
1134 x = (x-fLowerBound)/fScale;
1135 PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
1136 return;
1137 }
1138}
1139
1141{
1143}
1144
1146{
1148}
1149
1151{
1152 double fVal = GetDouble();
1153 if (std::abs(fVal) >= 1.0)
1155 else
1156 PushDouble(::atanh(fVal));
1157}
1158
1160{
1161 PushDouble( tanh( GetDouble()));
1162}
1163
1165{
1166 double nVal = GetDouble();
1167 if (nVal < 0.0)
1169 else
1170 PushDouble(Fakultaet(nVal));
1171}
1172
1174{
1175 if ( MustHaveParamCount( GetByte(), 2 ) )
1176 {
1177 double k = ::rtl::math::approxFloor(GetDouble());
1178 double n = ::rtl::math::approxFloor(GetDouble());
1179 if (k < 0.0 || n < 0.0 || k > n)
1181 else
1182 PushDouble(BinomKoeff(n, k));
1183 }
1184}
1185
1187{
1188 if ( MustHaveParamCount( GetByte(), 2 ) )
1189 {
1190 double k = ::rtl::math::approxFloor(GetDouble());
1191 double n = ::rtl::math::approxFloor(GetDouble());
1192 if (k < 0.0 || n < 0.0 || k > n)
1194 else
1195 PushDouble(BinomKoeff(n + k - 1, k));
1196 }
1197}
1198
1200{
1201 if ( !MustHaveParamCount( GetByte(), 2 ) )
1202 return;
1203
1204 double k = ::rtl::math::approxFloor(GetDouble());
1205 double n = ::rtl::math::approxFloor(GetDouble());
1206 if (n < 0.0 || k < 0.0 || k > n)
1208 else if (k == 0.0)
1209 PushInt(1); // (n! / (n - 0)!) == 1
1210 else
1211 {
1212 double nVal = n;
1213 for (sal_uLong i = static_cast<sal_uLong>(k)-1; i >= 1; i--)
1214 nVal *= n-static_cast<double>(i);
1215 PushDouble(nVal);
1216 }
1217}
1218
1220{
1221 if ( MustHaveParamCount( GetByte(), 2 ) )
1222 {
1223 double k = ::rtl::math::approxFloor(GetDouble());
1224 double n = ::rtl::math::approxFloor(GetDouble());
1225 if (n < 0.0 || k < 0.0)
1227 else
1228 PushDouble(pow(n,k));
1229 }
1230}
1231
1232double ScInterpreter::GetBinomDistPMF(double x, double n, double p)
1233// used in ScB and ScBinomDist
1234// preconditions: 0.0 <= x <= n, 0.0 < p < 1.0; x,n integral although double
1235{
1236 double q = (0.5 - p) + 0.5;
1237 double fFactor = pow(q, n);
1238 if (fFactor <=::std::numeric_limits<double>::min())
1239 {
1240 fFactor = pow(p, n);
1241 if (fFactor <= ::std::numeric_limits<double>::min())
1242 return GetBetaDistPDF(p, x+1.0, n-x+1.0)/(n+1.0);
1243 else
1244 {
1245 sal_uInt32 max = static_cast<sal_uInt32>(n - x);
1246 for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
1247 fFactor *= (n-i)/(i+1)*q/p;
1248 return fFactor;
1249 }
1250 }
1251 else
1252 {
1253 sal_uInt32 max = static_cast<sal_uInt32>(x);
1254 for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
1255 fFactor *= (n-i)/(i+1)*p/q;
1256 return fFactor;
1257 }
1258}
1259
1260static double lcl_GetBinomDistRange(double n, double xs,double xe,
1261 double fFactor /* q^n */, double p, double q)
1262//preconditions: 0.0 <= xs < xe <= n; xs,xe,n integral although double
1263{
1264 sal_uInt32 i;
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 KahanSum 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 std::min(fSum.get(), 1.0);
1277}
1278
1280{
1281 sal_uInt8 nParamCount = GetByte();
1282 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1283 return ;
1284 if (nParamCount == 3) // mass function
1285 {
1286 double x = ::rtl::math::approxFloor(GetDouble());
1287 double p = GetDouble();
1288 double n = ::rtl::math::approxFloor(GetDouble());
1289 if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1291 else if (p == 0.0)
1292 PushDouble( (x == 0.0) ? 1.0 : 0.0 );
1293 else if ( p == 1.0)
1294 PushDouble( (x == n) ? 1.0 : 0.0);
1295 else
1297 }
1298 else
1299 { // nParamCount == 4
1300 double xe = ::rtl::math::approxFloor(GetDouble());
1301 double xs = ::rtl::math::approxFloor(GetDouble());
1302 double p = GetDouble();
1303 double n = ::rtl::math::approxFloor(GetDouble());
1304 double q = (0.5 - p) + 0.5;
1305 bool bIsValidX = ( 0.0 <= xs && xs <= xe && xe <= n);
1306 if ( bIsValidX && 0.0 < p && p < 1.0)
1307 {
1308 if (xs == xe) // mass function
1310 else
1311 {
1312 double fFactor = pow(q, n);
1313 if (fFactor > ::std::numeric_limits<double>::min())
1314 PushDouble(lcl_GetBinomDistRange(n,xs,xe,fFactor,p,q));
1315 else
1316 {
1317 fFactor = pow(p, n);
1318 if (fFactor > ::std::numeric_limits<double>::min())
1319 {
1320 // sum from j=xs to xe {(n choose j) * p^j * q^(n-j)}
1321 // = sum from i = n-xe to n-xs { (n choose i) * q^i * p^(n-i)}
1322 PushDouble(lcl_GetBinomDistRange(n,n-xe,n-xs,fFactor,q,p));
1323 }
1324 else
1325 PushDouble(GetBetaDist(q,n-xe,xe+1.0)-GetBetaDist(q,n-xs+1,xs) );
1326 }
1327 }
1328 }
1329 else
1330 {
1331 if ( bIsValidX ) // not(0<p<1)
1332 {
1333 if ( p == 0.0 )
1334 PushDouble( (xs == 0.0) ? 1.0 : 0.0 );
1335 else if ( p == 1.0 )
1336 PushDouble( (xe == n) ? 1.0 : 0.0 );
1337 else
1339 }
1340 else
1342 }
1343 }
1344}
1345
1347{
1348 if ( !MustHaveParamCount( GetByte(), 4 ) )
1349 return;
1350
1351 bool bIsCum = GetBool(); // false=mass function; true=cumulative
1352 double p = GetDouble();
1353 double n = ::rtl::math::approxFloor(GetDouble());
1354 double x = ::rtl::math::approxFloor(GetDouble());
1355 double q = (0.5 - p) + 0.5; // get one bit more for p near 1.0
1356 if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1357 {
1359 return;
1360 }
1361 if ( p == 0.0)
1362 {
1363 PushDouble( (x==0.0 || bIsCum) ? 1.0 : 0.0 );
1364 return;
1365 }
1366 if ( p == 1.0)
1367 {
1368 PushDouble( (x==n) ? 1.0 : 0.0);
1369 return;
1370 }
1371 if (!bIsCum)
1373 else
1374 {
1375 if (x == n)
1376 PushDouble(1.0);
1377 else
1378 {
1379 double fFactor = pow(q, n);
1380 if (x == 0.0)
1381 PushDouble(fFactor);
1382 else if (fFactor <= ::std::numeric_limits<double>::min())
1383 {
1384 fFactor = pow(p, n);
1385 if (fFactor <= ::std::numeric_limits<double>::min())
1386 PushDouble(GetBetaDist(q,n-x,x+1.0));
1387 else
1388 {
1389 if (fFactor > fMachEps)
1390 {
1391 double fSum = 1.0 - fFactor;
1392 sal_uInt32 max = static_cast<sal_uInt32> (n - x) - 1;
1393 for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
1394 {
1395 fFactor *= (n-i)/(i+1)*q/p;
1396 fSum -= fFactor;
1397 }
1398 PushDouble( (fSum < 0.0) ? 0.0 : fSum );
1399 }
1400 else
1401 PushDouble(lcl_GetBinomDistRange(n,n-x,n,fFactor,q,p));
1402 }
1403 }
1404 else
1405 PushDouble( lcl_GetBinomDistRange(n,0.0,x,fFactor,p,q)) ;
1406 }
1407 }
1408}
1409
1411{
1412 if ( !MustHaveParamCount( GetByte(), 3 ) )
1413 return;
1414
1415 double alpha = GetDouble();
1416 double p = GetDouble();
1417 double n = ::rtl::math::approxFloor(GetDouble());
1418 if (n < 0.0 || alpha < 0.0 || alpha > 1.0 || p < 0.0 || p > 1.0)
1420 else if ( alpha == 0.0 )
1421 PushDouble( 0.0 );
1422 else if ( alpha == 1.0 )
1423 PushDouble( p == 0 ? 0.0 : n );
1424 else
1425 {
1426 double fFactor;
1427 double q = (0.5 - p) + 0.5; // get one bit more for p near 1.0
1428 if ( q > p ) // work from the side where the cumulative curve is
1429 {
1430 // work from 0 upwards
1431 fFactor = pow(q,n);
1432 if (fFactor > ::std::numeric_limits<double>::min())
1433 {
1434 KahanSum fSum = fFactor;
1435 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1436 for (i = 0; i < max && fSum < alpha; i++)
1437 {
1438 fFactor *= (n-i)/(i+1)*p/q;
1439 fSum += fFactor;
1440 }
1441 PushDouble(i);
1442 }
1443 else
1444 {
1445 // accumulate BinomDist until accumulated BinomDist reaches alpha
1446 KahanSum fSum = 0.0;
1447 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1448 for (i = 0; i < max && fSum < alpha; i++)
1449 {
1450 const double x = GetBetaDistPDF( p, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
1451 if ( nGlobalError == FormulaError::NONE )
1452 fSum += x;
1453 else
1454 {
1455 PushNoValue();
1456 return;
1457 }
1458 }
1459 PushDouble( i - 1 );
1460 }
1461 }
1462 else
1463 {
1464 // work from n backwards
1465 fFactor = pow(p, n);
1466 if (fFactor > ::std::numeric_limits<double>::min())
1467 {
1468 KahanSum fSum = 1.0 - fFactor;
1469 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1470 for (i = 0; i < max && fSum >= alpha; i++)
1471 {
1472 fFactor *= (n-i)/(i+1)*q/p;
1473 fSum -= fFactor;
1474 }
1475 PushDouble(n-i);
1476 }
1477 else
1478 {
1479 // accumulate BinomDist until accumulated BinomDist reaches alpha
1480 KahanSum fSum = 0.0;
1481 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
1482 alpha = 1 - alpha;
1483 for (i = 0; i < max && fSum < alpha; i++)
1484 {
1485 const double x = GetBetaDistPDF( q, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
1486 if ( nGlobalError == FormulaError::NONE )
1487 fSum += x;
1488 else
1489 {
1490 PushNoValue();
1491 return;
1492 }
1493 }
1494 PushDouble( n - i + 1 );
1495 }
1496 }
1497 }
1498}
1499
1501{
1502 if ( !MustHaveParamCount( GetByte(), 3 ) )
1503 return;
1504
1505 double p = GetDouble(); // probability
1506 double s = ::rtl::math::approxFloor(GetDouble()); // No of successes
1507 double f = ::rtl::math::approxFloor(GetDouble()); // No of failures
1508 if ((f + s) <= 1.0 || p < 0.0 || p > 1.0)
1510 else
1511 {
1512 double q = 1.0 - p;
1513 double fFactor = pow(p,s);
1514 for (double i = 0.0; i < f; i++)
1515 fFactor *= (i+s)/(i+1.0)*q;
1516 PushDouble(fFactor);
1517 }
1518}
1519
1521{
1522 if ( !MustHaveParamCount( GetByte(), 4 ) )
1523 return;
1524
1525 bool bCumulative = GetBool();
1526 double p = GetDouble(); // probability
1527 double s = ::rtl::math::approxFloor(GetDouble()); // No of successes
1528 double f = ::rtl::math::approxFloor(GetDouble()); // No of failures
1529 if ( s < 1.0 || f < 0.0 || p < 0.0 || p > 1.0 )
1531 else
1532 {
1533 double q = 1.0 - p;
1534 if ( bCumulative )
1535 PushDouble( 1.0 - GetBetaDist( q, f + 1, s ) );
1536 else
1537 {
1538 double fFactor = pow( p, s );
1539 for ( double i = 0.0; i < f; i++ )
1540 fFactor *= ( i + s ) / ( i + 1.0 ) * q;
1541 PushDouble( fFactor );
1542 }
1543 }
1544}
1545
1546void ScInterpreter::ScNormDist( int nMinParamCount )
1547{
1548 sal_uInt8 nParamCount = GetByte();
1549 if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
1550 return;
1551 bool bCumulative = nParamCount != 4 || GetBool();
1552 double sigma = GetDouble(); // standard deviation
1553 double mue = GetDouble(); // mean
1554 double x = GetDouble(); // x
1555 if (sigma <= 0.0)
1556 {
1558 return;
1559 }
1560 if (bCumulative)
1561 PushDouble(integralPhi((x-mue)/sigma));
1562 else
1563 PushDouble(phi((x-mue)/sigma)/sigma);
1564}
1565
1566void ScInterpreter::ScLogNormDist( int nMinParamCount ) //expanded, see #i100119# and fdo72158
1567{
1568 sal_uInt8 nParamCount = GetByte();
1569 if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
1570 return;
1571 bool bCumulative = nParamCount != 4 || GetBool(); // cumulative
1572 double sigma = nParamCount >= 3 ? GetDouble() : 1.0; // standard deviation
1573 double mue = nParamCount >= 2 ? GetDouble() : 0.0; // mean
1574 double x = GetDouble(); // x
1575 if (sigma <= 0.0)
1576 {
1578 return;
1579 }
1580 if (bCumulative)
1581 { // cumulative
1582 if (x <= 0.0)
1583 PushDouble(0.0);
1584 else
1585 PushDouble(integralPhi((log(x)-mue)/sigma));
1586 }
1587 else
1588 { // density
1589 if (x <= 0.0)
1591 else
1592 PushDouble(phi((log(x)-mue)/sigma)/sigma/x);
1593 }
1594}
1595
1597{
1599}
1600
1602{
1603 sal_uInt8 nParamCount = GetByte();
1604 if ( !MustHaveParamCount( nParamCount, 2 ) )
1605 return;
1606 bool bCumulative = GetBool(); // cumulative
1607 double x = GetDouble(); // x
1608
1609 if ( bCumulative )
1610 PushDouble( integralPhi( x ) );
1611 else
1612 PushDouble( exp( - pow( x, 2 ) / 2 ) / sqrt( 2 * M_PI ) );
1613}
1614
1616{
1617 if ( !MustHaveParamCount( GetByte(), 3 ) )
1618 return;
1619
1620 double kum = GetDouble(); // 0 or 1
1621 double lambda = GetDouble(); // lambda
1622 double x = GetDouble(); // x
1623 if (lambda <= 0.0)
1625 else if (kum == 0.0) // density
1626 {
1627 if (x >= 0.0)
1628 PushDouble(lambda * exp(-lambda*x));
1629 else
1630 PushInt(0);
1631 }
1632 else // distribution
1633 {
1634 if (x > 0.0)
1635 PushDouble(1.0 - exp(-lambda*x));
1636 else
1637 PushInt(0);
1638 }
1639}
1640
1642{
1643 if ( !MustHaveParamCount( GetByte(), 3 ) )
1644 return;
1645 double fFlag = ::rtl::math::approxFloor(GetDouble());
1646 double fDF = ::rtl::math::approxFloor(GetDouble());
1647 double T = GetDouble();
1648 if (fDF < 1.0 || T < 0.0 || (fFlag != 1.0 && fFlag != 2.0) )
1649 {
1651 return;
1652 }
1653 PushDouble( GetTDist( T, fDF, static_cast<int>(fFlag) ) );
1654}
1655
1657{
1658 if ( !MustHaveParamCount( GetByte(), 2 ) )
1659 return;
1660 double fDF = ::rtl::math::approxFloor( GetDouble() );
1661 double fT = GetDouble();
1662 if ( fDF < 1.0 || ( nTails == 2 && fT < 0.0 ) )
1663 {
1665 return;
1666 }
1667 double fRes = GetTDist( fT, fDF, nTails );
1668 if ( nTails == 1 && fT < 0.0 )
1669 PushDouble( 1.0 - fRes ); // tdf#105937, right tail, negative X
1670 else
1671 PushDouble( fRes );
1672}
1673
1675{
1676 if ( !MustHaveParamCount( GetByte(), 3 ) )
1677 return;
1678 bool bCumulative = GetBool();
1679 double fDF = ::rtl::math::approxFloor( GetDouble() );
1680 double T = GetDouble();
1681 if ( fDF < 1.0 )
1682 {
1684 return;
1685 }
1686 PushDouble( GetTDist( T, fDF, ( bCumulative ? 4 : 3 ) ) );
1687}
1688
1690{
1691 if ( !MustHaveParamCount( GetByte(), 3 ) )
1692 return;
1693 double fF2 = ::rtl::math::approxFloor(GetDouble());
1694 double fF1 = ::rtl::math::approxFloor(GetDouble());
1695 double fF = GetDouble();
1696 if (fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10)
1697 {
1699 return;
1700 }
1701 PushDouble(GetFDist(fF, fF1, fF2));
1702}
1703
1705{
1706 int nParamCount = GetByte();
1707 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1708 return;
1709 bool bCum;
1710 if ( nParamCount == 3 )
1711 bCum = true;
1712 else if ( IsMissing() )
1713 {
1714 bCum = true;
1715 Pop();
1716 }
1717 else
1718 bCum = GetBool();
1719 double fF2 = ::rtl::math::approxFloor( GetDouble() );
1720 double fF1 = ::rtl::math::approxFloor( GetDouble() );
1721 double fF = GetDouble();
1722 if ( fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 )
1723 {
1725 return;
1726 }
1727 if ( bCum )
1728 {
1729 // left tail cumulative distribution
1730 PushDouble( 1.0 - GetFDist( fF, fF1, fF2 ) );
1731 }
1732 else
1733 {
1734 // probability density function
1735 PushDouble( pow( fF1 / fF2, fF1 / 2 ) * pow( fF, ( fF1 / 2 ) - 1 ) /
1736 ( pow( ( 1 + ( fF * fF1 / fF2 ) ), ( fF1 + fF2 ) / 2 ) *
1737 GetBeta( fF1 / 2, fF2 / 2 ) ) );
1738 }
1739}
1740
1742{
1743 double fResult;
1744 if ( !MustHaveParamCount( GetByte(), 2 ) )
1745 return;
1746 double fDF = ::rtl::math::approxFloor(GetDouble());
1747 double fChi = GetDouble();
1748 if ( fDF < 1.0 // x<=0 returns 1, see ODFF1.2 6.18.11
1749 || ( !bODFF && fChi < 0 ) ) // Excel does not accept negative fChi
1750 {
1752 return;
1753 }
1754 fResult = GetChiDist( fChi, fDF);
1755 if (nGlobalError != FormulaError::NONE)
1756 {
1758 return;
1759 }
1760 PushDouble(fResult);
1761}
1762
1764{
1765 if ( !MustHaveParamCount( GetByte(), 4 ) )
1766 return;
1767
1768 double kum = GetDouble(); // 0 or 1
1769 double beta = GetDouble(); // beta
1770 double alpha = GetDouble(); // alpha
1771 double x = GetDouble(); // x
1772 if (alpha <= 0.0 || beta <= 0.0 || x < 0.0)
1774 else if (kum == 0.0) // Density
1775 PushDouble(alpha/pow(beta,alpha)*pow(x,alpha-1.0)*
1776 exp(-pow(x/beta,alpha)));
1777 else // Distribution
1778 PushDouble(1.0 - exp(-pow(x/beta,alpha)));
1779}
1780
1782{
1783 sal_uInt8 nParamCount = GetByte();
1784 if ( !MustHaveParamCount( nParamCount, ( bODFF ? 2 : 3 ), 3 ) )
1785 return;
1786
1787 bool bCumulative = nParamCount != 3 || GetBool(); // default cumulative
1788 double lambda = GetDouble(); // Mean
1789 double x = ::rtl::math::approxFloor(GetDouble()); // discrete distribution
1790 if (lambda <= 0.0 || x < 0.0)
1792 else if (!bCumulative) // Probability mass function
1793 {
1794 if (lambda >712.0) // underflow in exp(-lambda)
1795 { // accuracy 11 Digits
1796 PushDouble( exp(x*log(lambda)-lambda-GetLogGamma(x+1.0)));
1797 }
1798 else
1799 {
1800 double fPoissonVar = 1.0;
1801 for ( double f = 0.0; f < x; ++f )
1802 fPoissonVar *= lambda / ( f + 1.0 );
1803 PushDouble( fPoissonVar * exp( -lambda ) );
1804 }
1805 }
1806 else // Cumulative distribution function
1807 {
1808 if (lambda > 712.0) // underflow in exp(-lambda)
1809 { // accuracy 12 Digits
1810 PushDouble(GetUpRegIGamma(x+1.0,lambda));
1811 }
1812 else
1813 {
1814 if (x >= 936.0) // result is always indistinguishable from 1
1815 PushDouble (1.0);
1816 else
1817 {
1818 double fSummand = std::exp(-lambda);
1819 KahanSum fSum = fSummand;
1820 int nEnd = sal::static_int_cast<int>( x );
1821 for (int i = 1; i <= nEnd; i++)
1822 {
1823 fSummand = (fSummand * lambda)/static_cast<double>(i);
1824 fSum += fSummand;
1825 }
1826 PushDouble(fSum.get());
1827 }
1828 }
1829 }
1830}
1831
1834static void lcl_PutFactorialElements( ::std::vector< double >& cn, double fLower, double fUpper, double fBase )
1835{
1836 for ( double i = fLower; i <= fUpper; ++i )
1837 {
1838 double fVal = fBase - i;
1839 if ( fVal > 1.0 )
1840 cn.push_back( fVal );
1841 }
1842}
1843
1856void ScInterpreter::ScHypGeomDist( int nMinParamCount )
1857{
1858 sal_uInt8 nParamCount = GetByte();
1859 if ( !MustHaveParamCount( nParamCount, nMinParamCount, 5 ) )
1860 return;
1861
1862 bool bCumulative = ( nParamCount == 5 && GetBool() );
1863 double N = ::rtl::math::approxFloor(GetDouble());
1864 double M = ::rtl::math::approxFloor(GetDouble());
1865 double n = ::rtl::math::approxFloor(GetDouble());
1866 double x = ::rtl::math::approxFloor(GetDouble());
1867
1868 if ( (x < 0.0) || (n < x) || (N < n) || (N < M) || (M < 0.0) )
1869 {
1871 return;
1872 }
1873
1874 KahanSum fVal = 0.0;
1875
1876 for ( int i = ( bCumulative ? 0 : x ); i <= x && nGlobalError == FormulaError::NONE; i++ )
1877 {
1878 if ( (n - i <= N - M) && (i <= M) )
1879 fVal += GetHypGeomDist( i, n, M, N );
1880 }
1881
1882 PushDouble( fVal.get() );
1883}
1884
1895double ScInterpreter::GetHypGeomDist( double x, double n, double M, double N )
1896{
1897 const size_t nMaxArraySize = 500000; // arbitrary max array size
1898
1899 std::vector<double> cnNumer, cnDenom;
1900
1901 size_t nEstContainerSize = static_cast<size_t>( x + ::std::min( n, M ) );
1902 size_t nMaxSize = ::std::min( cnNumer.max_size(), nMaxArraySize );
1903 if ( nEstContainerSize > nMaxSize )
1904 {
1905 PushNoValue();
1906 return 0;
1907 }
1908 cnNumer.reserve( nEstContainerSize + 10 );
1909 cnDenom.reserve( nEstContainerSize + 10 );
1910
1911 // Trim coefficient C first
1912 double fCNumVarUpper = N - n - M + x - 1.0;
1913 double fCDenomVarLower = 1.0;
1914 if ( N - n - M + x >= M - x + 1.0 )
1915 {
1916 fCNumVarUpper = M - x - 1.0;
1917 fCDenomVarLower = N - n - 2.0*(M - x) + 1.0;
1918 }
1919
1920 double fCNumLower = N - n - fCNumVarUpper;
1921 double fCDenomUpper = N - n - M + x + 1.0 - fCDenomVarLower;
1922
1923 double fDNumVarLower = n - M;
1924
1925 if ( n >= M + 1.0 )
1926 {
1927 if ( N - M < n + 1.0 )
1928 {
1929 // Case 1
1930
1931 if ( N - n < n + 1.0 )
1932 {
1933 // no overlap
1934 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1935 lcl_PutFactorialElements( cnDenom, 0.0, N - n - 1.0, N );
1936 }
1937 else
1938 {
1939 // overlap
1940 OSL_ENSURE( fCNumLower < n + 1.0, "ScHypGeomDist: wrong assertion" );
1941 lcl_PutFactorialElements( cnNumer, N - 2.0*n, fCNumVarUpper, N - n );
1942 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1943 }
1944
1945 OSL_ENSURE( fCDenomUpper <= N - M, "ScHypGeomDist: wrong assertion" );
1946
1947 if ( fCDenomUpper < n - x + 1.0 )
1948 // no overlap
1949 lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
1950 else
1951 {
1952 // overlap
1953 lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1954
1955 fCDenomUpper = n - x;
1956 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1957 }
1958 }
1959 else
1960 {
1961 // Case 2
1962
1963 if ( n > M - 1.0 )
1964 {
1965 // no overlap
1966 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1967 lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
1968 }
1969 else
1970 {
1971 lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
1972 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1973 }
1974
1975 OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
1976
1977 if ( fCDenomUpper < n - x + 1.0 )
1978 // no overlap
1979 lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - n + x, N - M + 1.0 );
1980 else
1981 {
1982 lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1983 fCDenomUpper = n - x;
1984 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1985 }
1986 }
1987
1988 OSL_ENSURE( fCDenomUpper <= M, "ScHypGeomDist: wrong assertion" );
1989 }
1990 else
1991 {
1992 if ( N - M < M + 1.0 )
1993 {
1994 // Case 3
1995
1996 if ( N - n < M + 1.0 )
1997 {
1998 // No overlap
1999 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
2000 lcl_PutFactorialElements( cnDenom, 0.0, N - M - 1.0, N );
2001 }
2002 else
2003 {
2004 lcl_PutFactorialElements( cnNumer, N - n - M, fCNumVarUpper, N - n );
2005 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
2006 }
2007
2008 if ( n - x + 1.0 > fCDenomUpper )
2009 // No overlap
2010 lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
2011 else
2012 {
2013 // Overlap
2014 lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
2015
2016 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
2017 fCDenomUpper = n - x;
2018 }
2019 }
2020 else
2021 {
2022 // Case 4
2023
2024 OSL_ENSURE( M >= n - x, "ScHypGeomDist: wrong assertion" );
2025 OSL_ENSURE( M - x <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
2026
2027 if ( N - n < N - M + 1.0 )
2028 {
2029 // No overlap
2030 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
2031 lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
2032 }
2033 else
2034 {
2035 // Overlap
2036 OSL_ENSURE( fCNumLower <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
2037 lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
2038 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
2039 }
2040
2041 if ( n - x + 1.0 > fCDenomUpper )
2042 // No overlap
2043 lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - n + x, N - M + 1.0 );
2044 else if ( M >= fCDenomUpper )
2045 {
2046 lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
2047
2048 fCDenomUpper = n - x;
2049 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
2050 }
2051 else
2052 {
2053 OSL_ENSURE( M <= fCDenomUpper, "ScHypGeomDist: wrong assertion" );
2054 lcl_PutFactorialElements( cnDenom, fCDenomVarLower, N - n - 2.0*M + x,
2055 N - n - M + x + 1.0 );
2056
2057 fCDenomUpper = n - x;
2058 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
2059 }
2060 }
2061
2062 OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
2063
2064 fDNumVarLower = 0.0;
2065 }
2066
2067 double nDNumVarUpper = fCDenomUpper < x + 1.0 ? n - x - 1.0 : n - fCDenomUpper - 1.0;
2068 double nDDenomVarLower = fCDenomUpper < x + 1.0 ? fCDenomVarLower : N - n - M + 1.0;
2069 lcl_PutFactorialElements( cnNumer, fDNumVarLower, nDNumVarUpper, n );
2070 lcl_PutFactorialElements( cnDenom, nDDenomVarLower, N - n - M + x, N - n - M + x + 1.0 );
2071
2072 ::std::sort( cnNumer.begin(), cnNumer.end() );
2073 ::std::sort( cnDenom.begin(), cnDenom.end() );
2074 auto it1 = cnNumer.rbegin(), it1End = cnNumer.rend();
2075 auto it2 = cnDenom.rbegin(), it2End = cnDenom.rend();
2076
2077 double fFactor = 1.0;
2078 for ( ; it1 != it1End || it2 != it2End; )
2079 {
2080 double fEnum = 1.0, fDenom = 1.0;
2081 if ( it1 != it1End )
2082 fEnum = *it1++;
2083 if ( it2 != it2End )
2084 fDenom = *it2++;
2085 fFactor *= fEnum / fDenom;
2086 }
2087
2088 return fFactor;
2089}
2090
2092{
2093 sal_uInt8 nMinParamCount = ( bODFF ? 3 : 4 );
2094 sal_uInt8 nParamCount = GetByte();
2095 if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
2096 return;
2097 bool bCumulative;
2098 if (nParamCount == 4)
2099 bCumulative = GetBool();
2100 else
2101 bCumulative = true;
2102 double fBeta = GetDouble(); // scale
2103 double fAlpha = GetDouble(); // shape
2104 double fX = GetDouble(); // x
2105 if ((!bODFF && fX < 0) || fAlpha <= 0.0 || fBeta <= 0.0)
2107 else
2108 {
2109 if (bCumulative) // distribution
2110 PushDouble( GetGammaDist( fX, fAlpha, fBeta));
2111 else // density
2112 PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta));
2113 }
2114}
2115
2117{
2118 if ( MustHaveParamCount( GetByte(), 3 ) )
2119 {
2120 double sigma = GetDouble();
2121 double mue = GetDouble();
2122 double x = GetDouble();
2123 if (sigma <= 0.0 || x < 0.0 || x > 1.0)
2125 else if (x == 0.0 || x == 1.0)
2126 PushNoValue();
2127 else
2128 PushDouble(gaussinv(x)*sigma + mue);
2129 }
2130}
2131
2133{
2134 double x = GetDouble();
2135 if (x < 0.0 || x > 1.0)
2137 else if (x == 0.0 || x == 1.0)
2138 PushNoValue();
2139 else
2141}
2142
2144{
2145 sal_uInt8 nParamCount = GetByte();
2146 if ( MustHaveParamCount( nParamCount, 1, 3 ) )
2147 {
2148 double fSigma = ( nParamCount == 3 ? GetDouble() : 1.0 ); // Stddev
2149 double fMue = ( nParamCount >= 2 ? GetDouble() : 0.0 ); // Mean
2150 double fP = GetDouble(); // p
2151 if ( fSigma <= 0.0 || fP <= 0.0 || fP >= 1.0 )
2153 else
2154 PushDouble( exp( fMue + fSigma * gaussinv( fP ) ) );
2155 }
2156}
2157
2158class ScGammaDistFunction : public ScDistFunc
2159{
2161 double fp, fAlpha, fBeta;
2162
2163public:
2164 ScGammaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2165 rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2166
2168
2169 double GetValue( double x ) const override { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); }
2170};
2171
2173{
2174 if ( !MustHaveParamCount( GetByte(), 3 ) )
2175 return;
2176 double fBeta = GetDouble();
2177 double fAlpha = GetDouble();
2178 double fP = GetDouble();
2179 if (fAlpha <= 0.0 || fBeta <= 0.0 || fP < 0.0 || fP >= 1.0 )
2180 {
2182 return;
2183 }
2184 if (fP == 0.0)
2185 PushInt(0);
2186 else
2187 {
2188 bool bConvError;
2189 ScGammaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2190 double fStart = fAlpha * fBeta;
2191 double fVal = lcl_IterateInverse( aFunc, fStart*0.5, fStart, bConvError );
2192 if (bConvError)
2193 SetError(FormulaError::NoConvergence);
2194 PushDouble(fVal);
2195 }
2196}
2197
2198class ScBetaDistFunction : public ScDistFunc
2199{
2201 double fp, fAlpha, fBeta;
2202
2203public:
2204 ScBetaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2205 rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2206
2208
2209 double GetValue( double x ) const override { return fp - rInt.GetBetaDist(x, fAlpha, fBeta); }
2210};
2211
2213{
2214 sal_uInt8 nParamCount = GetByte();
2215 if ( !MustHaveParamCount( nParamCount, 3, 5 ) )
2216 return;
2217 double fP, fA, fB, fAlpha, fBeta;
2218 if (nParamCount == 5)
2219 fB = GetDouble();
2220 else
2221 fB = 1.0;
2222 if (nParamCount >= 4)
2223 fA = GetDouble();
2224 else
2225 fA = 0.0;
2226 fBeta = GetDouble();
2227 fAlpha = GetDouble();
2228 fP = GetDouble();
2229 if (fP < 0.0 || fP > 1.0 || fA >= fB || fAlpha <= 0.0 || fBeta <= 0.0)
2230 {
2232 return;
2233 }
2234
2235 bool bConvError;
2236 ScBetaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2237 // 0..1 as range for iteration so it isn't extended beyond the valid range
2238 double fVal = lcl_IterateInverse( aFunc, 0.0, 1.0, bConvError );
2239 if (bConvError)
2240 PushError( FormulaError::NoConvergence);
2241 else
2242 PushDouble(fA + fVal*(fB-fA)); // scale to (A,B)
2243}
2244
2245// Note: T, F, and Chi are
2246// monotonically decreasing,
2247// therefore 1-Dist as function
2248
2249class ScTDistFunction : public ScDistFunc
2250{
2252 double fp, fDF;
2253 int nT;
2254
2255public:
2256 ScTDistFunction( ScInterpreter& rI, double fpVal, double fDFVal, int nType ) :
2257 rInt( rI ), fp( fpVal ), fDF( fDFVal ), nT( nType ) {}
2258
2259 virtual ~ScTDistFunction() {}
2260
2261 double GetValue( double x ) const override { return fp - rInt.GetTDist( x, fDF, nT ); }
2262};
2263
2265{
2266 if ( !MustHaveParamCount( GetByte(), 2 ) )
2267 return;
2268 double fDF = ::rtl::math::approxFloor(GetDouble());
2269 double fP = GetDouble();
2270 if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
2271 {
2273 return;
2274 }
2275 if ( nType == 4 ) // left-tailed cumulative t-distribution
2276 {
2277 if ( fP == 1.0 )
2279 else if ( fP < 0.5 )
2280 PushDouble( -GetTInv( 1 - fP, fDF, nType ) );
2281 else
2282 PushDouble( GetTInv( fP, fDF, nType ) );
2283 }
2284 else
2285 PushDouble( GetTInv( fP, fDF, nType ) );
2286};
2287
2288double ScInterpreter::GetTInv( double fAlpha, double fSize, int nType )
2289{
2290 bool bConvError;
2291 ScTDistFunction aFunc( *this, fAlpha, fSize, nType );
2292 double fVal = lcl_IterateInverse( aFunc, fSize * 0.5, fSize, bConvError );
2293 if (bConvError)
2294 SetError(FormulaError::NoConvergence);
2295 return fVal;
2296}
2297
2298class ScFDistFunction : public ScDistFunc
2299{
2301 double fp, fF1, fF2;
2302
2303public:
2304 ScFDistFunction( ScInterpreter& rI, double fpVal, double fF1Val, double fF2Val ) :
2305 rInt(rI), fp(fpVal), fF1(fF1Val), fF2(fF2Val) {}
2306
2307 virtual ~ScFDistFunction() {}
2308
2309 double GetValue( double x ) const override { return fp - rInt.GetFDist(x, fF1, fF2); }
2310};
2311
2313{
2314 if ( !MustHaveParamCount( GetByte(), 3 ) )
2315 return;
2316 double fF2 = ::rtl::math::approxFloor(GetDouble());
2317 double fF1 = ::rtl::math::approxFloor(GetDouble());
2318 double fP = GetDouble();
2319 if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
2320 {
2322 return;
2323 }
2324
2325 bool bConvError;
2326 ScFDistFunction aFunc( *this, fP, fF1, fF2 );
2327 double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
2328 if (bConvError)
2329 SetError(FormulaError::NoConvergence);
2330 PushDouble(fVal);
2331}
2332
2334{
2335 if ( !MustHaveParamCount( GetByte(), 3 ) )
2336 return;
2337 double fF2 = ::rtl::math::approxFloor(GetDouble());
2338 double fF1 = ::rtl::math::approxFloor(GetDouble());
2339 double fP = GetDouble();
2340 if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
2341 {
2343 return;
2344 }
2345
2346 bool bConvError;
2347 ScFDistFunction aFunc( *this, ( 1.0 - fP ), fF1, fF2 );
2348 double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
2349 if (bConvError)
2350 SetError(FormulaError::NoConvergence);
2351 PushDouble(fVal);
2352}
2353
2354class ScChiDistFunction : public ScDistFunc
2355{
2357 double fp, fDF;
2358
2359public:
2360 ScChiDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2361 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2362
2364
2365 double GetValue( double x ) const override { return fp - rInt.GetChiDist(x, fDF); }
2366};
2367
2369{
2370 if ( !MustHaveParamCount( GetByte(), 2 ) )
2371 return;
2372 double fDF = ::rtl::math::approxFloor(GetDouble());
2373 double fP = GetDouble();
2374 if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
2375 {
2377 return;
2378 }
2379
2380 bool bConvError;
2381 ScChiDistFunction aFunc( *this, fP, fDF );
2382 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2383 if (bConvError)
2384 SetError(FormulaError::NoConvergence);
2385 PushDouble(fVal);
2386}
2387
2388/***********************************************/
2389class ScChiSqDistFunction : public ScDistFunc
2390{
2392 double fp, fDF;
2393
2394public:
2395 ScChiSqDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2396 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2397
2399
2400 double GetValue( double x ) const override { return fp - rInt.GetChiSqDistCDF(x, fDF); }
2401};
2402
2404{
2405 if ( !MustHaveParamCount( GetByte(), 2 ) )
2406 return;
2407 double fDF = ::rtl::math::approxFloor(GetDouble());
2408 double fP = GetDouble();
2409 if (fDF < 1.0 || fP < 0.0 || fP >= 1.0 )
2410 {
2412 return;
2413 }
2414
2415 bool bConvError;
2416 ScChiSqDistFunction aFunc( *this, fP, fDF );
2417 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2418 if (bConvError)
2419 SetError(FormulaError::NoConvergence);
2420 PushDouble(fVal);
2421}
2422
2424{
2425 if ( MustHaveParamCount( GetByte(), 3 ) )
2426 {
2427 double n = ::rtl::math::approxFloor(GetDouble());
2428 double sigma = GetDouble();
2429 double alpha = GetDouble();
2430 if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
2432 else
2433 PushDouble( gaussinv(1.0-alpha/2.0) * sigma/sqrt(n) );
2434 }
2435}
2436
2438{
2439 if ( MustHaveParamCount( GetByte(), 3 ) )
2440 {
2441 double n = ::rtl::math::approxFloor(GetDouble());
2442 double sigma = GetDouble();
2443 double alpha = GetDouble();
2444 if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
2446 else if (n == 1.0) // for interoperability with Excel
2447 PushError(FormulaError::DivisionByZero);
2448 else
2449 PushDouble( sigma * GetTInv( alpha, n - 1, 2 ) / sqrt( n ) );
2450 }
2451}
2452
2454{
2455 sal_uInt8 nParamCount = GetByte();
2456 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
2457 return;
2458 double sigma = 0.0, x;
2459 if (nParamCount == 3)
2460 {
2461 sigma = GetDouble();
2462 if (sigma <= 0.0)
2463 {
2465 return;
2466 }
2467 }
2468 x = GetDouble();
2469
2470 KahanSum fSum = 0.0;
2471 KahanSum fSumSqr = 0.0;
2472 double fVal;
2473 double rValCount = 0.0;
2474 switch (GetStackType())
2475 {
2476 case svDouble :
2477 {
2478 fVal = GetDouble();
2479 fSum += fVal;
2480 fSumSqr += fVal*fVal;
2481 rValCount++;
2482 }
2483 break;
2484 case svSingleRef :
2485 {
2486 ScAddress aAdr;
2487 PopSingleRef( aAdr );
2488 ScRefCellValue aCell(mrDoc, aAdr);
2489 if (aCell.hasNumeric())
2490 {
2491 fVal = GetCellValue(aAdr, aCell);
2492 fSum += fVal;
2493 fSumSqr += fVal*fVal;
2494 rValCount++;
2495 }
2496 }
2497 break;
2498 case svRefList :
2499 case svDoubleRef :
2500 {
2501 short nParam = 1;
2502 size_t nRefInList = 0;
2503 while (nParam-- > 0)
2504 {
2505 ScRange aRange;
2506 FormulaError nErr = FormulaError::NONE;
2507 PopDoubleRef( aRange, nParam, nRefInList);
2508 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
2509 if (aValIter.GetFirst(fVal, nErr))
2510 {
2511 fSum += fVal;
2512 fSumSqr += fVal*fVal;
2513 rValCount++;
2514 while ((nErr == FormulaError::NONE) && aValIter.GetNext(fVal, nErr))
2515 {
2516 fSum += fVal;
2517 fSumSqr += fVal*fVal;
2518 rValCount++;
2519 }
2520 SetError(nErr);
2521 }
2522 }
2523 }
2524 break;
2525 case svMatrix :
2528 {
2529 ScMatrixRef pMat = GetMatrix();
2530 if (pMat)
2531 {
2532 SCSIZE nCount = pMat->GetElementCount();
2533 if (pMat->IsNumeric())
2534 {
2535 for ( SCSIZE i = 0; i < nCount; i++ )
2536 {
2537 fVal= pMat->GetDouble(i);
2538 fSum += fVal;
2539 fSumSqr += fVal * fVal;
2540 rValCount++;
2541 }
2542 }
2543 else
2544 {
2545 for (SCSIZE i = 0; i < nCount; i++)
2546 if (!pMat->IsStringOrEmpty(i))
2547 {
2548 fVal= pMat->GetDouble(i);
2549 fSum += fVal;
2550 fSumSqr += fVal * fVal;
2551 rValCount++;
2552 }
2553 }
2554 }
2555 }
2556 break;
2557 default : SetError(FormulaError::IllegalParameter); break;
2558 }
2559 if (rValCount <= 1.0)
2560 PushError( FormulaError::DivisionByZero);
2561 else
2562 {
2563 double mue = fSum.get()/rValCount;
2564
2565 if (nParamCount != 3)
2566 {
2567 sigma = (fSumSqr - fSum*fSum/rValCount).get()/(rValCount-1.0);
2568 if (sigma == 0.0)
2569 {
2570 PushError(FormulaError::DivisionByZero);
2571 return;
2572 }
2573 PushDouble(0.5 - gauss((mue-x)/sqrt(sigma/rValCount)));
2574 }
2575 else
2576 PushDouble(0.5 - gauss((mue-x)*sqrt(rValCount)/sigma));
2577 }
2578}
2579
2581 ,const SCSIZE nC1, const SCSIZE nC2,const SCSIZE nR1,const SCSIZE nR2
2582 ,const ScMatrixRef& pMat1,const ScMatrixRef& pMat2
2583 ,double& fT,double& fF)
2584{
2585 double fCount1 = 0.0;
2586 double fCount2 = 0.0;
2587 KahanSum fSum1 = 0.0;
2588 KahanSum fSumSqr1 = 0.0;
2589 KahanSum fSum2 = 0.0;
2590 KahanSum fSumSqr2 = 0.0;
2591 double fVal;
2592 SCSIZE i,j;
2593 for (i = 0; i < nC1; i++)
2594 for (j = 0; j < nR1; j++)
2595 {
2596 if (!pMat1->IsStringOrEmpty(i,j))
2597 {
2598 fVal = pMat1->GetDouble(i,j);
2599 fSum1 += fVal;
2600 fSumSqr1 += fVal * fVal;
2601 fCount1++;
2602 }
2603 }
2604 for (i = 0; i < nC2; i++)
2605 for (j = 0; j < nR2; j++)
2606 {
2607 if (!pMat2->IsStringOrEmpty(i,j))
2608 {
2609 fVal = pMat2->GetDouble(i,j);
2610 fSum2 += fVal;
2611 fSumSqr2 += fVal * fVal;
2612 fCount2++;
2613 }
2614 }
2615 if (fCount1 < 2.0 || fCount2 < 2.0)
2616 {
2617 PushNoValue();
2618 return false;
2619 } // if (fCount1 < 2.0 || fCount2 < 2.0)
2620 if ( _bTemplin )
2621 {
2622 double fS1 = (fSumSqr1-fSum1*fSum1/fCount1).get() / (fCount1-1.0) / fCount1;
2623 double fS2 = (fSumSqr2-fSum2*fSum2/fCount2).get() / (fCount2-1.0) / fCount2;
2624 if (fS1 + fS2 == 0.0)
2625 {
2626 PushNoValue();
2627 return false;
2628 }
2629 fT = std::abs(( fSum1/fCount1 - fSum2/fCount2 ).get())/sqrt(fS1+fS2);
2630 double c = fS1/(fS1+fS2);
2631 // GetTDist is calculated via GetBetaDist and also works with non-integral
2632 // degrees of freedom. The result matches Excel
2633 fF = 1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0));
2634 }
2635 else
2636 {
2637 // according to Bronstein-Semendjajew
2638 double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1).get() / (fCount1 - 1.0); // Variance
2639 double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2).get() / (fCount2 - 1.0);
2640 fT = std::abs( fSum1.get()/fCount1 - fSum2.get()/fCount2 ) /
2641 sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
2642 sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
2643 fF = fCount1 + fCount2 - 2;
2644 }
2645 return true;
2646}
2648{
2649 if ( !MustHaveParamCount( GetByte(), 4 ) )
2650 return;
2651 double fTyp = ::rtl::math::approxFloor(GetDouble());
2652 double fTails = ::rtl::math::approxFloor(GetDouble());
2653 if (fTails != 1.0 && fTails != 2.0)
2654 {
2656 return;
2657 }
2658
2659 ScMatrixRef pMat2 = GetMatrix();
2660 ScMatrixRef pMat1 = GetMatrix();
2661 if (!pMat1 || !pMat2)
2662 {
2664 return;
2665 }
2666 double fT, fF;
2667 SCSIZE nC1, nC2;
2668 SCSIZE nR1, nR2;
2669 SCSIZE i, j;
2670 pMat1->GetDimensions(nC1, nR1);
2671 pMat2->GetDimensions(nC2, nR2);
2672 if (fTyp == 1.0)
2673 {
2674 if (nC1 != nC2 || nR1 != nR2)
2675 {
2677 return;
2678 }
2679 double fCount = 0.0;
2680 KahanSum fSum1 = 0.0;
2681 KahanSum fSum2 = 0.0;
2682 KahanSum fSumSqrD = 0.0;
2683 double fVal1, fVal2;
2684 for (i = 0; i < nC1; i++)
2685 for (j = 0; j < nR1; j++)
2686 {
2687 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
2688 {
2689 fVal1 = pMat1->GetDouble(i,j);
2690 fVal2 = pMat2->GetDouble(i,j);
2691 fSum1 += fVal1;
2692 fSum2 += fVal2;
2693 fSumSqrD += (fVal1 - fVal2)*(fVal1 - fVal2);
2694 fCount++;
2695 }
2696 }
2697 if (fCount < 1.0)
2698 {
2699 PushNoValue();
2700 return;
2701 }
2702 KahanSum fSumD = fSum1 - fSum2;
2703 double fDivider = ( fSumSqrD*fCount - fSumD*fSumD ).get();
2704 if ( fDivider == 0.0 )
2705 {
2706 PushError(FormulaError::DivisionByZero);
2707 return;
2708 }
2709 fT = std::abs(fSumD.get()) * sqrt((fCount-1.0) / fDivider);
2710 fF = fCount - 1.0;
2711 }
2712 else if (fTyp == 2.0)
2713 {
2714 if (!CalculateTest(false,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF))
2715 return; // error was pushed
2716 }
2717 else if (fTyp == 3.0)
2718 {
2719 if (!CalculateTest(true,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF))
2720 return; // error was pushed
2721 }
2722 else
2723 {
2725 return;
2726 }
2727 PushDouble( GetTDist( fT, fF, static_cast<int>(fTails) ) );
2728}
2729
2731{
2732 if ( !MustHaveParamCount( GetByte(), 2 ) )
2733 return;
2734 ScMatrixRef pMat2 = GetMatrix();
2735 ScMatrixRef pMat1 = GetMatrix();
2736 if (!pMat1 || !pMat2)
2737 {
2739 return;
2740 }
2741
2742 auto aVal1 = pMat1->CollectKahan(sc::op::kOpSumAndSumSquare);
2743 auto aVal2 = pMat2->CollectKahan(sc::op::kOpSumAndSumSquare);
2744 double fCount1 = aVal1.mnCount;
2745 double fCount2 = aVal2.mnCount;
2746 KahanSum fSum1 = aVal1.maAccumulator[0];
2747 KahanSum fSumSqr1 = aVal1.maAccumulator[1];
2748 KahanSum fSum2 = aVal2.maAccumulator[0];
2749 KahanSum fSumSqr2 = aVal2.maAccumulator[1];
2750
2751 if (fCount1 < 2.0 || fCount2 < 2.0)
2752 {
2753 PushNoValue();
2754 return;
2755 }
2756 double fS1 = (fSumSqr1-fSum1*fSum1/fCount1).get() / (fCount1-1.0);
2757 double fS2 = (fSumSqr2-fSum2*fSum2/fCount2).get() / (fCount2-1.0);
2758 if (fS1 == 0.0 || fS2 == 0.0)
2759 {
2760 PushNoValue();
2761 return;
2762 }
2763 double fF, fF1, fF2;
2764 if (fS1 > fS2)
2765 {
2766 fF = fS1/fS2;
2767 fF1 = fCount1-1.0;
2768 fF2 = fCount2-1.0;
2769 }
2770 else
2771 {
2772 fF = fS2/fS1;
2773 fF1 = fCount2-1.0;
2774 fF2 = fCount1-1.0;
2775 }
2776 double fFcdf = GetFDist(fF, fF1, fF2);
2777 PushDouble(2.0*std::min(fFcdf, 1.0 - fFcdf));
2778}
2779
2781{
2782 if ( !MustHaveParamCount( GetByte(), 2 ) )
2783 return;
2784 ScMatrixRef pMat2 = GetMatrix();
2785 ScMatrixRef pMat1 = GetMatrix();
2786 if (!pMat1 || !pMat2)
2787 {
2789 return;
2790 }
2791 SCSIZE nC1, nC2;
2792 SCSIZE nR1, nR2;
2793 pMat1->GetDimensions(nC1, nR1);
2794 pMat2->GetDimensions(nC2, nR2);
2795 if (nR1 != nR2 || nC1 != nC2)
2796 {
2798 return;
2799 }
2800 KahanSum fChi = 0.0;
2801 bool bEmpty = true;
2802 for (SCSIZE i = 0; i < nC1; i++)
2803 {
2804 for (SCSIZE j = 0; j < nR1; j++)
2805 {
2806 if (!(pMat1->IsEmpty(i,j) || pMat2->IsEmpty(i,j)))
2807 {
2808 bEmpty = false;
2809 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
2810 {
2811 double fValX = pMat1->GetDouble(i,j);
2812 double fValE = pMat2->GetDouble(i,j);
2813 if ( fValE == 0.0 )
2814 {
2815 PushError(FormulaError::DivisionByZero);
2816 return;
2817 }
2818 // These fTemp values guard against a failure when compiled
2819 // with optimization (using g++ 4.8.2 on tinderbox 71-TDF),
2820 // where ((fValX - fValE) * (fValX - fValE)) with
2821 // fValE==1e+308 should had produced Infinity but did
2822 // not, instead the result of divide() then was 1e+308.
2823 volatile double fTemp1 = (fValX - fValE) * (fValX - fValE);
2824 double fTemp2 = fTemp1;
2825 if (std::isinf(fTemp2))
2826 {
2827 PushError(FormulaError::NoConvergence);
2828 return;
2829 }
2830 fChi += sc::divide( fTemp2, fValE);
2831 }
2832 else
2833 {
2835 return;
2836 }
2837 }
2838 }
2839 }
2840 if ( bEmpty )
2841 {
2842 // not in ODFF1.2, but for interoperability with Excel
2844 return;
2845 }
2846 double fDF;
2847 if (nC1 == 1 || nR1 == 1)
2848 {
2849 fDF = static_cast<double>(nC1*nR1 - 1);
2850 if (fDF == 0.0)
2851 {
2852 PushNoValue();
2853 return;
2854 }
2855 }
2856 else
2857 fDF = static_cast<double>(nC1-1)*static_cast<double>(nR1-1);
2858 PushDouble(GetChiDist(fChi.get(), fDF));
2859}
2860
2862{
2863 KahanSum fSum;
2864 double fCount;
2865 std::vector<double> values;
2866 if ( !CalculateSkew(fSum, fCount, values) )
2867 return;
2868
2869 // ODF 1.2 constraints: # of numbers >= 4
2870 if (fCount < 4.0)
2871 {
2872 // for interoperability with Excel
2873 PushError( FormulaError::DivisionByZero);
2874 return;
2875 }
2876
2877 KahanSum vSum;
2878 double fMean = fSum.get() / fCount;
2879 for (double v : values)
2880 vSum += (v - fMean) * (v - fMean);
2881
2882 double fStdDev = sqrt(vSum.get() / (fCount - 1.0));
2883 if (fStdDev == 0.0)
2884 {
2885 PushError( FormulaError::DivisionByZero);
2886 return;
2887 }
2888
2889 KahanSum xpower4 = 0.0;
2890 for (double v : values)
2891 {
2892 double dx = (v - fMean) / fStdDev;
2893 xpower4 += (dx * dx) * (dx * dx);
2894 }
2895
2896 double k_d = (fCount - 2.0) * (fCount - 3.0);
2897 double k_l = fCount * (fCount + 1.0) / ((fCount - 1.0) * k_d);
2898 double k_t = 3.0 * (fCount - 1.0) * (fCount - 1.0) / k_d;
2899
2900 PushDouble(xpower4.get() * k_l - k_t);
2901}
2902
2904{
2905 short nParamCount = GetByte();
2906 KahanSum nVal = 0.0;
2907 double nValCount = 0.0;
2908 ScAddress aAdr;
2909 ScRange aRange;
2910 size_t nRefInList = 0;
2911 while ((nGlobalError == FormulaError::NONE) && (nParamCount-- > 0))
2912 {
2913 switch (GetStackType())
2914 {
2915 case svDouble :
2916 {
2917 double x = GetDouble();
2918 if (x > 0.0)
2919 {
2920 nVal += 1.0/x;
2921 nValCount++;
2922 }
2923 else
2924 SetError( FormulaError::IllegalArgument);
2925 break;
2926 }
2927 case svSingleRef :
2928 {
2929 PopSingleRef( aAdr );
2930 ScRefCellValue aCell(mrDoc, aAdr);
2931 if (aCell.hasNumeric())
2932 {
2933 double x = GetCellValue(aAdr, aCell);
2934 if (x > 0.0)
2935 {
2936 nVal += 1.0/x;
2937 nValCount++;
2938 }
2939 else
2940 SetError( FormulaError::IllegalArgument);
2941 }
2942 break;
2943 }
2944 case svDoubleRef :
2945 case svRefList :
2946 {
2947 FormulaError nErr = FormulaError::NONE;
2948 PopDoubleRef( aRange, nParamCount, nRefInList);
2949 double nCellVal;
2950 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
2951 if (aValIter.GetFirst(nCellVal, nErr))
2952 {
2953 if (nCellVal > 0.0)
2954 {
2955 nVal += 1.0/nCellVal;
2956 nValCount++;
2957 }
2958 else
2959 SetError( FormulaError::IllegalArgument);
2960 SetError(nErr);
2961 while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
2962 {
2963 if (nCellVal > 0.0)
2964 {
2965 nVal += 1.0/nCellVal;
2966 nValCount++;
2967 }
2968 else
2969 SetError( FormulaError::IllegalArgument);
2970 }
2971 SetError(nErr);
2972 }
2973 }
2974 break;
2975 case svMatrix :
2978 {
2979 ScMatrixRef pMat = GetMatrix();
2980 if (pMat)
2981 {
2982 SCSIZE nCount = pMat->GetElementCount();
2983 if (pMat->IsNumeric())
2984 {
2985 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
2986 {
2987 double x = pMat->GetDouble(nElem);
2988 if (x > 0.0)
2989 {
2990 nVal += 1.0/x;
2991 nValCount++;
2992 }
2993 else
2994 SetError( FormulaError::IllegalArgument);
2995 }
2996 }
2997 else
2998 {
2999 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3000 if (!pMat->IsStringOrEmpty(nElem))
3001 {
3002 double x = pMat->GetDouble(nElem);
3003 if (x > 0.0)
3004 {
3005 nVal += 1.0/x;
3006 nValCount++;
3007 }
3008 else
3009 SetError( FormulaError::IllegalArgument);
3010 }
3011 }
3012 }
3013 }
3014 break;
3015 default : SetError(FormulaError::IllegalParameter); break;
3016 }
3017 }
3018 if (nGlobalError == FormulaError::NONE)
3019 PushDouble( nValCount / nVal.get() );
3020 else
3022}
3023
3025{
3026 short nParamCount = GetByte();
3027 KahanSum nVal = 0.0;
3028 double nValCount = 0.0;
3029 ScAddress aAdr;
3030 ScRange aRange;
3031
3032 size_t nRefInList = 0;
3033 while ((nGlobalError == FormulaError::NONE) && (nParamCount-- > 0))
3034 {
3035 switch (GetStackType())
3036 {
3037 case svDouble :
3038 {
3039 double x = GetDouble();
3040 if (x > 0.0)
3041 {
3042 nVal += log(x);
3043 nValCount++;
3044 }
3045 else if ( x == 0.0 )
3046 {
3047 // value of 0 means that function result will be 0
3048 while ( nParamCount-- > 0 )
3049 PopError();
3050 PushDouble( 0.0 );
3051 return;
3052 }
3053 else
3054 SetError( FormulaError::IllegalArgument);
3055 break;
3056 }
3057 case svSingleRef :
3058 {
3059 PopSingleRef( aAdr );
3060 ScRefCellValue aCell(mrDoc, aAdr);
3061 if (aCell.hasNumeric())
3062 {
3063 double x = GetCellValue(aAdr, aCell);
3064 if (x > 0.0)
3065 {
3066 nVal += log(x);
3067 nValCount++;
3068 }
3069 else if ( x == 0.0 )
3070 {
3071 // value of 0 means that function result will be 0
3072 while ( nParamCount-- > 0 )
3073 PopError();
3074 PushDouble( 0.0 );
3075 return;
3076 }
3077 else
3078 SetError( FormulaError::IllegalArgument);
3079 }
3080 break;
3081 }
3082 case svDoubleRef :
3083 case svRefList :
3084 {
3085 FormulaError nErr = FormulaError::NONE;
3086 PopDoubleRef( aRange, nParamCount, nRefInList);
3087 double nCellVal;
3088 ScValueIterator aValIter(mrContext, aRange, mnSubTotalFlags);
3089 if (aValIter.GetFirst(nCellVal, nErr))
3090 {
3091 if (nCellVal > 0.0)
3092 {
3093 nVal += log(nCellVal);
3094 nValCount++;
3095 }
3096 else if ( nCellVal == 0.0 )
3097 {
3098 // value of 0 means that function result will be 0
3099 while ( nParamCount-- > 0 )
3100 PopError();
3101 PushDouble( 0.0 );
3102 return;
3103 }
3104 else
3105 SetError( FormulaError::IllegalArgument);
3106 SetError(nErr);
3107 while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
3108 {
3109 if (nCellVal > 0.0)
3110 {
3111 nVal += log(nCellVal);
3112 nValCount++;
3113 }
3114 else if ( nCellVal == 0.0 )
3115 {
3116 // value of 0 means that function result will be 0
3117 while ( nParamCount-- > 0 )
3118 PopError();
3119 PushDouble( 0.0 );
3120 return;
3121 }
3122 else
3123 SetError( FormulaError::IllegalArgument);
3124 }
3125 SetError(nErr);
3126 }
3127 }
3128 break;
3129 case svMatrix :
3132 {
3133 ScMatrixRef pMat = GetMatrix();
3134 if (pMat)
3135 {
3136 SCSIZE nCount = pMat->GetElementCount();
3137 if (pMat->IsNumeric())
3138 {
3139 for (SCSIZE ui = 0; ui < nCount; ui++)
3140 {
3141 double x = pMat->GetDouble(ui);
3142 if (x > 0.0)
3143 {
3144 nVal += log(x);
3145 nValCount++;
3146 }
3147 else if ( x == 0.0 )
3148 {
3149 // value of 0 means that function result will be 0
3150 while ( nParamCount-- > 0 )
3151 PopError();
3152 PushDouble( 0.0 );
3153 return;
3154 }
3155 else
3156 SetError( FormulaError::IllegalArgument);
3157 }
3158 }
3159 else
3160 {
3161 for (SCSIZE ui = 0; ui < nCount; ui++)
3162 {
3163 if (!pMat->IsStringOrEmpty(ui))
3164 {
3165 double x = pMat->GetDouble(ui);
3166 if (x > 0.0)
3167 {
3168 nVal += log(x);
3169 nValCount++;
3170 }
3171 else if ( x == 0.0 )
3172 {
3173 // value of 0 means that function result will be 0
3174 while ( nParamCount-- > 0 )
3175 PopError();
3176 PushDouble( 0.0 );
3177 return;
3178 }
3179 else
3180 SetError( FormulaError::IllegalArgument);
3181 }
3182 }
3183 }
3184 }
3185 }
3186 break;
3187 default : SetError(FormulaError::IllegalParameter); break;
3188 }
3189 }
3190 if (nGlobalError == FormulaError::NONE)
3191 PushDouble(exp(nVal.get() / nValCount));
3192 else
3194}
3195
3197{
3198 if ( MustHaveParamCount( GetByte(), 3 ) )
3199 {
3200 double sigma = GetDouble();
3201 double mue = GetDouble();
3202 double x = GetDouble();
3203 if (sigma < 0.0)
3204 PushError( FormulaError::IllegalArgument);
3205 else if (sigma == 0.0)
3206 PushError( FormulaError::DivisionByZero);
3207 else
3208 PushDouble((x-mue)/sigma);
3209 }
3210}
3211bool ScInterpreter::CalculateSkew(KahanSum& fSum, double& fCount, std::vector<double>& values)
3212{
3213 short nParamCount = GetByte();
3214 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3215 return false;
3216
3217 fSum = 0.0;
3218 fCount = 0.0;
3219 double fVal = 0.0;
3220 ScAddress aAdr;
3221 ScRange aRange;
3222 size_t nRefInList = 0;
3223 while (nParamCount-- > 0)
3224 {
3225 switch (GetStackType())
3226 {
3227 case svDouble :
3228 {
3229 fVal = GetDouble();
3230 fSum += fVal;
3231 values.push_back(fVal);
3232 fCount++;
3233 }
3234 break;
3235 case svSingleRef :
3236 {
3237 PopSingleRef( aAdr );
3238 ScRefCellValue aCell(mrDoc, aAdr);
3239 if (aCell.hasNumeric())
3240 {
3241 fVal = GetCellValue(aAdr, aCell);
3242 fSum += fVal;
3243 values.push_back(fVal);
3244 fCount++;
3245 }
3246 }
3247 break;
3248 case svDoubleRef :
3249 case svRefList :
3250 {
3251 PopDoubleRef( aRange, nParamCount, nRefInList);
3252 FormulaError nErr = FormulaError::NONE;
3253 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
3254 if (aValIter.GetFirst(fVal, nErr))
3255 {
3256 fSum += fVal;
3257 values.push_back(fVal);
3258 fCount++;
3259 SetError(nErr);
3260 while ((nErr == FormulaError::NONE) && aValIter.GetNext(fVal, nErr))
3261 {
3262 fSum += fVal;
3263 values.push_back(fVal);
3264 fCount++;
3265 }
3266 SetError(nErr);
3267 }
3268 }
3269 break;
3270 case svMatrix :
3273 {
3274 ScMatrixRef pMat = GetMatrix();
3275 if (pMat)
3276 {
3277 SCSIZE nCount = pMat->GetElementCount();
3278 if (pMat->IsNumeric())
3279 {
3280 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3281 {
3282 fVal = pMat->GetDouble(nElem);
3283 fSum += fVal;
3284 values.push_back(fVal);
3285 fCount++;
3286 }
3287 }
3288 else
3289 {
3290 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3291 if (!pMat->IsStringOrEmpty(nElem))
3292 {
3293 fVal = pMat->GetDouble(nElem);
3294 fSum += fVal;
3295 values.push_back(fVal);
3296 fCount++;
3297 }
3298 }
3299 }
3300 }
3301 break;
3302 default :
3303 SetError(FormulaError::IllegalParameter);
3304 break;
3305 }
3306 }
3307
3308 if (nGlobalError != FormulaError::NONE)
3309 {
3311 return false;
3312 } // if (nGlobalError != FormulaError::NONE)
3313 return true;
3314}
3315
3317{
3318 KahanSum fSum;
3319 double fCount;
3320 std::vector<double> values;
3321 if (!CalculateSkew( fSum, fCount, values))
3322 return;
3323 // SKEW/SKEWP's constraints: they require at least three numbers
3324 if (fCount < 3.0)
3325 {
3326 // for interoperability with Excel
3327 PushError(FormulaError::DivisionByZero);
3328 return;
3329 }
3330
3331 KahanSum vSum;
3332 double fMean = fSum.get() / fCount;
3333 for (double v : values)
3334 vSum += (v - fMean) * (v - fMean);
3335
3336 double fStdDev = sqrt( vSum.get() / (bSkewp ? fCount : (fCount - 1.0)));
3337 if (fStdDev == 0)
3338 {
3340 return;
3341 }
3342
3343 KahanSum xcube = 0.0;
3344 for (double v : values)
3345 {
3346 double dx = (v - fMean) / fStdDev;
3347 xcube += dx * dx * dx;
3348 }
3349
3350 if (bSkewp)
3351 PushDouble( xcube.get() / fCount );
3352 else
3353 PushDouble( ((xcube.get() * fCount) / (fCount - 1.0)) / (fCount - 2.0) );
3354}
3355
3357{
3358 CalculateSkewOrSkewp( false );
3359}
3360
3362{
3363 CalculateSkewOrSkewp( true );
3364}
3365
3366double ScInterpreter::GetMedian( vector<double> & rArray )
3367{
3368 size_t nSize = rArray.size();
3369 if (nSize == 0 || nGlobalError != FormulaError::NONE)
3370 {
3371 SetError( FormulaError::NoValue);
3372 return 0.0;
3373 }
3374
3375 // Upper median.
3376 size_t nMid = nSize / 2;
3377 vector<double>::iterator iMid = rArray.begin() + nMid;
3378 ::std::nth_element( rArray.begin(), iMid, rArray.end());
3379 if (nSize & 1)
3380 return *iMid; // Lower and upper median are equal.
3381 else
3382 {
3383 double fUp = *iMid;
3384 // Lower median.
3385 iMid = ::std::max_element( rArray.begin(), rArray.begin() + nMid);
3386 return (fUp + *iMid) / 2;
3387 }
3388}
3389
3391{
3392 sal_uInt8 nParamCount = GetByte();
3393 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3394 return;
3395 vector<double> aArray;
3396 GetNumberSequenceArray( nParamCount, aArray, false );
3397 PushDouble( GetMedian( aArray));
3398}
3399
3400double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile )
3401{
3402 size_t nSize = rArray.size();
3403 if (nSize == 1)
3404 return rArray[0];
3405 else
3406 {
3407 size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * (nSize-1)));
3408 double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1));
3409 OSL_ENSURE(nIndex < nSize, "GetPercentile: wrong index(1)");
3410 vector<double>::iterator iter = rArray.begin() + nIndex;
3411 ::std::nth_element( rArray.begin(), iter, rArray.end());
3412 if (fDiff == 0.0)
3413 return *iter;
3414 else
3415 {
3416 OSL_ENSURE(nIndex < nSize-1, "GetPercentile: wrong index(2)");
3417 double fVal = *iter;
3418 iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
3419 return fVal + fDiff * (*iter - fVal);
3420 }
3421 }
3422}
3423
3424double ScInterpreter::GetPercentileExclusive( vector<double> & rArray, double fPercentile )
3425{
3426 size_t nSize1 = rArray.size() + 1;
3427 if ( rArray.empty() || nSize1 == 1 || nGlobalError != FormulaError::NONE)
3428 {
3429 SetError( FormulaError::NoValue );
3430 return 0.0;
3431 }
3432 if ( fPercentile * nSize1 < 1.0 || fPercentile * nSize1 > static_cast<double>( nSize1 - 1 ) )
3433 {
3434 SetError( FormulaError::IllegalParameter );
3435 return 0.0;
3436 }
3437
3438 size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * nSize1 - 1 ));
3439 double fDiff = fPercentile * nSize1 - 1 - ::rtl::math::approxFloor( fPercentile * nSize1 - 1 );
3440 OSL_ENSURE(nIndex < ( nSize1 - 1 ), "GetPercentile: wrong index(1)");
3441 vector<double>::iterator iter = rArray.begin() + nIndex;
3442 ::std::nth_element( rArray.begin(), iter, rArray.end());
3443 if (fDiff == 0.0)
3444 return *iter;
3445 else
3446 {
3447 OSL_ENSURE(nIndex < nSize1, "GetPercentile: wrong index(2)");
3448 double fVal = *iter;
3449 iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
3450 return fVal + fDiff * (*iter - fVal);
3451 }
3452}
3453
3454void ScInterpreter::ScPercentile( bool bInclusive )
3455{
3456 if ( !MustHaveParamCount( GetByte(), 2 ) )
3457 return;
3458 double alpha = GetDouble();
3459 if ( bInclusive ? ( alpha < 0.0 || alpha > 1.0 ) : ( alpha <= 0.0 || alpha >= 1.0 ) )
3460 {
3462 return;
3463 }
3464 vector<double> aArray;
3465 GetNumberSequenceArray( 1, aArray, false );
3466 if ( aArray.empty() || nGlobalError != FormulaError::NONE )
3467 {
3468 SetError( FormulaError::NoValue );
3469 return;
3470 }
3471 if ( bInclusive )
3472 PushDouble( GetPercentile( aArray, alpha ));
3473 else
3474 PushDouble( GetPercentileExclusive( aArray, alpha ));
3475}
3476
3477void ScInterpreter::ScQuartile( bool bInclusive )
3478{
3479 if ( !MustHaveParamCount( GetByte(), 2 ) )
3480 return;
3481 double fFlag = ::rtl::math::approxFloor(GetDouble());
3482 if ( bInclusive ? ( fFlag < 0.0 || fFlag > 4.0 ) : ( fFlag <= 0.0 || fFlag >= 4.0 ) )
3483 {
3485 return;
3486 }
3487 vector<double> aArray;
3488 GetNumberSequenceArray( 1, aArray, false );
3489 if ( aArray.empty() || nGlobalError != FormulaError::NONE )
3490 {
3491 SetError( FormulaError::NoValue );
3492 return;
3493 }
3494 if ( bInclusive )
3495 PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentile( aArray, 0.25 * fFlag ) );
3496 else
3497 PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentileExclusive( aArray, 0.25 * fFlag ) );
3498}
3499
3501{
3502 sal_uInt8 nParamCount = GetByte();
3503 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3504 return;
3505 vector<double> aSortArray;
3506 GetSortArray( nParamCount, aSortArray, nullptr, false, false );
3507 SCSIZE nSize = aSortArray.size();
3508 if (nSize == 0 || nGlobalError != FormulaError::NONE)
3509 PushNoValue();
3510 else
3511 {
3512 SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
3513 double nOldVal = aSortArray[0];
3514 SCSIZE i;
3515 for ( i = 1; i < nSize; i++)
3516 {
3517 if (aSortArray[i] == nOldVal)
3518 nCount++;
3519 else
3520 {
3521 nOldVal = aSortArray[i];
3522 if (nCount > nMax)
3523 {
3524 nMax = nCount;
3525 nMaxIndex = i-1;
3526 }
3527 nCount = 1;
3528 }
3529 }
3530 if (nCount > nMax)
3531 {
3532 nMax = nCount;
3533 nMaxIndex = i-1;
3534 }
3535 if (nMax == 1 && nCount == 1)
3536 PushNoValue();
3537 else if (nMax == 1)
3538 PushDouble(nOldVal);
3539 else
3540 PushDouble(aSortArray[nMaxIndex]);
3541 }
3542}
3543
3545{
3546 sal_uInt8 nParamCount = GetByte();
3547 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3548 return;
3549 vector<double> aArray;
3550 GetNumberSequenceArray( nParamCount, aArray, false );
3551 vector< double > aSortArray( aArray );
3552 QuickSort( aSortArray, nullptr );
3553 SCSIZE nSize = aSortArray.size();
3554 if ( nSize == 0 || nGlobalError != FormulaError::NONE )
3555 PushNoValue();
3556 else
3557 {
3558 SCSIZE nMax = 1, nCount = 1;
3559 double nOldVal = aSortArray[ 0 ];
3560 vector< double > aResultArray( 1 );
3561 SCSIZE i;
3562 for ( i = 1; i < nSize; i++ )
3563 {
3564 if ( aSortArray[ i ] == nOldVal )
3565 nCount++;
3566 else
3567 {
3568 if ( nCount >= nMax && nCount > 1 )
3569 {
3570 if ( nCount > nMax )
3571 {
3572 nMax = nCount;
3573 if ( aResultArray.size() != 1 )
3574 vector< double >( 1 ).swap( aResultArray );
3575 aResultArray[ 0 ] = nOldVal;
3576 }
3577 else
3578 aResultArray.emplace_back( nOldVal );
3579 }
3580 nOldVal = aSortArray[ i ];
3581 nCount = 1;
3582 }
3583 }
3584 if ( nCount >= nMax && nCount > 1 )
3585 {
3586 if ( nCount > nMax )
3587 vector< double >().swap( aResultArray );
3588 aResultArray.emplace_back( nOldVal );
3589 }
3590 if ( nMax == 1 && nCount == 1 )
3591 PushNoValue();
3592 else if ( nMax == 1 )
3593 PushDouble( nOldVal ); // there is only 1 result, no reordering needed
3594 else
3595 {
3596 // sort resultArray according to ordering of aArray
3597 vector< vector< double > > aOrder;
3598 aOrder.resize( aResultArray.size(), vector< double >( 2 ) );
3599 for ( i = 0; i < aResultArray.size(); i++ )
3600 {
3601 for ( SCSIZE j = 0; j < nSize; j++ )
3602 {
3603 if ( aArray[ j ] == aResultArray[ i ] )
3604 {
3605 aOrder[ i ][ 0 ] = aResultArray[ i ];
3606 aOrder[ i ][ 1 ] = j;
3607 break;
3608 }
3609 }
3610 }
3611 sort( aOrder.begin(), aOrder.end(), []( const std::vector< double >& lhs,
3612 const std::vector< double >& rhs )
3613 { return lhs[ 1 ] < rhs[ 1 ]; } );
3614
3615 if ( bSingle )
3616 PushDouble( aOrder[ 0 ][ 0 ] );
3617 else
3618 {
3619 // put result in correct order in aResultArray
3620 for ( i = 0; i < aResultArray.size(); i++ )
3621 aResultArray[ i ] = aOrder[ i ][ 0 ];
3622 ScMatrixRef pResMatrix = GetNewMat( 1, aResultArray.size(), true );
3623 pResMatrix->PutDoubleVector( aResultArray, 0, 0 );
3624 PushMatrix( pResMatrix );
3625 }
3626 }
3627 }
3628}
3629
3631{
3632 if ( !MustHaveParamCount( GetByte(), 2 ) )
3633 return;
3634
3635 SCSIZE nCol = 0, nRow = 0;
3636 auto aArray = GetTopNumberArray(nCol, nRow);
3637 const auto nRankArraySize = aArray.size();
3638 if (nRankArraySize == 0 || nGlobalError != FormulaError::NONE)
3639 {
3640 PushNoValue();
3641 return;
3642 }
3643 assert(nRankArraySize == nCol * nRow);
3644
3645 std::vector<SCSIZE> aRankArray;
3646 aRankArray.reserve(nRankArraySize);
3647 std::transform(aArray.begin(), aArray.end(), std::back_inserter(aRankArray),
3648 [](double f) {
3649 f = rtl::math::approxFloor(f);
3650 // Valid ranks are >= 1.
3651 if (f < 1.0 || !o3tl::convertsToAtMost(f, std::numeric_limits<SCSIZE>::max()))
3652 return static_cast<SCSIZE>(0);
3653 return static_cast<SCSIZE>(f);
3654 });
3655
3656 vector<double> aSortArray;
3657 GetNumberSequenceArray(1, aSortArray, false );
3658 const SCSIZE nSize = aSortArray.size();
3659 if (nSize == 0 || nGlobalError != FormulaError::NONE)
3660 PushNoValue();
3661 else if (nRankArraySize == 1)
3662 {
3663 const SCSIZE k = aRankArray[0];
3664 if (k < 1 || nSize < k)
3665 PushNoValue();
3666 else
3667 {
3668 vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
3669 ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
3670 PushDouble( *iPos);
3671 }
3672 }
3673 else
3674 {
3675 std::set<SCSIZE> aIndices;
3676 for (SCSIZE n : aRankArray)
3677 {
3678 if (1 <= n && n <= nSize)
3679 aIndices.insert(bSmall ? n-1 : nSize-n);
3680 }
3681 // We can spare sorting when the total number of ranks is small enough.
3682 // Find only the elements at given indices if, arbitrarily, the index size is
3683 // smaller than 1/3 of the haystack array's size; just sort it squarely, otherwise.
3684 if (aIndices.size() < nSize/3)
3685 {
3686 auto itBegin = aSortArray.begin();
3687 for (SCSIZE i : aIndices)
3688 {
3689 auto it = aSortArray.begin() + i;
3690 std::nth_element(itBegin, it, aSortArray.end());
3691 itBegin = ++it;
3692 }
3693 }
3694 else
3695 std::sort(aSortArray.begin(), aSortArray.end());
3696
3697 aArray.clear();
3698 for (SCSIZE n : aRankArray)
3699 {
3700 if (1 <= n && n <= nSize)
3701 aArray.push_back( aSortArray[bSmall ? n-1 : nSize-n]);
3702 else
3703 aArray.push_back( CreateDoubleError( FormulaError::NoValue));
3704 }
3705 ScMatrixRef pResult = GetNewMat(nCol, nRow, aArray);
3706 PushMatrix(pResult);
3707 }
3708}
3709
3711{
3712 CalculateSmallLarge(false);
3713}
3714
3716{
3717 CalculateSmallLarge(true);
3718}
3719
3720void ScInterpreter::ScPercentrank( bool bInclusive )
3721{
3722 sal_uInt8 nParamCount = GetByte();
3723 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
3724 return;
3725 double fSignificance = ( nParamCount == 3 ? ::rtl::math::approxFloor( GetDouble() ) : 3.0 );
3726 if ( fSignificance < 1.0 )
3727 {
3729 return;
3730 }
3731 double fNum = GetDouble();
3732 vector<double> aSortArray;
3733 GetSortArray( 1, aSortArray, nullptr, false, false );
3734 SCSIZE nSize = aSortArray.size();
3735 if ( nSize == 0 || nGlobalError != FormulaError::NONE )
3736 PushNoValue();
3737 else
3738 {
3739 if ( fNum < aSortArray[ 0 ] || fNum > aSortArray[ nSize - 1 ] )
3740 PushNoValue();
3741 else
3742 {
3743 double fRes;
3744 if ( nSize == 1 )
3745 fRes = 1.0; // fNum == aSortArray[ 0 ], see test above
3746 else
3747 fRes = GetPercentrank( aSortArray, fNum, bInclusive );
3748 if ( fRes != 0.0 )
3749 {
3750 double fExp = ::rtl::math::approxFloor( log10( fRes ) ) + 1.0 - fSignificance;
3751 fRes = ::rtl::math::round( fRes * pow( 10, -fExp ) ) / pow( 10, -fExp );
3752 }
3753 PushDouble( fRes );
3754 }
3755 }
3756}
3757
3758double ScInterpreter::GetPercentrank( ::std::vector<double> & rArray, double fVal, bool bInclusive )
3759{
3760 SCSIZE nSize = rArray.size();
3761 double fRes;
3762 if ( fVal == rArray[ 0 ] )
3763 {
3764 if ( bInclusive )
3765 fRes = 0.0;
3766 else
3767 fRes = 1.0 / static_cast<double>( nSize + 1 );
3768 }
3769 else
3770 {
3771 SCSIZE nOldCount = 0;
3772 double fOldVal = rArray[ 0 ];
3773 SCSIZE i;
3774 for ( i = 1; i < nSize && rArray[ i ] < fVal; i++ )
3775 {
3776 if ( rArray[ i ] != fOldVal )
3777 {
3778 nOldCount = i;
3779 fOldVal = rArray[ i ];
3780 }
3781 }
3782 if ( rArray[ i ] != fOldVal )
3783 nOldCount = i;
3784 if ( fVal == rArray[ i ] )
3785 {
3786 if ( bInclusive )
3787 fRes = div( nOldCount, nSize - 1 );
3788 else
3789 fRes = static_cast<double>( i + 1 ) / static_cast<double>( nSize + 1 );
3790 }
3791 else
3792 {
3793 // nOldCount is the count of smaller entries
3794 // fVal is between rArray[ nOldCount - 1 ] and rArray[ nOldCount ]
3795 // use linear interpolation to find a position between the entries
3796 if ( nOldCount == 0 )
3797 {
3798 OSL_FAIL( "should not happen" );
3799 fRes = 0.0;
3800 }
3801 else
3802 {
3803 double fFract = ( fVal - rArray[ nOldCount - 1 ] ) /
3804 ( rArray[ nOldCount ] - rArray[ nOldCount - 1 ] );
3805 if ( bInclusive )
3806 fRes = div( static_cast<double>( nOldCount - 1 ) + fFract, nSize - 1 );
3807 else
3808 fRes = ( static_cast<double>(nOldCount) + fFract ) / static_cast<double>( nSize + 1 );
3809 }
3810 }
3811 }
3812 return fRes;
3813}
3814
3816{
3817 if ( !MustHaveParamCount( GetByte(), 2 ) )
3818 return;
3819 double alpha = GetDouble();
3820 if (alpha < 0.0 || alpha >= 1.0)
3821 {
3823 return;
3824 }
3825 vector<double> aSortArray;
3826 GetSortArray( 1, aSortArray, nullptr, false, false );
3827 SCSIZE nSize = aSortArray.size();
3828 if (nSize == 0 || nGlobalError != FormulaError::NONE)
3829 PushNoValue();
3830 else
3831 {
3832 sal_uLong nIndex = static_cast<sal_uLong>(::rtl::math::approxFloor(alpha*static_cast<double>(nSize)));
3833 if (nIndex % 2 != 0)
3834 nIndex--;
3835 nIndex /= 2;
3836 OSL_ENSURE(nIndex < nSize, "ScTrimMean: wrong index");
3837 KahanSum fSum = 0.0;
3838 for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
3839 fSum += aSortArray[i];
3840 PushDouble(fSum.get()/static_cast<double>(nSize-2*nIndex));
3841 }
3842}
3843
3844std::vector<double> ScInterpreter::GetTopNumberArray( SCSIZE& rCol, SCSIZE& rRow )
3845{
3846 std::vector<double> aArray;
3847 switch (GetStackType())
3848 {
3849 case svDouble:
3850 aArray.push_back(PopDouble());
3851 rCol = rRow = 1;
3852 break;
3853 case svSingleRef:
3854 {
3855 ScAddress aAdr;
3856 PopSingleRef(aAdr);
3857 ScRefCellValue aCell(mrDoc, aAdr);
3858 if (aCell.hasNumeric())
3859 {
3860 aArray.push_back(GetCellValue(aAdr, aCell));
3861 rCol = rRow = 1;
3862 }
3863 }
3864 break;
3865 case svDoubleRef:
3866 {
3867 ScRange aRange;
3868 PopDoubleRef(aRange, true);
3869 if (nGlobalError != FormulaError::NONE)
3870 break;
3871
3872 // give up unless the start and end are in the same sheet
3873 if (aRange.aStart.Tab() != aRange.aEnd.Tab())
3874 {
3875 SetError(FormulaError::IllegalParameter);
3876 break;
3877 }
3878
3879 // the range already is in order
3880 assert(aRange.aStart.Col() <= aRange.aEnd.Col());
3881 assert(aRange.aStart.Row() <= aRange.aEnd.Row());
3882 rCol = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3883 rRow = aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3884 aArray.reserve(rCol * rRow);
3885
3886 FormulaError nErr = FormulaError::NONE;
3887 double fCellVal;
3888 ScValueIterator aValIter(mrContext, aRange, mnSubTotalFlags);
3889 if (aValIter.GetFirst(fCellVal, nErr))
3890 {
3891 do
3892 aArray.push_back(fCellVal);
3893 while (aValIter.GetNext(fCellVal, nErr) && nErr == FormulaError::NONE);
3894 }
3895 if (aArray.size() != rCol * rRow)
3896 {
3897 aArray.clear();
3898 SetError(nErr);
3899 }
3900 }
3901 break;
3902 case svMatrix:
3905 {
3906 ScMatrixRef pMat = GetMatrix();
3907 if (!pMat)
3908 break;
3909
3910 if (pMat->IsNumeric())
3911 {
3912 SCSIZE nCount = pMat->GetElementCount();
3913 aArray.reserve(nCount);
3914 for (SCSIZE i = 0; i < nCount; ++i)
3915 aArray.push_back(pMat->GetDouble(i));
3916 pMat->GetDimensions(rCol, rRow);
3917 }
3918 else
3919 SetError(FormulaError::IllegalParameter);
3920 }
3921 break;
3922 default:
3923 SetError(FormulaError::IllegalParameter);
3924 break;
3925 }
3926 return aArray;
3927}
3928
3929void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray, bool bConvertTextInArray )
3930{
3931 ScAddress aAdr;
3932 ScRange aRange;
3933 const bool bIgnoreErrVal = bool(mnSubTotalFlags & SubtotalFlags::IgnoreErrVal);
3934 short nParam = nParamCount;
3935 size_t nRefInList = 0;
3936 ReverseStack( nParamCount );
3937 while (nParam-- > 0)
3938 {
3939 const StackVar eStackType = GetStackType();
3940 switch (eStackType)
3941 {
3942 case svDouble :
3943 rArray.push_back( PopDouble());
3944 break;
3945 case svSingleRef :
3946 {
3947 PopSingleRef( aAdr );
3948 ScRefCellValue aCell(mrDoc, aAdr);
3949 if (bIgnoreErrVal && aCell.hasError())
3950 ; // nothing
3951 else if (aCell.hasNumeric())
3952 rArray.push_back(GetCellValue(aAdr, aCell));
3953 }
3954 break;
3955 case svDoubleRef :
3956 case svRefList :
3957 {
3958 PopDoubleRef( aRange, nParam, nRefInList);
3959 if (nGlobalError != FormulaError::NONE)
3960 break;
3961
3962 aRange.PutInOrder();
3963 SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3964 nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3965 rArray.reserve( rArray.size() + nCellCount);
3966
3967 FormulaError nErr = FormulaError::NONE;
3968 double fCellVal;
3969 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
3970 if (aValIter.GetFirst( fCellVal, nErr))
3971 {
3972 if (bIgnoreErrVal)
3973 {
3974 if (nErr == FormulaError::NONE)
3975 rArray.push_back( fCellVal);
3976 while (aValIter.GetNext( fCellVal, nErr))
3977 {
3978 if (nErr == FormulaError::NONE)
3979 rArray.push_back( fCellVal);
3980 }
3981 }
3982 else
3983 {
3984 rArray.push_back( fCellVal);
3985 SetError(nErr);
3986 while ((nErr == FormulaError::NONE) && aValIter.GetNext( fCellVal, nErr))
3987 rArray.push_back( fCellVal);
3988 SetError(nErr);
3989 }
3990 }
3991 }
3992 break;
3993 case svMatrix :
3996 {
3997 ScMatrixRef pMat = GetMatrix();
3998 if (!pMat)
3999 break;
4000
4001 SCSIZE nCount = pMat->GetElementCount();
4002 rArray.reserve( rArray.size() + nCount);
4003 if (pMat->IsNumeric())
4004 {
4005 if (bIgnoreErrVal)
4006 {
4007 for (SCSIZE i = 0; i < nCount; ++i)
4008 {
4009 const double fVal = pMat->GetDouble(i);
4010 if (nGlobalError == FormulaError::NONE)
4011 rArray.push_back( fVal);
4012 else
4013 nGlobalError = FormulaError::NONE;
4014 }
4015 }
4016 else
4017 {
4018 for (SCSIZE i = 0; i < nCount; ++i)
4019 rArray.push_back( pMat->GetDouble(i));
4020 }
4021 }
4022 else if (bConvertTextInArray && eStackType == svMatrix)
4023 {
4024 for (SCSIZE i = 0; i < nCount; ++i)
4025 {
4026 if ( pMat->IsValue( i ) )
4027 {
4028 if (bIgnoreErrVal)
4029 {
4030 const double fVal = pMat->GetDouble(i);
4031 if (nGlobalError == FormulaError::NONE)
4032 rArray.push_back( fVal);
4033 else
4034 nGlobalError = FormulaError::NONE;
4035 }
4036 else
4037 rArray.push_back( pMat->GetDouble(i));
4038 }
4039 else
4040 {
4041 // tdf#88547 try to convert string to (date)value
4042 OUString aStr = pMat->GetString( i ).getString();
4043 if ( aStr.getLength() > 0 )
4044 {
4046 nGlobalError = FormulaError::NONE;
4047 double fVal = ConvertStringToValue( aStr );
4048 if ( nGlobalError == FormulaError::NONE )
4049 {
4050 rArray.push_back( fVal );
4051 nGlobalError = nErr;
4052 }
4053 else
4054 {
4055 if (!bIgnoreErrVal)
4056 rArray.push_back( CreateDoubleError( FormulaError::NoValue));
4057 // Propagate previous error if any, else
4058 // the current #VALUE! error, unless
4059 // ignoring error values.
4060 if (nErr != FormulaError::NONE)
4061 nGlobalError = nErr;
4062 else if (!bIgnoreErrVal)
4063 nGlobalError = FormulaError::NoValue;
4064 else
4065 nGlobalError = FormulaError::NONE;
4066 }
4067 }
4068 }
4069 }
4070 }
4071 else
4072 {
4073 if (bIgnoreErrVal)
4074 {
4075 for (SCSIZE i = 0; i < nCount; ++i)
4076 {
4077 if (pMat->IsValue(i))
4078 {
4079 const double fVal = pMat->GetDouble(i);
4080 if (nGlobalError == FormulaError::NONE)
4081 rArray.push_back( fVal);
4082 else
4083 nGlobalError = FormulaError::NONE;
4084 }
4085 }
4086 }
4087 else
4088 {
4089 for (SCSIZE i = 0; i < nCount; ++i)
4090 {
4091 if (pMat->IsValue(i))
4092 rArray.push_back( pMat->GetDouble(i));
4093 }
4094 }
4095 }
4096 }
4097 break;
4098 default :
4099 PopError();
4100 SetError( FormulaError::IllegalParameter);
4101 break;
4102 }
4103 if (nGlobalError != FormulaError::NONE)
4104 break; // while
4105 }
4106 // nParam > 0 in case of error, clean stack environment and obtain earlier
4107 // error if there was one.
4108 while (nParam-- > 0)
4109 PopError();
4110}
4111
4112void ScInterpreter::GetSortArray( sal_uInt8 nParamCount, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder, bool bConvertTextInArray, bool bAllowEmptyArray )
4113{
4114 GetNumberSequenceArray( nParamCount, rSortArray, bConvertTextInArray );
4115 if (rSortArray.size() > MAX_COUNT_DOUBLE_FOR_SORT(mrDoc.GetSheetLimits()))
4116 SetError( FormulaError::MatrixSize);
4117 else if ( rSortArray.empty() )
4118 {
4119 if ( bAllowEmptyArray )
4120 return;
4121 SetError( FormulaError::NoValue);
4122 }
4123
4124 if (nGlobalError == FormulaError::NONE)
4125 QuickSort( rSortArray, pIndexOrder);
4126}
4127
4128static void lcl_QuickSort( tools::Long nLo, tools::Long nHi, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
4129{
4130 // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().
4131
4132 using ::std::swap;
4133
4134 if (nHi - nLo == 1)
4135 {
4136 if (rSortArray[nLo] > rSortArray[nHi])
4137 {
4138 swap(rSortArray[nLo], rSortArray[nHi]);
4139 if (pIndexOrder)
4140 swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
4141 }
4142 return;
4143 }
4144
4145 tools::Long ni = nLo;
4146 tools::Long nj = nHi;
4147 do
4148 {
4149 double fLo = rSortArray[nLo];
4150 while (ni <= nHi && rSortArray[ni] < fLo) ni++;
4151 while (nj >= nLo && fLo < rSortArray[nj]) nj--;
4152 if (ni <= nj)
4153 {
4154 if (ni != nj)
4155 {
4156 swap(rSortArray[ni], rSortArray[nj]);
4157 if (pIndexOrder)
4158 swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
4159 }
4160
4161 ++ni;
4162 --nj;
4163 }
4164 }
4165 while (ni < nj);
4166
4167 if ((nj - nLo) < (nHi - ni))
4168 {
4169 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
4170 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
4171 }
4172 else
4173 {
4174 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
4175 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
4176 }
4177}
4178
4179void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
4180{
4181 tools::Long n = static_cast<tools::Long>(rSortArray.size());
4182
4183 if (pIndexOrder)
4184 {
4185 pIndexOrder->clear();
4186 pIndexOrder->reserve(n);
4187 for (tools::Long i = 0; i < n; ++i)
4188 pIndexOrder->push_back(i);
4189 }
4190
4191 if (n < 2)
4192 return;
4193
4194 size_t nValCount = rSortArray.size();
4195 for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
4196 {
4197 size_t nInd = comphelper::rng::uniform_size_distribution(0, nValCount-2);
4198 ::std::swap( rSortArray[i], rSortArray[nInd]);
4199 if (pIndexOrder)
4200 ::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
4201 }
4202
4203 lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
4204}
4205
4206void ScInterpreter::ScRank( bool bAverage )
4207{
4208 sal_uInt8 nParamCount = GetByte();
4209 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
4210 return;
4211 bool bAscending;
4212 if ( nParamCount == 3 )
4213 bAscending = GetBool();
4214 else
4215 bAscending = false;
4216
4217 vector<double> aSortArray;
4218 GetSortArray( 1, aSortArray, nullptr, false, false );
4219 double fVal = GetDouble();
4220 SCSIZE nSize = aSortArray.size();
4221 if ( nSize == 0 || nGlobalError != FormulaError::NONE )
4222 PushNoValue();
4223 else
4224 {
4225 if ( fVal < aSortArray[ 0 ] || fVal > aSortArray[ nSize - 1 ] )
4226 PushNoValue();
4227 else
4228 {
4229 double fLastPos = 0;
4230 double fFirstPos = -1.0;
4231 bool bFinished = false;
4232 SCSIZE i;
4233 for (i = 0; i < nSize && !bFinished; i++)
4234 {
4235 if ( aSortArray[ i ] == fVal )
4236 {
4237 if ( fFirstPos < 0 )
4238 fFirstPos = i + 1.0;
4239 }
4240 else
4241 {
4242 if ( aSortArray[ i ] > fVal )
4243 {
4244 fLastPos = i;
4245 bFinished = true;
4246 }
4247 }
4248 }
4249 if ( !bFinished )
4250 fLastPos = i;
4251 if ( !bAverage )
4252 {
4253 if ( bAscending )
4254 PushDouble( fFirstPos );
4255 else
4256 PushDouble( nSize + 1.0 - fLastPos );
4257 }
4258 else
4259 {
4260 if ( bAscending )
4261 PushDouble( ( fFirstPos + fLastPos ) / 2.0 );
4262 else
4263 PushDouble( nSize + 1.0 - ( fFirstPos + fLastPos ) / 2.0 );
4264 }
4265 }
4266 }
4267}
4268
4270{
4271 sal_uInt8 nParamCount = GetByte();
4272 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
4273 return;
4274 sal_uInt16 SaveSP = sp;
4275 double nMiddle = 0.0;
4276 KahanSum rVal = 0.0;
4277 double rValCount = 0.0;
4278 ScAddress aAdr;
4279 ScRange aRange;
4280 short nParam = nParamCount;
4281 size_t nRefInList = 0;
4282 while (nParam-- > 0)
4283 {
4284 switch (GetStackType())
4285 {
4286 case svDouble :
4287 rVal += GetDouble();
4288 rValCount++;
4289 break;
4290 case svSingleRef :
4291 {
4292 PopSingleRef( aAdr );
4293 ScRefCellValue aCell(mrDoc, aAdr);
4294 if (aCell.hasNumeric())
4295 {
4296 rVal += GetCellValue(aAdr, aCell);
4297 rValCount++;
4298 }
4299 }
4300 break;
4301 case svDoubleRef :
4302 case svRefList :
4303 {
4304 FormulaError nErr = FormulaError::NONE;
4305 double nCellVal;
4306 PopDoubleRef( aRange, nParam, nRefInList);
4307 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
4308 if (aValIter.GetFirst(nCellVal, nErr))
4309 {
4310 rVal += nCellVal;
4311 rValCount++;
4312 SetError(nErr);
4313 while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
4314 {
4315 rVal += nCellVal;
4316 rValCount++;
4317 }
4318 SetError(nErr);
4319 }
4320 }
4321 break;
4322 case svMatrix :
4325 {
4326 ScMatrixRef pMat = GetMatrix();
4327 if (pMat)
4328 {
4329 SCSIZE nCount = pMat->GetElementCount();
4330 if (pMat->IsNumeric())
4331 {
4332 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4333 {
4334 rVal += pMat->GetDouble(nElem);
4335 rValCount++;
4336 }
4337 }
4338 else
4339 {
4340 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4341 if (!pMat->IsStringOrEmpty(nElem))
4342 {
4343 rVal += pMat->GetDouble(nElem);
4344 rValCount++;
4345 }
4346 }
4347 }
4348 }
4349 break;
4350 default :
4351 SetError(FormulaError::IllegalParameter);
4352 break;
4353 }
4354 }
4355 if (nGlobalError != FormulaError::NONE)
4356 {
4358 return;
4359 }
4360 nMiddle = rVal.get() / rValCount;
4361 sp = SaveSP;
4362 rVal = 0.0;
4363 nParam = nParamCount;
4364 nRefInList = 0;
4365 while (nParam-- > 0)
4366 {
4367 switch (GetStackType())
4368 {
4369 case svDouble :
4370 rVal += std::abs(GetDouble() - nMiddle);
4371 break;
4372 case svSingleRef :
4373 {
4374 PopSingleRef( aAdr );
4375 ScRefCellValue aCell(mrDoc, aAdr);
4376 if (aCell.hasNumeric())
4377 rVal += std::abs(GetCellValue(aAdr, aCell) - nMiddle);
4378 }
4379 break;
4380 case svDoubleRef :
4381 case svRefList :
4382 {
4383 FormulaError nErr = FormulaError::NONE;
4384 double nCellVal;
4385 PopDoubleRef( aRange, nParam, nRefInList);
4386 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
4387 if (aValIter.GetFirst(nCellVal, nErr))
4388 {
4389 rVal += std::abs(nCellVal - nMiddle);
4390 while (aValIter.GetNext(nCellVal, nErr))
4391 rVal += std::abs(nCellVal - nMiddle);
4392 }
4393 }
4394 break;
4395 case svMatrix :
4398 {
4399 ScMatrixRef pMat = GetMatrix();
4400 if (pMat)
4401 {
4402 SCSIZE nCount = pMat->GetElementCount();
4403 if (pMat->IsNumeric())
4404 {
4405 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4406 {
4407 rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
4408 }
4409 }
4410 else
4411 {
4412 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4413 {
4414 if (!pMat->IsStringOrEmpty(nElem))
4415 rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
4416 }
4417 }
4418 }
4419 }
4420 break;
4421 default : SetError(FormulaError::IllegalParameter); break;
4422 }
4423 }
4424 PushDouble(rVal.get() / rValCount);
4425}
4426
4428{
4429 auto VarResult = []( double fVal, size_t /*nValCount*/ )
4430 {
4431 return fVal;
4432 };
4433 GetStVarParams( false /*bTextAsZero*/, VarResult);
4434}
4435
4437{
4438 sal_uInt8 nParamCount = GetByte();
4439 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
4440 return;
4441 double fUp, fLo;
4442 fUp = GetDouble();
4443 if (nParamCount == 4)
4444 fLo = GetDouble();
4445 else
4446 fLo = fUp;
4447 if (fLo > fUp)
4448 std::swap( fLo, fUp );
4449 ScMatrixRef pMatP = GetMatrix();
4450 ScMatrixRef pMatW = GetMatrix();
4451 if (!pMatP || !pMatW)
4453 else
4454 {
4455 SCSIZE nC1, nC2;
4456 SCSIZE nR1, nR2;
4457 pMatP->GetDimensions(nC1, nR1);
4458 pMatW->GetDimensions(nC2, nR2);
4459 if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
4460 nC2 == 0 || nR2 == 0)
4461 PushNA();
4462 else
4463 {
4464 KahanSum fSum = 0.0;
4465 KahanSum fRes = 0.0;
4466 bool bStop = false;
4467 double fP, fW;
4468 for ( SCSIZE i = 0; i < nC1 && !bStop; i++ )
4469 {
4470 for (SCSIZE j = 0; j < nR1 && !bStop; ++j )
4471 {
4472 if (pMatP->IsValue(i,j) && pMatW->IsValue(i,j))
4473 {
4474 fP = pMatP->GetDouble(i,j);
4475 fW = pMatW->GetDouble(i,j);
4476 if (fP < 0.0 || fP > 1.0)
4477 bStop = true;
4478 else
4479 {
4480 fSum += fP;
4481 if (fW >= fLo && fW <= fUp)
4482 fRes += fP;
4483 }
4484 }
4485 else
4486 SetError( FormulaError::IllegalArgument);
4487 }
4488 }
4489 if (bStop || std::abs((fSum -1.0).get()) > 1.0E-7)
4490 PushNoValue();
4491 else
4492 PushDouble(fRes.get());
4493 }
4494 }
4495}
4496
4498{
4499 // This is identical to ScPearson()
4500 ScPearson();
4501}
4502
4504{
4505 CalculatePearsonCovar( false, false, false );
4506}
4507
4509{
4510 CalculatePearsonCovar( false, false, true );
4511}
4512
4514{
4515 CalculatePearsonCovar( true, false, false );
4516}
4517
4518void ScInterpreter::CalculatePearsonCovar( bool _bPearson, bool _bStexy, bool _bSample )
4519{
4520 if ( !MustHaveParamCount( GetByte(), 2 ) )
4521 return;
4522 ScMatrixRef pMat1 = GetMatrix();
4523 ScMatrixRef pMat2 = GetMatrix();
4524 if (!pMat1 || !pMat2)
4525 {
4527 return;
4528 }
4529 SCSIZE nC1, nC2;
4530 SCSIZE nR1, nR2;
4531 pMat1->GetDimensions(nC1, nR1);
4532 pMat2->GetDimensions(nC2, nR2);
4533 if (nR1 != nR2 || nC1 != nC2)
4534 {
4536 return;
4537 }
4538 /* #i78250#
4539 * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
4540 * but the latter produces wrong results if the absolute values are high,
4541 * for example above 10^8
4542 */
4543 double fCount = 0.0;
4544 KahanSum fSumX = 0.0;
4545 KahanSum fSumY = 0.0;
4546
4547 for (SCSIZE i = 0; i < nC1; i++)
4548 {
4549 for (SCSIZE j = 0; j < nR1; j++)
4550 {
4551 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4552 {
4553 fSumX += pMat1->GetDouble(i,j);
4554 fSumY += pMat2->GetDouble(i,j);
4555 fCount++;
4556 }
4557 }
4558 }
4559 if (fCount < (_bStexy ? 3.0 : (_bSample ? 2.0 : 1.0)))
4560 PushNoValue();
4561 else
4562 {
4563 KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4564 KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4565 KahanSum fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2
4566 const double fMeanX = fSumX.get() / fCount;
4567 const double fMeanY = fSumY.get() / fCount;
4568 for (SCSIZE i = 0; i < nC1; i++)
4569 {
4570 for (SCSIZE j = 0; j < nR1; j++)
4571 {
4572 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4573 {
4574 const double fValX = pMat1->GetDouble(i,j);
4575 const double fValY = pMat2->GetDouble(i,j);
4576 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4577 if ( _bPearson )
4578 {
4579 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4580 fSumSqrDeltaY += (fValY - fMeanY) * (fValY - fMeanY);
4581 }
4582 }
4583 }
4584 }
4585 if ( _bPearson )
4586 {
4587 // tdf#94962 - Values below the numerical limit lead to unexpected results
4588 if (fSumSqrDeltaX < ::std::numeric_limits<double>::min()
4589 || (!_bStexy && fSumSqrDeltaY < ::std::numeric_limits<double>::min()))
4590 PushError( FormulaError::DivisionByZero);
4591 else if ( _bStexy )
4592 PushDouble( sqrt( ( fSumSqrDeltaY - fSumDeltaXDeltaY *
4593 fSumDeltaXDeltaY / fSumSqrDeltaX ).get() / (fCount-2)));
4594 else
4595 PushDouble( fSumDeltaXDeltaY.get() / sqrt( fSumSqrDeltaX.get() * fSumSqrDeltaY.get() ));
4596 }
4597 else
4598 {
4599 if ( _bSample )
4600 PushDouble( fSumDeltaXDeltaY.get() / ( fCount - 1 ) );
4601 else
4602 PushDouble( fSumDeltaXDeltaY.get() / fCount);
4603 }
4604 }
4605}
4606
4608{
4609 // Same as ScPearson()*ScPearson()
4610 ScPearson();
4611 if (nGlobalError != FormulaError::NONE)
4612 return;
4613
4614 switch (GetStackType())
4615 {
4616 case svDouble:
4617 {
4618 double fVal = PopDouble();
4619 PushDouble( fVal * fVal);
4620 }
4621 break;
4622 default:
4623 PopError();
4624 PushNoValue();
4625 }
4626}
4627
4629{
4630 CalculatePearsonCovar( true, true, false );
4631}
4633{
4634 if ( !MustHaveParamCount( GetByte(), 2 ) )
4635 return;
4636 ScMatrixRef pMat1 = GetMatrix();
4637 ScMatrixRef pMat2 = GetMatrix();
4638 if (!pMat1 || !pMat2)
4639 {
4641 return;
4642 }
4643 SCSIZE nC1, nC2;
4644 SCSIZE nR1, nR2;
4645 pMat1->GetDimensions(nC1, nR1);
4646 pMat2->GetDimensions(nC2, nR2);
4647 if (nR1 != nR2 || nC1 != nC2)
4648 {
4650 return;
4651 }
4652 // #i78250# numerical stability improved
4653 double fCount = 0.0;
4654 KahanSum fSumX = 0.0;
4655 KahanSum fSumY = 0.0;
4656
4657 for (SCSIZE i = 0; i < nC1; i++)
4658 {
4659 for (SCSIZE j = 0; j < nR1; j++)
4660 {
4661 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4662 {
4663 fSumX += pMat1->GetDouble(i,j);
4664 fSumY += pMat2->GetDouble(i,j);
4665 fCount++;
4666 }
4667 }
4668 }
4669 if (fCount < 1.0)
4670 PushNoValue();
4671 else
4672 {
4673 KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4674 KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4675 double fMeanX = fSumX.get() / fCount;
4676 double fMeanY = fSumY.get() / fCount;
4677 for (SCSIZE i = 0; i < nC1; i++)
4678 {
4679 for (SCSIZE j = 0; j < nR1; j++)
4680 {
4681 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4682 {
4683 double fValX = pMat1->GetDouble(i,j);
4684 double fValY = pMat2->GetDouble(i,j);
4685 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4686 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4687 }
4688 }
4689 }
4690 if (fSumSqrDeltaX == 0.0)
4691 PushError( FormulaError::DivisionByZero);
4692 else
4693 {
4694 if ( bSlope )
4695 PushDouble( fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get());
4696 else
4697 PushDouble( fMeanY - fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * fMeanX);
4698 }
4699 }
4700}
4701
4703{
4705}
4706
4708{
4710}
4711
4713{
4714 if ( !MustHaveParamCount( GetByte(), 3 ) )
4715 return;
4716 ScMatrixRef pMat1 = GetMatrix();
4717 ScMatrixRef pMat2 = GetMatrix();
4718 if (!pMat1 || !pMat2)
4719 {
4721 return;
4722 }
4723 SCSIZE nC1, nC2;
4724 SCSIZE nR1, nR2;
4725 pMat1->GetDimensions(nC1, nR1);
4726 pMat2->GetDimensions(nC2, nR2);
4727 if (nR1 != nR2 || nC1 != nC2)
4728 {
4730 return;
4731 }
4732 double fVal = GetDouble();
4733 // #i78250# numerical stability improved
4734 double fCount = 0.0;
4735 KahanSum fSumX = 0.0;
4736 KahanSum fSumY = 0.0;
4737
4738 for (SCSIZE i = 0; i < nC1; i++)
4739 {
4740 for (SCSIZE j = 0; j < nR1; j++)
4741 {
4742 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4743 {
4744 fSumX += pMat1->GetDouble(i,j);
4745 fSumY += pMat2->GetDouble(i,j);
4746 fCount++;
4747 }
4748 }
4749 }
4750 if (fCount < 1.0)
4751 PushNoValue();
4752 else
4753 {
4754 KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4755 KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4756 double fMeanX = fSumX.get() / fCount;
4757 double fMeanY = fSumY.get() / fCount;
4758 for (SCSIZE i = 0; i < nC1; i++)
4759 {
4760 for (SCSIZE j = 0; j < nR1; j++)
4761 {
4762 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4763 {
4764 double fValX = pMat1->GetDouble(i,j);
4765 double fValY = pMat2->GetDouble(i,j);
4766 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4767 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4768 }
4769 }
4770 }
4771 if (fSumSqrDeltaX == 0.0)
4772 PushError( FormulaError::DivisionByZero);
4773 else
4774 PushDouble( fMeanY + fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * (fVal - fMeanX));
4775 }
4776}
4777
4778static void lcl_roundUpNearestPow2(SCSIZE& nNum, SCSIZE& nNumBits)
4779{
4780 // Find the least power of 2 that is less than or equal to nNum.
4781 SCSIZE nPow2(1);
4782 nNumBits = std::numeric_limits<SCSIZE>::digits;
4783 nPow2 <<= (nNumBits - 1);
4784 while (nPow2 >= 1)
4785 {
4786 if (nNum & nPow2)
4787 break;
4788
4789 --nNumBits;
4790 nPow2 >>= 1;
4791 }
4792
4793 if (nPow2 != nNum)
4794 nNum = nPow2 ? (nPow2 << 1) : 1;
4795 else
4796 --nNumBits;
4797}
4798
4800{
4801 SCSIZE nOut = 0;
4802 for (SCSIZE nMask = 1; nMask < nBound; nMask <<= 1)
4803 {
4804 nOut <<= 1;
4805
4806 if (nIn & nMask)
4807 nOut |= 1;
4808 }
4809
4810 return nOut;
4811}
4812
4813namespace {
4814
4815// Computes and stores twiddle factors for computing DFT later.
4816struct ScTwiddleFactors
4817{
4818 ScTwiddleFactors(SCSIZE nN, bool bInverse) :
4819 mfWReal(nN),
4820 mfWImag(nN),
4821 mnN(nN),
4822 mbInverse(bInverse)
4823 {}
4824
4825 void Compute();
4826
4827 void Conjugate()
4828 {
4829 mbInverse = !mbInverse;
4830 for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
4831 mfWImag[nIdx] = -mfWImag[nIdx];
4832 }
4833
4834 std::vector<double> mfWReal;
4835 std::vector<double> mfWImag;
4836 SCSIZE mnN;
4837 bool mbInverse;
4838};
4839
4840}
4841
4842void ScTwiddleFactors::Compute()
4843{
4844 mfWReal.resize(mnN);
4845 mfWImag.resize(mnN);
4846
4847 double nW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnN);
4848
4849 if (mnN == 1)
4850 {
4851 mfWReal[0] = 1.0;
4852 mfWImag[0] = 0.0;
4853 }
4854 else if (mnN == 2)
4855 {
4856 mfWReal[0] = 1;
4857 mfWImag[0] = 0;
4858
4859 mfWReal[1] = -1;
4860 mfWImag[1] = 0;
4861 }
4862 else if (mnN == 4)
4863 {
4864 mfWReal[0] = 1;
4865 mfWImag[0] = 0;
4866
4867 mfWReal[1] = 0;
4868 mfWImag[1] = (mbInverse ? 1.0 : -1.0);
4869
4870 mfWReal[2] = -1;
4871 mfWImag[2] = 0;
4872
4873 mfWReal[3] = 0;
4874 mfWImag[3] = (mbInverse ? -1.0 : 1.0);
4875 }
4876 else if ((mnN % 4) == 0)
4877 {
4878 const SCSIZE nQSize = mnN >> 2;
4879 // Compute cos of the start quadrant.
4880 // This is the first quadrant if mbInverse == true, else it is the fourth quadrant.
4881 for (SCSIZE nIdx = 0; nIdx <= nQSize; ++nIdx)
4882 mfWReal[nIdx] = cos(nW*static_cast<double>(nIdx));
4883
4884 if (mbInverse)
4885 {
4886 const SCSIZE nQ1End = nQSize;
4887 // First quadrant
4888 for (SCSIZE nIdx = 0; nIdx <= nQ1End; ++nIdx)
4889 mfWImag[nIdx] = mfWReal[nQ1End-nIdx];
4890
4891 // Second quadrant
4892 const SCSIZE nQ2End = nQ1End << 1;
4893 for (SCSIZE nIdx = nQ1End+1; nIdx <= nQ2End; ++nIdx)
4894 {
4895 mfWReal[nIdx] = -mfWReal[nQ2End - nIdx];
4896 mfWImag[nIdx] = mfWImag[nQ2End - nIdx];
4897 }
4898
4899 // Third quadrant
4900 const SCSIZE nQ3End = nQ2End + nQ1End;
4901 for (SCSIZE nIdx = nQ2End+1; nIdx <= nQ3End; ++nIdx)
4902 {
4903 mfWReal[nIdx] = -mfWReal[nIdx - nQ2End];
4904 mfWImag[nIdx] = -mfWImag[nIdx - nQ2End];
4905 }
4906
4907 // Fourth Quadrant
4908 for (SCSIZE nIdx = nQ3End+1; nIdx < mnN; ++nIdx)
4909 {
4910 mfWReal[nIdx] = mfWReal[mnN - nIdx];
4911 mfWImag[nIdx] = -mfWImag[mnN - nIdx];
4912 }
4913 }
4914 else
4915 {
4916 const SCSIZE nQ4End = nQSize;
4917 const SCSIZE nQ3End = nQSize << 1;
4918 const SCSIZE nQ2End = nQ3End + nQSize;
4919
4920 // Fourth quadrant.
4921 for (SCSIZE nIdx = 0; nIdx <= nQ4End; ++nIdx)
4922 mfWImag[nIdx] = -mfWReal[nQ4End - nIdx];
4923
4924 // Third quadrant.
4925 for (SCSIZE nIdx = nQ4End+1; nIdx <= nQ3End; ++nIdx)
4926 {
4927 mfWReal[nIdx] = -mfWReal[nQ3End - nIdx];
4928 mfWImag[nIdx] = mfWImag[nQ3End - nIdx];
4929 }
4930
4931 // Second quadrant.
4932 for (SCSIZE nIdx = nQ3End+1; nIdx <= nQ2End; ++nIdx)
4933 {
4934 mfWReal[nIdx] = -mfWReal[nIdx - nQ3End];
4935 mfWImag[nIdx] = -mfWImag[nIdx - nQ3End];
4936 }
4937
4938 // First quadrant.
4939 for (SCSIZE nIdx = nQ2End+1; nIdx < mnN; ++nIdx)
4940 {
4941 mfWReal[nIdx] = mfWReal[mnN - nIdx];
4942 mfWImag[nIdx] = -mfWImag[mnN - nIdx];
4943 }
4944 }
4945 }
4946 else
4947 {
4948 for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
4949 {
4950 double fAngle = nW*static_cast<double>(nIdx);
4951 mfWReal[nIdx] = cos(fAngle);
4952 mfWImag[nIdx] = sin(fAngle);
4953 }
4954 }
4955}
4956
4957namespace {
4958
4959// A radix-2 decimation in time FFT algorithm for complex valued input.
4960class ScComplexFFT2
4961{
4962public:
4963 // rfArray.size() would always be even and a power of 2. (asserted in prepare())
4964 // rfArray's first half contains the real parts and the later half contains the imaginary parts.
4965 ScComplexFFT2(std::vector<double>& raArray, bool bInverse, bool bPolar, double fMinMag,
4966 ScTwiddleFactors& rTF, bool bSubSampleTFs = false, bool bDisableNormalize = false) :
4967 mrArray(raArray),
4968 mfWReal(rTF.mfWReal),
4969 mfWImag(rTF.mfWImag),
4970 mnPoints(raArray.size()/2),
4971 mnStages(0),
4972 mfMinMag(fMinMag),
4973 mbInverse(bInverse),
4974 mbPolar(bPolar),
4975 mbDisableNormalize(bDisableNormalize),
4976 mbSubSampleTFs(bSubSampleTFs)
4977 {}
4978
4979 void Compute();
4980
4981private:
4982
4983 void prepare();
4984
4985 double getReal(SCSIZE nIdx)
4986 {
4987 return mrArray[nIdx];
4988 }
4989
4990 void setReal(double fVal, SCSIZE nIdx)
4991 {
4992 mrArray[nIdx] = fVal;
4993 }
4994
4995 double getImag(SCSIZE nIdx)
4996 {
4997 return mrArray[mnPoints + nIdx];
4998 }
4999
5000 void setImag(double fVal, SCSIZE nIdx)
5001 {
5002 mrArray[mnPoints + nIdx] = fVal;
5003 }
5004
5005 SCSIZE getTFactorIndex(SCSIZE nPtIndex, SCSIZE nTfIdxScaleBits)
5006 {
5007 return ( ( nPtIndex << nTfIdxScaleBits ) & ( mnPoints - 1 ) ); // (x & (N-1)) is same as (x % N) but faster.
5008 }
5009
5010 void computeFly(SCSIZE nTopIdx, SCSIZE nBottomIdx, SCSIZE nWIdx1, SCSIZE nWIdx2)
5011 {
5012 if (mbSubSampleTFs)
5013 {
5014 nWIdx1 <<= 1;
5015 nWIdx2 <<= 1;
5016 }
5017
5018 const double x1r = getReal(nTopIdx);
5019 const double x2r = getReal(nBottomIdx);
5020
5021 const double& w1r = mfWReal[nWIdx1];
5022 const double& w1i = mfWImag[nWIdx1];
5023
5024 const double& w2r = mfWReal[nWIdx2];
5025 const double& w2i = mfWImag[nWIdx2];
5026
5027 const double x1i = getImag(nTopIdx);
5028 const double x2i = getImag(nBottomIdx);
5029
5030 setReal(x1r + x2r*w1r - x2i*w1i, nTopIdx);
5031 setImag(x1i + x2i*w1r + x2r*w1i, nTopIdx);
5032
5033 setReal(x1r + x2r*w2r - x2i*w2i, nBottomIdx);
5034 setImag(x1i + x2i*w2r + x2r*w2i, nBottomIdx);
5035 }
5036
5037 std::vector<double>& mrArray;
5038 std::vector<double>& mfWReal;
5039 std::vector<double>& mfWImag;
5040 SCSIZE mnPoints;
5041 SCSIZE mnStages;
5042 double mfMinMag;
5043 bool mbInverse:1;
5044 bool mbPolar:1;
5045 bool mbDisableNormalize:1;
5046 bool mbSubSampleTFs:1;
5047};
5048
5049}
5050
5051void ScComplexFFT2::prepare()
5052{
5053 SCSIZE nPoints = mnPoints;
5054 lcl_roundUpNearestPow2(nPoints, mnStages);
5055 assert(nPoints == mnPoints);
5056
5057 // Reorder array by bit-reversed indices.
5058 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5059 {
5060 SCSIZE nRevIdx = lcl_bitReverse(nIdx, mnPoints);
5061 if (nIdx < nRevIdx)
5062 {
5063 double fTmp = getReal(nIdx);
5064 setReal(getReal(nRevIdx), nIdx);
5065 setReal(fTmp, nRevIdx);
5066
5067 fTmp = getImag(nIdx);
5068 setImag(getImag(nRevIdx), nIdx);
5069 setImag(fTmp, nRevIdx);
5070 }
5071 }
5072}
5073
5074static void lcl_normalize(std::vector<double>& rCmplxArray, bool bScaleOnlyReal)
5075{
5076 const SCSIZE nPoints = rCmplxArray.size()/2;
5077 const double fScale = 1.0/static_cast<double>(nPoints);
5078
5079 // Scale the real part
5080 for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
5081 rCmplxArray[nIdx] *= fScale;
5082
5083 if (!bScaleOnlyReal)
5084 {
5085 const SCSIZE nLen = nPoints*2;
5086 for (SCSIZE nIdx = nPoints; nIdx < nLen; ++nIdx)
5087 rCmplxArray[nIdx] *= fScale;
5088 }
5089}
5090
5091static void lcl_convertToPolar(std::vector<double>& rCmplxArray, double fMinMag)
5092{
5093 const SCSIZE nPoints = rCmplxArray.size()/2;
5094 double fMag, fPhase, fR, fI;
5095 for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
5096 {
5097 fR = rCmplxArray[nIdx];
5098 fI = rCmplxArray[nPoints+nIdx];
5099 fMag = sqrt(fR*fR + fI*fI);
5100 if (fMag < fMinMag)
5101 {
5102 fMag = 0.0;
5103 fPhase = 0.0;
5104 }
5105 else
5106 {
5107 fPhase = atan2(fI, fR);
5108 }
5109
5110 rCmplxArray[nIdx] = fMag;
5111 rCmplxArray[nPoints+nIdx] = fPhase;
5112 }
5113}
5114
5115void ScComplexFFT2::Compute()
5116{
5117 prepare();
5118
5119 const SCSIZE nFliesInStage = mnPoints/2;
5120 for (SCSIZE nStage = 0; nStage < mnStages; ++nStage)
5121 {
5122 const SCSIZE nTFIdxScaleBits = mnStages - nStage - 1; // Twiddle factor index's scale factor in bits.
5123 const SCSIZE nFliesInGroup = SCSIZE(1) << nStage;
5124 const SCSIZE nGroups = nFliesInStage/nFliesInGroup;
5125 const SCSIZE nFlyWidth = nFliesInGroup;
5126 for (SCSIZE nGroup = 0, nFlyTopIdx = 0; nGroup < nGroups; ++nGroup)
5127 {
5128 for (SCSIZE nFly = 0; nFly < nFliesInGroup; ++nFly, ++nFlyTopIdx)
5129 {
5130 SCSIZE nFlyBottomIdx = nFlyTopIdx + nFlyWidth;
5131 SCSIZE nWIdx1 = getTFactorIndex(nFlyTopIdx, nTFIdxScaleBits);
5132 SCSIZE nWIdx2 = getTFactorIndex(nFlyBottomIdx, nTFIdxScaleBits);
5133
5134 computeFly(nFlyTopIdx, nFlyBottomIdx, nWIdx1, nWIdx2);
5135 }
5136
5137 nFlyTopIdx += nFlyWidth;
5138 }
5139 }
5140
5141 if (mbPolar)
5142 lcl_convertToPolar(mrArray, mfMinMag);
5143
5144 // Normalize after converting to polar, so we have a chance to
5145 // save O(mnPoints) flops.
5146 if (mbInverse && !mbDisableNormalize)
5147 lcl_normalize(mrArray, mbPolar);
5148}
5149
5150namespace {
5151
5152// Bluestein's algorithm or chirp z-transform algorithm that can be used to
5153// compute DFT of a complex valued input of any length N in O(N lgN) time.
5154class ScComplexBluesteinFFT
5155{
5156public:
5157
5158 ScComplexBluesteinFFT(std::vector<double>& rArray, bool bReal, bool bInverse,
5159 bool bPolar, double fMinMag, bool bDisableNormalize = false) :
5160 mrArray(rArray),
5161 mnPoints(rArray.size()/2), // rArray should have space for imaginary parts even if real input.
5162 mfMinMag(fMinMag),
5163 mbReal(bReal),
5164 mbInverse(bInverse),
5165 mbPolar(bPolar),
5166 mbDisableNormalize(bDisableNormalize)
5167 {}
5168
5169 void Compute();
5170
5171private:
5172 std::vector<double>& mrArray;
5173 const SCSIZE mnPoints;
5174 double mfMinMag;
5175 bool mbReal:1;
5176 bool mbInverse:1;
5177 bool mbPolar:1;
5178 bool mbDisableNormalize:1;
5179};
5180
5181}
5182
5183void ScComplexBluesteinFFT::Compute()
5184{
5185 std::vector<double> aRealScalars(mnPoints);
5186 std::vector<double> aImagScalars(mnPoints);
5187 double fW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnPoints);
5188 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5189 {
5190 double fAngle = 0.5*fW*static_cast<double>(nIdx*nIdx);
5191 aRealScalars[nIdx] = cos(fAngle);
5192 aImagScalars[nIdx] = sin(fAngle);
5193 }
5194
5195 SCSIZE nMinSize = mnPoints*2 - 1;
5196 SCSIZE nExtendedLength = nMinSize, nTmp = 0;
5197 lcl_roundUpNearestPow2(nExtendedLength, nTmp);
5198 std::vector<double> aASignal(nExtendedLength*2); // complex valued
5199 std::vector<double> aBSignal(nExtendedLength*2); // complex valued
5200
5201 double fReal, fImag;
5202 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5203 {
5204 // Real part of A signal.
5205 aASignal[nIdx] = mrArray[nIdx]*aRealScalars[nIdx] + (mbReal ? 0.0 : -mrArray[mnPoints+nIdx]*aImagScalars[nIdx]);
5206 // Imaginary part of A signal.
5207 aASignal[nExtendedLength + nIdx] = mrArray[nIdx]*aImagScalars[nIdx] + (mbReal ? 0.0 : mrArray[mnPoints+nIdx]*aRealScalars[nIdx]);
5208
5209 // Real part of B signal.
5210 aBSignal[nIdx] = fReal = aRealScalars[nIdx];
5211 // Imaginary part of B signal.
5212 aBSignal[nExtendedLength + nIdx] = fImag = -aImagScalars[nIdx]; // negative sign because B signal is the conjugation of the scalars.
5213
5214 if (nIdx)
5215 {
5216 // B signal needs a mirror of its part in 0 < n < mnPoints at the tail end.
5217 aBSignal[nExtendedLength - nIdx] = fReal;
5218 aBSignal[(nExtendedLength<<1) - nIdx] = fImag;
5219 }
5220 }
5221
5222 {
5223 ScTwiddleFactors aTF(nExtendedLength, false /*not inverse*/);
5224 aTF.Compute();
5225
5226 // Do complex-FFT2 of both A and B signal.
5227 ScComplexFFT2 aFFT2A(aASignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
5228 aTF, false /*no subsample*/, true /* disable normalize */);
5229 aFFT2A.Compute();
5230
5231 ScComplexFFT2 aFFT2B(aBSignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
5232 aTF, false /*no subsample*/, true /* disable normalize */);
5233 aFFT2B.Compute();
5234
5235 double fAR, fAI, fBR, fBI;
5236 for (SCSIZE nIdx = 0; nIdx < nExtendedLength; ++nIdx)
5237 {
5238 fAR = aASignal[nIdx];
5239 fAI = aASignal[nExtendedLength + nIdx];
5240 fBR = aBSignal[nIdx];
5241 fBI = aBSignal[nExtendedLength + nIdx];
5242
5243 // Do point-wise product inplace in A signal.
5244 aASignal[nIdx] = fAR*fBR - fAI*fBI;
5245 aASignal[nExtendedLength + nIdx] = fAR*fBI + fAI*fBR;
5246 }
5247
5248 // Do complex-inverse-FFT2 of aASignal.
5249 aTF.Conjugate();
5250 ScComplexFFT2 aFFT2AI(aASignal, true /*inverse*/, false /*no polar*/, 0.0 /* no clipping */, aTF); // Need normalization here.
5251 aFFT2AI.Compute();
5252 }
5253
5254 // Point-wise multiply with scalars.
5255 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5256 {
5257 fReal = aASignal[nIdx];
5258 fImag = aASignal[nExtendedLength + nIdx];
5259 mrArray[nIdx] = fReal*aRealScalars[nIdx] - fImag*aImagScalars[nIdx]; // no conjugation needed here.
5260 mrArray[mnPoints + nIdx] = fReal*aImagScalars[nIdx] + fImag*aRealScalars[nIdx];
5261 }
5262
5263 // Normalize/Polar operations
5264 if (mbPolar)
5265 lcl_convertToPolar(mrArray, mfMinMag);
5266
5267 // Normalize after converting to polar, so we have a chance to
5268 // save O(mnPoints) flops.
5269 if (mbInverse && !mbDisableNormalize)
5270 lcl_normalize(mrArray, mbPolar);
5271}
5272
5273namespace {
5274
5275// Computes DFT of an even length(N) real-valued input by using a
5276// ScComplexFFT2 if N == 2^k for some k or else by using a ScComplexBluesteinFFT
5277// with a complex valued input of length = N/2.
5278class ScRealFFT
5279{
5280public:
5281
5282 ScRealFFT(std::vector<double>& rInArray, std::vector<double>& rOutArray, bool bInverse,
5283 bool bPolar, double fMinMag) :
5284 mrInArray(rInArray),
5285 mrOutArray(rOutArray),
5286 mfMinMag(fMinMag),
5287 mbInverse(bInverse),
5288 mbPolar(bPolar)
5289 {}
5290
5291 void Compute();
5292
5293private:
5294 std::vector<double>& mrInArray;
5295 std::vector<double>& mrOutArray;
5296 double mfMinMag;
5297 bool mbInverse:1;
5298 bool mbPolar:1;
5299};
5300
5301}
5302
5303void ScRealFFT::Compute()
5304{
5305 // input length has to be even to do this optimization.
5306 assert(mrInArray.size() % 2 == 0);
5307 assert(mrInArray.size()*2 == mrOutArray.size());
5308 // nN is the number of points in the complex-fft input
5309 // which will be half of the number of points in real array.
5310 const SCSIZE nN = mrInArray.size()/2;
5311 if (nN == 0)
5312 {
5313 mrOutArray[0] = mrInArray[0];
5314 mrOutArray[1] = 0.0;
5315 return;
5316 }
5317
5318 // work array should be the same length as mrInArray
5319 std::vector<double> aWorkArray(nN*2);
5320 for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
5321 {
5322 SCSIZE nDoubleIdx = 2*nIdx;
5323 // Use even elements as real part
5324 aWorkArray[nIdx] = mrInArray[nDoubleIdx];
5325 // and odd elements as imaginary part of the contrived complex sequence.
5326 aWorkArray[nN+nIdx] = mrInArray[nDoubleIdx+1];
5327 }
5328
5329 ScTwiddleFactors aTFs(nN*2, mbInverse);
5330 aTFs.Compute();
5331 SCSIZE nNextPow2 = nN, nTmp = 0;
5332 lcl_roundUpNearestPow2(nNextPow2, nTmp);
5333
5334 if (nNextPow2 == nN)
5335 {
5336 ScComplexFFT2 aFFT2(aWorkArray, mbInverse, false /*disable polar*/, 0.0 /* no clipping */,
5337 aTFs, true /*subsample tf*/, true /*disable normalize*/);
5338 aFFT2.Compute();
5339 }
5340 else
5341 {
5342 ScComplexBluesteinFFT aFFT(aWorkArray, false /*complex input*/, mbInverse, false /*disable polar*/,
5343 0.0 /* no clipping */, true /*disable normalize*/);
5344 aFFT.Compute();
5345 }
5346
5347 // Post process aWorkArray to populate mrOutArray
5348
5349 const SCSIZE nTwoN = 2*nN, nThreeN = 3*nN;
5350 double fY1R, fY2R, fY1I, fY2I, fResR, fResI, fWR, fWI;
5351 for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
5352 {
5353 const SCSIZE nIdxRev = nIdx ? (nN - nIdx) : 0;
5354 fY1R = aWorkArray[nIdx];
5355 fY2R = aWorkArray[nIdxRev];
5356 fY1I = aWorkArray[nN + nIdx];
5357 fY2I = aWorkArray[nN + nIdxRev];
5358 fWR = aTFs.mfWReal[nIdx];
5359 fWI = aTFs.mfWImag[nIdx];
5360
5361 // mrOutArray has length = 4*nN
5362 // Real part of the final output (only half of the symmetry around Nyquist frequency)
5363 // Fills the first quarter.
5364 mrOutArray[nIdx] = fResR = 0.5*(
5365 fY1R + fY2R +
5366 fWR * (fY1I + fY2I) +
5367 fWI * (fY1R - fY2R) );
5368 // Imaginary part of the final output (only half of the symmetry around Nyquist frequency)
5369 // Fills the third quarter.
5370 mrOutArray[nTwoN + nIdx] = fResI = 0.5*(
5371 fY1I - fY2I +
5372 fWI * (fY1I + fY2I) -
5373 fWR * (fY1R - fY2R) );
5374
5375 // Fill the missing 2 quarters using symmetry argument.
5376 if (nIdx)
5377 {
5378 // Fills the 2nd quarter.
5379 mrOutArray[nN + nIdxRev] = fResR;
5380 // Fills the 4th quarter.
5381 mrOutArray[nThreeN + nIdxRev] = -fResI;
5382 }
5383 else
5384 {
5385 mrOutArray[nN] = fY1R - fY1I;
5386 mrOutArray[nThreeN] = 0.0;
5387 }
5388 }
5389
5390 // Normalize/Polar operations
5391 if (mbPolar)
5392 lcl_convertToPolar(mrOutArray, mfMinMag);
5393
5394 // Normalize after converting to polar, so we have a chance to
5395 // save O(mnPoints) flops.
5396 if (mbInverse)
5397 lcl_normalize(mrOutArray, mbPolar);
5398}
5399
5400using ScMatrixGenerator = ScMatrixRef(SCSIZE, SCSIZE, std::vector<double>&);
5401
5402namespace {
5403
5404// Generic FFT class that decides which FFT implementation to use.
5405class ScFFT
5406{
5407public:
5408
5409 ScFFT(ScMatrixRef& pMat, bool bReal, bool bInverse, bool bPolar, double fMinMag) :
5410 mpInputMat(pMat),
5411 mfMinMag(fMinMag),
5412 mbReal(bReal),
5413 mbInverse(bInverse),
5414 mbPolar(bPolar)
5415 {}
5416
5417 ScMatrixRef Compute(const std::function<ScMatrixGenerator>& rMatGenFunc);
5418
5419private:
5420 ScMatrixRef& mpInputMat;
5421 double mfMinMag;
5422 bool mbReal:1;
5423 bool mbInverse:1;
5424 bool mbPolar:1;
5425};
5426
5427}
5428
5429ScMatrixRef ScFFT::Compute(const std::function<ScMatrixGenerator>& rMatGenFunc)
5430{
5431 std::vector<double> aArray;
5432 mpInputMat->GetDoubleArray(aArray);
5433 SCSIZE nPoints = mbReal ? aArray.size() : (aArray.size()/2);
5434 if (nPoints == 1)
5435 {
5436 std::vector<double> aOutArray(2);
5437 aOutArray[0] = aArray[0];
5438 aOutArray[1] = mbReal ? 0.0 : aArray[1];
5439 if (mbPolar)
5440 lcl_convertToPolar(aOutArray, mfMinMag);
5441 return rMatGenFunc(2, 1, aOutArray);
5442 }
5443
5444 if (mbReal && (nPoints % 2) == 0)
5445 {
5446 std::vector<double> aOutArray(nPoints*2);
5447 ScRealFFT aFFT(aArray, aOutArray, mbInverse, mbPolar, mfMinMag);
5448 aFFT.Compute();
5449 return rMatGenFunc(2, nPoints, aOutArray);
5450 }
5451
5452 SCSIZE nNextPow2 = nPoints, nTmp = 0;
5453 lcl_roundUpNearestPow2(nNextPow2, nTmp);
5454 if (nNextPow2 == nPoints && !mbReal)
5455 {
5456 ScTwiddleFactors aTF(nPoints, mbInverse);
5457 aTF.Compute();
5458 ScComplexFFT2 aFFT2(aArray, mbInverse, mbPolar, mfMinMag, aTF);
5459 aFFT2.Compute();
5460 return rMatGenFunc(2, nPoints, aArray);
5461 }
5462
5463 if (mbReal)
5464 aArray.resize(nPoints*2, 0.0);
5465 ScComplexBluesteinFFT aFFT(aArray, mbReal, mbInverse, mbPolar, mfMinMag);
5466 aFFT.Compute();
5467 return rMatGenFunc(2, nPoints, aArray);
5468}
5469
5471{
5472 sal_uInt8 nParamCount = GetByte();
5473 if ( !MustHaveParamCount( nParamCount, 2, 5 ) )
5474 return;
5475
5476 bool bInverse = false;
5477 bool bPolar = false;
5478 double fMinMag = 0.0;
5479
5480 if (nParamCount == 5)
5481 {
5482 if (IsMissing())
5483 Pop();
5484 else
5485 fMinMag = GetDouble();
5486 }
5487
5488 if (nParamCount >= 4)
5489 {
5490 if (IsMissing())
5491 Pop();
5492 else