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