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