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 * M_SQRT1_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 {
3414 // Note: neg fDiff seen with forum-mso-en4-719754.xlsx with
3415 // fPercentile of near 1 where approxFloor gave nIndex of nSize-1
3416 // resulting in a non-zero tiny negative fDiff.
3417 return *iter;
3418 }
3419 else
3420 {
3421 OSL_ENSURE(nIndex < nSize-1, "GetPercentile: wrong index(2)");
3422 double fVal = *iter;
3423 iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
3424 return fVal + fDiff * (*iter - fVal);
3425 }
3426 }
3427}
3428
3429double ScInterpreter::GetPercentileExclusive( vector<double> & rArray, double fPercentile )
3430{
3431 size_t nSize1 = rArray.size() + 1;
3432 if ( rArray.empty() || nSize1 == 1 || nGlobalError != FormulaError::NONE)
3433 {
3434 SetError( FormulaError::NoValue );
3435 return 0.0;
3436 }
3437 if ( fPercentile * nSize1 < 1.0 || fPercentile * nSize1 > static_cast<double>( nSize1 - 1 ) )
3438 {
3439 SetError( FormulaError::IllegalParameter );
3440 return 0.0;
3441 }
3442
3443 size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * nSize1 - 1 ));
3444 double fDiff = fPercentile * nSize1 - 1 - ::rtl::math::approxFloor( fPercentile * nSize1 - 1 );
3445 OSL_ENSURE(nIndex < ( nSize1 - 1 ), "GetPercentile: wrong index(1)");
3446 vector<double>::iterator iter = rArray.begin() + nIndex;
3447 ::std::nth_element( rArray.begin(), iter, rArray.end());
3448 if (fDiff == 0.0)
3449 return *iter;
3450 else
3451 {
3452 OSL_ENSURE(nIndex < nSize1, "GetPercentile: wrong index(2)");
3453 double fVal = *iter;
3454 iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
3455 return fVal + fDiff * (*iter - fVal);
3456 }
3457}
3458
3459void ScInterpreter::ScPercentile( bool bInclusive )
3460{
3461 if ( !MustHaveParamCount( GetByte(), 2 ) )
3462 return;
3463 double alpha = GetDouble();
3464 if ( bInclusive ? ( alpha < 0.0 || alpha > 1.0 ) : ( alpha <= 0.0 || alpha >= 1.0 ) )
3465 {
3467 return;
3468 }
3469 vector<double> aArray;
3470 GetNumberSequenceArray( 1, aArray, false );
3471 if ( aArray.empty() || nGlobalError != FormulaError::NONE )
3472 {
3473 PushNoValue();
3474 return;
3475 }
3476 if ( bInclusive )
3477 PushDouble( GetPercentile( aArray, alpha ));
3478 else
3480}
3481
3482void ScInterpreter::ScQuartile( bool bInclusive )
3483{
3484 if ( !MustHaveParamCount( GetByte(), 2 ) )
3485 return;
3486 double fFlag = ::rtl::math::approxFloor(GetDouble());
3487 if ( bInclusive ? ( fFlag < 0.0 || fFlag > 4.0 ) : ( fFlag <= 0.0 || fFlag >= 4.0 ) )
3488 {
3490 return;
3491 }
3492 vector<double> aArray;
3493 GetNumberSequenceArray( 1, aArray, false );
3494 if ( aArray.empty() || nGlobalError != FormulaError::NONE )
3495 {
3496 PushNoValue();
3497 return;
3498 }
3499 if ( bInclusive )
3500 PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentile( aArray, 0.25 * fFlag ) );
3501 else
3502 PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentileExclusive( aArray, 0.25 * fFlag ) );
3503}
3504
3506{
3507 sal_uInt8 nParamCount = GetByte();
3508 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3509 return;
3510 vector<double> aSortArray;
3511 GetSortArray( nParamCount, aSortArray, nullptr, false, false );
3512 SCSIZE nSize = aSortArray.size();
3513 if (nSize == 0 || nGlobalError != FormulaError::NONE)
3514 PushNoValue();
3515 else
3516 {
3517 SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
3518 double nOldVal = aSortArray[0];
3519 SCSIZE i;
3520 for ( i = 1; i < nSize; i++)
3521 {
3522 if (aSortArray[i] == nOldVal)
3523 nCount++;
3524 else
3525 {
3526 nOldVal = aSortArray[i];
3527 if (nCount > nMax)
3528 {
3529 nMax = nCount;
3530 nMaxIndex = i-1;
3531 }
3532 nCount = 1;
3533 }
3534 }
3535 if (nCount > nMax)
3536 {
3537 nMax = nCount;
3538 nMaxIndex = i-1;
3539 }
3540 if (nMax == 1 && nCount == 1)
3541 PushNoValue();
3542 else if (nMax == 1)
3543 PushDouble(nOldVal);
3544 else
3545 PushDouble(aSortArray[nMaxIndex]);
3546 }
3547}
3548
3550{
3551 sal_uInt8 nParamCount = GetByte();
3552 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3553 return;
3554 vector<double> aArray;
3555 GetNumberSequenceArray( nParamCount, aArray, false );
3556 vector< double > aSortArray( aArray );
3557 QuickSort( aSortArray, nullptr );
3558 SCSIZE nSize = aSortArray.size();
3559 if ( nSize == 0 || nGlobalError != FormulaError::NONE )
3560 PushNoValue();
3561 else
3562 {
3563 SCSIZE nMax = 1, nCount = 1;
3564 double nOldVal = aSortArray[ 0 ];
3565 vector< double > aResultArray( 1 );
3566 SCSIZE i;
3567 for ( i = 1; i < nSize; i++ )
3568 {
3569 if ( aSortArray[ i ] == nOldVal )
3570 nCount++;
3571 else
3572 {
3573 if ( nCount >= nMax && nCount > 1 )
3574 {
3575 if ( nCount > nMax )
3576 {
3577 nMax = nCount;
3578 if ( aResultArray.size() != 1 )
3579 vector< double >( 1 ).swap( aResultArray );
3580 aResultArray[ 0 ] = nOldVal;
3581 }
3582 else
3583 aResultArray.emplace_back( nOldVal );
3584 }
3585 nOldVal = aSortArray[ i ];
3586 nCount = 1;
3587 }
3588 }
3589 if ( nCount >= nMax && nCount > 1 )
3590 {
3591 if ( nCount > nMax )
3592 vector< double >().swap( aResultArray );
3593 aResultArray.emplace_back( nOldVal );
3594 }
3595 if ( nMax == 1 && nCount == 1 )
3596 PushNoValue();
3597 else if ( nMax == 1 )
3598 PushDouble( nOldVal ); // there is only 1 result, no reordering needed
3599 else
3600 {
3601 // sort resultArray according to ordering of aArray
3602 vector< vector< double > > aOrder;
3603 aOrder.resize( aResultArray.size(), vector< double >( 2 ) );
3604 for ( i = 0; i < aResultArray.size(); i++ )
3605 {
3606 for ( SCSIZE j = 0; j < nSize; j++ )
3607 {
3608 if ( aArray[ j ] == aResultArray[ i ] )
3609 {
3610 aOrder[ i ][ 0 ] = aResultArray[ i ];
3611 aOrder[ i ][ 1 ] = j;
3612 break;
3613 }
3614 }
3615 }
3616 sort( aOrder.begin(), aOrder.end(), []( const std::vector< double >& lhs,
3617 const std::vector< double >& rhs )
3618 { return lhs[ 1 ] < rhs[ 1 ]; } );
3619
3620 if ( bSingle )
3621 PushDouble( aOrder[ 0 ][ 0 ] );
3622 else
3623 {
3624 // put result in correct order in aResultArray
3625 for ( i = 0; i < aResultArray.size(); i++ )
3626 aResultArray[ i ] = aOrder[ i ][ 0 ];
3627 ScMatrixRef pResMatrix = GetNewMat( 1, aResultArray.size(), true );
3628 pResMatrix->PutDoubleVector( aResultArray, 0, 0 );
3629 PushMatrix( pResMatrix );
3630 }
3631 }
3632 }
3633}
3634
3636{
3637 if ( !MustHaveParamCount( GetByte(), 2 ) )
3638 return;
3639
3640 SCSIZE nCol = 0, nRow = 0;
3641 const auto aArray = GetRankNumberArray(nCol, nRow);
3642 const size_t nRankArraySize = aArray.size();
3643 if (nRankArraySize == 0 || nGlobalError != FormulaError::NONE)
3644 {
3645 PushNoValue();
3646 return;
3647 }
3648 assert(nRankArraySize == nCol * nRow);
3649
3650 std::vector<SCSIZE> aRankArray;
3651 aRankArray.reserve(nRankArraySize);
3652 std::transform(aArray.begin(), aArray.end(), std::back_inserter(aRankArray),
3653 [bSmall](double f) {
3654 f = (bSmall ? rtl::math::approxFloor(f) : rtl::math::approxCeil(f));
3655 // Valid ranks are >= 1.
3656 if (f < 1.0 || !o3tl::convertsToAtMost(f, std::numeric_limits<SCSIZE>::max()))
3657 return static_cast<SCSIZE>(0);
3658 return static_cast<SCSIZE>(f);
3659 });
3660
3661 vector<double> aSortArray;
3662 GetNumberSequenceArray(1, aSortArray, false );
3663 const SCSIZE nSize = aSortArray.size();
3664 if (nSize == 0 || nGlobalError != FormulaError::NONE)
3665 PushNoValue();
3666 else if (nRankArraySize == 1)
3667 {
3668 const SCSIZE k = aRankArray[0];
3669 if (k < 1 || nSize < k)
3670 {
3671 if (!std::isfinite(aArray[0]))
3672 PushDouble(aArray[0]); // propagates error
3673 else
3674 PushNoValue();
3675 }
3676 else
3677 {
3678 vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
3679 ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
3680 PushDouble( *iPos);
3681 }
3682 }
3683 else
3684 {
3685 std::set<SCSIZE> aIndices;
3686 for (SCSIZE n : aRankArray)
3687 {
3688 if (1 <= n && n <= nSize)
3689 aIndices.insert(bSmall ? n-1 : nSize-n);
3690 }
3691 // We can spare sorting when the total number of ranks is small enough.
3692 // Find only the elements at given indices if, arbitrarily, the index size is
3693 // smaller than 1/3 of the haystack array's size; just sort it squarely, otherwise.
3694 if (aIndices.size() < nSize/3)
3695 {
3696 auto itBegin = aSortArray.begin();
3697 for (SCSIZE i : aIndices)
3698 {
3699 auto it = aSortArray.begin() + i;
3700 std::nth_element(itBegin, it, aSortArray.end());
3701 itBegin = ++it;
3702 }
3703 }
3704 else
3705 std::sort(aSortArray.begin(), aSortArray.end());
3706
3707 std::vector<double> aResultArray;
3708 aResultArray.reserve(nRankArraySize);
3709 for (size_t i = 0; i < nRankArraySize; ++i)
3710 {
3711 const SCSIZE n = aRankArray[i];
3712 if (1 <= n && n <= nSize)
3713 aResultArray.push_back( aSortArray[bSmall ? n-1 : nSize-n]);
3714 else if (!std::isfinite( aArray[i]))
3715 aResultArray.push_back( aArray[i]); // propagate error
3716 else
3717 aResultArray.push_back( CreateDoubleError( FormulaError::IllegalArgument));
3718 }
3719 ScMatrixRef pResult = GetNewMat(nCol, nRow, aResultArray);
3720 PushMatrix(pResult);
3721 }
3722}
3723
3725{
3726 CalculateSmallLarge(false);
3727}
3728
3730{
3731 CalculateSmallLarge(true);
3732}
3733
3734void ScInterpreter::ScPercentrank( bool bInclusive )
3735{
3736 sal_uInt8 nParamCount = GetByte();
3737 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
3738 return;
3739 double fSignificance = ( nParamCount == 3 ? ::rtl::math::approxFloor( GetDouble() ) : 3.0 );
3740 if ( fSignificance < 1.0 )
3741 {
3743 return;
3744 }
3745 double fNum = GetDouble();
3746 vector<double> aSortArray;
3747 GetSortArray( 1, aSortArray, nullptr, false, false );
3748 SCSIZE nSize = aSortArray.size();
3749 if ( nSize == 0 || nGlobalError != FormulaError::NONE )
3750 PushNoValue();
3751 else
3752 {
3753 if ( fNum < aSortArray[ 0 ] || fNum > aSortArray[ nSize - 1 ] )
3754 PushNoValue();
3755 else
3756 {
3757 double fRes;
3758 if ( nSize == 1 )
3759 fRes = 1.0; // fNum == aSortArray[ 0 ], see test above
3760 else
3761 fRes = GetPercentrank( aSortArray, fNum, bInclusive );
3762 if ( fRes != 0.0 )
3763 {
3764 double fExp = ::rtl::math::approxFloor( log10( fRes ) ) + 1.0 - fSignificance;
3765 fRes = ::rtl::math::round( fRes * pow( 10, -fExp ) ) / pow( 10, -fExp );
3766 }
3767 PushDouble( fRes );
3768 }
3769 }
3770}
3771
3772double ScInterpreter::GetPercentrank( ::std::vector<double> & rArray, double fVal, bool bInclusive )
3773{
3774 SCSIZE nSize = rArray.size();
3775 double fRes;
3776 if ( fVal == rArray[ 0 ] )
3777 {
3778 if ( bInclusive )
3779 fRes = 0.0;
3780 else
3781 fRes = 1.0 / static_cast<double>( nSize + 1 );
3782 }
3783 else
3784 {
3785 SCSIZE nOldCount = 0;
3786 double fOldVal = rArray[ 0 ];
3787 SCSIZE i;
3788 for ( i = 1; i < nSize && rArray[ i ] < fVal; i++ )
3789 {
3790 if ( rArray[ i ] != fOldVal )
3791 {
3792 nOldCount = i;
3793 fOldVal = rArray[ i ];
3794 }
3795 }
3796 if ( rArray[ i ] != fOldVal )
3797 nOldCount = i;
3798 if ( fVal == rArray[ i ] )
3799 {
3800 if ( bInclusive )
3801 fRes = div( nOldCount, nSize - 1 );
3802 else
3803 fRes = static_cast<double>( i + 1 ) / static_cast<double>( nSize + 1 );
3804 }
3805 else
3806 {
3807 // nOldCount is the count of smaller entries
3808 // fVal is between rArray[ nOldCount - 1 ] and rArray[ nOldCount ]
3809 // use linear interpolation to find a position between the entries
3810 if ( nOldCount == 0 )
3811 {
3812 OSL_FAIL( "should not happen" );
3813 fRes = 0.0;
3814 }
3815 else
3816 {
3817 double fFract = ( fVal - rArray[ nOldCount - 1 ] ) /
3818 ( rArray[ nOldCount ] - rArray[ nOldCount - 1 ] );
3819 if ( bInclusive )
3820 fRes = div( static_cast<double>( nOldCount - 1 ) + fFract, nSize - 1 );
3821 else
3822 fRes = ( static_cast<double>(nOldCount) + fFract ) / static_cast<double>( nSize + 1 );
3823 }
3824 }
3825 }
3826 return fRes;
3827}
3828
3830{
3831 if ( !MustHaveParamCount( GetByte(), 2 ) )
3832 return;
3833 double alpha = GetDouble();
3834 if (alpha < 0.0 || alpha >= 1.0)
3835 {
3837 return;
3838 }
3839 vector<double> aSortArray;
3840 GetSortArray( 1, aSortArray, nullptr, false, false );
3841 SCSIZE nSize = aSortArray.size();
3842 if (nSize == 0 || nGlobalError != FormulaError::NONE)
3843 PushNoValue();
3844 else
3845 {
3846 sal_uLong nIndex = static_cast<sal_uLong>(::rtl::math::approxFloor(alpha*static_cast<double>(nSize)));
3847 if (nIndex % 2 != 0)
3848 nIndex--;
3849 nIndex /= 2;
3850 OSL_ENSURE(nIndex < nSize, "ScTrimMean: wrong index");
3851 KahanSum fSum = 0.0;
3852 for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
3853 fSum += aSortArray[i];
3854 PushDouble(fSum.get()/static_cast<double>(nSize-2*nIndex));
3855 }
3856}
3857
3858std::vector<double> ScInterpreter::GetRankNumberArray( SCSIZE& rCol, SCSIZE& rRow )
3859{
3860 std::vector<double> aArray;
3861 switch (GetStackType())
3862 {
3863 case svDouble:
3864 aArray.push_back(PopDouble());
3865 rCol = rRow = 1;
3866 break;
3867 case svSingleRef:
3868 {
3869 ScAddress aAdr;
3870 PopSingleRef(aAdr);
3871 ScRefCellValue aCell(mrDoc, aAdr);
3872 if (aCell.hasNumeric())
3873 {
3874 aArray.push_back(GetCellValue(aAdr, aCell));
3875 rCol = rRow = 1;
3876 }
3877 }
3878 break;
3879 case svDoubleRef:
3880 {
3881 ScRange aRange;
3882 PopDoubleRef(aRange, true);
3883 if (nGlobalError != FormulaError::NONE)
3884 break;
3885
3886 // give up unless the start and end are in the same sheet
3887 if (aRange.aStart.Tab() != aRange.aEnd.Tab())
3888 {
3889 SetError(FormulaError::IllegalParameter);
3890 break;
3891 }
3892
3893 // the range already is in order
3894 assert(aRange.aStart.Col() <= aRange.aEnd.Col());
3895 assert(aRange.aStart.Row() <= aRange.aEnd.Row());
3896 rCol = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3897 rRow = aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3898 aArray.reserve(rCol * rRow);
3899
3900 FormulaError nErr = FormulaError::NONE;
3901 double fCellVal;
3902 ScValueIterator aValIter(mrContext, aRange, mnSubTotalFlags);
3903 if (aValIter.GetFirst(fCellVal, nErr))
3904 {
3905 do
3906 aArray.push_back(fCellVal);
3907 while (aValIter.GetNext(fCellVal, nErr) && nErr == FormulaError::NONE);
3908 }
3909 // Note that SMALL() and LARGE() rank parameters (2nd) have
3910 // ParamClass::Value, so in array mode this is never hit and
3911 // argument was converted to matrix instead, but for normal
3912 // evaluation any non-numeric value including empty cell will
3913 // result in error anyway, so just clear and propagate an existing
3914 // error here already.
3915 if (aArray.size() != rCol * rRow)
3916 {
3917 aArray.clear();
3918 SetError(nErr);
3919 }
3920 }
3921 break;
3922 case svMatrix:
3925 {
3926 ScMatrixRef pMat = GetMatrix();
3927 if (!pMat)
3928 break;
3929
3930 const SCSIZE nCount = pMat->GetElementCount();
3931 aArray.reserve(nCount);
3932 // Do not propagate errors from matrix elements as global error.
3933 pMat->SetErrorInterpreter(nullptr);
3934 if (pMat->IsNumeric())
3935 {
3936 for (SCSIZE i = 0; i < nCount; ++i)
3937 aArray.push_back(pMat->GetDouble(i));
3938 }
3939 else
3940 {
3941 for (SCSIZE i = 0; i < nCount; ++i)
3942 {
3943 if (pMat->IsValue(i))
3944 aArray.push_back( pMat->GetDouble(i));
3945 else
3946 aArray.push_back( CreateDoubleError( FormulaError::NoValue));
3947 }
3948 }
3949 pMat->GetDimensions(rCol, rRow);
3950 }
3951 break;
3952 default:
3953 PopError();
3954 SetError(FormulaError::IllegalParameter);
3955 break;
3956 }
3957 return aArray;
3958}
3959
3960void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray, bool bConvertTextInArray )
3961{
3962 ScAddress aAdr;
3963 ScRange aRange;
3964 const bool bIgnoreErrVal = bool(mnSubTotalFlags & SubtotalFlags::IgnoreErrVal);
3965 short nParam = nParamCount;
3966 size_t nRefInList = 0;
3967 ReverseStack( nParamCount );
3968 while (nParam-- > 0)
3969 {
3970 const StackVar eStackType = GetStackType();
3971 switch (eStackType)
3972 {
3973 case svDouble :
3974 rArray.push_back( PopDouble());
3975 break;
3976 case svSingleRef :
3977 {
3978 PopSingleRef( aAdr );
3979 ScRefCellValue aCell(mrDoc, aAdr);
3980 if (bIgnoreErrVal && aCell.hasError())
3981 ; // nothing
3982 else if (aCell.hasNumeric())
3983 rArray.push_back(GetCellValue(aAdr, aCell));
3984 }
3985 break;
3986 case svDoubleRef :
3987 case svRefList :
3988 {
3989 PopDoubleRef( aRange, nParam, nRefInList);
3990 if (nGlobalError != FormulaError::NONE)
3991 break;
3992
3993 aRange.PutInOrder();
3994 SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3995 nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3996 rArray.reserve( rArray.size() + nCellCount);
3997
3998 FormulaError nErr = FormulaError::NONE;
3999 double fCellVal;
4000 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
4001 if (aValIter.GetFirst( fCellVal, nErr))
4002 {
4003 if (bIgnoreErrVal)
4004 {
4005 if (nErr == FormulaError::NONE)
4006 rArray.push_back( fCellVal);
4007 while (aValIter.GetNext( fCellVal, nErr))
4008 {
4009 if (nErr == FormulaError::NONE)
4010 rArray.push_back( fCellVal);
4011 }
4012 }
4013 else
4014 {
4015 rArray.push_back( fCellVal);
4016 SetError(nErr);
4017 while ((nErr == FormulaError::NONE) && aValIter.GetNext( fCellVal, nErr))
4018 rArray.push_back( fCellVal);
4019 SetError(nErr);
4020 }
4021 }
4022 }
4023 break;
4024 case svMatrix :
4027 {
4028 ScMatrixRef pMat = GetMatrix();
4029 if (!pMat)
4030 break;
4031
4032 SCSIZE nCount = pMat->GetElementCount();
4033 rArray.reserve( rArray.size() + nCount);
4034 if (pMat->IsNumeric())
4035 {
4036 if (bIgnoreErrVal)
4037 {
4038 for (SCSIZE i = 0; i < nCount; ++i)
4039 {
4040 const double fVal = pMat->GetDouble(i);
4041 if (nGlobalError == FormulaError::NONE)
4042 rArray.push_back( fVal);
4043 else
4044 nGlobalError = FormulaError::NONE;
4045 }
4046 }
4047 else
4048 {
4049 for (SCSIZE i = 0; i < nCount; ++i)
4050 rArray.push_back( pMat->GetDouble(i));
4051 }
4052 }
4053 else if (bConvertTextInArray && eStackType == svMatrix)
4054 {
4055 for (SCSIZE i = 0; i < nCount; ++i)
4056 {
4057 if ( pMat->IsValue( i ) )
4058 {
4059 if (bIgnoreErrVal)
4060 {
4061 const double fVal = pMat->GetDouble(i);
4062 if (nGlobalError == FormulaError::NONE)
4063 rArray.push_back( fVal);
4064 else
4065 nGlobalError = FormulaError::NONE;
4066 }
4067 else
4068 rArray.push_back( pMat->GetDouble(i));
4069 }
4070 else
4071 {
4072 // tdf#88547 try to convert string to (date)value
4073 OUString aStr = pMat->GetString( i ).getString();
4074 if ( aStr.getLength() > 0 )
4075 {
4077 nGlobalError = FormulaError::NONE;
4078 double fVal = ConvertStringToValue( aStr );
4079 if ( nGlobalError == FormulaError::NONE )
4080 {
4081 rArray.push_back( fVal );
4082 nGlobalError = nErr;
4083 }
4084 else
4085 {
4086 if (!bIgnoreErrVal)
4087 rArray.push_back( CreateDoubleError( FormulaError::NoValue));
4088 // Propagate previous error if any, else
4089 // the current #VALUE! error, unless
4090 // ignoring error values.
4091 if (nErr != FormulaError::NONE)
4092 nGlobalError = nErr;
4093 else if (!bIgnoreErrVal)
4094 nGlobalError = FormulaError::NoValue;
4095 else
4096 nGlobalError = FormulaError::NONE;
4097 }
4098 }
4099 }
4100 }
4101 }
4102 else
4103 {
4104 if (bIgnoreErrVal)
4105 {
4106 for (SCSIZE i = 0; i < nCount; ++i)
4107 {
4108 if (pMat->IsValue(i))
4109 {
4110 const double fVal = pMat->GetDouble(i);
4111 if (nGlobalError == FormulaError::NONE)
4112 rArray.push_back( fVal);
4113 else
4114 nGlobalError = FormulaError::NONE;
4115 }
4116 }
4117 }
4118 else
4119 {
4120 for (SCSIZE i = 0; i < nCount; ++i)
4121 {
4122 if (pMat->IsValue(i))
4123 rArray.push_back( pMat->GetDouble(i));
4124 }
4125 }
4126 }
4127 }
4128 break;
4129 default :
4130 PopError();
4131 SetError( FormulaError::IllegalParameter);
4132 break;
4133 }
4134 if (nGlobalError != FormulaError::NONE)
4135 break; // while
4136 }
4137 // nParam > 0 in case of error, clean stack environment and obtain earlier
4138 // error if there was one.
4139 while (nParam-- > 0)
4140 PopError();
4141}
4142
4143void ScInterpreter::GetSortArray( sal_uInt8 nParamCount, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder, bool bConvertTextInArray, bool bAllowEmptyArray )
4144{
4145 GetNumberSequenceArray( nParamCount, rSortArray, bConvertTextInArray );
4146 if (rSortArray.size() > MAX_COUNT_DOUBLE_FOR_SORT(mrDoc.GetSheetLimits()))
4147 SetError( FormulaError::MatrixSize);
4148 else if ( rSortArray.empty() )
4149 {
4150 if ( bAllowEmptyArray )
4151 return;
4152 SetError( FormulaError::NoValue);
4153 }
4154
4155 if (nGlobalError == FormulaError::NONE)
4156 QuickSort( rSortArray, pIndexOrder);
4157}
4158
4159static void lcl_QuickSort( tools::Long nLo, tools::Long nHi, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
4160{
4161 // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().
4162
4163 using ::std::swap;
4164
4165 if (nHi - nLo == 1)
4166 {
4167 if (rSortArray[nLo] > rSortArray[nHi])
4168 {
4169 swap(rSortArray[nLo], rSortArray[nHi]);
4170 if (pIndexOrder)
4171 swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
4172 }
4173 return;
4174 }
4175
4176 tools::Long ni = nLo;
4177 tools::Long nj = nHi;
4178 do
4179 {
4180 double fLo = rSortArray[nLo];
4181 while (ni <= nHi && rSortArray[ni] < fLo) ni++;
4182 while (nj >= nLo && fLo < rSortArray[nj]) nj--;
4183 if (ni <= nj)
4184 {
4185 if (ni != nj)
4186 {
4187 swap(rSortArray[ni], rSortArray[nj]);
4188 if (pIndexOrder)
4189 swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
4190 }
4191
4192 ++ni;
4193 --nj;
4194 }
4195 }
4196 while (ni < nj);
4197
4198 if ((nj - nLo) < (nHi - ni))
4199 {
4200 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
4201 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
4202 }
4203 else
4204 {
4205 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
4206 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
4207 }
4208}
4209
4210void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
4211{
4212 tools::Long n = static_cast<tools::Long>(rSortArray.size());
4213
4214 if (pIndexOrder)
4215 {
4216 pIndexOrder->clear();
4217 pIndexOrder->reserve(n);
4218 for (tools::Long i = 0; i < n; ++i)
4219 pIndexOrder->push_back(i);
4220 }
4221
4222 if (n < 2)
4223 return;
4224
4225 size_t nValCount = rSortArray.size();
4226 for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
4227 {
4228 size_t nInd = comphelper::rng::uniform_size_distribution(0, nValCount-2);
4229 ::std::swap( rSortArray[i], rSortArray[nInd]);
4230 if (pIndexOrder)
4231 ::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
4232 }
4233
4234 lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
4235}
4236
4237void ScInterpreter::ScRank( bool bAverage )
4238{
4239 sal_uInt8 nParamCount = GetByte();
4240 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
4241 return;
4242 bool bAscending;
4243 if ( nParamCount == 3 )
4244 bAscending = GetBool();
4245 else
4246 bAscending = false;
4247
4248 vector<double> aSortArray;
4249 GetSortArray( 1, aSortArray, nullptr, false, false );
4250 double fVal = GetDouble();
4251 SCSIZE nSize = aSortArray.size();
4252 if ( nSize == 0 || nGlobalError != FormulaError::NONE )
4253 PushNoValue();
4254 else
4255 {
4256 if ( fVal < aSortArray[ 0 ] || fVal > aSortArray[ nSize - 1 ] )
4257 PushError( FormulaError::NotAvailable);
4258 else
4259 {
4260 double fLastPos = 0;
4261 double fFirstPos = -1.0;
4262 bool bFinished = false;
4263 SCSIZE i;
4264 for (i = 0; i < nSize && !bFinished; i++)
4265 {
4266 if ( aSortArray[ i ] == fVal )
4267 {
4268 if ( fFirstPos < 0 )
4269 fFirstPos = i + 1.0;
4270 }
4271 else
4272 {
4273 if ( aSortArray[ i ] > fVal )
4274 {
4275 fLastPos = i;
4276 bFinished = true;
4277 }
4278 }
4279 }
4280 if ( !bFinished )
4281 fLastPos = i;
4282 if (fFirstPos <= 0)
4283 {
4284 PushError( FormulaError::NotAvailable);
4285 }
4286 else if ( !bAverage )
4287 {
4288 if ( bAscending )
4289 PushDouble( fFirstPos );
4290 else
4291 PushDouble( nSize + 1.0 - fLastPos );
4292 }
4293 else
4294 {
4295 if ( bAscending )
4296 PushDouble( ( fFirstPos + fLastPos ) / 2.0 );
4297 else
4298 PushDouble( nSize + 1.0 - ( fFirstPos + fLastPos ) / 2.0 );
4299 }
4300 }
4301 }
4302}
4303
4305{
4306 sal_uInt8 nParamCount = GetByte();
4307 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
4308 return;
4309 sal_uInt16 SaveSP = sp;
4310 double nMiddle = 0.0;
4311 KahanSum rVal = 0.0;
4312 double rValCount = 0.0;
4313 ScAddress aAdr;
4314 ScRange aRange;
4315 short nParam = nParamCount;
4316 size_t nRefInList = 0;
4317 while (nParam-- > 0)
4318 {
4319 switch (GetStackType())
4320 {
4321 case svDouble :
4322 rVal += GetDouble();
4323 rValCount++;
4324 break;
4325 case svSingleRef :
4326 {
4327 PopSingleRef( aAdr );
4328 ScRefCellValue aCell(mrDoc, aAdr);
4329 if (aCell.hasNumeric())
4330 {
4331 rVal += GetCellValue(aAdr, aCell);
4332 rValCount++;
4333 }
4334 }
4335 break;
4336 case svDoubleRef :
4337 case svRefList :
4338 {
4339 FormulaError nErr = FormulaError::NONE;
4340 double nCellVal;
4341 PopDoubleRef( aRange, nParam, nRefInList);
4342 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
4343 if (aValIter.GetFirst(nCellVal, nErr))
4344 {
4345 rVal += nCellVal;
4346 rValCount++;
4347 SetError(nErr);
4348 while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
4349 {
4350 rVal += nCellVal;
4351 rValCount++;
4352 }
4353 SetError(nErr);
4354 }
4355 }
4356 break;
4357 case svMatrix :
4360 {
4361 ScMatrixRef pMat = GetMatrix();
4362 if (pMat)
4363 {
4364 SCSIZE nCount = pMat->GetElementCount();
4365 if (pMat->IsNumeric())
4366 {
4367 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4368 {
4369 rVal += pMat->GetDouble(nElem);
4370 rValCount++;
4371 }
4372 }
4373 else
4374 {
4375 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4376 if (!pMat->IsStringOrEmpty(nElem))
4377 {
4378 rVal += pMat->GetDouble(nElem);
4379 rValCount++;
4380 }
4381 }
4382 }
4383 }
4384 break;
4385 default :
4386 SetError(FormulaError::IllegalParameter);
4387 break;
4388 }
4389 }
4390 if (nGlobalError != FormulaError::NONE)
4391 {
4393 return;
4394 }
4395 nMiddle = rVal.get() / rValCount;
4396 sp = SaveSP;
4397 rVal = 0.0;
4398 nParam = nParamCount;
4399 nRefInList = 0;
4400 while (nParam-- > 0)
4401 {
4402 switch (GetStackType())
4403 {
4404 case svDouble :
4405 rVal += std::abs(GetDouble() - nMiddle);
4406 break;
4407 case svSingleRef :
4408 {
4409 PopSingleRef( aAdr );
4410 ScRefCellValue aCell(mrDoc, aAdr);
4411 if (aCell.hasNumeric())
4412 rVal += std::abs(GetCellValue(aAdr, aCell) - nMiddle);
4413 }
4414 break;
4415 case svDoubleRef :
4416 case svRefList :
4417 {
4418 FormulaError nErr = FormulaError::NONE;
4419 double nCellVal;
4420 PopDoubleRef( aRange, nParam, nRefInList);
4421 ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
4422 if (aValIter.GetFirst(nCellVal, nErr))
4423 {
4424 rVal += std::abs(nCellVal - nMiddle);
4425 while (aValIter.GetNext(nCellVal, nErr))
4426 rVal += std::abs(nCellVal - nMiddle);
4427 }
4428 }
4429 break;
4430 case svMatrix :
4433 {
4434 ScMatrixRef pMat = GetMatrix();
4435 if (pMat)
4436 {
4437 SCSIZE nCount = pMat->GetElementCount();
4438 if (pMat->IsNumeric())
4439 {
4440 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4441 {
4442 rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
4443 }
4444 }
4445 else
4446 {
4447 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4448 {
4449 if (!pMat->IsStringOrEmpty(nElem))
4450 rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
4451 }
4452 }
4453 }
4454 }
4455 break;
4456 default : SetError(FormulaError::IllegalParameter); break;
4457 }
4458 }
4459 PushDouble(rVal.get() / rValCount);
4460}
4461
4463{
4464 auto VarResult = []( double fVal, size_t /*nValCount*/ )
4465 {
4466 return fVal;
4467 };
4468 GetStVarParams( false /*bTextAsZero*/, VarResult);
4469}
4470
4472{
4473 sal_uInt8 nParamCount = GetByte();
4474 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
4475 return;
4476 double fUp, fLo;
4477 fUp = GetDouble();
4478 if (nParamCount == 4)
4479 fLo = GetDouble();
4480 else
4481 fLo = fUp;
4482 if (fLo > fUp)
4483 std::swap( fLo, fUp );
4484 ScMatrixRef pMatP = GetMatrix();
4485 ScMatrixRef pMatW = GetMatrix();
4486 if (!pMatP || !pMatW)
4488 else
4489 {
4490 SCSIZE nC1, nC2;
4491 SCSIZE nR1, nR2;
4492 pMatP->GetDimensions(nC1, nR1);
4493 pMatW->GetDimensions(nC2, nR2);
4494 if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
4495 nC2 == 0 || nR2 == 0)
4496 PushNA();
4497 else
4498 {
4499 KahanSum fSum = 0.0;
4500 KahanSum fRes = 0.0;
4501 bool bStop = false;
4502 double fP, fW;
4503 for ( SCSIZE i = 0; i < nC1 && !bStop; i++ )
4504 {
4505 for (SCSIZE j = 0; j < nR1 && !bStop; ++j )
4506 {
4507 if (pMatP->IsValue(i,j) && pMatW->IsValue(i,j))
4508 {
4509 fP = pMatP->GetDouble(i,j);
4510 fW = pMatW->GetDouble(i,j);
4511 if (fP < 0.0 || fP > 1.0)
4512 bStop = true;
4513 else
4514 {
4515 fSum += fP;
4516 if (fW >= fLo && fW <= fUp)
4517 fRes += fP;
4518 }
4519 }
4520 else
4521 SetError( FormulaError::IllegalArgument);
4522 }
4523 }
4524 if (bStop || std::abs((fSum -1.0).get()) > 1.0E-7)
4525 PushNoValue();
4526 else
4527 PushDouble(fRes.get());
4528 }
4529 }
4530}
4531
4533{
4534 // This is identical to ScPearson()
4535 ScPearson();
4536}
4537
4539{
4540 CalculatePearsonCovar( false, false, false );
4541}
4542
4544{
4545 CalculatePearsonCovar( false, false, true );
4546}
4547
4549{
4550 CalculatePearsonCovar( true, false, false );
4551}
4552
4553void ScInterpreter::CalculatePearsonCovar( bool _bPearson, bool _bStexy, bool _bSample )
4554{
4555 if ( !MustHaveParamCount( GetByte(), 2 ) )
4556 return;
4557 ScMatrixRef pMat1 = GetMatrix();
4558 ScMatrixRef pMat2 = GetMatrix();
4559 if (!pMat1 || !pMat2)
4560 {
4562 return;
4563 }
4564 SCSIZE nC1, nC2;
4565 SCSIZE nR1, nR2;
4566 pMat1->GetDimensions(nC1, nR1);
4567 pMat2->GetDimensions(nC2, nR2);
4568 if (nR1 != nR2 || nC1 != nC2)
4569 {
4571 return;
4572 }
4573 /* #i78250#
4574 * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
4575 * but the latter produces wrong results if the absolute values are high,
4576 * for example above 10^8
4577 */
4578 double fCount = 0.0;
4579 KahanSum fSumX = 0.0;
4580 KahanSum fSumY = 0.0;
4581
4582 for (SCSIZE i = 0; i < nC1; i++)
4583 {
4584 for (SCSIZE j = 0; j < nR1; j++)
4585 {
4586 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4587 {
4588 fSumX += pMat1->GetDouble(i,j);
4589 fSumY += pMat2->GetDouble(i,j);
4590 fCount++;
4591 }
4592 }
4593 }
4594 if (fCount < (_bStexy ? 3.0 : (_bSample ? 2.0 : 1.0)))
4595 PushNoValue();
4596 else
4597 {
4598 KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4599 KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4600 KahanSum fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2
4601 const double fMeanX = fSumX.get() / fCount;
4602 const double fMeanY = fSumY.get() / fCount;
4603 for (SCSIZE i = 0; i < nC1; i++)
4604 {
4605 for (SCSIZE j = 0; j < nR1; j++)
4606 {
4607 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4608 {
4609 const double fValX = pMat1->GetDouble(i,j);
4610 const double fValY = pMat2->GetDouble(i,j);
4611 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4612 if ( _bPearson )
4613 {
4614 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4615 fSumSqrDeltaY += (fValY - fMeanY) * (fValY - fMeanY);
4616 }
4617 }
4618 }
4619 }
4620 if ( _bPearson )
4621 {
4622 // tdf#94962 - Values below the numerical limit lead to unexpected results
4623 if (fSumSqrDeltaX < ::std::numeric_limits<double>::min()
4624 || (!_bStexy && fSumSqrDeltaY < ::std::numeric_limits<double>::min()))
4625 PushError( FormulaError::DivisionByZero);
4626 else if ( _bStexy )
4627 PushDouble( sqrt( ( fSumSqrDeltaY - fSumDeltaXDeltaY *
4628 fSumDeltaXDeltaY / fSumSqrDeltaX ).get() / (fCount-2)));
4629 else
4630 PushDouble( fSumDeltaXDeltaY.get() / sqrt( fSumSqrDeltaX.get() * fSumSqrDeltaY.get() ));
4631 }
4632 else
4633 {
4634 if ( _bSample )
4635 PushDouble( fSumDeltaXDeltaY.get() / ( fCount - 1 ) );
4636 else
4637 PushDouble( fSumDeltaXDeltaY.get() / fCount);
4638 }
4639 }
4640}
4641
4643{
4644 // Same as ScPearson()*ScPearson()
4645 ScPearson();
4646 if (nGlobalError != FormulaError::NONE)
4647 return;
4648
4649 switch (GetStackType())
4650 {
4651 case svDouble:
4652 {
4653 double fVal = PopDouble();
4654 PushDouble( fVal * fVal);
4655 }
4656 break;
4657 default:
4658 PopError();
4659 PushNoValue();
4660 }
4661}
4662
4664{
4665 CalculatePearsonCovar( true, true, false );
4666}
4668{
4669 if ( !MustHaveParamCount( GetByte(), 2 ) )
4670 return;
4671 ScMatrixRef pMat1 = GetMatrix();
4672 ScMatrixRef pMat2 = GetMatrix();
4673 if (!pMat1 || !pMat2)
4674 {
4676 return;
4677 }
4678 SCSIZE nC1, nC2;
4679 SCSIZE nR1, nR2;
4680 pMat1->GetDimensions(nC1, nR1);
4681 pMat2->GetDimensions(nC2, nR2);
4682 if (nR1 != nR2 || nC1 != nC2)
4683 {
4685 return;
4686 }
4687 // #i78250# numerical stability improved
4688 double fCount = 0.0;
4689 KahanSum fSumX = 0.0;
4690 KahanSum fSumY = 0.0;
4691
4692 for (SCSIZE i = 0; i < nC1; i++)
4693 {
4694 for (SCSIZE j = 0; j < nR1; j++)
4695 {
4696 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4697 {
4698 fSumX += pMat1->GetDouble(i,j);
4699 fSumY += pMat2->GetDouble(i,j);
4700 fCount++;
4701 }
4702 }
4703 }
4704 if (fCount < 1.0)
4705 PushNoValue();
4706 else
4707 {
4708 KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4709 KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4710 double fMeanX = fSumX.get() / fCount;
4711 double fMeanY = fSumY.get() / fCount;
4712 for (SCSIZE i = 0; i < nC1; i++)
4713 {
4714 for (SCSIZE j = 0; j < nR1; j++)
4715 {
4716 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4717 {
4718 double fValX = pMat1->GetDouble(i,j);
4719 double fValY = pMat2->GetDouble(i,j);
4720 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4721 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4722 }
4723 }
4724 }
4725 if (fSumSqrDeltaX == 0.0)
4726 PushError( FormulaError::DivisionByZero);
4727 else
4728 {
4729 if ( bSlope )
4730 PushDouble( fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get());
4731 else
4732 PushDouble( fMeanY - fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * fMeanX);
4733 }
4734 }
4735}
4736
4738{
4740}
4741
4743{
4745}
4746
4748{
4749 if ( !MustHaveParamCount( GetByte(), 3 ) )
4750 return;
4751 ScMatrixRef pMat1 = GetMatrix();
4752 ScMatrixRef pMat2 = GetMatrix();
4753 if (!pMat1 || !pMat2)
4754 {
4756 return;
4757 }
4758 SCSIZE nC1, nC2;
4759 SCSIZE nR1, nR2;
4760 pMat1->GetDimensions(nC1, nR1);
4761 pMat2->GetDimensions(nC2, nR2);
4762 if (nR1 != nR2 || nC1 != nC2)
4763 {
4765 return;
4766 }
4767 double fVal = GetDouble();
4768 // #i78250# numerical stability improved
4769 double fCount = 0.0;
4770 KahanSum fSumX = 0.0;
4771 KahanSum fSumY = 0.0;
4772
4773 for (SCSIZE i = 0; i < nC1; i++)
4774 {
4775 for (SCSIZE j = 0; j < nR1; j++)
4776 {
4777 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4778 {
4779 fSumX += pMat1->GetDouble(i,j);
4780 fSumY += pMat2->GetDouble(i,j);
4781 fCount++;
4782 }
4783 }
4784 }
4785 if (fCount < 1.0)
4786 PushNoValue();
4787 else
4788 {
4789 KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4790 KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4791 double fMeanX = fSumX.get() / fCount;
4792 double fMeanY = fSumY.get() / fCount;
4793 for (SCSIZE i = 0; i < nC1; i++)
4794 {
4795 for (SCSIZE j = 0; j < nR1; j++)
4796 {
4797 if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4798 {
4799 double fValX = pMat1->GetDouble(i,j);
4800 double fValY = pMat2->GetDouble(i,j);
4801 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4802 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4803 }
4804 }
4805 }
4806 if (fSumSqrDeltaX == 0.0)
4807 PushError( FormulaError::DivisionByZero);
4808 else
4809 PushDouble( fMeanY + fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * (fVal - fMeanX));
4810 }
4811}
4812
4813static void lcl_roundUpNearestPow2(SCSIZE& nNum, SCSIZE& nNumBits)
4814{
4815 // Find the least power of 2 that is less than or equal to nNum.
4816 SCSIZE nPow2(1);
4817 nNumBits = std::numeric_limits<SCSIZE>::digits;
4818 nPow2 <<= (nNumBits - 1);
4819 while (nPow2 >= 1)
4820 {
4821 if (nNum & nPow2)
4822 break;
4823
4824 --nNumBits;
4825 nPow2 >>= 1;
4826 }
4827
4828 if (nPow2 != nNum)
4829 nNum = nPow2 ? (nPow2 << 1) : 1;
4830 else
4831 --nNumBits;
4832}
4833
4835{
4836 SCSIZE nOut = 0;
4837 for (SCSIZE nMask = 1; nMask < nBound; nMask <<= 1)
4838 {
4839 nOut <<= 1;
4840
4841 if (nIn & nMask)
4842 nOut |= 1;
4843 }
4844
4845 return nOut;
4846}
4847
4848namespace {
4849
4850// Computes and stores twiddle factors for computing DFT later.
4851struct ScTwiddleFactors
4852{
4853 ScTwiddleFactors(SCSIZE nN, bool bInverse) :
4854 mfWReal(nN),
4855 mfWImag(nN),
4856 mnN(nN),
4857 mbInverse(bInverse)
4858 {}
4859
4860 void Compute();
4861
4862 void Conjugate()
4863 {
4864 mbInverse = !mbInverse;
4865 for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
4866 mfWImag[nIdx] = -mfWImag[nIdx];
4867 }
4868
4869 std::vector<double> mfWReal;
4870 std::vector<double> mfWImag;
4871 SCSIZE mnN;
4872 bool mbInverse;
4873};
4874
4875}
4876
4877void ScTwiddleFactors::Compute()
4878{
4879 mfWReal.resize(mnN);
4880 mfWImag.resize(mnN);
4881
4882 double nW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnN);
4883
4884 if (mnN == 1)
4885 {
4886 mfWReal[0] = 1.0;
4887 mfWImag[0] = 0.0;
4888 }
4889 else if (mnN == 2)
4890 {
4891 mfWReal[0] = 1;
4892 mfWImag[0] = 0;
4893
4894 mfWReal[1] = -1;
4895 mfWImag[1] = 0;
4896 }
4897 else if (mnN == 4)
4898 {
4899 mfWReal[0] = 1;
4900 mfWImag[0] = 0;
4901
4902 mfWReal[1] = 0;
4903 mfWImag[1] = (mbInverse ? 1.0 : -1.0);
4904
4905 mfWReal[2] = -1;
4906 mfWImag[2] = 0;
4907
4908 mfWReal[3] = 0;
4909 mfWImag[3] = (mbInverse ? -1.0 : 1.0);
4910 }
4911 else if ((mnN % 4) == 0)
4912 {
4913 const SCSIZE nQSize = mnN >> 2;
4914 // Compute cos of the start quadrant.
4915 // This is the first quadrant if mbInverse == true, else it is the fourth quadrant.
4916 for (SCSIZE nIdx = 0; nIdx <= nQSize; ++nIdx)
4917 mfWReal[nIdx] = cos(nW*static_cast<double>(nIdx));
4918
4919 if (mbInverse)
4920 {
4921 const SCSIZE nQ1End = nQSize;
4922 // First quadrant
4923 for (SCSIZE nIdx = 0; nIdx <= nQ1End; ++nIdx)
4924 mfWImag[nIdx] = mfWReal[nQ1End-nIdx];
4925
4926 // Second quadrant
4927 const SCSIZE nQ2End = nQ1End << 1;
4928 for (SCSIZE nIdx = nQ1End+1; nIdx <= nQ2End; ++nIdx)
4929 {
4930 mfWReal[nIdx] = -mfWReal[nQ2End - nIdx];
4931 mfWImag[nIdx] = mfWImag[nQ2End - nIdx];
4932 }
4933
4934 // Third quadrant
4935 const SCSIZE nQ3End = nQ2End + nQ1End;
4936 for (SCSIZE nIdx = nQ2End+1; nIdx <= nQ3End; ++nIdx)
4937 {
4938 mfWReal[nIdx] = -mfWReal[nIdx - nQ2End];
4939 mfWImag[nIdx] = -mfWImag[nIdx - nQ2End];
4940 }
4941
4942 // Fourth Quadrant
4943 for (SCSIZE nIdx = nQ3End+1; nIdx < mnN; ++nIdx)
4944 {
4945 mfWReal[nIdx] = mfWReal[mnN - nIdx];
4946 mfWImag[nIdx] = -mfWImag[mnN - nIdx];
4947 }
4948 }
4949 else
4950 {
4951 const SCSIZE nQ4End = nQSize;
4952 const SCSIZE nQ3End = nQSize << 1;
4953 const SCSIZE nQ2End = nQ3End + nQSize;
4954
4955 // Fourth quadrant.
4956 for (SCSIZE nIdx = 0; nIdx <= nQ4End; ++nIdx)
4957 mfWImag[nIdx] = -mfWReal[nQ4End - nIdx];
4958
4959 // Third quadrant.
4960 for (SCSIZE nIdx = nQ4End+1; nIdx <= nQ3End; ++nIdx)
4961 {
4962 mfWReal[nIdx] = -mfWReal[nQ3End - nIdx];
4963 mfWImag[nIdx] = mfWImag[nQ3End - nIdx];
4964 }
4965
4966 // Second quadrant.
4967 for (SCSIZE nIdx = nQ3End+1; nIdx <= nQ2End; ++nIdx)
4968 {
4969 mfWReal[nIdx] = -mfWReal[nIdx - nQ3End];
4970 mfWImag[nIdx] = -mfWImag[nIdx - nQ3End];
4971 }
4972
4973 // First quadrant.
4974 for (SCSIZE nIdx = nQ2End+1; nIdx < mnN; ++nIdx)
4975 {
4976 mfWReal[nIdx] = mfWReal[mnN - nIdx];
4977 mfWImag[nIdx] = -mfWImag[mnN - nIdx];
4978 }
4979 }
4980 }
4981 else
4982 {
4983 for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
4984 {
4985 double fAngle = nW*static_cast<double>(nIdx);
4986 mfWReal[nIdx] = cos(fAngle);
4987 mfWImag[nIdx] = sin(fAngle);
4988 }
4989 }
4990}
4991
4992namespace {
4993
4994// A radix-2 decimation in time FFT algorithm for complex valued input.
4995class ScComplexFFT2
4996{
4997public:
4998 // rfArray.size() would always be even and a power of 2. (asserted in prepare())
4999 // rfArray's first half contains the real parts and the later half contains the imaginary parts.
5000 ScComplexFFT2(std::vector<double>& raArray, bool bInverse, bool bPolar, double fMinMag,
5001 ScTwiddleFactors& rTF, bool bSubSampleTFs = false, bool bDisableNormalize = false) :
5002 mrArray(raArray),
5003 mfWReal(rTF.mfWReal),
5004 mfWImag(rTF.mfWImag),
5005 mnPoints(raArray.size()/2),
5006 mnStages(0),
5007 mfMinMag(fMinMag),
5008 mbInverse(bInverse),
5009 mbPolar(bPolar),
5010 mbDisableNormalize(bDisableNormalize),
5011 mbSubSampleTFs(bSubSampleTFs)
5012 {}
5013
5014 void Compute();
5015
5016private:
5017
5018 void prepare();
5019
5020 double getReal(SCSIZE nIdx)
5021 {
5022 return mrArray[nIdx];
5023 }
5024
5025 void setReal(double fVal, SCSIZE nIdx)
5026 {
5027 mrArray[nIdx] = fVal;
5028 }
5029
5030 double getImag(SCSIZE nIdx)
5031 {
5032 return mrArray[mnPoints + nIdx];
5033 }
5034
5035 void setImag(double fVal, SCSIZE nIdx)
5036 {
5037 mrArray[mnPoints + nIdx] = fVal;
5038 }
5039
5040 SCSIZE getTFactorIndex(SCSIZE nPtIndex, SCSIZE nTfIdxScaleBits)
5041 {
5042 return ( ( nPtIndex << nTfIdxScaleBits ) & ( mnPoints - 1 ) ); // (x & (N-1)) is same as (x % N) but faster.
5043 }
5044
5045 void computeFly(SCSIZE nTopIdx, SCSIZE nBottomIdx, SCSIZE nWIdx1, SCSIZE nWIdx2)
5046 {
5047 if (mbSubSampleTFs)
5048 {
5049 nWIdx1 <<= 1;
5050 nWIdx2 <<= 1;
5051 }
5052
5053 const double x1r = getReal(nTopIdx);
5054 const double x2r = getReal(nBottomIdx);
5055
5056 const double& w1r = mfWReal[nWIdx1];
5057 const double& w1i = mfWImag[nWIdx1];
5058
5059 const double& w2r = mfWReal[nWIdx2];
5060 const double& w2i = mfWImag[nWIdx2];
5061
5062 const double x1i = getImag(nTopIdx);
5063 const double x2i = getImag(nBottomIdx);
5064
5065 setReal(x1r + x2r*w1r - x2i*w1i, nTopIdx);
5066 setImag(x1i + x2i*w1r + x2r*w1i, nTopIdx);
5067
5068 setReal(x1r + x2r*w2r - x2i*w2i, nBottomIdx);
5069 setImag(x1i + x2i*w2r + x2r*w2i, nBottomIdx);
5070 }
5071
5072 std::vector<double>& mrArray;
5073 std::vector<double>& mfWReal;
5074 std::vector<double>& mfWImag;
5075 SCSIZE mnPoints;
5076 SCSIZE mnStages;
5077 double mfMinMag;
5078 bool mbInverse:1;
5079 bool mbPolar:1;
5080 bool mbDisableNormalize:1;
5081 bool mbSubSampleTFs:1;
5082};
5083
5084}
5085
5086void ScComplexFFT2::prepare()
5087{
5088 SCSIZE nPoints = mnPoints;
5089 lcl_roundUpNearestPow2(nPoints, mnStages);
5090 assert(nPoints == mnPoints);
5091
5092 // Reorder array by bit-reversed indices.
5093 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5094 {
5095 SCSIZE nRevIdx = lcl_bitReverse(nIdx, mnPoints);
5096 if (nIdx < nRevIdx)
5097 {
5098 double fTmp = getReal(nIdx);
5099 setReal(getReal(nRevIdx), nIdx);
5100 setReal(fTmp, nRevIdx);
5101
5102 fTmp = getImag(nIdx);
5103 setImag(getImag(nRevIdx), nIdx);
5104 setImag(fTmp, nRevIdx);
5105 }
5106 }
5107}
5108
5109static void lcl_normalize(std::vector<double>& rCmplxArray, bool bScaleOnlyReal)
5110{
5111 const SCSIZE nPoints = rCmplxArray.size()/2;
5112 const double fScale = 1.0/static_cast<double>(nPoints);
5113
5114 // Scale the real part
5115 for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
5116 rCmplxArray[nIdx] *= fScale;
5117
5118 if (!bScaleOnlyReal)
5119 {
5120 const SCSIZE nLen = nPoints*2;
5121 for (SCSIZE nIdx = nPoints; nIdx < nLen; ++nIdx)
5122 rCmplxArray[nIdx] *= fScale;
5123 }
5124}
5125
5126static void lcl_convertToPolar(std::vector<double>& rCmplxArray, double fMinMag)
5127{
5128 const SCSIZE nPoints = rCmplxArray.size()/2;
5129 double fMag, fPhase, fR, fI;
5130 for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
5131 {
5132 fR = rCmplxArray[nIdx];
5133 fI = rCmplxArray[nPoints+nIdx];
5134 fMag = std::hypot(fR, fI);
5135 if (fMag < fMinMag)
5136 {
5137 fMag = 0.0;
5138 fPhase = 0.0;
5139 }
5140 else
5141 {
5142 fPhase = atan2(fI, fR);
5143 }
5144
5145 rCmplxArray[nIdx] = fMag;
5146 rCmplxArray[nPoints+nIdx] = fPhase;
5147 }
5148}
5149
5150void ScComplexFFT2::Compute()
5151{
5152 prepare();
5153
5154 const SCSIZE nFliesInStage = mnPoints/2;
5155 for (SCSIZE nStage = 0; nStage < mnStages; ++nStage)
5156 {
5157 const SCSIZE nTFIdxScaleBits = mnStages - nStage - 1; // Twiddle factor index's scale factor in bits.
5158 const SCSIZE nFliesInGroup = SCSIZE(1) << nStage;
5159 const SCSIZE nGroups = nFliesInStage/nFliesInGroup;
5160 const SCSIZE nFlyWidth = nFliesInGroup;
5161 for (SCSIZE nGroup = 0, nFlyTopIdx = 0; nGroup < nGroups; ++nGroup)
5162 {
5163 for (SCSIZE nFly = 0; nFly < nFliesInGroup; ++nFly, ++nFlyTopIdx)
5164 {
5165 SCSIZE nFlyBottomIdx = nFlyTopIdx + nFlyWidth;
5166 SCSIZE nWIdx1 = getTFactorIndex(nFlyTopIdx, nTFIdxScaleBits);
5167 SCSIZE nWIdx2 = getTFactorIndex(nFlyBottomIdx, nTFIdxScaleBits);
5168
5169 computeFly(nFlyTopIdx, nFlyBottomIdx, nWIdx1, nWIdx2);
5170 }
5171
5172 nFlyTopIdx += nFlyWidth;
5173 }
5174 }
5175
5176 if (mbPolar)
5177 lcl_convertToPolar(mrArray, mfMinMag);
5178
5179 // Normalize after converting to polar, so we have a chance to
5180 // save O(mnPoints) flops.
5181 if (mbInverse && !mbDisableNormalize)
5182 lcl_normalize(mrArray, mbPolar);
5183}
5184
5185namespace {
5186
5187// Bluestein's algorithm or chirp z-transform algorithm that can be used to
5188// compute DFT of a complex valued input of any length N in O(N lgN) time.
5189class ScComplexBluesteinFFT
5190{
5191public:
5192
5193 ScComplexBluesteinFFT(std::vector<double>& rArray, bool bReal, bool bInverse,
5194 bool bPolar, double fMinMag, bool bDisableNormalize = false) :
5195 mrArray(rArray),
5196 mnPoints(rArray.size()/2), // rArray should have space for imaginary parts even if real input.
5197 mfMinMag(fMinMag),
5198 mbReal(bReal),
5199 mbInverse(bInverse),
5200 mbPolar(bPolar),
5201 mbDisableNormalize(bDisableNormalize)
5202 {}
5203
5204 void Compute();
5205
5206private:
5207 std::vector<double>& mrArray;
5208 const SCSIZE mnPoints;
5209 double mfMinMag;
5210 bool mbReal:1;
5211 bool mbInverse:1;
5212 bool mbPolar:1;
5213 bool mbDisableNormalize:1;
5214};
5215
5216}
5217
5218void ScComplexBluesteinFFT::Compute()
5219{
5220 std::vector<double> aRealScalars(mnPoints);
5221 std::vector<double> aImagScalars(mnPoints);
5222 double fW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnPoints);
5223 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5224 {
5225 double fAngle = 0.5*fW*static_cast<double>(nIdx*nIdx);
5226 aRealScalars[nIdx] = cos(fAngle);
5227 aImagScalars[nIdx] = sin(fAngle);
5228 }
5229
5230 SCSIZE nMinSize = mnPoints*2 - 1;
5231 SCSIZE nExtendedLength = nMinSize, nTmp = 0;
5232 lcl_roundUpNearestPow2(nExtendedLength, nTmp);
5233 std::vector<double> aASignal(nExtendedLength*2); // complex valued
5234 std::vector<double> aBSignal(nExtendedLength*2); // complex valued
5235
5236 double fReal, fImag;
5237 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5238 {
5239 // Real part of A signal.
5240 aASignal[nIdx] = mrArray[nIdx]*aRealScalars[nIdx] + (mbReal ? 0.0 : -mrArray[mnPoints+nIdx]*aImagScalars[nIdx]);
5241 // Imaginary part of A signal.
5242 aASignal[nExtendedLength + nIdx] = mrArray[nIdx]*aImagScalars[nIdx] + (mbReal ? 0.0 : mrArray[mnPoints+nIdx]*aRealScalars[nIdx]);
5243
5244 // Real part of B signal.
5245 aBSignal[nIdx] = fReal = aRealScalars[nIdx];
5246 // Imaginary part of B signal.
5247 aBSignal[nExtendedLength + nIdx] = fImag = -aImagScalars[nIdx]; // negative sign because B signal is the conjugation of the scalars.
5248
5249 if (nIdx)
5250 {
5251 // B signal needs a mirror of its part in 0 < n < mnPoints at the tail end.
5252 aBSignal[nExtendedLength - nIdx] = fReal;
5253 aBSignal[(nExtendedLength<<1) - nIdx] = fImag;
5254 }
5255 }
5256
5257 {
5258 ScTwiddleFactors aTF(nExtendedLength, false /*not inverse*/);
5259 aTF.Compute();
5260
5261 // Do complex-FFT2 of both A and B signal.
5262 ScComplexFFT2 aFFT2A(aASignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
5263 aTF, false /*no subsample*/, true /* disable normalize */);
5264 aFFT2A.Compute();
5265
5266 ScComplexFFT2 aFFT2B(aBSignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
5267 aTF, false /*no subsample*/, true /* disable normalize */);
5268 aFFT2B.Compute();
5269
5270 double fAR, fAI, fBR, fBI;
5271 for (SCSIZE nIdx = 0; nIdx < nExtendedLength; ++nIdx)
5272 {
5273 fAR = aASignal[nIdx];
5274 fAI = aASignal[nExtendedLength + nIdx];
5275 fBR = aBSignal[nIdx];
5276 fBI = aBSignal[nExtendedLength + nIdx];
5277
5278 // Do point-wise product inplace in A signal.
5279 aASignal[nIdx] = fAR*fBR - fAI*fBI;
5280 aASignal[nExtendedLength + nIdx] = fAR*fBI + fAI*fBR;
5281 }
5282
5283 // Do complex-inverse-FFT2 of aASignal.
5284 aTF.Conjugate();
5285 ScComplexFFT2 aFFT2AI(aASignal, true /*inverse*/, false /*no polar*/, 0.0 /* no clipping */, aTF); // Need normalization here.
5286 aFFT2AI.Compute();
5287 }
5288
5289 // Point-wise multiply with scalars.
5290 for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5291 {
5292 fReal = aASignal[nIdx];
5293 fImag = aASignal[nExtendedLength + nIdx];
5294 mrArray[nIdx] = fReal*aRealScalars[nIdx] - fImag*aImagScalars[nIdx]; // no conjugation needed here.
5295 mrArray[mnPoints + nIdx] = fReal*aImagScalars[nIdx] + fImag*aRealScalars[nIdx];
5296 }
5297
5298 // Normalize/Polar operations
5299 if (mbPolar)
5300 lcl_convertToPolar(mrArray, mfMinMag);
5301
5302 // Normalize after converting to polar, so we have a chance to
5303 // save O(mnPoints) flops.
5304 if (mbInverse && !mbDisableNormalize)
5305 lcl_normalize(mrArray, mbPolar);
5306}
5307
5308namespace {
5309
5310// Computes DFT of an even length(N) real-valued input by using a
5311// ScComplexFFT2 if N == 2^k for some k or else by using a ScComplexBluesteinFFT
5312// with a complex valued input of length = N/2.
5313class ScRealFFT
5314{
5315public:
5316
5317 ScRealFFT(std::vector<double>& rInArray, std::vector<double>& rOutArray, bool bInverse,
5318 bool bPolar, double fMinMag) :
5319 mrInArray(rInArray),
5320 mrOutArray(rOutArray),
5321 mfMinMag(fMinMag),
5322 mbInverse(bInverse),
5323 mbPolar(bPolar)
5324 {}
5325
5326 void Compute();
5327
5328private:
5329 std::vector<double>& mrInArray;
5330 std::vector<double>& mrOutArray;
5331 double mfMinMag;
5332 bool mbInverse:1;
5333 bool mbPolar:1;
5334};
5335
5336}
5337
5338void ScRealFFT::Compute()
5339{
5340 // input length has to be even to do this optimization.
5341 assert(mrInArray.size() % 2 == 0);
5342 assert(mrInArray.size()*2 == mrOutArray.size());
5343 // nN is the number of points in the complex-fft input
5344 // which will be half of the number of points in real array.
5345 const SCSIZE nN = mrInArray.size()/2;
5346 if (nN == 0)
5347 {
5348 mrOutArray[0] = mrInArray[0];
5349 mrOutArray[1] = 0.0;
5350 return;
5351 }
5352
5353 // work array should be the same length as mrInArray
5354 std::vector<double> aWorkArray(nN*2);
5355 for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
5356 {
5357 SCSIZE nDoubleIdx = 2*nIdx;
5358 // Use even elements as real part
5359 aWorkArray[nIdx] = mrInArray[nDoubleIdx];
5360 // and odd elements as imaginary part of the contrived complex sequence.
5361 aWorkArray[nN+nIdx] = mrInArray[nDoubleIdx+1];
5362 }
5363
5364 ScTwiddleFactors aTFs(nN*2, mbInverse);
5365 aTFs.Compute();
5366 SCSIZE nNextPow2 = nN, nTmp = 0;
5367 lcl_roundUpNearestPow2(nNextPow2, nTmp);
5368
5369 if (nNextPow2 == nN)
5370 {
5371 ScComplexFFT2 aFFT2(aWorkArray, mbInverse, false /*disable polar*/, 0.0 /* no clipping */,
5372 aTFs, true /*subsample tf*/, true /*disable normalize*/);
5373 aFFT2.Compute();
5374 }
5375 else
5376 {
5377 ScComplexBluesteinFFT aFFT(aWorkArray, false /*complex input*/, mbInverse, false /*disable polar*/,
5378 0.0 /* no clipping */, true /*disable normalize*/);
5379 aFFT.Compute();
5380 }
5381
5382 // Post process aWorkArray to populate mrOutArray
5383
5384 const SCSIZE nTwoN = 2*nN, nThreeN = 3*nN;
5385 double fY1R, fY2R, fY1I, fY2I, fResR, fResI, fWR, fWI;
5386 for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
5387 {
5388 const SCSIZE nIdxRev = nIdx ? (nN - nIdx) : 0;
5389 fY1R = aWorkArray[nIdx];
5390 fY2R = aWorkArray[nIdxRev];
5391 fY1I = aWorkArray[nN + nIdx];
5392 fY2I = aWorkArray[nN + nIdxRev];
5393 fWR = aTFs.mfWReal[nIdx];
5394 fWI = aTFs.mfWImag[nIdx];
5395
5396 // mrOutArray has length = 4*nN
5397 // Real part of the final output (only half of the symmetry around Nyquist frequency)
5398 // Fills the first quarter.
5399 mrOutArray[nIdx] = fResR = 0.5*(
5400 fY1R + fY2R +
5401 fWR * (fY1I + fY2I) +
5402 fWI * (fY1R - fY2R) );
5403 // Imaginary part of the final output (only half of the symmetry around Nyquist frequency)
5404 // Fills the third quarter.
5405 mrOutArray[nTwoN + nIdx] = fResI = 0.5*(
5406 fY1I - fY2I +
5407 fWI * (fY1I + fY2I) -
5408 fWR * (fY1R - fY2R) );
5409
5410 // Fill the missing 2 quarters using symmetry argument.
5411 if (nIdx)
5412 {
5413 // Fills the 2nd quarter.
5414 mrOutArray[nN + nIdxRev] = fResR;
5415 // Fills the 4th quarter.
5416 mrOutArray[nThreeN + nIdxRev] = -fResI;
5417 }
5418 else
5419 {
5420 mrOutArray[nN] = fY1R - fY1I;
5421 mrOutArray[nThreeN] = 0.0;
5422 }
5423 }
5424
5425 // Normalize/Polar operations
5426 if (mbPolar)
5427 lcl_convertToPolar(mrOutArray, mfMinMag);
5428
5429 // Normalize after converting to polar, so we have a chance to
5430 // save O(mnPoints) flops.
5431 if (mbInverse)
5432 lcl_normalize(mrOutArray, mbPolar);
5433}
5434
5435using ScMatrixGenerator = ScMatrixRef(SCSIZE, SCSIZE, std::vector<double>&);
5436
5437namespace {
5438
5439// Generic FFT class that decides which FFT implementation to use.
5440class ScFFT
5441{
5442public:
5443
5444 ScFFT(ScMatrixRef& pMat, bool bReal, bool bInverse, bool bPolar, double fMinMag) :
5445 mpInputMat(pMat),
5446 mfMinMag(fMinMag),
5447 mbReal(bReal),
5448 mbInverse(bInverse),
5449 mbPolar(bPolar)
5450 {}
5451
5452 ScMatrixRef Compute(const std::function<ScMatrixGenerator>& rMatGenFunc);
5453
5454private:
5455 ScMatrixRef& mpInputMat;
5456 double mfMinMag;<