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 
40 using ::std::vector;
41 using namespace formula;
42 
44 // This is an arbitrary limit.
45 static size_t MAX_COUNT_DOUBLE_FOR_SORT(const ScSheetLimits& rSheetLimits)
46 {
47  return rSheetLimits.GetMaxRowCount() * 2;
48 }
49 
50 const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental
51 const double fMachEps = ::std::numeric_limits<double>::epsilon();
52 
53 namespace {
54 
55 class ScDistFunc
56 {
57 public:
58  virtual double GetValue(double x) const = 0;
59 
60 protected:
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 
71 static bool lcl_HasChangeOfSign( double u, double w )
72 {
73  return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);
74 }
75 
76 static 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 
195 double ScInterpreter::phi(double x)
196 {
197  return 0.39894228040143268 * exp(-(x * x) / 2.0);
198 }
199 
201 { // Using gauss(x)+0.5 has severe cancellation errors for x<-4
202  return 0.5 * std::erfc(-x * 0.7071067811865475); // * 1/sqrt(2)
203 }
204 
205 double 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 
215 double 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 
268 double ScInterpreter::gaussinv(double x)
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 
439 double ScInterpreter::Fakultaet(double x)
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 
460 double 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
491 static 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
559 static 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
577 static 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 
585 double 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) - rtl::math::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 
629 double ScInterpreter::GetLogGamma(double fZ)
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) - rtl::math::log1p(fZ) - log(fZ);
638 }
639 
640 double 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 
648 double 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
669 double 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
681 double 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 
689 double 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)
739  PushIllegalArgument();
740  else
741  {
742  double fX = GetDouble();
743  if (bCumulative)
744  PushDouble(GetChiSqDistCDF(fX,fDF));
745  else
746  PushDouble(GetChiSqDistPDF(fX,fDF));
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 )
758  PushIllegalArgument();
759  else
760  {
761  double fX = GetDouble();
762  if ( fX < 0 )
763  PushIllegalArgument();
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))
778  PushIllegalArgument();
779  else
780  {
781  double fResult = GetGamma(x);
782  if (nGlobalError != FormulaError::NONE)
783  {
784  PushError( nGlobalError);
785  return;
786  }
787  PushDouble(fResult);
788  }
789 }
790 
792 {
793  double x = GetDouble();
794  if (x > 0.0) // constraint from ODFF
795  PushDouble( GetLogGamma(x));
796  else
797  PushIllegalArgument();
798 }
799 
800 double 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 * ::rtl::math::log1p(fTempA)
827  -fB * ::rtl::math::log1p(fTempB)-fgm);
828  fResult *= fLanczos;
829  return fResult;
830 }
831 
832 // Same as GetBeta but with logarithm
833 double 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 * ::rtl::math::log1p(fTempA)
856  -fB * ::rtl::math::log1p(fTempB)-fgm;
857  fResult += fLogLanczos;
858  return fResult;
859 }
860 
861 // beta distribution probability density function
862 double 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) * ::rtl::math::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) ? ::rtl::math::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 */
937 static 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
985 double 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 * ::rtl::math::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 = ::rtl::math::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  {
1062  PushIllegalArgument();
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  {
1122  PushIllegalArgument();
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 {
1142  PushDouble(phi(GetDouble()));
1143 }
1144 
1146 {
1147  PushDouble(gauss(GetDouble()));
1148 }
1149 
1151 {
1152  double fVal = GetDouble();
1153  if (std::abs(fVal) >= 1.0)
1154  PushIllegalArgument();
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)
1168  PushIllegalArgument();
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)
1180  PushIllegalArgument();
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)
1193  PushIllegalArgument();
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)
1207  PushIllegalArgument();
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)
1226  PushIllegalArgument();
1227  else
1228  PushDouble(pow(n,k));
1229  }
1230 }
1231 
1232 double 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 
1260 static 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)
1290  PushIllegalArgument();
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
1296  PushDouble(GetBinomDistPMF(x,n,p));
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
1309  PushDouble(GetBinomDistPMF(xs,n,p));
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
1338  PushIllegalArgument();
1339  }
1340  else
1341  PushIllegalArgument();
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  {
1358  PushIllegalArgument();
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)
1372  PushDouble( GetBinomDistPMF(x,n,p));
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)
1419  PushIllegalArgument();
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)
1509  PushIllegalArgument();
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 )
1530  PushIllegalArgument();
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 
1546 void 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  {
1557  PushIllegalArgument();
1558  return;
1559  }
1560  if (bCumulative)
1561  PushDouble(integralPhi((x-mue)/sigma));
1562  else
1563  PushDouble(phi((x-mue)/sigma)/sigma);
1564 }
1565 
1566 void 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  {
1577  PushIllegalArgument();
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)
1590  PushIllegalArgument();
1591  else
1592  PushDouble(phi((log(x)-mue)/sigma)/sigma/x);
1593  }
1594 }
1595 
1597 {
1598  PushDouble(integralPhi(GetDouble()));
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)
1624  PushIllegalArgument();
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  {
1650  PushIllegalArgument();
1651  return;
1652  }
1653  PushDouble( GetTDist( T, fDF, static_cast<int>(fFlag) ) );
1654 }
1655 
1656 void ScInterpreter::ScTDist_T( int nTails )
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  {
1664  PushIllegalArgument();
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  {
1683  PushIllegalArgument();
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  {
1698  PushIllegalArgument();
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  {
1724  PushIllegalArgument();
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 
1741 void ScInterpreter::ScChiDist( bool bODFF )
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  {
1751  PushIllegalArgument();
1752  return;
1753  }
1754  fResult = GetChiDist( fChi, fDF);
1755  if (nGlobalError != FormulaError::NONE)
1756  {
1757  PushError( nGlobalError);
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)
1773  PushIllegalArgument();
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)
1791  PushIllegalArgument();
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 
1834 static 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 
1856 void 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  {
1870  PushIllegalArgument();
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 
1895 double 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 
2091 void ScInterpreter::ScGammaDist( bool bODFF )
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)
2106  PushIllegalArgument();
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)
2124  PushIllegalArgument();
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)
2136  PushIllegalArgument();
2137  else if (x == 0.0 || x == 1.0)
2138  PushNoValue();
2139  else
2140  PushDouble(gaussinv(x));
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 )
2152  PushIllegalArgument();
2153  else
2154  PushDouble( exp( fMue + fSigma * gaussinv( fP ) ) );
2155  }
2156 }
2157 
2158 class ScGammaDistFunction : public ScDistFunc
2159 {
2161  double fp, fAlpha, fBeta;
2162 
2163 public:
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  {
2181  PushIllegalArgument();
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 
2198 class ScBetaDistFunction : public ScDistFunc
2199 {
2201  double fp, fAlpha, fBeta;
2202 
2203 public:
2204  ScBetaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2205  rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2206 
2207  virtual ~ScBetaDistFunction() {}
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  {
2231  PushIllegalArgument();
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 
2249 class ScTDistFunction : public ScDistFunc
2250 {
2252  double fp, fDF;
2253  int nT;
2254 
2255 public:
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  {
2272  PushIllegalArgument();
2273  return;
2274  }
2275  if ( nType == 4 ) // left-tailed cumulative t-distribution
2276  {
2277  if ( fP == 1.0 )
2278  PushIllegalArgument();
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 
2288 double 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 
2298 class ScFDistFunction : public ScDistFunc
2299 {
2301  double fp, fF1, fF2;
2302 
2303 public:
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  {
2321  PushIllegalArgument();
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  {
2342  PushIllegalArgument();
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 
2354 class ScChiDistFunction : public ScDistFunc
2355 {
2357  double fp, fDF;
2358 
2359 public:
2360  ScChiDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2361  rInt(rI), fp(fpVal), fDF(fDFVal) {}
2362 
2363  virtual ~ScChiDistFunction() {}
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  {
2376  PushIllegalArgument();
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 /***********************************************/
2389 class ScChiSqDistFunction : public ScDistFunc
2390 {
2392  double fp, fDF;
2393 
2394 public:
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  {
2411  PushIllegalArgument();
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)
2431  PushIllegalArgument();
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)
2445  PushIllegalArgument();
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  {
2464  PushIllegalArgument();
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( mrDoc, 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 :
2526  case svExternalSingleRef:
2527  case svExternalDoubleRef:
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 
2580 bool ScInterpreter::CalculateTest(bool _bTemplin
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  {
2655  PushIllegalArgument();
2656  return;
2657  }
2658 
2659  ScMatrixRef pMat2 = GetMatrix();
2660  ScMatrixRef pMat1 = GetMatrix();
2661  if (!pMat1 || !pMat2)
2662  {
2663  PushIllegalParameter();
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  {
2676  PushIllegalArgument();
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  {
2724  PushIllegalArgument();
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  {
2738  PushIllegalParameter();
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  {
2788  PushIllegalParameter();
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  {
2797  PushIllegalArgument();
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  {
2834  PushIllegalArgument();
2835  return;
2836  }
2837  }
2838  }
2839  }
2840  if ( bEmpty )
2841  {
2842  // not in ODFF1.2, but for interoperability with Excel
2843  PushIllegalArgument();
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( mrDoc, 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 :
2976  case svExternalSingleRef:
2977  case svExternalDoubleRef:
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
3021  PushError( nGlobalError);
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(mrDoc, 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 :
3130  case svExternalSingleRef:
3131  case svExternalDoubleRef:
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
3193  PushError( nGlobalError);
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 }
3211 bool 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( mrDoc, 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 :
3271  case svExternalSingleRef:
3272  case svExternalDoubleRef:
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  {
3310  PushError( nGlobalError);
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  {
3339  PushIllegalArgument();
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 
3366 double 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 
3400 double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile )
3401 {
3402  size_t nSize = rArray.size();
3403  if (nSize == 1)
3404  return rArray[0];
3405  else
3406  {
3407  size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * (nSize-1)));
3408  double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1));
3409  OSL_ENSURE(nIndex < nSize, "GetPercentile: wrong index(1)");
3410  vector<double>::iterator iter = rArray.begin() + nIndex;
3411  ::std::nth_element( rArray.begin(), iter, rArray.end());
3412  if (fDiff == 0.0)
3413  return *iter;
3414  else
3415  {
3416  OSL_ENSURE(nIndex < nSize-1, "GetPercentile: wrong index(2)");
3417  double fVal = *iter;
3418  iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
3419  return fVal + fDiff * (*iter - fVal);
3420  }
3421  }
3422 }
3423 
3424 double ScInterpreter::GetPercentileExclusive( vector<double> & rArray, double fPercentile )
3425 {
3426  size_t nSize1 = rArray.size() + 1;
3427  if ( rArray.empty() || nSize1 == 1 || nGlobalError != FormulaError::NONE)
3428  {
3429  SetError( FormulaError::NoValue );
3430  return 0.0;
3431  }
3432  if ( fPercentile * nSize1 < 1.0 || fPercentile * nSize1 > static_cast<double>( nSize1 - 1 ) )
3433  {
3434  SetError( FormulaError::IllegalParameter );
3435  return 0.0;
3436  }
3437 
3438  size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * nSize1 - 1 ));
3439  double fDiff = fPercentile * nSize1 - 1 - ::rtl::math::approxFloor( fPercentile * nSize1 - 1 );
3440  OSL_ENSURE(nIndex < ( nSize1 - 1 ), "GetPercentile: wrong index(1)");
3441  vector<double>::iterator iter = rArray.begin() + nIndex;
3442  ::std::nth_element( rArray.begin(), iter, rArray.end());
3443  if (fDiff == 0.0)
3444  return *iter;
3445  else
3446  {
3447  OSL_ENSURE(nIndex < nSize1, "GetPercentile: wrong index(2)");
3448  double fVal = *iter;
3449  iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
3450  return fVal + fDiff * (*iter - fVal);
3451  }
3452 }
3453 
3454 void ScInterpreter::ScPercentile( bool bInclusive )
3455 {
3456  if ( !MustHaveParamCount( GetByte(), 2 ) )
3457  return;
3458  double alpha = GetDouble();
3459  if ( bInclusive ? ( alpha < 0.0 || alpha > 1.0 ) : ( alpha <= 0.0 || alpha >= 1.0 ) )
3460  {
3461  PushIllegalArgument();
3462  return;
3463  }
3464  vector<double> aArray;
3465  GetNumberSequenceArray( 1, aArray, false );
3466  if ( aArray.empty() || nGlobalError != FormulaError::NONE )
3467  {
3468  SetError( FormulaError::NoValue );
3469  return;
3470  }
3471  if ( bInclusive )
3472  PushDouble( GetPercentile( aArray, alpha ));
3473  else
3474  PushDouble( GetPercentileExclusive( aArray, alpha ));
3475 }
3476 
3477 void ScInterpreter::ScQuartile( bool bInclusive )
3478 {
3479  if ( !MustHaveParamCount( GetByte(), 2 ) )
3480  return;
3481  double fFlag = ::rtl::math::approxFloor(GetDouble());
3482  if ( bInclusive ? ( fFlag < 0.0 || fFlag > 4.0 ) : ( fFlag <= 0.0 || fFlag >= 4.0 ) )
3483  {
3484  PushIllegalArgument();
3485  return;
3486  }
3487  vector<double> aArray;
3488  GetNumberSequenceArray( 1, aArray, false );
3489  if ( aArray.empty() || nGlobalError != FormulaError::NONE )
3490  {
3491  SetError( FormulaError::NoValue );
3492  return;
3493  }
3494  if ( bInclusive )
3495  PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentile( aArray, 0.25 * fFlag ) );
3496  else
3497  PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentileExclusive( aArray, 0.25 * fFlag ) );
3498 }
3499 
3501 {
3502  sal_uInt8 nParamCount = GetByte();
3503  if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3504  return;
3505  vector<double> aSortArray;
3506  GetSortArray( nParamCount, aSortArray, nullptr, false, false );
3507  SCSIZE nSize = aSortArray.size();
3508  if (nSize == 0 || nGlobalError != FormulaError::NONE)
3509  PushNoValue();
3510  else
3511  {
3512  SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
3513  double nOldVal = aSortArray[0];
3514  SCSIZE i;
3515  for ( i = 1; i < nSize; i++)
3516  {
3517  if (aSortArray[i] == nOldVal)
3518  nCount++;
3519  else
3520  {
3521  nOldVal = aSortArray[i];
3522  if (nCount > nMax)
3523  {
3524  nMax = nCount;
3525  nMaxIndex = i-1;
3526  }
3527  nCount = 1;
3528  }
3529  }
3530  if (nCount > nMax)
3531  {
3532  nMax = nCount;
3533  nMaxIndex = i-1;
3534  }
3535  if (nMax == 1 && nCount == 1)
3536  PushNoValue();
3537  else if (nMax == 1)
3538  PushDouble(nOldVal);
3539  else
3540  PushDouble(aSortArray[nMaxIndex]);
3541  }
3542 }
3543 
3545 {
3546  sal_uInt8 nParamCount = GetByte();
3547  if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3548  return;
3549  vector<double> aArray;
3550  GetNumberSequenceArray( nParamCount, aArray, false );
3551  vector< double > aSortArray( aArray );
3552  QuickSort( aSortArray, nullptr );
3553  SCSIZE nSize = aSortArray.size();
3554  if ( nSize == 0 || nGlobalError != FormulaError::NONE )
3555  PushNoValue();
3556  else
3557  {
3558  SCSIZE nMax = 1, nCount = 1;
3559  double nOldVal = aSortArray[ 0 ];
3560  vector< double > aResultArray( 1 );
3561  SCSIZE i;
3562  for ( i = 1; i < nSize; i++ )
3563  {
3564  if ( aSortArray[ i ] == nOldVal )
3565  nCount++;
3566  else
3567  {
3568  if ( nCount >= nMax && nCount > 1 )
3569  {
3570  if ( nCount > nMax )
3571  {
3572  nMax = nCount;
3573  if ( aResultArray.size() != 1 )
3574  vector< double >( 1 ).swap( aResultArray );
3575  aResultArray[ 0 ] = nOldVal;
3576  }
3577  else
3578  aResultArray.emplace_back( nOldVal );
3579  }
3580  nOldVal = aSortArray[ i ];
3581  nCount = 1;
3582  }
3583  }
3584  if ( nCount >= nMax && nCount > 1 )
3585  {
3586  if ( nCount > nMax )
3587  vector< double >().swap( aResultArray );
3588  aResultArray.emplace_back( nOldVal );
3589  }
3590  if ( nMax == 1 && nCount == 1 )
3591  PushNoValue();
3592  else if ( nMax == 1 )
3593  PushDouble( nOldVal ); // there is only 1 result, no reordering needed
3594  else
3595  {
3596  // sort resultArray according to ordering of aArray
3597  vector< vector< double > > aOrder;
3598  aOrder.resize( aResultArray.size(), vector< double >( 2 ) );
3599  for ( i = 0; i < aResultArray.size(); i++ )
3600  {
3601  for ( SCSIZE j = 0; j < nSize; j++ )
3602  {
3603  if ( aArray[ j ] == aResultArray[ i ] )
3604  {
3605  aOrder[ i ][ 0 ] = aResultArray[ i ];
3606  aOrder[ i ][ 1 ] = j;
3607  break;
3608  }
3609  }
3610  }
3611  sort( aOrder.begin(), aOrder.end(), []( const std::vector< double >& lhs,
3612  const std::vector< double >& rhs )
3613  { return lhs[ 1 ] < rhs[ 1 ]; } );
3614 
3615  if ( bSingle )
3616  PushDouble( aOrder[ 0 ][ 0 ] );
3617  else
3618  {
3619  // put result in correct order in aResultArray
3620  for ( i = 0; i < aResultArray.size(); i++ )
3621  aResultArray[ i ] = aOrder[ i ][ 0 ];
3622  ScMatrixRef pResMatrix = GetNewMat( 1, aResultArray.size(), true );
3623  pResMatrix->PutDoubleVector( aResultArray, 0, 0 );
3624  PushMatrix( pResMatrix );
3625  }
3626  }
3627  }
3628 }
3629 
3631 {
3632  if ( !MustHaveParamCount( GetByte(), 2 ) )
3633  return;
3634 
3635  SCSIZE nCol = 0, nRow = 0;
3636  auto aArray = GetTopNumberArray(nCol, nRow);
3637  const auto nRankArraySize = aArray.size();
3638  if (nRankArraySize == 0 || nGlobalError != FormulaError::NONE)
3639  {
3640  PushNoValue();
3641  return;
3642  }
3643  assert(nRankArraySize == nCol * nRow);
3644 
3645  std::vector<SCSIZE> aRankArray;
3646  aRankArray.reserve(nRankArraySize);
3647  std::transform(aArray.begin(), aArray.end(), std::back_inserter(aRankArray),
3648  [](double f) {
3649  f = rtl::math::approxFloor(f);
3650  // Valid ranks are >= 1.
3651  if (f < 1.0 || !o3tl::convertsToAtMost(f, std::numeric_limits<SCSIZE>::max()))
3652  return static_cast<SCSIZE>(0);
3653  return static_cast<SCSIZE>(f);
3654  });
3655 
3656  vector<double> aSortArray;
3657  GetNumberSequenceArray(1, aSortArray, false );
3658  const SCSIZE nSize = aSortArray.size();
3659  if (nSize == 0 || nGlobalError != FormulaError::NONE)
3660  PushNoValue();
3661  else if (nRankArraySize == 1)
3662  {
3663  const SCSIZE k = aRankArray[0];
3664  if (k < 1 || nSize < k)
3665  PushNoValue();
3666  else
3667  {
3668  vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
3669  ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
3670  PushDouble( *iPos);
3671  }
3672  }
3673  else
3674  {
3675  std::set<SCSIZE> aIndices;
3676  for (SCSIZE n : aRankArray)
3677  {
3678  if (1 <= n && n <= nSize)
3679  aIndices.insert(bSmall ? n-1 : nSize-n);
3680  }
3681  // We can spare sorting when the total number of ranks is small enough.
3682  // Find only the elements at given indices if, arbitrarily, the index size is
3683  // smaller than 1/3 of the haystack array's size; just sort it squarely, otherwise.
3684  if (aIndices.size() < nSize/3)
3685  {
3686  auto itBegin = aSortArray.begin();
3687  for (SCSIZE i : aIndices)
3688  {
3689  auto it = aSortArray.begin() + i;
3690  std::nth_element(itBegin, it, aSortArray.end());
3691  itBegin = ++it;
3692  }
3693  }
3694  else
3695  std::sort(aSortArray.begin(), aSortArray.end());
3696 
3697  aArray.clear();
3698  for (SCSIZE n : aRankArray)
3699  {
3700  if (1 <= n && n <= nSize)
3701  aArray.push_back( aSortArray[bSmall ? n-1 : nSize-n]);
3702  else
3703  aArray.push_back( CreateDoubleError( FormulaError::NoValue));
3704  }
3705  ScMatrixRef pResult = GetNewMat(nCol, nRow, aArray);
3706  PushMatrix(pResult);
3707  }
3708 }
3709 
3711 {
3712  CalculateSmallLarge(false);
3713 }
3714 
3716 {
3717  CalculateSmallLarge(true);
3718 }
3719 
3720 void ScInterpreter::ScPercentrank( bool bInclusive )
3721 {
3722  sal_uInt8 nParamCount = GetByte();
3723  if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
3724  return;
3725  double fSignificance = ( nParamCount == 3 ? ::rtl::math::approxFloor( GetDouble() ) : 3.0 );
3726  if ( fSignificance < 1.0 )
3727  {
3728  PushIllegalArgument();
3729  return;
3730  }
3731  double fNum = GetDouble();
3732  vector<double> aSortArray;
3733  GetSortArray( 1, aSortArray, nullptr, false, false );
3734  SCSIZE nSize = aSortArray.size();
3735  if ( nSize == 0 || nGlobalError != FormulaError::NONE )
3736  PushNoValue();
3737  else
3738  {
3739  if ( fNum < aSortArray[ 0 ] || fNum > aSortArray[ nSize - 1 ] )
3740  PushNoValue();
3741  else
3742  {
3743  double fRes;
3744  if ( nSize == 1 )
3745  fRes = 1.0; // fNum == aSortArray[ 0 ], see test above
3746  else
3747  fRes = GetPercentrank( aSortArray, fNum, bInclusive );
3748  if ( fRes != 0.0 )
3749  {
3750  double fExp = ::rtl::math::approxFloor( log10( fRes ) ) + 1.0 - fSignificance;
3751  fRes = ::rtl::math::round( fRes * pow( 10, -fExp ) ) / pow( 10, -fExp );
3752  }
3753  PushDouble( fRes );
3754  }
3755  }
3756 }
3757 
3758 double ScInterpreter::GetPercentrank( ::std::vector<double> & rArray, double fVal, bool bInclusive )
3759 {
3760  SCSIZE nSize = rArray.size();
3761  double fRes;
3762  if ( fVal == rArray[ 0 ] )
3763  {
3764  if ( bInclusive )
3765  fRes = 0.0;
3766  else
3767  fRes = 1.0 / static_cast<double>( nSize + 1 );
3768  }
3769  else
3770  {
3771  SCSIZE nOldCount = 0;
3772  double fOldVal = rArray[ 0 ];
3773  SCSIZE i;
3774  for ( i = 1; i < nSize && rArray[ i ] < fVal; i++ )
3775  {
3776  if ( rArray[ i ] != fOldVal )
3777  {
3778  nOldCount = i;
3779  fOldVal = rArray[ i ];
3780  }
3781  }
3782  if ( rArray[ i ] != fOldVal )
3783  nOldCount = i;
3784  if ( fVal == rArray[ i ] )
3785  {
3786  if ( bInclusive )
3787  fRes = div( nOldCount, nSize - 1 );
3788  else
3789  fRes = static_cast<double>( i + 1 ) / static_cast<double>( nSize + 1 );
3790  }
3791  else
3792  {
3793  // nOldCount is the count of smaller entries
3794  // fVal is between rArray[ nOldCount - 1 ] and rArray[ nOldCount ]
3795  // use linear interpolation to find a position between the entries
3796  if ( nOldCount == 0 )
3797  {
3798  OSL_FAIL( "should not happen" );
3799  fRes = 0.0;
3800  }
3801  else
3802  {
3803  double fFract = ( fVal - rArray[ nOldCount - 1 ] ) /
3804  ( rArray[ nOldCount ] - rArray[ nOldCount - 1 ] );
3805  if ( bInclusive )
3806  fRes = div( static_cast<double>( nOldCount - 1 ) + fFract, nSize - 1 );
3807  else
3808  fRes = ( static_cast<double>(nOldCount) + fFract ) / static_cast<double>( nSize + 1 );
3809  }
3810  }
3811  }
3812  return fRes;
3813 }
3814 
3816 {
3817  if ( !MustHaveParamCount( GetByte(), 2 ) )
3818  return;
3819  double alpha = GetDouble();
3820  if (alpha < 0.0 || alpha >= 1.0)
3821  {
3822  PushIllegalArgument();
3823  return;
3824  }
3825  vector<double> aSortArray;
3826  GetSortArray( 1, aSortArray, nullptr, false, false );
3827  SCSIZE nSize = aSortArray.size();
3828  if (nSize == 0 || nGlobalError != FormulaError::NONE)
3829  PushNoValue();
3830  else
3831  {
3832  sal_uLong nIndex = static_cast<sal_uLong>(::rtl::math::approxFloor(alpha*static_cast<double>(nSize)));
3833  if (nIndex % 2 != 0)
3834  nIndex--;
3835  nIndex /= 2;
3836  OSL_ENSURE(nIndex < nSize, "ScTrimMean: wrong index");
3837  KahanSum fSum = 0.0;
3838  for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
3839  fSum += aSortArray[i];
3840  PushDouble(fSum.get()/static_cast<double>(nSize-2*nIndex));
3841  }
3842 }
3843 
3844 std::vector<double> ScInterpreter::GetTopNumberArray( SCSIZE& rCol, SCSIZE& rRow )
3845 {
3846  std::vector<double> aArray;
3847  switch (GetStackType())
3848  {
3849  case svDouble:
3850  aArray.push_back(PopDouble());
3851  rCol = rRow = 1;
3852  break;
3853  case svSingleRef:
3854  {
3855  ScAddress aAdr;
3856  PopSingleRef(aAdr);
3857  ScRefCellValue aCell(mrDoc, aAdr);
3858  if (aCell.hasNumeric())
3859  {
3860  aArray.push_back(GetCellValue(aAdr, aCell));
3861  rCol = rRow = 1;
3862  }
3863  }
3864  break;
3865  case svDoubleRef:
3866  {
3867  ScRange aRange;
3868  PopDoubleRef(aRange, true);
3869  if (nGlobalError != FormulaError::NONE)
3870  break;
3871 
3872  // give up unless the start and end are in the same sheet
3873  if (aRange.aStart.Tab() != aRange.aEnd.Tab())
3874  {
3875  SetError(FormulaError::IllegalParameter);
3876  break;
3877  }
3878 
3879  // the range already is in order
3880  assert(aRange.aStart.Col() <= aRange.aEnd.Col());
3881  assert(aRange.aStart.Row() <= aRange.aEnd.Row());
3882  rCol = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3883  rRow = aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3884  aArray.reserve(rCol * rRow);
3885 
3886  FormulaError nErr = FormulaError::NONE;
3887  double fCellVal;
3888  ScValueIterator aValIter(mrDoc, aRange, mnSubTotalFlags);
3889  if (aValIter.GetFirst(fCellVal, nErr))
3890  {
3891  do
3892  aArray.push_back(fCellVal);
3893  while (aValIter.GetNext(fCellVal, nErr) && nErr == FormulaError::NONE);
3894  }
3895  if (aArray.size() != rCol * rRow)
3896  {
3897  aArray.clear();
3898  SetError(nErr);
3899  }
3900  }
3901  break;
3902  case svMatrix:
3903  case svExternalSingleRef:
3904  case svExternalDoubleRef:
3905  {
3906  ScMatrixRef pMat = GetMatrix();
3907  if (!pMat)
3908  break;
3909 
3910  if (pMat->IsNumeric())
3911  {
3912  SCSIZE nCount = pMat->GetElementCount();
3913  aArray.reserve(nCount);
3914  for (SCSIZE i = 0; i < nCount; ++i)
3915  aArray.push_back(pMat->GetDouble(i));
3916  pMat->GetDimensions(rCol, rRow);
3917  }
3918  else
3919  SetError(FormulaError::IllegalParameter);
3920  }
3921  break;
3922  default:
3923  SetError(FormulaError::IllegalParameter);
3924  break;
3925  }
3926  return aArray;
3927 }
3928 
3929 void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray, bool bConvertTextInArray )
3930 {
3931  ScAddress aAdr;
3932  ScRange aRange;
3933  const bool bIgnoreErrVal = bool(mnSubTotalFlags & SubtotalFlags::IgnoreErrVal);
3934  short nParam = nParamCount;
3935  size_t nRefInList = 0;
3936  ReverseStack( nParamCount );
3937  while (nParam-- > 0)
3938  {
3939  const StackVar eStackType = GetStackType();
3940  switch (eStackType)
3941  {
3942  case svDouble :
3943  rArray.push_back( PopDouble());
3944  break;
3945  case svSingleRef :
3946  {
3947  PopSingleRef( aAdr );
3948  ScRefCellValue aCell(mrDoc, aAdr);
3949  if (bIgnoreErrVal && aCell.hasError())
3950  ; // nothing
3951  else if (aCell.hasNumeric())
3952  rArray.push_back(GetCellValue(aAdr, aCell));
3953  }
3954  break;
3955  case svDoubleRef :
3956  case svRefList :
3957  {
3958  PopDoubleRef( aRange, nParam, nRefInList);
3959  if (nGlobalError != FormulaError::NONE)
3960  break;
3961 
3962  aRange.PutInOrder();
3963  SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3964  nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3965  rArray.reserve( rArray.size() + nCellCount);
3966 
3967  FormulaError nErr = FormulaError::NONE;
3968  double fCellVal;
3969  ScValueIterator aValIter( mrDoc, aRange, mnSubTotalFlags );
3970  if (aValIter.GetFirst( fCellVal, nErr))
3971  {
3972  if (bIgnoreErrVal)
3973  {
3974  if (nErr == FormulaError::NONE)
3975  rArray.push_back( fCellVal);
3976  while (aValIter.GetNext( fCellVal, nErr))
3977  {
3978  if (nErr == FormulaError::NONE)
3979  rArray.push_back( fCellVal);
3980  }
3981  }
3982  else
3983  {
3984  rArray.push_back( fCellVal);
3985  SetError(nErr);
3986  while ((nErr == FormulaError::NONE) && aValIter.GetNext( fCellVal, nErr))
3987  rArray.push_back( fCellVal);
3988  SetError(nErr);
3989  }
3990  }
3991  }
3992  break;
3993  case svMatrix :
3994  case svExternalSingleRef:
3995  case svExternalDoubleRef:
3996  {
3997  ScMatrixRef pMat = GetMatrix();
3998  if (!pMat)
3999  break;
4000 
4001  SCSIZE nCount = pMat->GetElementCount();
4002  rArray.reserve( rArray.size() + nCount);
4003  if (pMat->IsNumeric())
4004  {
4005  if (bIgnoreErrVal)
4006  {
4007  for (SCSIZE i = 0; i < nCount; ++i)
4008  {
4009  const double fVal = pMat->GetDouble(i);
4010  if (nGlobalError == FormulaError::NONE)
4011  rArray.push_back( fVal);
4012  else
4013  nGlobalError = FormulaError::NONE;
4014  }
4015  }
4016  else
4017  {
4018  for (SCSIZE i = 0; i < nCount; ++i)
4019  rArray.push_back( pMat->GetDouble(i));
4020  }
4021  }
4022  else if (bConvertTextInArray && eStackType == svMatrix)
4023  {
4024  for (SCSIZE i = 0; i < nCount; ++i)
4025  {
4026  if ( pMat->IsValue( i ) )
4027  {
4028  if (bIgnoreErrVal)
4029  {
4030  const double fVal = pMat->GetDouble(i);
4031  if (nGlobalError == FormulaError::NONE)
4032  rArray.push_back( fVal);
4033  else
4034  nGlobalError = FormulaError::NONE;
4035  }
4036  else
4037  rArray.push_back( pMat->GetDouble(i));
4038  }
4039  else
4040  {
4041  // tdf#88547 try to convert string to (date)value
4042  OUString aStr = pMat->GetString( i ).getString();
4043  if ( aStr.getLength() > 0 )
4044  {
4045  FormulaError nErr = nGlobalError;
4046  nGlobalError = FormulaError::NONE;
4047  double fVal = ConvertStringToValue( aStr );
4048  if ( nGlobalError == FormulaError::NONE )
4049  {
4050  rArray.push_back( fVal );
4051  nGlobalError = nErr;
4052  }
4053  else
4054  {
4055  if (!bIgnoreErrVal)
4056  rArray.push_back( CreateDoubleError( FormulaError::NoValue));
4057  // Propagate previous error if any, else
4058  // the current #VALUE! error, unless
4059  // ignoring error values.
4060  if (nErr != FormulaError::NONE)
4061  nGlobalError = nErr;
4062  else if (!bIgnoreErrVal)
4063  nGlobalError = FormulaError::NoValue;
4064  else
4065  nGlobalError = FormulaError::NONE;
4066  }
4067  }
4068  }
4069  }
4070  }
4071  else
4072  {
4073  if (bIgnoreErrVal)
4074  {
4075  for (SCSIZE i = 0; i < nCount; ++i)
4076  {
4077  if (pMat->IsValue(i))
4078  {
4079  const double fVal = pMat->GetDouble(i);
4080  if (nGlobalError == FormulaError::NONE)
4081  rArray.push_back( fVal);
4082  else
4083  nGlobalError = FormulaError::NONE;
4084  }
4085  }
4086  }
4087  else
4088  {
4089  for (SCSIZE i = 0; i < nCount; ++i)
4090  {
4091  if (pMat->IsValue(i))
4092  rArray.push_back( pMat->GetDouble(i));
4093  }
4094  }
4095  }
4096  }
4097  break;
4098  default :
4099  PopError();
4100  SetError( FormulaError::IllegalParameter);
4101  break;
4102  }
4103  if (nGlobalError != FormulaError::NONE)
4104  break; // while
4105  }
4106  // nParam > 0 in case of error, clean stack environment and obtain earlier
4107  // error if there was one.
4108  while (nParam-- > 0)
4109  PopError();
4110 }
4111 
4112 void ScInterpreter::GetSortArray( sal_uInt8 nParamCount, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder, bool bConvertTextInArray, bool bAllowEmptyArray )
4113 {
4114  GetNumberSequenceArray( nParamCount, rSortArray, bConvertTextInArray );
4115  if (rSortArray.size() > MAX_COUNT_DOUBLE_FOR_SORT(mrDoc.GetSheetLimits()))
4116  SetError( FormulaError::MatrixSize);
4117  else if ( rSortArray.empty() )
4118  {
4119  if ( bAllowEmptyArray )
4120  return;
4121  SetError( FormulaError::NoValue);
4122  }
4123 
4124  if (nGlobalError == FormulaError::NONE)
4125  QuickSort( rSortArray, pIndexOrder);
4126 }
4127 
4128 static void lcl_QuickSort( tools::Long nLo, tools::Long nHi, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
4129 {
4130  // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().
4131 
4132  using ::std::swap;
4133 
4134  if (nHi - nLo == 1)
4135  {
4136  if (rSortArray[nLo] > rSortArray[nHi])
4137  {
4138  swap(rSortArray[nLo], rSortArray[nHi]);
4139  if (pIndexOrder)
4140  swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
4141  }
4142  return;
4143  }
4144 
4145  tools::Long ni = nLo;
4146  tools::Long nj = nHi;
4147  do
4148  {
4149  double fLo = rSortArray[nLo];
4150  while (ni <= nHi && rSortArray[ni] < fLo) ni++;
4151  while (nj >= nLo && fLo < rSortArray[nj]) nj--;
4152  if (ni <= nj)
4153  {
4154  if (ni != nj)
4155  {
4156  swap(rSortArray[ni], rSortArray[nj]);
4157  if (pIndexOrder)
4158  swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
4159  }
4160 
4161  ++ni;
4162  --nj;
4163  }
4164  }
4165  while (ni < nj);
4166 
4167  if ((nj - nLo) < (nHi - ni))
4168  {
4169  if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
4170  if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
4171  }
4172  else
4173  {
4174  if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
4175  if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
4176  }
4177 }
4178 
4179 void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
4180 {
4181  tools::Long n = static_cast<tools::Long>(rSortArray.size());
4182 
4183  if (pIndexOrder)
4184  {
4185  pIndexOrder->clear();
4186  pIndexOrder->reserve(n);
4187  for (tools::Long i = 0; i < n; ++i)
4188  pIndexOrder->push_back(i);
4189  }
4190 
4191  if (n < 2)
4192  return;
4193 
4194  size_t nValCount = rSortArray.size();
4195  for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
4196  {
4197  size_t nInd = comphelper::rng::uniform_size_distribution(0, nValCount-2);
4198  ::std::swap( rSortArray[i], rSortArray[nInd]);
4199  if (pIndexOrder)
4200  ::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
4201  }
4202 
4203  lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
4204 }
4205 
4206 void ScInterpreter::ScRank( bool bAverage )
4207 {
4208  sal_uInt8 nParamCount = GetByte();
4209  if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
4210  return;
4211  bool bAscending;
4212  if ( nParamCount == 3 )
4213  bAscending = GetBool();
4214  else
4215  bAscending = false;
4216 
4217  vector<double> aSortArray;
4218  GetSortArray( 1, aSortArray, nullptr, false, false );
4219  double fVal = GetDouble();
4220  SCSIZE nSize = aSortArray.size();
4221  if ( nSize == 0 || nGlobalError != FormulaError::NONE )
4222  PushNoValue();
4223  else
4224  {
4225  if ( fVal < aSortArray[ 0 ] || fVal > aSortArray[ nSize - 1 ] )
4226  PushNoValue();
4227  else
4228  {
4229  double fLastPos = 0;
4230  double fFirstPos = -1.0;
4231  bool bFinished = false;
4232  SCSIZE i;
4233  for (i = 0; i < nSize && !bFinished; i++)
4234  {
4235  if ( aSortArray[ i ] == fVal )
4236  {
4237  if ( fFirstPos < 0 )
4238  fFirstPos = i + 1.0;
4239  }
4240  else
4241  {
4242  if ( aSortArray[ i ] > fVal )
4243  {
4244  fLastPos = i;
4245  bFinished = true;
4246  }
4247  }
4248  }
4249  if ( !bFinished )
4250  fLastPos = i;
4251  if ( !bAverage )
4252  {
4253  if ( bAscending )
4254  PushDouble( fFirstPos );
4255  else
4256  PushDouble( nSize + 1.0 - fLastPos );
4257  }
4258  else
4259  {
4260  if ( bAscending )
4261  PushDouble( ( fFirstPos + fLastPos ) / 2.0 );
4262  else
4263  PushDouble( nSize + 1.0 - ( fFirstPos + fLastPos ) / 2.0 );
4264  }
4265  }
4266  }
4267 }
4268 
4270 {
4271  sal_uInt8 nParamCount = GetByte();
4272  if ( !MustHaveParamCountMin( nParamCount, 1 ) )
4273  return;
4274  sal_uInt16 SaveSP = sp;
4275  double nMiddle = 0.0;
4276  KahanSum rVal = 0.0;
4277  double rValCount = 0.0;
4278  ScAddress aAdr;
4279  ScRange aRange;
4280  short nParam = nParamCount;
4281  size_t nRefInList = 0;
4282  while (nParam-- > 0)
4283  {
4284  switch (GetStackType())
4285  {
4286  case svDouble :
4287  rVal += GetDouble();
4288  rValCount++;
4289  break;
4290  case svSingleRef :
4291  {
4292  PopSingleRef( aAdr );
4293  ScRefCellValue aCell(mrDoc, aAdr);
4294  if (aCell.hasNumeric())
4295  {
4296  rVal += GetCellValue(aAdr, aCell);
4297  rValCount++;
4298  }
4299  }
4300  break;
4301  case svDoubleRef :
4302  case svRefList :
4303  {
4304  FormulaError nErr = FormulaError::NONE;
4305  double nCellVal;
4306  PopDoubleRef( aRange, nParam, nRefInList);
4307  ScValueIterator aValIter( mrDoc, aRange, mnSubTotalFlags );
4308  if (aValIter.GetFirst(nCellVal, nErr))
4309  {
4310  rVal += nCellVal;
4311  rValCount++;
4312  SetError(nErr);
4313  while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
4314  {
4315  rVal += nCellVal;
4316  rValCount++;
4317  }
4318  SetError(nErr);
4319  }
4320  }
4321  break;
4322  case svMatrix :
4323  case svExternalSingleRef:
4324  case svExternalDoubleRef:
4325  {
4326  ScMatrixRef pMat = GetMatrix();
4327  if (pMat)
4328  {
4329  SCSIZE nCount = pMat->GetElementCount();
4330  if (pMat->IsNumeric())
4331  {
4332  for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4333  {
4334  rVal += pMat->GetDouble(nElem);
4335  rValCount++;
4336  }
4337  }
4338  else
4339  {
4340  for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4341  if (!pMat->IsStringOrEmpty(nElem))
4342  {
4343  rVal += pMat->GetDouble(nElem);
4344  rValCount++;
4345  }
4346  }
4347  }
4348  }
4349  break;
4350  default :
4351  SetError(FormulaError::IllegalParameter);
4352  break;
4353  }
4354  }
4355  if (nGlobalError != FormulaError::NONE)
4356  {
4357  PushError( nGlobalError);
4358  return;
4359  }
4360  nMiddle = rVal.get() / rValCount;
4361  sp = SaveSP;
4362  rVal = 0.0;
4363  nParam = nParamCount;
4364  nRefInList = 0;
4365  while (nParam-- > 0)
4366  {
4367  switch (GetStackType())
4368  {
4369  case svDouble :
4370  rVal += std::abs(GetDouble() - nMiddle);
4371  break;
4372  case svSingleRef :
4373  {
4374  PopSingleRef( aAdr );
4375  ScRefCellValue aCell(mrDoc, aAdr);
4376  if (aCell.hasNumeric())
4377  rVal += std::abs(GetCellValue(aAdr, aCell) - nMiddle);
4378  }
4379  break;
4380  case svDoubleRef :
4381  case svRefList :
4382  {
4383  FormulaError nErr = FormulaError::NONE;
4384  double nCellVal;
4385  PopDoubleRef( aRange, nParam, nRefInList);
4386  ScValueIterator aValIter( mrDoc, aRange, mnSubTotalFlags );
4387  if (aValIter.GetFirst(nCellVal, nErr))
4388  {
4389  rVal += std::abs(nCellVal - nMiddle);
4390  while (aValIter.GetNext(nCellVal, nErr))
4391  rVal += std::abs(nCellVal - nMiddle);
4392  }
4393  }
4394  break;
4395  case svMatrix :
4396  case svExternalSingleRef:
4397  case svExternalDoubleRef:
4398  {
4399  ScMatrixRef pMat = GetMatrix();
4400  if (pMat)
4401  {
4402  SCSIZE nCount = pMat->GetElementCount();
4403  if (pMat->IsNumeric())
4404  {
4405  for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4406  {
4407  rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
4408  }
4409  }
4410  else
4411  {
4412  for (SCSIZE nElem = 0; nElem < nCount; nElem++)
4413  {
4414  if (!pMat->IsStringOrEmpty(nElem))
4415  rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
4416  }
4417  }
4418  }
4419  }
4420  break;
4421  default : SetError(FormulaError::IllegalParameter); break;
4422  }
4423  }
4424  PushDouble(rVal.get() / rValCount);
4425 }
4426 
4428 {
4429  auto VarResult = []( double fVal, size_t /*nValCount*/ )
4430  {
4431  return fVal;
4432  };
4433  GetStVarParams( false /*bTextAsZero*/, VarResult);
4434 }
4435 
4437 {
4438  sal_uInt8 nParamCount = GetByte();
4439  if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
4440  return;
4441  double fUp, fLo;
4442  fUp = GetDouble();
4443  if (nParamCount == 4)
4444  fLo = GetDouble();
4445  else
4446  fLo = fUp;
4447  if (fLo > fUp)
4448  {
4449  double fTemp = fLo;
4450  fLo = fUp;
4451  fUp = fTemp;
4452  }
4453  ScMatrixRef pMatP = GetMatrix();
4454  ScMatrixRef pMatW = GetMatrix();
4455  if (!pMatP || !pMatW)
4456  PushIllegalParameter();
4457  else
4458  {
4459  SCSIZE nC1, nC2;
4460  SCSIZE nR1, nR2;
4461  pMatP->GetDimensions(nC1, nR1);
4462  pMatW->GetDimensions(nC2, nR2);
4463  if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
4464  nC2 == 0 || nR2 == 0)
4465  PushNA();
4466  else
4467  {
4468  KahanSum fSum = 0.0;
4469  KahanSum fRes = 0.0;
4470  bool bStop = false;
4471  double fP, fW;
4472  for ( SCSIZE i = 0; i < nC1 && !bStop; i++ )
4473  {
4474  for (SCSIZE j = 0; j < nR1 && !bStop; ++j )
4475  {
4476  if (pMatP->IsValue(i,j) && pMatW->IsValue(i,j))
4477  {
4478  fP = pMatP->GetDouble(i,j);
4479  fW = pMatW->GetDouble(i,j);
4480  if (fP < 0.0 || fP > 1.0)
4481  bStop = true;
4482  else
4483  {
4484  fSum += fP;
4485  if (fW >= fLo && fW <= fUp)
4486  fRes += fP;
4487  }
4488  }
4489  else
4490  SetError( FormulaError::IllegalArgument);
4491  }
4492  }
4493  if (bStop || std::abs((fSum -1.0).get()) > 1.0E-7)
4494  PushNoValue();
4495  else
4496  PushDouble(fRes.get());
4497  }
4498  }
4499 }
4500 
4502 {
4503  // This is identical to ScPearson()
4504  ScPearson();
4505 }
4506 
4508 {
4509  CalculatePearsonCovar( false, false, false );
4510 }
4511 
4513 {
4514  CalculatePearsonCovar( false, false, true );
4515 }
4516 
4518 {
4519  CalculatePearsonCovar( true, false, false );
4520 }
4521 
4522 void ScInterpreter::CalculatePearsonCovar( bool _bPearson, bool _bStexy, bool _bSample )
4523 {
4524  if ( !MustHaveParamCount( GetByte(), 2 ) )
4525  return;
4526  ScMatrixRef pMat1 = GetMatrix();
4527  ScMatrixRef pMat2 = GetMatrix();
4528  if (!pMat1 || !pMat2)
4529  {
4530  PushIllegalParameter();
4531  return;
4532  }
4533  SCSIZE nC1, nC2;
4534  SCSIZE nR1, nR2;
4535  pMat1->GetDimensions(nC1, nR1);
4536  pMat2->GetDimensions(nC2, nR2);
4537  if (nR1 != nR2 || nC1 != nC2)
4538  {
4539  PushIllegalArgument();
4540  return;
4541  }
4542  /* #i78250#
4543  * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
4544  * but the latter produces wrong results if the absolute values are high,
4545  * for example above 10^8
4546  */
4547  double fCount = 0.0;
4548  KahanSum fSumX = 0.0;
4549  KahanSum fSumY = 0.0;
4550 
4551  for (SCSIZE i = 0; i < nC1; i++)
4552  {
4553  for (SCSIZE j = 0; j < nR1; j++)
4554  {
4555  if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4556  {
4557  fSumX += pMat1->GetDouble(i,j);
4558  fSumY += pMat2->GetDouble(i,j);
4559  fCount++;
4560  }
4561  }
4562  }
4563  if (fCount < (_bStexy ? 3.0 : (_bSample ? 2.0 : 1.0)))
4564  PushNoValue();
4565  else
4566  {
4567  KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4568  KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4569  KahanSum fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2
4570  const double fMeanX = fSumX.get() / fCount;
4571  const double fMeanY = fSumY.get() / fCount;
4572  for (SCSIZE i = 0; i < nC1; i++)
4573  {
4574  for (SCSIZE j = 0; j < nR1; j++)
4575  {
4576  if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4577  {
4578  const double fValX = pMat1->GetDouble(i,j);
4579  const double fValY = pMat2->GetDouble(i,j);
4580  fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4581  if ( _bPearson )
4582  {
4583  fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4584  fSumSqrDeltaY += (fValY - fMeanY) * (fValY - fMeanY);
4585  }
4586  }
4587  }
4588  }
4589  if ( _bPearson )
4590  {
4591  // tdf#94962 - Values below the numerical limit lead to unexpected results
4592  if (fSumSqrDeltaX < ::std::numeric_limits<double>::min()
4593  || (!_bStexy && fSumSqrDeltaY < ::std::numeric_limits<double>::min()))
4594  PushError( FormulaError::DivisionByZero);
4595  else if ( _bStexy )
4596  PushDouble( sqrt( ( fSumSqrDeltaY - fSumDeltaXDeltaY *
4597  fSumDeltaXDeltaY / fSumSqrDeltaX ).get() / (fCount-2)));
4598  else
4599  PushDouble( fSumDeltaXDeltaY.get() / sqrt( fSumSqrDeltaX.get() * fSumSqrDeltaY.get() ));
4600  }
4601  else
4602  {
4603  if ( _bSample )
4604  PushDouble( fSumDeltaXDeltaY.get() / ( fCount - 1 ) );
4605  else
4606  PushDouble( fSumDeltaXDeltaY.get() / fCount);
4607  }
4608  }
4609 }
4610 
4612 {
4613  // Same as ScPearson()*ScPearson()
4614  ScPearson();
4615  if (nGlobalError != FormulaError::NONE)
4616  return;
4617 
4618  switch (GetStackType())
4619  {
4620  case svDouble:
4621  {
4622  double fVal = PopDouble();
4623  PushDouble( fVal * fVal);
4624  }
4625  break;
4626  default:
4627  PopError();
4628  PushNoValue();
4629  }
4630 }
4631 
4633 {
4634  CalculatePearsonCovar( true, true, false );
4635 }
4637 {
4638  if ( !MustHaveParamCount( GetByte(), 2 ) )
4639  return;
4640  ScMatrixRef pMat1 = GetMatrix();
4641  ScMatrixRef pMat2 = GetMatrix();
4642  if (!pMat1 || !pMat2)
4643  {
4644  PushIllegalParameter();
4645  return;
4646  }
4647  SCSIZE nC1, nC2;
4648  SCSIZE nR1, nR2;
4649  pMat1->GetDimensions(nC1, nR1);
4650  pMat2->GetDimensions(nC2, nR2);
4651  if (nR1 != nR2 || nC1 != nC2)
4652  {
4653  PushIllegalArgument();
4654  return;
4655  }
4656  // #i78250# numerical stability improved
4657  double fCount = 0.0;
4658  KahanSum fSumX = 0.0;
4659  KahanSum fSumY = 0.0;
4660 
4661  for (SCSIZE i = 0; i < nC1; i++)
4662  {
4663  for (SCSIZE j = 0; j < nR1; j++)
4664  {
4665  if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4666  {
4667  fSumX += pMat1->GetDouble(i,j);
4668  fSumY += pMat2->GetDouble(i,j);
4669  fCount++;
4670  }
4671  }
4672  }
4673  if (fCount < 1.0)
4674  PushNoValue();
4675  else
4676  {
4677  KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4678  KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4679  double fMeanX = fSumX.get() / fCount;
4680  double fMeanY = fSumY.get() / fCount;
4681  for (SCSIZE i = 0; i < nC1; i++)
4682  {
4683  for (SCSIZE j = 0; j < nR1; j++)
4684  {
4685  if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4686  {
4687  double fValX = pMat1->GetDouble(i,j);
4688  double fValY = pMat2->GetDouble(i,j);
4689  fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4690  fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4691  }
4692  }
4693  }
4694  if (fSumSqrDeltaX == 0.0)
4695  PushError( FormulaError::DivisionByZero);
4696  else
4697  {
4698  if ( bSlope )
4699  PushDouble( fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get());
4700  else
4701  PushDouble( fMeanY - fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * fMeanX);
4702  }
4703  }
4704 }
4705 
4707 {
4708  CalculateSlopeIntercept(true);
4709 }
4710 
4712 {
4713  CalculateSlopeIntercept(false);
4714 }
4715 
4717 {
4718  if ( !MustHaveParamCount( GetByte(), 3 ) )
4719  return;
4720  ScMatrixRef pMat1 = GetMatrix();
4721  ScMatrixRef pMat2 = GetMatrix();
4722  if (!pMat1 || !pMat2)
4723  {
4724  PushIllegalParameter();
4725  return;
4726  }
4727  SCSIZE nC1, nC2;
4728  SCSIZE nR1, nR2;
4729  pMat1->GetDimensions(nC1, nR1);
4730  pMat2->GetDimensions(nC2, nR2);
4731  if (nR1 != nR2 || nC1 != nC2)
4732  {
4733  PushIllegalArgument();
4734  return;
4735  }
4736  double fVal = GetDouble();
4737  // #i78250# numerical stability improved
4738  double fCount = 0.0;
4739  KahanSum fSumX = 0.0;
4740  KahanSum fSumY = 0.0;
4741 
4742  for (SCSIZE i = 0; i < nC1; i++)
4743  {
4744  for (SCSIZE j = 0; j < nR1; j++)
4745  {
4746  if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4747  {
4748  fSumX += pMat1->GetDouble(i,j);
4749  fSumY += pMat2->GetDouble(i,j);
4750  fCount++;
4751  }
4752  }
4753  }
4754  if (fCount < 1.0)
4755  PushNoValue();
4756  else
4757  {
4758  KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4759  KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4760  double fMeanX = fSumX.get() / fCount;
4761  double fMeanY = fSumY.get() / fCount;
4762  for (SCSIZE i = 0; i < nC1; i++)
4763  {
4764  for (SCSIZE j = 0; j < nR1; j++)
4765  {
4766  if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
4767  {
4768  double fValX = pMat1->GetDouble(i,j);
4769  double fValY = pMat2->GetDouble(i,j);
4770  fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4771  fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4772  }
4773  }
4774  }
4775  if (fSumSqrDeltaX == 0.0)
4776  PushError( FormulaError::DivisionByZero);
4777  else
4778  PushDouble( fMeanY + fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * (fVal - fMeanX));
4779  }
4780 }
4781 
4782 static void lcl_roundUpNearestPow2(SCSIZE& nNum, SCSIZE& nNumBits)
4783 {
4784  // Find the least power of 2 that is less than or equal to nNum.
4785  SCSIZE nPow2(1);
4786  nNumBits = std::numeric_limits<SCSIZE>::digits;
4787  nPow2 <<= (nNumBits - 1);
4788  while (nPow2 >= 1)
4789  {
4790  if (nNum & nPow2)
4791  break;
4792 
4793  --nNumBits;
4794  nPow2 >>= 1;
4795  }
4796 
4797  if (nPow2 != nNum)
4798  nNum = nPow2 ? (nPow2 << 1) : 1;
4799  else
4800  --nNumBits;
4801 }
4802 
4803 static SCSIZE lcl_bitReverse(SCSIZE nIn, SCSIZE nBound)
4804 {
4805  SCSIZE nOut = 0;
4806  for (SCSIZE nMask = 1; nMask < nBound; nMask <<= 1)
4807  {
4808  nOut <<= 1;
4809 
4810  if (nIn & nMask)
4811  nOut |= 1;
4812  }
4813 
4814  return nOut;
4815 }
4816 
4817 namespace {
4818 
4819 // Computes and stores twiddle factors for computing DFT later.
4820 struct ScTwiddleFactors
4821 {
4822  ScTwiddleFactors(SCSIZE nN, bool bInverse) :
4823  mfWReal(nN),
4824  mfWImag(nN),
4825  mnN(nN),
4826  mbInverse(bInverse)
4827  {}
4828 
4829  void Compute();
4830 
4831  void Conjugate()
4832  {
4833  mbInverse = !mbInverse;
4834  for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
4835  mfWImag[nIdx] = -mfWImag[nIdx];
4836  }
4837 
4838  std::vector<double> mfWReal;
4839  std::vector<double> mfWImag;
4840  SCSIZE mnN;
4841  bool mbInverse;
4842 };
4843 
4844 }
4845 
4846 void ScTwiddleFactors::Compute()
4847 {
4848  mfWReal.resize(mnN);
4849  mfWImag.resize(mnN);
4850 
4851  double nW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnN);
4852 
4853  if (mnN == 1)
4854  {
4855  mfWReal[0] = 1.0;
4856  mfWImag[0] = 0.0;
4857  }
4858  else if (mnN == 2)
4859  {
4860  mfWReal[0] = 1;
4861  mfWImag[0] = 0;
4862 
4863  mfWReal[1] = -1;
4864  mfWImag[1] = 0;
4865  }
4866  else if (mnN == 4)
4867  {
4868  mfWReal[0] = 1;
4869  mfWImag[0] = 0;
4870 
4871  mfWReal[1] = 0;
4872  mfWImag[1] = (mbInverse ? 1.0 : -1.0);
4873 
4874  mfWReal[2] = -1;
4875  mfWImag[2] = 0;
4876 
4877  mfWReal[3] = 0;
4878  mfWImag[3] = (mbInverse ? -1.0 : 1.0);
4879  }
4880  else if ((mnN % 4) == 0)
4881  {
4882  const SCSIZE nQSize = mnN >> 2;
4883  // Compute cos of the start quadrant.
4884  // This is the first quadrant if mbInverse == true, else it is the fourth quadrant.
4885  for (SCSIZE nIdx = 0; nIdx <= nQSize; ++nIdx)
4886  mfWReal[nIdx] = cos(nW*static_cast<double>(nIdx));
4887 
4888  if (mbInverse)
4889  {
4890  const SCSIZE nQ1End = nQSize;
4891  // First quadrant
4892  for (SCSIZE nIdx = 0; nIdx <= nQ1End; ++nIdx)
4893  mfWImag[nIdx] = mfWReal[nQ1End-nIdx];
4894 
4895  // Second quadrant
4896  const SCSIZE nQ2End = nQ1End << 1;
4897  for (SCSIZE nIdx = nQ1End+1; nIdx <= nQ2End; ++nIdx)
4898  {
4899  mfWReal[nIdx] = -mfWReal[nQ2End - nIdx];
4900  mfWImag[nIdx] = mfWImag[nQ2End - nIdx];
4901  }
4902 
4903  // Third quadrant
4904  const SCSIZE nQ3End = nQ2End + nQ1End;
4905  for (SCSIZE nIdx = nQ2End+1; nIdx <= nQ3End; ++nIdx)
4906  {
4907  mfWReal[nIdx] = -mfWReal[nIdx - nQ2End];
4908  mfWImag[nIdx] = -mfWImag[nIdx - nQ2End];
4909  }
4910 
4911  // Fourth Quadrant
4912  for (SCSIZE nIdx = nQ3End+1; nIdx < mnN; ++nIdx)
4913  {
4914  mfWReal[nIdx] = mfWReal[mnN - nIdx];
4915  mfWImag[nIdx] = -mfWImag[mnN - nIdx];
4916  }
4917  }
4918  else
4919  {
4920  const SCSIZE nQ4End = nQSize;
4921  const SCSIZE nQ3End = nQSize << 1;
4922  const SCSIZE nQ2End = nQ3End + nQSize;
4923 
4924  // Fourth quadrant.
4925  for (SCSIZE nIdx = 0; nIdx <= nQ4End; ++nIdx)
4926  mfWImag[nIdx] = -mfWReal[nQ4End - nIdx];
4927 
4928  // Third quadrant.
4929  for (SCSIZE nIdx = nQ4End+1; nIdx <= nQ3End; ++nIdx)
4930  {
4931  mfWReal[nIdx] = -mfWReal[nQ3End - nIdx];
4932  mfWImag[nIdx] = mfWImag[nQ3End - nIdx];
4933  }
4934 
4935  // Second quadrant.
4936  for (SCSIZE nIdx = nQ3End+1; nIdx <= nQ2End; ++nIdx)
4937  {
4938  mfWReal[nIdx] = -mfWReal[nIdx - nQ3End];
4939  mfWImag[nIdx] = -mfWImag[nIdx - nQ3End];
4940  }
4941 
4942  // First quadrant.
4943  for (SCSIZE nIdx = nQ2End+1; nIdx < mnN; ++nIdx)
4944  {
4945  mfWReal[nIdx] = mfWReal[mnN - nIdx];
4946  mfWImag[nIdx] = -mfWImag[mnN - nIdx];
4947  }
4948  }
4949  }
4950  else
4951  {
4952  for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
4953  {
4954  double fAngle = nW*static_cast<double>(nIdx);
4955  mfWReal[nIdx] = cos(fAngle);
4956  mfWImag[nIdx] = sin(fAngle);
4957  }
4958  }
4959 }
4960 
4961 namespace {
4962 
4963 // A radix-2 decimation in time FFT algorithm for complex valued input.
4964 class ScComplexFFT2
4965 {
4966 public:
4967  // rfArray.size() would always be even and a power of 2. (asserted in prepare())
4968  // rfArray's first half contains the real parts and the later half contains the imaginary parts.
4969  ScComplexFFT2(std::vector<double>& raArray, bool bInverse, bool bPolar, double fMinMag,
4970  ScTwiddleFactors& rTF, bool bSubSampleTFs = false, bool bDisableNormalize = false) :
4971  mrArray(raArray),
4972  mfWReal(rTF.mfWReal),
4973  mfWImag(rTF.mfWImag),
4974  mnPoints(raArray.size()/2),
4975  mnStages(0),
4976  mfMinMag(fMinMag),
4977  mbInverse(bInverse),
4978  mbPolar(bPolar),
4979  mbDisableNormalize(bDisableNormalize),
4980  mbSubSampleTFs(bSubSampleTFs)
4981  {}
4982 
4983  void Compute();
4984 
4985 private:
4986 
4987  void prepare();
4988 
4989  double getReal(SCSIZE nIdx)
4990  {
4991  return mrArray[nIdx];
4992  }
4993 
4994  void setReal(double fVal, SCSIZE nIdx)
4995  {
4996  mrArray[nIdx] = fVal;
4997  }
4998 
4999  double getImag(SCSIZE nIdx)
5000  {
5001  return mrArray[mnPoints + nIdx];
5002  }
5003 
5004  void setImag(double fVal, SCSIZE nIdx)
5005  {
5006  mrArray[mnPoints + nIdx] = fVal;
5007  }
5008 
5009  SCSIZE getTFactorIndex(SCSIZE nPtIndex, SCSIZE nTfIdxScaleBits)
5010  {
5011  return ( ( nPtIndex << nTfIdxScaleBits ) & ( mnPoints - 1 ) ); // (x & (N-1)) is same as (x % N) but faster.
5012  }
5013 
5014  void computeFly(SCSIZE nTopIdx, SCSIZE nBottomIdx, SCSIZE nWIdx1, SCSIZE nWIdx2)
5015  {
5016  if (mbSubSampleTFs)
5017  {
5018  nWIdx1 <<= 1;
5019  nWIdx2 <<= 1;
5020  }
5021 
5022  const double x1r = getReal(nTopIdx);
5023  const double x2r = getReal(nBottomIdx);
5024 
5025  const double& w1r = mfWReal[nWIdx1];
5026  const double& w1i = mfWImag[nWIdx1];
5027 
5028  const double& w2r = mfWReal[nWIdx2];
5029  const double& w2i = mfWImag[nWIdx2];
5030 
5031  const double x1i = getImag(nTopIdx);
5032  const double x2i = getImag(nBottomIdx);
5033 
5034  setReal(x1r + x2r*w1r - x2i*w1i, nTopIdx);
5035  setImag(x1i + x2i*w1r + x2r*w1i, nTopIdx);
5036 
5037  setReal(x1r + x2r*w2r - x2i*w2i, nBottomIdx);
5038  setImag(x1i + x2i*w2r + x2r*w2i, nBottomIdx);
5039  }
5040 
5041  std::vector<double>& mrArray;
5042  std::vector<double>& mfWReal;
5043  std::vector<double>& mfWImag;
5044  SCSIZE mnPoints;
5045  SCSIZE mnStages;
5046  double mfMinMag;
5047  bool mbInverse:1;
5048  bool mbPolar:1;
5049  bool mbDisableNormalize:1;
5050  bool mbSubSampleTFs:1;
5051 };
5052 
5053 }
5054 
5055 void ScComplexFFT2::prepare()
5056 {
5057  SCSIZE nPoints = mnPoints;
5058  lcl_roundUpNearestPow2(nPoints, mnStages);
5059  assert(nPoints == mnPoints);
5060 
5061  // Reorder array by bit-reversed indices.
5062  for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5063  {
5064  SCSIZE nRevIdx = lcl_bitReverse(nIdx, mnPoints);
5065  if (nIdx < nRevIdx)
5066  {
5067  double fTmp = getReal(nIdx);
5068  setReal(getReal(nRevIdx), nIdx);
5069  setReal(fTmp, nRevIdx);
5070 
5071  fTmp = getImag(nIdx);
5072  setImag(getImag(nRevIdx), nIdx);
5073  setImag(fTmp, nRevIdx);
5074  }
5075  }
5076 }
5077 
5078 static void lcl_normalize(std::vector<double>& rCmplxArray, bool bScaleOnlyReal)
5079 {
5080  const SCSIZE nPoints = rCmplxArray.size()/2;
5081  const double fScale = 1.0/static_cast<double>(nPoints);
5082 
5083  // Scale the real part
5084  for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
5085  rCmplxArray[nIdx] *= fScale;
5086 
5087  if (!bScaleOnlyReal)
5088  {
5089  const SCSIZE nLen = nPoints*2;
5090  for (SCSIZE nIdx = nPoints; nIdx < nLen; ++nIdx)
5091  rCmplxArray[nIdx] *= fScale;
5092  }
5093 }
5094 
5095 static void lcl_convertToPolar(std::vector<double>& rCmplxArray, double fMinMag)
5096 {
5097  const SCSIZE nPoints = rCmplxArray.size()/2;
5098  double fMag, fPhase, fR, fI;
5099  for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
5100  {
5101  fR = rCmplxArray[nIdx];
5102  fI = rCmplxArray[nPoints+nIdx];
5103  fMag = sqrt(fR*fR + fI*fI);
5104  if (fMag < fMinMag)
5105  {
5106  fMag = 0.0;
5107  fPhase = 0.0;
5108  }
5109  else
5110  {
5111  fPhase = atan2(fI, fR);
5112  }
5113 
5114  rCmplxArray[nIdx] = fMag;
5115  rCmplxArray[nPoints+nIdx] = fPhase;
5116  }
5117 }
5118 
5119 void ScComplexFFT2::Compute()
5120 {
5121  prepare();
5122 
5123  const SCSIZE nFliesInStage = mnPoints/2;
5124  for (SCSIZE nStage = 0; nStage < mnStages; ++nStage)
5125  {
5126  const SCSIZE nTFIdxScaleBits = mnStages - nStage - 1; // Twiddle factor index's scale factor in bits.
5127  const SCSIZE nFliesInGroup = SCSIZE(1) << nStage;
5128  const SCSIZE nGroups = nFliesInStage/nFliesInGroup;
5129  const SCSIZE nFlyWidth = nFliesInGroup;
5130  for (SCSIZE nGroup = 0, nFlyTopIdx = 0; nGroup < nGroups; ++nGroup)
5131  {
5132  for (SCSIZE nFly = 0; nFly < nFliesInGroup; ++nFly, ++nFlyTopIdx)
5133  {
5134  SCSIZE nFlyBottomIdx = nFlyTopIdx + nFlyWidth;
5135  SCSIZE nWIdx1 = getTFactorIndex(nFlyTopIdx, nTFIdxScaleBits);
5136  SCSIZE nWIdx2 = getTFactorIndex(nFlyBottomIdx, nTFIdxScaleBits);
5137 
5138  computeFly(nFlyTopIdx, nFlyBottomIdx, nWIdx1, nWIdx2);
5139  }
5140 
5141  nFlyTopIdx += nFlyWidth;
5142  }
5143  }
5144 
5145  if (mbPolar)
5146  lcl_convertToPolar(mrArray, mfMinMag);
5147 
5148  // Normalize after converting to polar, so we have a chance to
5149  // save O(mnPoints) flops.
5150  if (mbInverse && !mbDisableNormalize)
5151  lcl_normalize(mrArray, mbPolar);
5152 }
5153 
5154 namespace {
5155 
5156 // Bluestein's algorithm or chirp z-transform algorithm that can be used to
5157 // compute DFT of a complex valued input of any length N in O(N lgN) time.
5158 class ScComplexBluesteinFFT
5159 {
5160 public:
5161 
5162  ScComplexBluesteinFFT(std::vector<double>& rArray, bool bReal, bool bInverse,
5163  bool bPolar, double fMinMag, bool bDisableNormalize = false) :
5164  mrArray(rArray),
5165  mnPoints(rArray.size()/2), // rArray should have space for imaginary parts even if real input.
5166  mfMinMag(fMinMag),
5167  mbReal(bReal),
5168  mbInverse(bInverse),
5169  mbPolar(bPolar),
5170  mbDisableNormalize(bDisableNormalize)
5171  {}
5172 
5173  void Compute();
5174 
5175 private:
5176  std::vector<double>& mrArray;
5177  const SCSIZE mnPoints;
5178  double mfMinMag;
5179  bool mbReal:1;
5180  bool mbInverse:1;
5181  bool mbPolar:1;
5182  bool mbDisableNormalize:1;
5183 };
5184 
5185 }
5186 
5187 void ScComplexBluesteinFFT::Compute()
5188 {
5189  std::vector<double> aRealScalars(mnPoints);
5190  std::vector<double> aImagScalars(mnPoints);
5191  double fW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnPoints);
5192  for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5193  {
5194  double fAngle = 0.5*fW*static_cast<double>(nIdx*nIdx);
5195  aRealScalars[nIdx] = cos(fAngle);
5196  aImagScalars[nIdx] = sin(fAngle);
5197  }
5198 
5199  SCSIZE nMinSize = mnPoints*2 - 1;
5200  SCSIZE nExtendedLength = nMinSize, nTmp = 0;
5201  lcl_roundUpNearestPow2(nExtendedLength, nTmp);
5202  std::vector<double> aASignal(nExtendedLength*2); // complex valued
5203  std::vector<double> aBSignal(nExtendedLength*2); // complex valued
5204 
5205  double fReal, fImag;
5206  for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5207  {
5208  // Real part of A signal.
5209  aASignal[nIdx] = mrArray[nIdx]*aRealScalars[nIdx] + (mbReal ? 0.0 : -mrArray[mnPoints+nIdx]*aImagScalars[nIdx]);
5210  // Imaginary part of A signal.
5211  aASignal[nExtendedLength + nIdx] = mrArray[nIdx]*aImagScalars[nIdx] + (mbReal ? 0.0 : mrArray[mnPoints+nIdx]*aRealScalars[nIdx]);
5212 
5213  // Real part of B signal.
5214  aBSignal[nIdx] = fReal = aRealScalars[nIdx];
5215  // Imaginary part of B signal.
5216  aBSignal[nExtendedLength + nIdx] = fImag = -aImagScalars[nIdx]; // negative sign because B signal is the conjugation of the scalars.
5217 
5218  if (nIdx)
5219  {
5220  // B signal needs a mirror of its part in 0 < n < mnPoints at the tail end.
5221  aBSignal[nExtendedLength - nIdx] = fReal;
5222  aBSignal[(nExtendedLength<<1) - nIdx] = fImag;
5223  }
5224  }
5225 
5226  {
5227  ScTwiddleFactors aTF(nExtendedLength, false /*not inverse*/);
5228  aTF.Compute();
5229 
5230  // Do complex-FFT2 of both A and B signal.
5231  ScComplexFFT2 aFFT2A(aASignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
5232  aTF, false /*no subsample*/, true /* disable normalize */);
5233  aFFT2A.Compute();
5234 
5235  ScComplexFFT2 aFFT2B(aBSignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
5236  aTF, false /*no subsample*/, true /* disable normalize */);
5237  aFFT2B.Compute();
5238 
5239  double fAR, fAI, fBR, fBI;
5240  for (SCSIZE nIdx = 0; nIdx < nExtendedLength; ++nIdx)
5241  {
5242  fAR = aASignal[nIdx];
5243  fAI = aASignal[nExtendedLength + nIdx];
5244  fBR = aBSignal[nIdx];
5245  fBI = aBSignal[nExtendedLength + nIdx];
5246 
5247  // Do point-wise product inplace in A signal.
5248  aASignal[nIdx] = fAR*fBR - fAI*fBI;
5249  aASignal[nExtendedLength + nIdx] = fAR*fBI + fAI*fBR;
5250  }
5251 
5252  // Do complex-inverse-FFT2 of aASignal.
5253  aTF.Conjugate();
5254  ScComplexFFT2 aFFT2AI(aASignal, true /*inverse*/, false /*no polar*/, 0.0 /* no clipping */, aTF); // Need normalization here.
5255  aFFT2AI.Compute();
5256  }
5257 
5258  // Point-wise multiply with scalars.
5259  for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
5260  {
5261  fReal = aASignal[nIdx];
5262  fImag = aASignal[nExtendedLength + nIdx];
5263  mrArray[nIdx] = fReal*aRealScalars[nIdx] - fImag*aImagScalars[nIdx]; // no conjugation needed here.
5264  mrArray[mnPoints + nIdx] = fReal*aImagScalars[nIdx] + fImag*aRealScalars[nIdx];
5265  }
5266 
5267  // Normalize/Polar operations
5268  if (mbPolar)
5269  lcl_convertToPolar(mrArray, mfMinMag);
5270 
5271  // Normalize after converting to polar, so we have a chance to
5272  // save O(mnPoints) flops.
5273  if (mbInverse && !mbDisableNormalize)
5274  lcl_normalize(mrArray, mbPolar);
5275 }
5276 
5277 namespace {
5278 
5279 // Computes DFT of an even length(N) real-valued input by using a
5280 // ScComplexFFT2 if N == 2^k for some k or else by using a ScComplexBluesteinFFT
5281 // with a complex valued input of length = N/2.
5282 class ScRealFFT
5283 {
5284 public:
5285 
5286  ScRealFFT(std::vector<double>& rInArray, std::vector<double>& rOutArray, bool bInverse,
5287  bool bPolar, double fMinMag) :
5288  mrInArray(rInArray),
5289  mrOutArray(rOutArray),
5290  mfMinMag(fMinMag),
5291  mbInverse(bInverse),
5292  mbPolar(bPolar)
5293  {}
5294 
5295  void Compute();
5296 
5297 private:
5298  std::vector<double>& mrInArray;
5299  std::vector<double>& mrOutArray;
5300  double mfMinMag;
5301  bool mbInverse:1;
5302  bool mbPolar:1;
5303 };
5304 
5305 }
5306 
5307 void ScRealFFT::Compute()
5308 {
5309  // input length has to be even to do this optimization.
5310  assert(mrInArray.size() % 2 == 0);
5311  assert(mrInArray.size()*2 == mrOutArray.size());
5312  // nN is the number of points in the complex-fft input
5313  // which will be half of the number of points in real array.
5314  const SCSIZE nN = mrInArray.size()/2;
5315  if (nN == 0)
5316  {
5317  mrOutArray[0] = mrInArray[0];
5318  mrOutArray[1] = 0.0;
5319  return;
5320  }
5321 
5322  // work array should be the same length as mrInArray
5323  std::vector<double> aWorkArray(nN*2);
5324  for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
5325  {
5326  SCSIZE nDoubleIdx = 2*nIdx;
5327  // Use even elements as real part
5328  aWorkArray[nIdx] = mrInArray[nDoubleIdx];
5329  // and odd elements as imaginary part of the contrived complex sequence.
5330  aWorkArray[nN+nIdx] = mrInArray[nDoubleIdx+1];
5331  }
5332 
5333  ScTwiddleFactors aTFs(nN*2, mbInverse);
5334  aTFs.Compute();
5335  SCSIZE nNextPow2 = nN, nTmp = 0;
5336  lcl_roundUpNearestPow2(nNextPow2, nTmp);
5337 
5338  if (nNextPow2 == nN)
5339  {
5340  ScComplexFFT2 aFFT2(aWorkArray, mbInverse, false /*disable polar*/, 0.0 /* no clipping */,
5341  aTFs, true /*subsample tf*/, true /*disable normalize*/);
5342  aFFT2.Compute();
5343  }
5344  else
5345  {
5346  ScComplexBluesteinFFT aFFT(aWorkArray, false /*complex input*/, mbInverse, false /*disable polar*/,
5347  0.0 /* no clipping */, true /*disable normalize*/);
5348  aFFT.Compute();
5349  }
5350 
5351  // Post process aWorkArray to populate mrOutArray
5352 
5353  const SCSIZE nTwoN = 2*nN, nThreeN = 3*nN;
5354  double fY1R, fY2R, fY1I, fY2I, fResR, fResI, fWR, fWI;
5355  for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
5356  {
5357  const SCSIZE nIdxRev = nIdx ? (nN - nIdx) : 0;
5358  fY1R = aWorkArray[nIdx];
5359  fY2R = aWorkArray[nIdxRev];
5360  fY1I = aWorkArray[nN + nIdx];
5361  fY2I = aWorkArray[nN + nIdxRev];
5362  fWR = aTFs.mfWReal[nIdx];
5363  fWI = aTFs.mfWImag[nIdx];
5364 
5365  // mrOutArray has length = 4*nN
5366  // Real part of the final output (only half of the symmetry around Nyquist frequency)
5367  // Fills the first quarter.
5368  mrOutArray[nIdx] = fResR = 0.5*(
5369  fY1R + fY2R +
5370  fWR * (fY1I + fY2I) +
5371  fWI * (fY1R - fY2R) );
5372  // Imaginary part of the final output (only half of the symmetry around Nyquist frequency)
5373  // Fills the third quarter.
5374  mrOutArray[nTwoN + nIdx] = fResI = 0.5*(
5375  fY1I - fY2I +
5376  fWI * (fY1I + fY2I) -
5377  fWR * (fY1R - fY2R) );
5378 
5379  // Fill the missing 2 quarters using symmetry argument.
5380  if (nIdx)
5381  {
5382  // Fills the 2nd quarter.
5383  mrOutArray[nN + nIdxRev] = fResR;
5384  // Fills the 4th quarter.
5385  mrOutArray[nThreeN + nIdxRev] = -fResI;
5386  }
5387  else
5388  {
5389  mrOutArray[nN] = fY1R - fY1I;
5390  mrOutArray[nThreeN] = 0.0;
5391  }
5392  }
5393 
5394  // Normalize/Polar operations
5395  if (mbPolar)
5396  lcl_convertToPolar(mrOutArray, mfMinMag);
5397 
5398  // Normalize after converting to polar, so we have a chance to
5399  // save O(mnPoints) flops.
5400  if (mbInverse)
5401  lcl_normalize(mrOutArray, mbPolar);
5402 }
5403 
5404 using ScMatrixGenerator = ScMatrixRef(SCSIZE, SCSIZE, std::vector<double>&);
5405 
5406 namespace {
5407 
5408 // Generic FFT class that decides which FFT implementation to use.
5409 class ScFFT
5410 {
5411 public:
5412 
5413  ScFFT(ScMatrixRef& pMat, bool bReal, bool bInverse, bool bPolar, double fMinMag) :
5414  mpInputMat(pMat),
5415  mfMinMag(fMinMag),
5416  mbReal(bReal),
5417  mbInverse(bInverse),
5418  mbPolar(bPolar)
5419  {}
5420 
5421  ScMatrixRef Compute(const std::function<ScMatrixGenerator>& rMatGenFunc);
5422 
5423 private:
5424  ScMatrixRef& mpInputMat;
5425  double mfMinMag;
5426  bool mbReal:1;
5427  bool mbInverse:1;
5428  bool mbPolar:1;
5429 };
5430 
5431 }
5432 
5433 ScMatrixRef ScFFT::Compute(const std::function<ScMatrixGenerator>& rMatGenFunc)
5434 {
5435  std::vector<double> aArray;
5436  mpInputMat->GetDoubleArray(aArray);
5437  SCSIZE nPoints = mbReal ? aArray.size() : (aArray.size()/2);
5438  if (nPoints == 1)
5439  {
5440  std::vector<double> aOutArray(2);
5441  aOutArray[0] = aArray[0];
5442  aOutArray[1] = mbReal ? 0.0 : aArray[1];
5443  if (mbPolar)
5444  lcl_convertToPolar(aOutArray, mfMinMag);
5445  return rMatGenFunc(2, 1, aOutArray);
5446  }
5447 
5448  if (mbReal && (nPoints % 2) == 0)
5449  {
5450  std::vector<double> aOutArray(nPoints*2);
5451  ScRealFFT aFFT(aArray, aOutArray, mbInverse, mbPolar, mfMinMag);
5452  aFFT.Compute();
5453  return rMatGenFunc(2, nPoints, aOutArray);
5454  }
5455 
5456  SCSIZE nNextPow2 = nPoints, nTmp = 0;
5457  lcl_roundUpNearestPow2(nNextPow2, nTmp);
5458  if (nNextPow2 == nPoints && !mbReal)
5459  {
5460  ScTwiddleFactors aTF(nPoints, mbInverse);
5461  aTF.Compute();
5462  ScComplexFFT2 aFFT2(aArray, mbInverse, mbPolar, mfMinMag, aTF);
5463  aFFT2.Compute();
5464  return rMatGenFunc(2, nPoints, aArray);
5465  }
5466 
5467  if (mbReal)
5468  aArray.resize(nPoints*2, 0.0);
5469  ScComplexBluesteinFFT aFFT(aArray, mbReal, mbInverse, mbPolar, mfMinMag);
5470  aFFT.Compute();
5471  return rMatGenFunc(2, nPoints, aArray);
5472 }
5473 
5475 {
5476  sal_uInt8 nParamCount = GetByte();
5477  if ( !MustHaveParamCount( nParamCount, 2, 5 ) )
5478  return;
5479 
5480  bool bInverse = false;
5481  bool bPolar = false;
5482  double fMinMag = 0.0;
5483 
5484  if (nParamCount == 5)
5485  {
5486  if (IsMissing())
5487  Pop();
5488  else
5489  fMinMag = GetDouble();
5490  }
5491 
5492  if (nParamCount >= 4)
5493  {
5494  if (IsMissing())
5495  Pop();
5496  else
5497  bPolar = GetBool();
5498  }
5499 
5500  if (nParamCount >= 3)
5501  {
5502  if (IsMissing())
5503  Pop();
5504  else
5505  bInverse = GetBool();
5506  }
5507 
5508  bool bGroupedByColumn = GetBool();
5509 
5510  ScMatrixRef pInputMat = GetMatrix();
5511  if (!pInputMat)
5512  {
5513  PushIllegalParameter();
5514  return;
5515  }
5516 
5517  SCSIZE nC, nR;
5518  pInputMat->GetDimensions(nC, nR);
5519 
5520  if ((bGroupedByColumn && nC > 2) || (!bGroupedByColumn && nR > 2))
5521  {
5522  // There can be no more than 2 columns (real, imaginary) if data grouped by columns.
5523  // and no more than 2 rows if data is grouped by rows.
5524  PushIllegalArgument();
5525  return;
5526  }
5527 
5528  if (!pInputMat->IsNumeric())
5529  {
5530  PushNoValue();
5531  return;
5532  }
5533 
5534  bool bRealInput = true;
5535  if (!bGroupedByColumn)
5536  {
5537  pInputMat->MatTrans(*pInputMat);
5538  bRealInput = (nR == 1);
5539  }
5540  else
5541  {
5542  bRealInput = (nC == 1);
5543  }
5544 
5545  ScFFT aFFT(pInputMat, bRealInput, bInverse, bPolar, fMinMag);
5546  std::function<ScMatrixGenerator> aFunc = [this](SCSIZE nCol, SCSIZE nRow, std::vector<double>& rVec) -> ScMatrixRef
5547  {
5548  return this->GetNewMat(nCol, nRow, rVec);
5549  };
5550  ScMatrixRef pOut = aFFT.Compute(aFunc);
5551  PushMatrix(pOut);
5552 }
5553 
5554 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
void ScFDist_LT()
Definition: interpr3.cxx:1704
StackVar
void ScBadName()
Definition: interpr3.cxx:185
void ScPearson()
Definition: interpr3.cxx:4517
static void lcl_QuickSort(tools::Long nLo, tools::Long nHi, vector< double > &rSortArray, vector< tools::Long > *pIndexOrder)
Definition: interpr3.cxx:4128
sal_Int32 nIndex
ScInterpreter & rInt
Definition: interpr3.cxx:2160
bool CalculateSkew(KahanSum &fSum, double &fCount, std::vector< double > &values)
Definition: interpr3.cxx:3211
ScAddress aStart
Definition: address.hxx:497
std::vector< kOp > kOpSumAndSumSquare
void ScBetaInv()
Definition: interpr3.cxx:2212
bool hasError() const
Definition: cellvalue.cxx:627
void ScRank(bool bAverage)
Definition: interpr3.cxx:4206
static double GetPercentile(::std::vector< double > &rArray, double fPercentile)
Definition: interpr3.cxx:3400
SCROW Row() const
Definition: address.hxx:274
void ScGauss()
Definition: interpr3.cxx:1145
void ScCritBinom()
Definition: interpr3.cxx:1410
std::vector< char * > values
static double phi(double x)
Definition: interpr3.cxx:195
static double GetChiSqDistPDF(double fX, double fDF)
Definition: interpr3.cxx:689
std::string GetValue
double divide(const double &fNumerator, const double &fDenominator)
Return fNumerator/fDenominator if fDenominator!=0 else +-Infinity if fNumerator!=0 or NaN if fNumerat...
Definition: math.hxx:49
std::string GetChiSqDistPDF
int SetError()
void ScStdNormDist_MS()
Definition: interpr3.cxx:1601
ScChiDistFunction(ScInterpreter &rI, double fpVal, double fDFVal)
Definition: interpr3.cxx:2360
double GetFDist(double x, double fF1, double fF2)
Definition: interpr3.cxx:640
void ScChiTest()
Definition: interpr3.cxx:2780
void ScTDist_MS()
Definition: interpr3.cxx:1674
ScFDistFunction(ScInterpreter &rI, double fpVal, double fF1Val, double fF2Val)
Definition: interpr3.cxx:2304
void ScStandard()
Definition: interpr3.cxx:3196
void ScHypGeomDist(int nMinParamCount)
Calculates a value of the hypergeometric distribution.
Definition: interpr3.cxx:1856
void GetSortArray(sal_uInt8 nParamCount,::std::vector< double > &rSortArray,::std::vector< tools::Long > *pIndexOrder, bool bConvertTextInArray, bool bAllowEmptyArray)
Definition: interpr3.cxx:4112
void ScLogNormInv()
Definition: interpr3.cxx:2143
void ScFourier()
Definition: interpr3.cxx:5474
sal_uIntPtr sal_uLong
double GetPercentileExclusive(::std::vector< double > &rArray, double fPercentile)
Definition: interpr3.cxx:3424
long Long
static void lcl_normalize(std::vector< double > &rCmplxArray, bool bScaleOnlyReal)
Definition: interpr3.cxx:5078
sal_Int64 n
static void lcl_roundUpNearestPow2(SCSIZE &nNum, SCSIZE &nNumBits)
Definition: interpr3.cxx:4782
#define max(a, b)
double GetGamma(double x)
You must ensure non integer arguments for fZ<1.
Definition: interpr3.cxx:585
#define N
void GetNumberSequenceArray(sal_uInt8 nParamCount,::std::vector< double > &rArray, bool bConvertTextInArray)
Definition: interpr3.cxx:3929
ScAddress aEnd
Definition: address.hxx:498
void ScDevSq()
Definition: interpr3.cxx:4427
double GetChiDist(double fChi, double fDF)
You must ensure fDF>0.0.
Definition: interpr3.cxx:669
std::string GetGammaDistPDF
void ScChiSqDist_MS()
Definition: interpr3.cxx:750
void ScIntercept()
Definition: interpr3.cxx:4711
This is very similar to ScCellValue, except that it references the original value instead of copying ...
Definition: cellvalue.hxx:103
void ScChiSqDist()
Definition: interpr3.cxx:727
void ScStdNormDist()
Definition: interpr3.cxx:1596
static double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)
Definition: interpr3.cxx:937
SCROW GetMaxRowCount() const
Definition: sheetlimits.hxx:66
void ScModalValue_MS(bool bSingle)
Definition: interpr3.cxx:3544
static size_t MAX_COUNT_DOUBLE_FOR_SORT(const ScSheetLimits &rSheetLimits)
Two columns of data should be sortable with GetSortArray() and QuickSort()
Definition: interpr3.cxx:45
float x
std::string gaussinv
void ScPermut()
Definition: interpr3.cxx:1199
void ScGammaInv()
Definition: interpr3.cxx:2172
void ScPoissonDist(bool bODFF)
Definition: interpr3.cxx:1781
std::vector< double > GetTopNumberArray(SCSIZE &rCol, SCSIZE &rRow)
Definition: interpr3.cxx:3844
std::string GetBetaDistPDF
void ScGammaDist(bool bODFF)
Definition: interpr3.cxx:2091
void ScNoName()
Definition: interpr3.cxx:180
tuple log
double GetValue(double x) const override
Definition: interpr3.cxx:2209
std::string GetBetaDist
ScInterpreter & rInt
Definition: interpr3.cxx:2356
double GetBeta(double fAlpha, double fBeta)
Definition: interpr3.cxx:800
char sal_uInt16 & nParamCount
Definition: callform.cxx:53
double Fakultaet(double x)
Definition: interpr3.cxx:439
This class provides LO with Kahan summation algorithm About this algorithm: https://en.wikipedia.org/wiki/Kahan_summation_algorithm For general purpose software we assume first order error is enough.
Definition: kahan.hxx:25
void ScCorrel()
Definition: interpr3.cxx:4501
size_t SCSIZE
size_t typedef to be able to find places where code was changed from USHORT to size_t and is used to ...
Definition: address.hxx:44
void ScNegBinomDist_MS()
Definition: interpr3.cxx:1520
exports com.sun.star.deployment. ui
void ScExpDist()
Definition: interpr3.cxx:1615
static double lcl_IterateInverse(const ScDistFunc &rFunction, double fAx, double fBx, bool &rConvError)
Definition: interpr3.cxx:76
void ScFisher()
Definition: interpr3.cxx:1150
void CalculateSlopeIntercept(bool bSlope)
Definition: interpr3.cxx:4636
void ScNormInv()
Definition: interpr3.cxx:2116
int nCount
static double gaussinv(double x)
Definition: interpr3.cxx:268
std::string GetBinomDistPMF
oslFileHandle & pOut
static const double fMaxGammaArgument
Definition: interpre.hxx:983
::boost::intrusive_ptr< ScMatrix > ScMatrixRef
Definition: types.hxx:25
SCTAB Tab() const
Definition: address.hxx:283
virtual ~ScGammaDistFunction()
Definition: interpr3.cxx:2167
ScInterpreter & rInt
Definition: interpr3.cxx:2200
void ScChiSqInv()
Definition: interpr3.cxx:2403
ScInterpreter & rInt
Definition: interpr3.cxx:2300
double GetTDist(double T, double fDF, int nType)
Definition: interpr3.cxx:648
static double lcl_GetGammaHelper(double fZ)
You must ensure fZ>0; fZ>171.624376956302 will overflow.
Definition: interpr3.cxx:559
std::string gauss
void ScFisherInv()
Definition: interpr3.cxx:1159
static double gauss(double x)
Definition: interpr3.cxx:215
constexpr std::enable_if_t< std::is_floating_point_v< F > &&std::is_integral_v< I >, bool > convertsToAtMost(F value, I max)
std::string taylor
ScInterpreter & rInt
Definition: interpr3.cxx:2251
void ScSkewp()
Definition: interpr3.cxx:3361
void swap(cow_wrapper< T, P > &a, cow_wrapper< T, P > &b)
ScMatrixRef(SCSIZE, SCSIZE, std::vector< double > &) ScMatrixGenerator
Definition: interpr3.cxx:5404
void ScGeoMean()
Definition: interpr3.cxx:3024
void CalculateSmallLarge(bool bSmall)
Definition: interpr3.cxx:3630
static SCSIZE lcl_bitReverse(SCSIZE nIn, SCSIZE nBound)
Definition: interpr3.cxx:4803
ScGammaDistFunction(ScInterpreter &rI, double fpVal, double fAlphaVal, double fBetaVal)
Definition: interpr3.cxx:2164
void ScTInv(int nType)
Definition: interpr3.cxx:2264
void ScTrimMean()
Definition: interpr3.cxx:3815
int i
static double lcl_GetLogGammaHelper(double fZ)
You must ensure fZ>0.
Definition: interpr3.cxx:577
static bool lcl_HasChangeOfSign(double u, double w)
u*w<0.0 fails for values near zero
Definition: interpr3.cxx:71
std::string GetUpRegIGamma
void ScLogGamma()
Definition: interpr3.cxx:791
void CalculatePearsonCovar(bool _bPearson, bool _bStexy, bool _bSample)
Definition: interpr3.cxx:4522
std::string GetLowRegIGamma
double GetValue(double x) const override
Definition: interpr3.cxx:2400
void ScBinomDist()
Definition: interpr3.cxx:1346
virtual ~ScTDistFunction()
Definition: interpr3.cxx:2259
virtual ~ScChiSqDistFunction()
Definition: interpr3.cxx:2398
void ScPercentile(bool bInclusive)
Definition: interpr3.cxx:3454
double GetGammaDist(double fX, double fAlpha, double fLambda)
Gamma distribution, cumulative distribution function.
Definition: interpr6.cxx:198
static double GetLogGamma(double x)
You must ensure fZ>0.
Definition: interpr3.cxx:629
std::string GetLogGamma
virtual ~ScChiDistFunction()
Definition: interpr3.cxx:2363
sal_uInt16 & nParam
Definition: callform.cxx:57
size_t uniform_size_distribution(size_t a, size_t b)
void ScMedian()
Definition: interpr3.cxx:3390
std::string GetChiDist
void ScSlope()
Definition: interpr3.cxx:4706
std::string GetBeta
static void QuickSort(::std::vector< double > &rSortArray,::std::vector< tools::Long > *pIndexOrder)
Definition: interpr3.cxx:4179
void ScTDist_T(int nTails)
Definition: interpr3.cxx:1656
std::string GetChiSqDistCDF
svExternalDoubleRef
void ScConfidence()
Definition: interpr3.cxx:2423
size
svExternalSingleRef
void ScFTest()
Definition: interpr3.cxx:2730
static double taylor(const double *pPolynom, sal_uInt16 nMax, double x)
Definition: interpr3.cxx:205
bool hasNumeric() const
Definition: cellvalue.cxx:622
void ScChiDist(bool bODFF)
Definition: interpr3.cxx:1741
double GetValue(double x) const override
Definition: interpr3.cxx:2261
ScBetaDistFunction(ScInterpreter &rI, double fpVal, double fAlphaVal, double fBetaVal)
Definition: interpr3.cxx:2204
std::string GetLogBeta
static void lcl_PutFactorialElements(::std::vector< double > &cn, double fLower, double fUpper, double fBase)
Local function used in the calculation of the hypergeometric distribution.
Definition: interpr3.cxx:1834
void ScConfidenceT()
Definition: interpr3.cxx:2437
FormulaError
SCCOL Col() const
Definition: address.hxx:279
XPropertyListType t
double GetBinomDistPMF(double x, double n, double p)
Definition: interpr3.cxx:1232
void CalculateSkewOrSkewp(bool bSkewp)
Definition: interpr3.cxx:3316
void ScSmall()
Definition: interpr3.cxx:3715
float v
std::string GetGammaDist
void PutInOrder()
Definition: address.hxx:622
std::string GetTDist
double CreateDoubleError(FormulaError nErr)
void ScModalValue()
Definition: interpr3.cxx:3500
static double lcl_GetBinomDistRange(double n, double xs, double xe, double fFactor, double p, double q)
Definition: interpr3.cxx:1260
static double GetLogBeta(double fAlpha, double fBeta)
Definition: interpr3.cxx:833
double GetValue(double x) const override
Definition: interpr3.cxx:2309
virtual ~ScBetaDistFunction()
Definition: interpr3.cxx:2207
void ScLarge()
Definition: interpr3.cxx:3710
static double BinomKoeff(double n, double k)
Definition: interpr3.cxx:460
double GetChiSqDistCDF(double fX, double fDF)
You must ensure fDF>0.0.
Definition: interpr3.cxx:681
ScInterpreter & rInt
Definition: interpr3.cxx:2391
void ScBetaDist()
Definition: interpr3.cxx:1036
double GetValue(double x) const override
Definition: interpr3.cxx:2169
unsigned char sal_uInt8
float z
bool GetFirst(double &rValue, FormulaError &rErr)
Does NOT reset rValue if no value found!
Definition: dociter.cxx:270
static double lcl_getLanczosSum(double fZ)
you must ensure fZ>0 Uses a variant of the Lanczos sum with a rational function.
Definition: interpr3.cxx:491
std::string GetFDist
void ScNegBinomDist()
Definition: interpr3.cxx:1500
void ScLogNormDist(int nMinParamCount)
Definition: interpr3.cxx:1566
double GetTInv(double fAlpha, double fSize, int nType)
Definition: interpr3.cxx:2288
void ScBetaDist_MS()
Microsoft version has parameters in different order Also, upper and lowerbound are optional and have ...
Definition: interpr3.cxx:1100
void ScNormDist(int nMinParamCount)
Definition: interpr3.cxx:1546
void * p
QPRO_FUNC_TYPE nType
Definition: qproform.cxx:398
bool GetNext(double &rValue, FormulaError &rErr)
Does NOT reset rValue if no value found!
Definition: dociter.cxx:294
void ScPermutationA()
Definition: interpr3.cxx:1219
void ScSTEYX()
Definition: interpr3.cxx:4632
virtual ~ScFDistFunction()
Definition: interpr3.cxx:2307
void ScChiInv()
Definition: interpr3.cxx:2368
bool CalculateTest(bool _bTemplin, const SCSIZE nC1, const SCSIZE nC2, const SCSIZE nR1, const SCSIZE nR2, const ScMatrixRef &pMat1, const ScMatrixRef &pMat2, double &fT, double &fF)
Definition: interpr3.cxx:2580
double GetHypGeomDist(double x, double n, double M, double N)
Calculates a value of the hypergeometric distribution.
Definition: interpr3.cxx:1895
void ScCovarianceS()
Definition: interpr3.cxx:4512
void ScWeibull()
Definition: interpr3.cxx:1763
void ScTDist()
Definition: interpr3.cxx:1641
void ScFInv_LT()
Definition: interpr3.cxx:2333
double GetValue(double x) const override
Definition: interpr3.cxx:2365
void ScProbability()
Definition: interpr3.cxx:4436
void(* f)(TrueTypeTable *)
void ScPercentrank(bool bInclusive)
Definition: interpr3.cxx:3720
double get() const
Returns the final sum.
Definition: kahan.hxx:211
void ScFDist()
Definition: interpr3.cxx:1689
void ScForecast()
Definition: interpr3.cxx:4716
void ScCombin()
Definition: interpr3.cxx:1173
ScTDistFunction(ScInterpreter &rI, double fpVal, double fDFVal, int nType)
Definition: interpr3.cxx:2256
ScChiSqDistFunction(ScInterpreter &rI, double fpVal, double fDFVal)
Definition: interpr3.cxx:2395
std::string phi
double GetBetaDistPDF(double fX, double fA, double fB)
Definition: interpr3.cxx:862
static void lcl_convertToPolar(std::vector< double > &rCmplxArray, double fMinMag)
Definition: interpr3.cxx:5095
void ScQuartile(bool bInclusive)
Definition: interpr3.cxx:3477
M
void ScZTest()
Definition: interpr3.cxx:2453
void ScCombinA()
Definition: interpr3.cxx:1186
aStr
const double fMachEps
Definition: interpr3.cxx:51
double GetMedian(::std::vector< double > &rArray)
Definition: interpr3.cxx:3366
double div(const double &fNumerator, const double &fDenominator)
Return fNumerator/fDenominator if fDenominator!=0 else #DIV/0! error coded into double.
Definition: math.hxx:29
static double GetPercentrank(::std::vector< double > &rArray, double fVal, bool bInclusive)
Definition: interpr3.cxx:3758
void ScTTest()
Definition: interpr3.cxx:2647
void ScHarMean()
Definition: interpr3.cxx:2903
void ScSNormInv()
Definition: interpr3.cxx:2132
double GetBetaDist(double x, double alpha, double beta)
Definition: interpr3.cxx:985
const sal_uInt8 R
Definition: xlformula.cxx:50
void ScAveDev()
Definition: interpr3.cxx:4269
static double integralPhi(double x)
Definition: interpr3.cxx:200
void ScGamma()
Definition: interpr3.cxx:774
void ScCovarianceP()
Definition: interpr3.cxx:4507