?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
update ooo310-m15
[ooovba.git] / sc / source / core / tool / interpr3.cxx
blob372ef420d66228ac5073598ef65b5f7610cf6d08
1 /*************************************************************************
2 *
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * Copyright 2008 by Sun Microsystems, Inc.
6 *
7 * OpenOffice.org - a multi-platform office productivity suite
8 *
9 * $RCSfile: interpr3.cxx,v $
10 * $Revision: 1.25 $
11 *
12 * This file is part of OpenOffice.org.
13 *
14 * OpenOffice.org is free software: you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License version 3
16 * only, as published by the Free Software Foundation.
17 *
18 * OpenOffice.org is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU Lesser General Public License version 3 for more details
22 * (a copy is included in the LICENSE file that accompanied this code).
23 *
24 * You should have received a copy of the GNU Lesser General Public License
25 * version 3 along with OpenOffice.org. If not, see
26 * <http://www.openoffice.org/license.html>
27 * for a copy of the LGPLv3 License.
28 *
29 ************************************************************************/
30
31 // MARKER(update_precomp.py): autogen include statement, do not remove
32 #include "precompiled_sc.hxx"
33
34 // INCLUDE ---------------------------------------------------------------
35
36 #include <tools/solar.h>
37 #include <stdlib.h>
38 #include <string.h>
39 #include <rtl/logfile.hxx>
40
41 #include "interpre.hxx"
42 #include "global.hxx"
43 #include "compiler.hxx"
44 #include "cell.hxx"
45 #include "document.hxx"
46 #include "dociter.hxx"
47 #include "scmatrix.hxx"
48 #include "globstr.hrc"
49
50 #include <math.h>
51 #include <vector>
52 #include <algorithm>
53
54 using ::std::vector;
55 using namespace formula;
56
57 // STATIC DATA -----------------------------------------------------------
58
59 #define SCdEpsilon 1.0E-7
60 #define SC_MAX_ITERATION_COUNT 20
61 #define MAX_ANZ_DOUBLE_FOR_SORT 100000
62 // PI jetzt als F_PI aus solar.h
63 //#define PI 3.1415926535897932
64
65 const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental
66 const double fMachEps = ::std::numeric_limits<double>::epsilon();
67
68 //-----------------------------------------------------------------------------
69
70 class ScDistFunc
71 {
72 public:
73 virtual double GetValue(double x) const = 0;
74 };
75
76 // iteration for inverse distributions
77
78 //template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, bool& rConvError )
79
80 /** u*w<0.0 fails for values near zero */
81 inline bool lcl_HasChangeOfSign( double u, double w )
82 {
83 return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);
84 }
85
86 double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError )
87 {
88 rConvError = false;
89 const double fYEps = 1.0E-307;
90 const double fXEps = ::std::numeric_limits<double>::epsilon();
91
92 DBG_ASSERT(fAx<fBx, "IterateInverse: wrong interval");
93
94 // find enclosing interval
95
96 double fAy = rFunction.GetValue(fAx);
97 double fBy = rFunction.GetValue(fBx);
98 double fTemp;
99 unsigned short nCount;
100 for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++)
101 {
102 if (fabs(fAy) <= fabs(fBy))
103 {
104 fTemp = fAx;
105 fAx += 2.0 * (fAx - fBx);
106 if (fAx < 0.0)
107 fAx = 0.0;
108 fBx = fTemp;
109 fBy = fAy;
110 fAy = rFunction.GetValue(fAx);
111 }
112 else
113 {
114 fTemp = fBx;
115 fBx += 2.0 * (fBx - fAx);
116 fAx = fTemp;
117 fAy = fBy;
118 fBy = rFunction.GetValue(fBx);
119 }
120 }
121
122 if (fAy == 0.0)
123 return fAx;
124 if (fBy == 0.0)
125 return fBx;
126 if (!lcl_HasChangeOfSign( fAy, fBy))
127 {
128 rConvError = true;
129 return 0.0;
130 }
131 // inverse quadric interpolation with additional brackets
132 // set three points
133 double fPx = fAx;
134 double fPy = fAy;
135 double fQx = fBx;
136 double fQy = fBy;
137 double fRx = fAx;
138 double fRy = fAy;
139 double fSx = 0.5 * (fAx + fBx); // potential next point
140 bool bHasToInterpolate = true;
141 nCount = 0;
142 while ( nCount < 500 && fabs(fRy) > fYEps &&
143 (fBx-fAx) > ::std::max( fabs(fAx), fabs(fBx)) * fXEps )
144 {
145 if (bHasToInterpolate)
146 {
147 if (fPy!=fQy && fQy!=fRy && fRy!=fPy)
148 {
149 fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)
150 + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)
151 + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);
152 bHasToInterpolate = (fAx < fSx) && (fSx < fBx); // inside the brackets?
153 }
154 else
155 bHasToInterpolate = false;
156 }
157 if(!bHasToInterpolate)
158 {
159 fSx = 0.5 * (fAx + fBx);
160 // reset points
161 fPx = fAx; fPy = fAy;
162 fQx = fBx; fQy = fBy;
163 bHasToInterpolate = true;
164 }
165 // shift points for next interpolation
166 fPx = fQx; fQx = fRx; fRx = fSx;
167 fPy = fQy; fQy = fRy; fRy = rFunction.GetValue(fSx);
168 // update brackets
169 if (lcl_HasChangeOfSign( fAy, fRy))
170 {
171 fBx = fRx; fBy = fRy;
172 }
173 else
174 {
175 fAx = fRx; fAy = fRy;
176 }
177 // if last interration brought to small advance, then do bisection next
178 // time, for safety
179 bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));
180 ++nCount;
181 }
182 return fRx;
183 }
184
185 //-----------------------------------------------------------------------------
186 // Allgemeine Funktionen
187 //-----------------------------------------------------------------------------
188
189 void ScInterpreter::ScNoName()
190 {
191 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScNoName" );
192 PushError(errNoName);
193 }
194
195 void ScInterpreter::ScBadName()
196 {
197 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScBadName" );
198 short nParamCount = GetByte();
199 while (nParamCount-- > 0)
200 {
201 PopError();
202 }
203 PushError( errNoName);
204 }
205
206 double ScInterpreter::phi(double x)
207 {
208 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::phi" );
209 return 0.39894228040143268 * exp(-(x * x) / 2.0);
210 }
211
212 double ScInterpreter::taylor(double* pPolynom, USHORT nMax, double x)
213 {
214 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::taylor" );
215 double nVal = pPolynom[nMax];
216 for (short i = nMax-1; i >= 0; i--)
217 {
218 nVal = pPolynom[i] + (nVal * x);
219 }
220 return nVal;
221 }
222
223 double ScInterpreter::gauss(double x)
224 {
225 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::gauss" );
226 double t0[] =
227 { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,
228 -0.00118732821548045, 0.00011543468761616, -0.00000944465625950,
229 0.00000066596935163, -0.00000004122667415, 0.00000000227352982,
230 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };
231 double t2[] =
232 { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,
233 0.02699548325659403, -0.00449924720943234, -0.00224962360471617,
234 0.00134977416282970, -0.00011783742691370, -0.00011515930357476,
235 0.00003704737285544, 0.00000282690796889, -0.00000354513195524,
236 0.00000037669563126, 0.00000019202407921, -0.00000005226908590,
237 -0.00000000491799345, 0.00000000366377919, -0.00000000015981997,
238 -0.00000000017381238, 0.00000000002624031, 0.00000000000560919,
239 -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };
240 double t4[] =
241 { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,
242 0.00033457556441221, -0.00028996548915725, 0.00018178605666397,
243 -0.00008252863922168, 0.00002551802519049, -0.00000391665839292,
244 -0.00000074018205222, 0.00000064422023359, -0.00000017370155340,
245 0.00000000909595465, 0.00000000944943118, -0.00000000329957075,
246 0.00000000029492075, 0.00000000011874477, -0.00000000004420396,
247 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };
248 double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };
249
250 double xAbs = fabs(x);
251 USHORT xShort = (USHORT)::rtl::math::approxFloor(xAbs);
252 double nVal = 0.0;
253 if (xShort == 0)
254 nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;
255 else if ((xShort >= 1) && (xShort <= 2))
256 nVal = taylor(t2, 23, (xAbs - 2.0));
257 else if ((xShort >= 3) && (xShort <= 4))
258 nVal = taylor(t4, 20, (xAbs - 4.0));
259 else
260 nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs;
261 if (x < 0.0)
262 return -nVal;
263 else
264 return nVal;
265 }
266
267 //
268 // #i26836# new gaussinv implementation by Martin Eitzenberger <m.eitzenberger@unix.net>
269 //
270
271 double ScInterpreter::gaussinv(double x)
272 {
273 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::gaussinv" );
274 double q,t,z;
275
276 q=x-0.5;
277
278 if(fabs(q)<=.425)
279 {
280 t=0.180625-q*q;
281
282 z=
283 q*
284 (
285 (
286 (
287 (
288 (
289 (
290 (
291 t*2509.0809287301226727+33430.575583588128105
292 )
293 *t+67265.770927008700853
294 )
295 *t+45921.953931549871457
296 )
297 *t+13731.693765509461125
298 )
299 *t+1971.5909503065514427
300 )
301 *t+133.14166789178437745
302 )
303 *t+3.387132872796366608
304 )
305 /
306 (
307 (
308 (
309 (
310 (
311 (
312 (
313 t*5226.495278852854561+28729.085735721942674
314 )
315 *t+39307.89580009271061
316 )
317 *t+21213.794301586595867
318 )
319 *t+5394.1960214247511077
320 )
321 *t+687.1870074920579083
322 )
323 *t+42.313330701600911252
324 )
325 *t+1.0
326 );
327
328 }
329 else
330 {
331 if(q>0) t=1-x;
332 else t=x;
333
334 t=sqrt(-log(t));
335
336 if(t<=5.0)
337 {
338 t+=-1.6;
339
340 z=
341 (
342 (
343 (
344 (
345 (
346 (
347 (
348 t*7.7454501427834140764e-4+0.0227238449892691845833
349 )
350 *t+0.24178072517745061177
351 )
352 *t+1.27045825245236838258
353 )
354 *t+3.64784832476320460504
355 )
356 *t+5.7694972214606914055
357 )
358 *t+4.6303378461565452959
359 )
360 *t+1.42343711074968357734
361 )
362 /
363 (
364 (
365 (
366 (
367 (
368 (
369 (
370 t*1.05075007164441684324e-9+5.475938084995344946e-4
371 )
372 *t+0.0151986665636164571966
373 )
374 *t+0.14810397642748007459
375 )
376 *t+0.68976733498510000455
377 )
378 *t+1.6763848301838038494
379 )
380 *t+2.05319162663775882187
381 )
382 *t+1.0
383 );
384
385 }
386 else
387 {
388 t+=-5.0;
389
390 z=
391 (
392 (
393 (
394 (
395 (
396 (
397 (
398 t*2.01033439929228813265e-7+2.71155556874348757815e-5
399 )
400 *t+0.0012426609473880784386
401 )
402 *t+0.026532189526576123093
403 )
404 *t+0.29656057182850489123
405 )
406 *t+1.7848265399172913358
407 )
408 *t+5.4637849111641143699
409 )
410 *t+6.6579046435011037772
411 )
412 /
413 (
414 (
415 (
416 (
417 (
418 (
419 (
420 t*2.04426310338993978564e-15+1.4215117583164458887e-7
421 )
422 *t+1.8463183175100546818e-5
423 )
424 *t+7.868691311456132591e-4
425 )
426 *t+0.0148753612908506148525
427 )
428 *t+0.13692988092273580531
429 )
430 *t+0.59983220655588793769
431 )
432 *t+1.0
433 );
434
435 }
436
437 if(q<0.0) z=-z;
438 }
439
440 return z;
441 }
442
443 double ScInterpreter::Fakultaet(double x)
444 {
445 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::Fakultaet" );
446 x = ::rtl::math::approxFloor(x);
447 if (x < 0.0)
448 return 0.0;
449 else if (x == 0.0)
450 return 1.0;
451 else if (x <= 170.0)
452 {
453 double fTemp = x;
454 while (fTemp > 2.0)
455 {
456 fTemp--;
457 x *= fTemp;
458 }
459 }
460 else
461 SetError(errNoValue);
462 /* // Stirlingsche Naeherung zu ungenau
463 else
464 x = pow(x/exp(1), x) * sqrt(x) * SQRT_2_PI * (1.0 + 1.0 / (12.0 * x));
465 */
466 return x;
467 }
468
469 double ScInterpreter::BinomKoeff(double n, double k)
470 {
471 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::BinomKoeff" );
472 double nVal = 0.0;
473 k = ::rtl::math::approxFloor(k);
474 if (n < k)
475 nVal = 0.0;
476 else if (k == 0.0)
477 nVal = 1.0;
478 else
479 {
480 nVal = n/k;
481 n--;
482 k--;
483 while (k > 0.0)
484 {
485 nVal *= n/k;
486 k--;
487 n--;
488 }
489 /*
490 double f1 = n; // Zaehler
491 double f2 = k; // Nenner
492 n--;
493 k--;
494 while (k > 0.0)
495 {
496 f2 *= k;
497 f1 *= n;
498 k--;
499 n--;
500 }
501 nVal = f1 / f2;
502 */
503 }
504 return nVal;
505 }
506
507
508 // The algorithm is based on lanczos13m53 in lanczos.hpp
509 // in math library from http://www.boost.org
510 /** you must ensure fZ>0
511 Uses a variant of the Lanczos sum with a rational function. */
512 double lcl_getLanczosSum(double fZ)
513 {
514 const double fNum[13] ={
515 23531376880.41075968857200767445163675473,
516 42919803642.64909876895789904700198885093,
517 35711959237.35566804944018545154716670596,
518 17921034426.03720969991975575445893111267,
519 6039542586.35202800506429164430729792107,
520 1439720407.311721673663223072794912393972,
521 248874557.8620541565114603864132294232163,
522 31426415.58540019438061423162831820536287,
523 2876370.628935372441225409051620849613599,
524 186056.2653952234950402949897160456992822,
525 8071.672002365816210638002902272250613822,
526 210.8242777515793458725097339207133627117,
527 2.506628274631000270164908177133837338626
528 };
529 const double fDenom[13] = {
530 0,
531 39916800,
532 120543840,
533 150917976,
534 105258076,
535 45995730,
536 13339535,
537 2637558,
538 357423,
539 32670,
540 1925,
541 66,
542 1
543 };
544 // Horner scheme
545 double fSumNum;
546 double fSumDenom;
547 int nI;
548 double fZInv;
549 if (fZ<=1.0)
550 {
551 fSumNum = fNum[12];
552 fSumDenom = fDenom[12];
553 for (nI = 11; nI >= 0; --nI)
554 {
555 fSumNum *= fZ;
556 fSumNum += fNum[nI];
557 fSumDenom *= fZ;
558 fSumDenom += fDenom[nI];
559 }
560 }
561 else
562 // Cancel down with fZ^12; Horner scheme with reverse coefficients
563 {
564 fZInv = 1/fZ;
565 fSumNum = fNum[0];
566 fSumDenom = fDenom[0];
567 for (nI = 1; nI <=12; ++nI)
568 {
569 fSumNum *= fZInv;
570 fSumNum += fNum[nI];
571 fSumDenom *= fZInv;
572 fSumDenom += fDenom[nI];
573 }
574 }
575 return fSumNum/fSumDenom;
576 }
577
578 // The algorithm is based on tgamma in gamma.hpp
579 // in math library from http://www.boost.org
580 /** You must ensure fZ>0; fZ>171.624376956302 will overflow. */
581 double lcl_GetGammaHelper(double fZ)
582 {
583 double fGamma = lcl_getLanczosSum(fZ);
584 const double fg = 6.024680040776729583740234375;
585 double fZgHelp = fZ + fg - 0.5;
586 // avoid intermediate overflow
587 double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25);
588 fGamma *= fHalfpower;
589 fGamma /= exp(fZgHelp);
590 fGamma *= fHalfpower;
591 if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ))
592 fGamma = ::rtl::math::round(fGamma);
593 return fGamma;
594 }
595
596 // The algorithm is based on tgamma in gamma.hpp
597 // in math library from http://www.boost.org
598 /** You must ensure fZ>0 */
599 double lcl_GetLogGammaHelper(double fZ)
600 {
601 const double fg = 6.024680040776729583740234375;
602 double fZgHelp = fZ + fg - 0.5;
603 return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;
604 }
605
606 /** You must ensure non integer arguments for fZ<1 */
607 double ScInterpreter::GetGamma(double fZ)
608 {
609 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetGamma" );
610 const double fLogPi = log(F_PI);
611 const double fLogDblMax = log( ::std::numeric_limits<double>::max());
612
613 if (fZ > fMaxGammaArgument)
614 {
615 SetError(errIllegalFPOperation);
616 return HUGE_VAL;
617 }
618
619 if (fZ >= 1.0)
620 return lcl_GetGammaHelper(fZ);
621
622 if (fZ >= 0.5) // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
623 return lcl_GetGammaHelper(fZ+1) / fZ;
624
625 if (fZ >= -0.5) // shift to x>=1, might overflow
626 {
627 double fLogTest = lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log( fabs(fZ));
628 if (fLogTest >= fLogDblMax)
629 {
630 SetError( errIllegalFPOperation);
631 return HUGE_VAL;
632 }
633 return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
634 }
635 // fZ<-0.5
636 // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) )
637 double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( fabs( ::rtl::math::sin( F_PI*fZ)));
638 if (fLogDivisor - fLogPi >= fLogDblMax) // underflow
639 return 0.0;
640
641 if (fLogDivisor<0.0)
642 if (fLogPi - fLogDivisor > fLogDblMax) // overflow
643 {
644 SetError(errIllegalFPOperation);
645 return HUGE_VAL;
646 }
647
648 return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( F_PI*fZ) < 0.0) ? -1.0 : 1.0);
649 }
650
651
652 /** You must ensure fZ>0 */
653 double ScInterpreter::GetLogGamma(double fZ)
654 {
655 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetLogGamma" );
656 if (fZ >= fMaxGammaArgument)
657 return lcl_GetLogGammaHelper(fZ);
658 if (fZ >= 1.0)
659 return log(lcl_GetGammaHelper(fZ));
660 if (fZ >= 0.5)
661 return log( lcl_GetGammaHelper(fZ+1) / fZ);
662 return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);
663 }
664
665 double ScInterpreter::GetFDist(double x, double fF1, double fF2)
666 {
667 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetFDist" );
668 double arg = fF2/(fF2+fF1*x);
669 double alpha = fF2/2.0;
670 double beta = fF1/2.0;
671 return (GetBetaDist(arg, alpha, beta));
672 /*
673 double Z = (pow(fF,1.0/3.0)*(1.0-2.0/(9.0*fF2)) - (1.0-2.0/(9.0*fF1))) /
674 sqrt(2.0/(9.0*fF1) + pow(fF,2.0/3.0)*2.0/(9.0*fF2));
675 return (0.5-gauss(Z));
676 */
677 }
678
679 double ScInterpreter::GetTDist(double T, double fDF)
680 {
681 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetTDist" );
682 return 0.5 * GetBetaDist(fDF/(fDF+T*T), fDF/2.0, 0.5);
683 /*
684 USHORT DF = (USHORT) fDF;
685 double A = T / sqrt(DF);
686 double B = 1.0 + A*A;
687 double R;
688 if (DF == 1)
689 R = 0.5 + atan(A)/F_PI;
690 else if (DF % 2 == 0)
691 {
692 double S0 = A/(2.0 * sqrt(B));
693 double C0 = S0;
694 for (USHORT i = 2; i <= DF-2; i+=2)
695 {
696 C0 *= (1.0 - 1.0/(double)i)/B;
697 S0 += C0;
698 }
699 R = 0.5 + S0;
700 }
701 else
702 {
703 double S1 = A / (B * F_PI);
704 double C1 = S1;
705 for (USHORT i = 3; i <= DF-2; i+=2)
706 {
707 C1 *= (1.0 - 1.0/(double)i)/B;
708 S1 += C1;
709 }
710 R = 0.5 + atan(A)/F_PI + S1;
711 }
712 return 1.0 - R;
713 */
714 }
715
716 // for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom
717 /** You must ensure fDF>0.0 */
718 double ScInterpreter::GetChiDist(double fX, double fDF)
719 {
720 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetChiDist" );
721 if (fX <= 0.0)
722 return 1.0; // see ODFF
723 else
724 return GetUpRegIGamma( fDF/2.0, fX/2.0);
725 }
726
727 // ready for ODF 1.2
728 // for ODF CHISQDIST; cumulative distribution function, fDF=degrees of freedom
729 // returns left tail
730 /** You must ensure fDF>0.0 */
731 double ScInterpreter::GetChiSqDistCDF(double fX, double fDF)
732 {
733 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetChiSqDistCDF" );
734 if (fX <= 0.0)
735 return 0.0; // see ODFF
736 else
737 return GetLowRegIGamma( fDF/2.0, fX/2.0);
738 }
739
740 double ScInterpreter::GetChiSqDistPDF(double fX, double fDF)
741 {
742 // you must ensure fDF is positive integer
743 double fValue;
744 double fCount;
745 if (fX <= 0.0)
746 return 0.0; // see ODFF
747 if (fDF*fX > 1391000.0)
748 {
749 // intermediate invalid values, use log
750 fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0) - GetLogGamma(0.5*fDF));
751 }
752 else // fDF is small in most cases, we can iterate
753 {
754 if (fmod(fDF,2.0)<0.5)
755 {
756 // even
757 fValue = 0.5;
758 fCount = 2.0;
759 }
760 else
761 {
762 fValue = 1/sqrt(fX*2*F_PI);
763 fCount = 1.0;
764 }
765 while ( fCount < fDF)
766 {
767 fValue *= (fX / fCount);
768 fCount += 2.0;
769 }
770 if (fX>=1425.0) // underflow in e^(-x/2)
771 fValue = exp(log(fValue)-fX/2);
772 else
773 fValue *= exp(-fX/2);
774 }
775 return fValue;
776 }
777
778 void ScInterpreter::ScChiSqDist()
779 {
780 BYTE nParamCount = GetByte();
781 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
782 return;
783 double fX;
784 bool bCumulative;
785 if (nParamCount == 3)
786 bCumulative = GetBool();
787 else
788 bCumulative = true;
789 double fDF = ::rtl::math::approxFloor(GetDouble());
790 if (fDF < 1.0)
791 PushIllegalArgument();
792 else
793 {
794 fX = GetDouble();
795 if (bCumulative)
796 PushDouble(GetChiSqDistCDF(fX,fDF));
797 else
798 PushDouble(GetChiSqDistPDF(fX,fDF));
799 }
800 }
801
802 void ScInterpreter::ScGamma()
803 {
804 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGamma" );
805 double x = GetDouble();
806 double fResult;
807 if (x <= 0.0 && x == ::rtl::math::approxFloor(x))
808 PushIllegalArgument();
809 else
810 {
811 fResult = GetGamma(x);
812 if (nGlobalError)
813 {
814 PushError( nGlobalError);
815 return;
816 }
817 PushDouble(fResult);
818 }
819 }
820
821
822 void ScInterpreter::ScLogGamma()
823 {
824 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScLogGamma" );
825 double x = GetDouble();
826 if (x > 0.0) // constraint from ODFF
827 PushDouble( GetLogGamma(x));
828 else
829 PushIllegalArgument();
830 }
831
832 double ScInterpreter::GetBeta(double fAlpha, double fBeta)
833 {
834 double fA;
835 double fB;
836 if (fAlpha > fBeta)
837 {
838 fA = fAlpha; fB = fBeta;
839 }
840 else
841 {
842 fA = fBeta; fB = fAlpha;
843 }
844 if (fA+fB < fMaxGammaArgument) // simple case
845 return GetGamma(fA)/GetGamma(fA+fB)*GetGamma(fB);
846 // need logarithm
847 // GetLogGamma is not accurate enough, back to Lanczos for all three
848 // GetGamma and arrange factors newly.
849 const double fg = 6.024680040776729583740234375; //see GetGamma
850 double fgm = fg - 0.5;
851 double fLanczos = lcl_getLanczosSum(fA);
852 fLanczos /= lcl_getLanczosSum(fA+fB);
853 fLanczos *= lcl_getLanczosSum(fB);
854 double fABgm = fA+fB+fgm;
855 fLanczos *= sqrt((fABgm/(fA+fgm))/(fB+fgm));
856 double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
857 double fTempB = fA/(fB+fgm);
858 double fResult = exp(-fA * ::rtl::math::log1p(fTempA)
859 -fB * ::rtl::math::log1p(fTempB)-fgm);
860 fResult *= fLanczos;
861 return fResult;
862 }
863
864 // Same as GetBeta but with logarithm
865 double ScInterpreter::GetLogBeta(double fAlpha, double fBeta)
866 {
867 double fA;
868 double fB;
869 if (fAlpha > fBeta)
870 {
871 fA = fAlpha; fB = fBeta;
872 }
873 else
874 {
875 fA = fBeta; fB = fAlpha;
876 }
877 const double fg = 6.024680040776729583740234375; //see GetGamma
878 double fgm = fg - 0.5;
879 double fLanczos = lcl_getLanczosSum(fA);
880 fLanczos /= lcl_getLanczosSum(fA+fB);
881 fLanczos *= lcl_getLanczosSum(fB);
882 double fLogLanczos = log(fLanczos);
883 double fABgm = fA+fB+fgm;
884 fLogLanczos += 0.5*(log(fABgm)-log(fA+fgm)-log(fB+fgm));
885 double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
886 double fTempB = fA/(fB+fgm);
887 double fResult = -fA * ::rtl::math::log1p(fTempA)
888 -fB * ::rtl::math::log1p(fTempB)-fgm;
889 fResult += fLogLanczos;
890 return fResult;
891 }
892
893 // beta distribution probability density function
894 double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB)
895 {
896 // special cases
897 if (fA == 1.0) // result b*(1-x)^(b-1)
898 {
899 if (fB == 1.0)
900 return 1.0;
901 if (fB == 2.0)
902 return -2.0*fX + 2.0;
903 if (fX == 1.0 && fB < 1.0)
904 {
905 SetError(errIllegalArgument);
906 return HUGE_VAL;
907 }
908 if (fX <= 0.01)
909 return fB + fB * ::rtl::math::expm1((fB-1.0) * ::rtl::math::log1p(-fX));
910 else
911 return fB * pow(0.5-fX+0.5,fB-1.0);
912 }
913 if (fB == 1.0) // result a*x^(a-1)
914 {
915 if (fA == 2.0)
916 return fA * fX;
917 if (fX == 0.0 && fA < 1.0)
918 {
919 SetError(errIllegalArgument);
920 return HUGE_VAL;
921 }
922 return fA * pow(fX,fA-1);
923 }
924 if (fX <= 0.0)
925 {
926 if (fA < 1.0 && fX == 0.0)
927 {
928 SetError(errIllegalArgument);
929 return HUGE_VAL;
930 }
931 else
932 return 0.0;
933 }
934 if (fX >= 1.0)
935 if (fB < 1.0 && fX == 1.0)
936 {
937 SetError(errIllegalArgument);
938 return HUGE_VAL;
939 }
940 else
941 return 0.0;
942
943 // normal cases; result x^(a-1)*(1-x)^(b-1)/Beta(a,b)
944 const double fLogDblMax = log( ::std::numeric_limits<double>::max());
945 const double fLogDblMin = log( ::std::numeric_limits<double>::min());
946 double fLogY = (fX < 0.1) ? ::rtl::math::log1p(-fX) : log(0.5-fX+0.5);
947 double fLogX = log(fX);
948 double fAm1 = fA-1.0;
949 double fBm1 = fB-1.0;
950 double fLogBeta = GetLogBeta(fA,fB);
951 // check whether parts over- or underflow
952 if ( fAm1 * fLogX < fLogDblMax && fAm1 * fLogX > fLogDblMin
953 && fBm1 * fLogY < fLogDblMax && fBm1* fLogY > fLogDblMin
954 && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin )
955 return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);
956 else // need logarithm;
957 // might overflow as a whole, but seldom, not worth to pre-detect it
958 return exp((fA-1.0)*fLogX + (fB-1.0)* fLogY - fLogBeta);
959 }
960
961
962 /*
963 x^a * (1-x)^b
964 I_x(a,b) = ---------------- * result of ContFrac
965 a * Beta(a,b)
966 */
967 double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)
968 { // like old version
969 double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, cfnew, cf;
970 a1 = 1.0; b1 = 1.0;
971 b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;
972 if (b2 == 0.0)
973 {
974 a2 = 0.0;
975 fnorm = 1.0;
976 cf = 1.0;
977 }
978 else
979 {
980 a2 = 1.0;
981 fnorm = 1.0/b2;
982 cf = a2*fnorm;
983 }
984 cfnew = 1.0;
985 double rm = 1.0;
986
987 const double fMaxIter = 50000.0;
988 // loop security, normal cases converge in less than 100 iterations.
989 // FIXME: You will get so much iteratons for fX near mean,
990 // I do not know a better algorithm.
991 bool bfinished = false;
992 do
993 {
994 apl2m = fA + 2.0*rm;
995 d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m);
996 d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0));
997 a1 = (a2+d2m*a1)*fnorm;
998 b1 = (b2+d2m*b1)*fnorm;
999 a2 = a1 + d2m1*a2*fnorm;
1000 b2 = b1 + d2m1*b2*fnorm;
1001 if (b2 != 0.0)
1002 {
1003 fnorm = 1.0/b2;
1004 cfnew = a2*fnorm;
1005 bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);
1006 }
1007 cf = cfnew;
1008 rm += 1.0;
1009 }
1010 while (rm < fMaxIter && !bfinished);
1011 return cf;
1012 }
1013
1014 // cumulative distribution function, normalized
1015 double ScInterpreter::GetBetaDist(double fXin, double fAlpha, double fBeta)
1016 {
1017 // special cases
1018 if (fXin <= 0.0) // values are valid, see spec
1019 return 0.0;
1020 if (fXin >= 1.0) // values are valid, see spec
1021 return 1.0;
1022 if (fBeta == 1.0)
1023 return pow(fXin, fAlpha);
1024 if (fAlpha == 1.0)
1025 // 1.0 - pow(1.0-fX,fBeta) is not accurate enough
1026 return -::rtl::math::expm1(fBeta * ::rtl::math::log1p(-fXin));
1027 //FIXME: need special algorithm for fX near fP for large fA,fB
1028 double fResult;
1029 // I use always continued fraction, power series are neither
1030 // faster nor more accurate.
1031 double fY = (0.5-fXin)+0.5;
1032 double flnY = ::rtl::math::log1p(-fXin);
1033 double fX = fXin;
1034 double flnX = log(fXin);
1035 double fA = fAlpha;
1036 double fB = fBeta;
1037 bool bReflect = fXin > fAlpha/(fAlpha+fBeta);
1038 if (bReflect)
1039 {
1040 fA = fBeta;
1041 fB = fAlpha;
1042 fX = fY;
1043 fY = fXin;
1044 flnX = flnY;
1045 flnY = log(fXin);
1046 }
1047 fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);
1048 fResult = fResult/fA;
1049 double fP = fA/(fA+fB);
1050 double fQ = fB/(fA+fB);
1051 double fTemp;
1052 if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental
1053 fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;
1054 else
1055 fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));
1056 fResult *= fTemp;
1057 if (bReflect)
1058 fResult = 0.5 - fResult + 0.5;
1059 if (fResult > 1.0) // ensure valid range
1060 fResult = 1.0;
1061 if (fResult < 0.0)
1062 fResult = 0.0;
1063 return fResult;
1064 }
1065
1066 void ScInterpreter::ScBetaDist()
1067 {
1068 BYTE nParamCount = GetByte();
1069 if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547#
1070 return;
1071 double fLowerBound, fUpperBound;
1072 double alpha, beta, x;
1073 bool bIsCumulative;
1074 if (nParamCount == 6)
1075 bIsCumulative = GetBool();
1076 else
1077 bIsCumulative = true;
1078 if (nParamCount >= 5)
1079 fUpperBound = GetDouble();
1080 else
1081 fUpperBound = 1.0;
1082 if (nParamCount >= 4)
1083 fLowerBound = GetDouble();
1084 else
1085 fLowerBound = 0.0;
1086 beta = GetDouble();
1087 alpha = GetDouble();
1088 x = GetDouble();
1089 double fScale = fUpperBound - fLowerBound;
1090 if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0)
1091 {
1092 PushIllegalArgument();
1093 return;
1094 }
1095 if (bIsCumulative) // cumulative distribution function
1096 {
1097 // special cases
1098 if (x < fLowerBound)
1099 {
1100 PushDouble(0.0); return; //see spec
1101 }
1102 if (x > fUpperBound)
1103 {
1104 PushDouble(1.0); return; //see spec
1105 }
1106 // normal cases
1107 x = (x-fLowerBound)/fScale; // convert to standard form
1108 PushDouble(GetBetaDist(x, alpha, beta));
1109 return;
1110 }
1111 else // probability density function
1112 {
1113 if (x < fLowerBound || x > fUpperBound)
1114 {
1115 PushDouble(0.0);
1116 return;
1117 }
1118 x = (x-fLowerBound)/fScale;
1119 PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
1120 return;
1121 }
1122 }
1123
1124 void ScInterpreter::ScPhi()
1125 {
1126 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPhi" );
1127 PushDouble(phi(GetDouble()));
1128 }
1129
1130 void ScInterpreter::ScGauss()
1131 {
1132 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGauss" );
1133 PushDouble(gauss(GetDouble()));
1134 }
1135
1136 void ScInterpreter::ScFisher()
1137 {
1138 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFisher" );
1139 double fVal = GetDouble();
1140 if (fabs(fVal) >= 1.0)
1141 PushIllegalArgument();
1142 else
1143 PushDouble( ::rtl::math::atanh( fVal));
1144 }
1145
1146 void ScInterpreter::ScFisherInv()
1147 {
1148 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFisherInv" );
1149 PushDouble( tanh( GetDouble()));
1150 }
1151
1152 void ScInterpreter::ScFact()
1153 {
1154 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFact" );
1155 double nVal = GetDouble();
1156 if (nVal < 0.0)
1157 PushIllegalArgument();
1158 else
1159 PushDouble(Fakultaet(nVal));
1160 }
1161
1162 void ScInterpreter::ScKombin()
1163 {
1164 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScKombin" );
1165 if ( MustHaveParamCount( GetByte(), 2 ) )
1166 {
1167 double k = ::rtl::math::approxFloor(GetDouble());
1168 double n = ::rtl::math::approxFloor(GetDouble());
1169 if (k < 0.0 || n < 0.0 || k > n)
1170 PushIllegalArgument();
1171 else
1172 PushDouble(BinomKoeff(n, k));
1173 }
1174 }
1175
1176 void ScInterpreter::ScKombin2()
1177 {
1178 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScKombin2" );
1179 if ( MustHaveParamCount( GetByte(), 2 ) )
1180 {
1181 double k = ::rtl::math::approxFloor(GetDouble());
1182 double n = ::rtl::math::approxFloor(GetDouble());
1183 if (k < 0.0 || n < 0.0 || k > n)
1184 PushIllegalArgument();
1185 else
1186 PushDouble(BinomKoeff(n + k - 1, k));
1187 }
1188 }
1189
1190 void ScInterpreter::ScVariationen()
1191 {
1192 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScVariationen" );
1193 if ( MustHaveParamCount( GetByte(), 2 ) )
1194 {
1195 double k = ::rtl::math::approxFloor(GetDouble());
1196 double n = ::rtl::math::approxFloor(GetDouble());
1197 if (n < 0.0 || k < 0.0 || k > n)
1198 PushIllegalArgument();
1199 else if (k == 0.0)
1200 PushInt(1); // (n! / (n - 0)!) == 1
1201 else
1202 {
1203 double nVal = n;
1204 for (ULONG i = (ULONG)k-1; i >= 1; i--)
1205 nVal *= n-(double)i;
1206 PushDouble(nVal);
1207 }
1208 }
1209 }
1210
1211 void ScInterpreter::ScVariationen2()
1212 {
1213 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScVariationen2" );
1214 if ( MustHaveParamCount( GetByte(), 2 ) )
1215 {
1216 double k = ::rtl::math::approxFloor(GetDouble());
1217 double n = ::rtl::math::approxFloor(GetDouble());
1218 if (n < 0.0 || k < 0.0 || k > n)
1219 PushIllegalArgument();
1220 else
1221 PushDouble(pow(n,k));
1222 }
1223 }
1224
1225 void ScInterpreter::ScB()
1226 {
1227 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScB" );
1228 BYTE nParamCount = GetByte();
1229 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1230 return ;
1231 if (nParamCount == 3)
1232 {
1233 double x = ::rtl::math::approxFloor(GetDouble());
1234 double p = GetDouble();
1235 double n = ::rtl::math::approxFloor(GetDouble());
1236 if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1237 PushIllegalArgument();
1238 else
1239 {
1240 double q = 1.0 - p;
1241 double fFactor = pow(q, n);
1242 if (fFactor == 0.0)
1243 {
1244 fFactor = pow(p, n);
1245 if (fFactor == 0.0)
1246 PushNoValue();
1247 else
1248 {
1249 ULONG max = (ULONG) (n - x);
1250 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1251 fFactor *= (n-i)/(i+1)*q/p;
1252 PushDouble(fFactor);
1253 }
1254 }
1255 else
1256 {
1257 ULONG max = (ULONG) x;
1258 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1259 fFactor *= (n-i)/(i+1)*p/q;
1260 PushDouble(fFactor);
1261 }
1262 }
1263 }
1264 else if (nParamCount == 4)
1265 {
1266 double xe = GetDouble();
1267 double xs = GetDouble();
1268 double p = GetDouble();
1269 double n = GetDouble();
1270 // alter Stand 300-SC
1271 // if ((xs < n) && (xe < n) && (p < 1.0))
1272 // {
1273 // double Varianz = sqrt(n * p * (1.0 - p));
1274 // xs = fabs(xs - (n * p /* / 2.0 STE */ ));
1275 // xe = fabs(xe - (n * p /* / 2.0 STE */ ));
1276 //// STE double nVal = gauss((xs + 0.5) / Varianz) + gauss((xe + 0.5) / Varianz);
1277 // double nVal = fabs(gauss(xs / Varianz) - gauss(xe / Varianz));
1278 // PushDouble(nVal);
1279 // }
1280 if (xe <= n && xs <= xe &&
1281 p < 1.0 && p > 0.0 && n >= 0.0 && xs >= 0.0 )
1282 {
1283 double q = 1.0 - p;
1284 double fFactor = pow(q, n);
1285 if (fFactor == 0.0)
1286 {
1287 fFactor = pow(p, n);
1288 if (fFactor == 0.0)
1289 PushNoValue();
1290 else
1291 {
1292 double fSum = 0.0;
1293 ULONG max;
1294 if (xe < (ULONG) n)
1295 max = (ULONG) (n-xe)-1;
1296 else
1297 max = 0;
1298 ULONG i;
1299 for (i = 0; i < max && fFactor > 0.0; i++)
1300 fFactor *= (n-i)/(i+1)*q/p;
1301 if (xs < (ULONG) n)
1302 max = (ULONG) (n-xs);
1303 else
1304 fSum = fFactor;
1305 for (; i < max && fFactor > 0.0; i++)
1306 {
1307 fFactor *= (n-i)/(i+1)*q/p;
1308 fSum += fFactor;
1309 }
1310 PushDouble(fSum);
1311 }
1312 }
1313 else
1314 {
1315 ULONG max;
1316 double fSum;
1317 if ( (ULONG) xs == 0)
1318 {
1319 fSum = fFactor;
1320 max = 0;
1321 }
1322 else
1323 {
1324 max = (ULONG) xs-1;
1325 fSum = 0.0;
1326 }
1327 ULONG i;
1328 for (i = 0; i < max && fFactor > 0.0; i++)
1329 fFactor *= (n-i)/(i+1)*p/q;
1330 if ((ULONG)xe == 0) // beide 0
1331 fSum = fFactor;
1332 else
1333 max = (ULONG) xe;
1334 for (; i < max && fFactor > 0.0; i++)
1335 {
1336 fFactor *= (n-i)/(i+1)*p/q;
1337 fSum += fFactor;
1338 }
1339 PushDouble(fSum);
1340 }
1341 }
1342 else
1343 PushIllegalArgument();
1344 }
1345 }
1346
1347 void ScInterpreter::ScBinomDist()
1348 {
1349 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScBinomDist" );
1350 if ( MustHaveParamCount( GetByte(), 4 ) )
1351 {
1352 double kum = GetDouble(); // 0 oder 1
1353 double p = GetDouble(); // p
1354 double n = ::rtl::math::approxFloor(GetDouble()); // n
1355 double x = ::rtl::math::approxFloor(GetDouble()); // x
1356 double fFactor, q, fSum;
1357 if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1358 PushIllegalArgument();
1359 else if (kum == 0.0) // Dichte
1360 {
1361 q = 1.0 - p;
1362 fFactor = pow(q, n);
1363 if (fFactor == 0.0)
1364 {
1365 fFactor = pow(p, n);
1366 if (fFactor == 0.0)
1367 PushNoValue();
1368 else
1369 {
1370 ULONG max = (ULONG) (n - x);
1371 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1372 fFactor *= (n-i)/(i+1)*q/p;
1373 PushDouble(fFactor);
1374 }
1375 }
1376 else
1377 {
1378 ULONG max = (ULONG) x;
1379 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1380 fFactor *= (n-i)/(i+1)*p/q;
1381 PushDouble(fFactor);
1382 }
1383 }
1384 else // Verteilung
1385 {
1386 if (n == x)
1387 PushDouble(1.0);
1388 else
1389 {
1390 q = 1.0 - p;
1391 fFactor = pow(q, n);
1392 if (fFactor == 0.0)
1393 {
1394 fFactor = pow(p, n);
1395 if (fFactor == 0.0)
1396 PushNoValue();
1397 else
1398 {
1399 fSum = 1.0 - fFactor;
1400 ULONG max = (ULONG) (n - x) - 1;
1401 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1402 {
1403 fFactor *= (n-i)/(i+1)*q/p;
1404 fSum -= fFactor;
1405 }
1406 if (fSum < 0.0)
1407 PushDouble(0.0);
1408 else
1409 PushDouble(fSum);
1410 }
1411 }
1412 else
1413 {
1414 fSum = fFactor;
1415 ULONG max = (ULONG) x;
1416 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1417 {
1418 fFactor *= (n-i)/(i+1)*p/q;
1419 fSum += fFactor;
1420 }
1421 PushDouble(fSum);
1422 }
1423 }
1424 }
1425 }
1426 }
1427
1428 void ScInterpreter::ScCritBinom()
1429 {
1430 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScCritBinom" );
1431 if ( MustHaveParamCount( GetByte(), 3 ) )
1432 {
1433 double alpha = GetDouble(); // alpha
1434 double p = GetDouble(); // p
1435 double n = ::rtl::math::approxFloor(GetDouble());
1436 if (n < 0.0 || alpha <= 0.0 || alpha >= 1.0 || p < 0.0 || p > 1.0)
1437 PushIllegalArgument();
1438 else
1439 {
1440 double q = 1.0 - p;
1441 double fFactor = pow(q,n);
1442 if (fFactor == 0.0)
1443 {
1444 fFactor = pow(p, n);
1445 if (fFactor == 0.0)
1446 PushNoValue();
1447 else
1448 {
1449 double fSum = 1.0 - fFactor; ULONG max = (ULONG) n;
1450 ULONG i;
1451
1452 for ( i = 0; i < max && fSum >= alpha; i++)
1453 {
1454 fFactor *= (n-i)/(i+1)*q/p;
1455 fSum -= fFactor;
1456 }
1457 PushDouble(n-i);
1458 }
1459 }
1460 else
1461 {
1462 double fSum = fFactor; ULONG max = (ULONG) n;
1463 ULONG i;
1464
1465 for ( i = 0; i < max && fSum < alpha; i++)
1466 {
1467 fFactor *= (n-i)/(i+1)*p/q;
1468 fSum += fFactor;
1469 }
1470 PushDouble(i);
1471 }
1472 }
1473 }
1474 }
1475
1476 void ScInterpreter::ScNegBinomDist()
1477 {
1478 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScNegBinomDist" );
1479 if ( MustHaveParamCount( GetByte(), 3 ) )
1480 {
1481 double p = GetDouble(); // p
1482 double r = GetDouble(); // r
1483 double x = GetDouble(); // x
1484 if (r < 0.0 || x < 0.0 || p < 0.0 || p > 1.0)
1485 PushIllegalArgument();
1486 else
1487 {
1488 double q = 1.0 - p;
1489 double fFactor = pow(p,r);
1490 for (double i = 0.0; i < x; i++)
1491 fFactor *= (i+r)/(i+1.0)*q;
1492 PushDouble(fFactor);
1493 }
1494 }
1495 }
1496
1497 void ScInterpreter::ScNormDist()
1498 {
1499 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScNormDist" );
1500 if ( MustHaveParamCount( GetByte(), 4 ) )
1501 {
1502 double kum = GetDouble(); // 0 oder 1
1503 double sigma = GetDouble(); // Stdabw
1504 double mue = GetDouble(); // Mittelwert
1505 double x = GetDouble(); // x
1506 if (sigma < 0.0)
1507 PushError( errIllegalArgument);
1508 else if (sigma == 0.0)
1509 PushError( errDivisionByZero);
1510 else if (kum == 0.0) // Dichte
1511 PushDouble(phi((x-mue)/sigma)/sigma);
1512 else // Verteilung
1513 PushDouble(0.5 + gauss((x-mue)/sigma));
1514 }
1515 }
1516
1517 void ScInterpreter::ScLogNormDist()
1518 {
1519 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScLogNormDist" );
1520 if ( MustHaveParamCount( GetByte(), 3 ) )
1521 {
1522 double sigma = GetDouble(); // Stdabw
1523 double mue = GetDouble(); // Mittelwert
1524 double x = GetDouble(); // x
1525 if (sigma < 0.0)
1526 PushError( errIllegalArgument);
1527 else if (sigma == 0.0)
1528 PushError( errDivisionByZero);
1529 else if (x <= 0.0)
1530 PushIllegalArgument();
1531 else
1532 PushDouble(0.5 + gauss((log(x)-mue)/sigma));
1533 }
1534 }
1535
1536 void ScInterpreter::ScStdNormDist()
1537 {
1538 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScStdNormDist" );
1539 PushDouble(0.5 + gauss(GetDouble()));
1540 }
1541
1542 void ScInterpreter::ScExpDist()
1543 {
1544 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScExpDist" );
1545 if ( MustHaveParamCount( GetByte(), 3 ) )
1546 {
1547 double kum = GetDouble(); // 0 oder 1
1548 double lambda = GetDouble(); // lambda
1549 double x = GetDouble(); // x
1550 if (lambda <= 0.0)
1551 PushIllegalArgument();
1552 else if (kum == 0.0) // Dichte
1553 {
1554 if (x >= 0.0)
1555 PushDouble(lambda * exp(-lambda*x));
1556 else
1557 PushInt(0);
1558 }
1559 else // Verteilung
1560 {
1561 if (x > 0.0)
1562 PushDouble(1.0 - exp(-lambda*x));
1563 else
1564 PushInt(0);
1565 }
1566 }
1567 }
1568
1569 void ScInterpreter::ScTDist()
1570 {
1571 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTDist" );
1572 if ( !MustHaveParamCount( GetByte(), 3 ) )
1573 return;
1574 double fFlag = ::rtl::math::approxFloor(GetDouble());
1575 double fDF = ::rtl::math::approxFloor(GetDouble());
1576 double T = GetDouble();
1577 if (fDF < 1.0 || T < 0.0 || (fFlag != 1.0 && fFlag != 2.0) )
1578 {
1579 PushIllegalArgument();
1580 return;
1581 }
1582 double R = GetTDist(T, fDF);
1583 if (fFlag == 1.0)
1584 PushDouble(R);
1585 else
1586 PushDouble(2.0*R);
1587 }
1588
1589 void ScInterpreter::ScFDist()
1590 {
1591 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFDist" );
1592 if ( !MustHaveParamCount( GetByte(), 3 ) )
1593 return;
1594 double fF2 = ::rtl::math::approxFloor(GetDouble());
1595 double fF1 = ::rtl::math::approxFloor(GetDouble());
1596 double fF = GetDouble();
1597 if (fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10)
1598 {
1599 PushIllegalArgument();
1600 return;
1601 }
1602 PushDouble(GetFDist(fF, fF1, fF2));
1603 }
1604
1605 void ScInterpreter::ScChiDist()
1606 {
1607 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScChiDist" );
1608 double fResult;
1609 if ( !MustHaveParamCount( GetByte(), 2 ) )
1610 return;
1611 double fDF = ::rtl::math::approxFloor(GetDouble());
1612 double fChi = GetDouble();
1613 if (fDF < 1.0) // x<=0 returns 1, see ODFF 6.17.10
1614 {
1615 PushIllegalArgument();
1616 return;
1617 }
1618 fResult = GetChiDist( fChi, fDF);
1619 if (nGlobalError)
1620 {
1621 PushError( nGlobalError);
1622 return;
1623 }
1624 PushDouble(fResult);
1625 }
1626
1627 void ScInterpreter::ScWeibull()
1628 {
1629 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScWeibull" );
1630 if ( MustHaveParamCount( GetByte(), 4 ) )
1631 {
1632 double kum = GetDouble(); // 0 oder 1
1633 double beta = GetDouble(); // beta
1634 double alpha = GetDouble(); // alpha
1635 double x = GetDouble(); // x
1636 if (alpha <= 0.0 || beta <= 0.0 || x < 0.0)
1637 PushIllegalArgument();
1638 else if (kum == 0.0) // Dichte
1639 PushDouble(alpha/pow(beta,alpha)*pow(x,alpha-1.0)*
1640 exp(-pow(x/beta,alpha)));
1641 else // Verteilung
1642 PushDouble(1.0 - exp(-pow(x/beta,alpha)));
1643 }
1644 }
1645
1646 void ScInterpreter::ScPoissonDist()
1647 {
1648 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPoissonDist" );
1649 BYTE nParamCount = GetByte();
1650 if ( MustHaveParamCount( nParamCount, 2, 3 ) )
1651 {
1652 double kum = (nParamCount == 3 ? GetDouble() : 1); // 0 oder 1
1653 double lambda = GetDouble(); // Mittelwert
1654 double x = ::rtl::math::approxFloor(GetDouble()); // x
1655 if (lambda < 0.0 || x < 0.0)
1656 PushIllegalArgument();
1657 else if (kum == 0.0) // Dichte
1658 {
1659 if (lambda == 0.0)
1660 PushInt(0);
1661 else
1662 {
1663 double fPoissonVar = 1.0;
1664 for ( double f = 0.0; f < x; ++f )
1665 fPoissonVar *= lambda / ( f + 1.0 );
1666 PushDouble( fPoissonVar*exp( -lambda ) );
1667 }
1668 }
1669 else // Verteilung
1670 {
1671 if (lambda == 0.0)
1672 PushInt(1);
1673 else
1674 {
1675 double sum = 1.0;
1676 double fFak = 1.0;
1677 ULONG nEnd = (ULONG) x;
1678 for (ULONG i = 1; i <= nEnd; i++)
1679 {
1680 fFak *= (double)i;
1681 sum += pow( lambda, (double)i ) / fFak;
1682 }
1683 sum *= exp(-lambda);
1684 PushDouble(sum);
1685 }
1686 }
1687 }
1688 }
1689
1690 /** Local function used in the calculation of the hypergeometric distribution.
1691 */
1692 void lcl_PutFactorialElements( ::std::vector< double >& cn, double fLower, double fUpper, double fBase )
1693 {
1694 for ( double i = fLower; i <= fUpper; ++i )
1695 {
1696 double fVal = fBase - i;
1697 if ( fVal > 1.0 )
1698 cn.push_back( fVal );
1699 }
1700 }
1701
1702 /** Calculates a value of the hypergeometric distribution.
1703
1704 The algorithm is designed to avoid unnecessary multiplications and division
1705 by expanding all factorial elements (9 of them). It is done by excluding
1706 those ranges that overlap in the numerator and the denominator. This allows
1707 for a fast calculation for large values which would otherwise cause an overflow
1708 in the intermediate values.
1709
1710 @author Kohei Yoshida <kohei@openoffice.org>
1711
1712 @see #i47296#
1713
1714 */
1715 void ScInterpreter::ScHypGeomDist()
1716 {
1717 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScHypGeomDist" );
1718 const size_t nMaxArraySize = 500000; // arbitrary max array size
1719
1720 if ( !MustHaveParamCount( GetByte(), 4 ) )
1721 return;
1722
1723 double N = ::rtl::math::approxFloor(GetDouble());
1724 double M = ::rtl::math::approxFloor(GetDouble());
1725 double n = ::rtl::math::approxFloor(GetDouble());
1726 double x = ::rtl::math::approxFloor(GetDouble());
1727
1728 if( (x < 0.0) || (n < x) || (M < x) || (N < n) || (N < M) || (x < n - N + M) )
1729 {
1730 PushIllegalArgument();
1731 return;
1732 }
1733
1734 typedef ::std::vector< double > HypContainer;
1735 HypContainer cnNumer, cnDenom;
1736
1737 size_t nEstContainerSize = static_cast<size_t>( x + ::std::min( n, M ) );
1738 size_t nMaxSize = ::std::min( cnNumer.max_size(), nMaxArraySize );
1739 if ( nEstContainerSize > nMaxSize )
1740 {
1741 PushNoValue();
1742 return;
1743 }
1744 cnNumer.reserve( nEstContainerSize + 10 );
1745 cnDenom.reserve( nEstContainerSize + 10 );
1746
1747 // Trim coefficient C first
1748 double fCNumVarUpper = N - n - M + x - 1.0;
1749 double fCDenomVarLower = 1.0;
1750 if ( N - n - M + x >= M - x + 1.0 )
1751 {
1752 fCNumVarUpper = M - x - 1.0;
1753 fCDenomVarLower = N - n - 2.0*(M - x) + 1.0;
1754 }
1755
1756 #ifdef DBG_UTIL
1757 double fCNumLower = N - n - fCNumVarUpper;
1758 #endif
1759 double fCDenomUpper = N - n - M + x + 1.0 - fCDenomVarLower;
1760
1761 double fDNumVarLower = n - M;
1762
1763 if ( n >= M + 1.0 )
1764 {
1765 if ( N - M < n + 1.0 )
1766 {
1767 // Case 1
1768
1769 if ( N - n < n + 1.0 )
1770 {
1771 // no overlap
1772 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1773 lcl_PutFactorialElements( cnDenom, 0.0, N - n - 1.0, N );
1774 }
1775 else
1776 {
1777 // overlap
1778 DBG_ASSERT( fCNumLower < n + 1.0, "ScHypGeomDist: wrong assertion" );
1779 lcl_PutFactorialElements( cnNumer, N - 2.0*n, fCNumVarUpper, N - n );
1780 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1781 }
1782
1783 DBG_ASSERT( fCDenomUpper <= N - M, "ScHypGeomDist: wrong assertion" );
1784
1785 if ( fCDenomUpper < n - x + 1.0 )
1786 // no overlap
1787 lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
1788 else
1789 {
1790 // overlap
1791 lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1792
1793 fCDenomUpper = n - x;
1794 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1795 }
1796 }
1797 else
1798 {
1799 // Case 2
1800
1801 if ( n > M - 1.0 )
1802 {
1803 // no overlap
1804 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1805 lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
1806 }
1807 else
1808 {
1809 lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
1810 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1811 }
1812
1813 DBG_ASSERT( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
1814
1815 if ( fCDenomUpper < n - x + 1.0 )
1816 // no overlap
1817 lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - n + x, N - M + 1.0 );
1818 else
1819 {
1820 lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1821 fCDenomUpper = n - x;
1822 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1823 }
1824 }
1825
1826 DBG_ASSERT( fCDenomUpper <= M, "ScHypGeomDist: wrong assertion" );
1827 }
1828 else
1829 {
1830 if ( N - M < M + 1.0 )
1831 {
1832 // Case 3
1833
1834 if ( N - n < M + 1.0 )
1835 {
1836 // No overlap
1837 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1838 lcl_PutFactorialElements( cnDenom, 0.0, N - M - 1.0, N );
1839 }
1840 else
1841 {
1842 lcl_PutFactorialElements( cnNumer, N - n - M, fCNumVarUpper, N - n );
1843 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1844 }
1845
1846 if ( n - x + 1.0 > fCDenomUpper )
1847 // No overlap
1848 lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
1849 else
1850 {
1851 // Overlap
1852 lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1853
1854 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1855 fCDenomUpper = n - x;
1856 }
1857 }
1858 else
1859 {
1860 // Case 4
1861
1862 DBG_ASSERT( M >= n - x, "ScHypGeomDist: wrong assertion" );
1863 DBG_ASSERT( M - x <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
1864
1865 if ( N - n < N - M + 1.0 )
1866 {
1867 // No overlap
1868 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1869 lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
1870 }
1871 else
1872 {
1873 // Overlap
1874 DBG_ASSERT( fCNumLower <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
1875
1876 lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
1877 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1878 }
1879
1880 if ( n - x + 1.0 > fCDenomUpper )
1881 // No overlap
1882 lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - n + x, N - M + 1.0 );
1883 else if ( M >= fCDenomUpper )
1884 {
1885 lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1886
1887 fCDenomUpper = n - x;
1888 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1889 }
1890 else
1891 {
1892 DBG_ASSERT( M <= fCDenomUpper, "ScHypGeomDist: wrong assertion" );
1893 lcl_PutFactorialElements( cnDenom, fCDenomVarLower, N - n - 2.0*M + x,
1894 N - n - M + x + 1.0 );
1895
1896 fCDenomUpper = n - x;
1897 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1898 }
1899 }
1900
1901 DBG_ASSERT( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
1902
1903 fDNumVarLower = 0.0;
1904 }
1905
1906 double nDNumVarUpper = fCDenomUpper < x + 1.0 ? n - x - 1.0 : n - fCDenomUpper - 1.0;
1907 double nDDenomVarLower = fCDenomUpper < x + 1.0 ? fCDenomVarLower : N - n - M + 1.0;
1908 lcl_PutFactorialElements( cnNumer, fDNumVarLower, nDNumVarUpper, n );
1909 lcl_PutFactorialElements( cnDenom, nDDenomVarLower, N - n - M + x, N - n - M + x + 1.0 );
1910
1911 ::std::sort( cnNumer.begin(), cnNumer.end() );
1912 ::std::sort( cnDenom.begin(), cnDenom.end() );
1913 HypContainer::reverse_iterator it1 = cnNumer.rbegin(), it1End = cnNumer.rend();
1914 HypContainer::reverse_iterator it2 = cnDenom.rbegin(), it2End = cnDenom.rend();
1915
1916 double fFactor = 1.0;
1917 for ( ; it1 != it1End || it2 != it2End; )
1918 {
1919 double fEnum = 1.0, fDenom = 1.0;
1920 if ( it1 != it1End )
1921 fEnum = *it1++;
1922 if ( it2 != it2End )
1923 fDenom = *it2++;
1924 fFactor *= fEnum / fDenom;
1925 }
1926
1927 PushDouble(fFactor);
1928 }
1929
1930 void ScInterpreter::ScGammaDist()
1931 {
1932 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGammaDist" );
1933 BYTE nParamCount = GetByte();
1934 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1935 return;
1936 double bCumulative;
1937 if (nParamCount == 4)
1938 bCumulative = GetBool();
1939 else
1940 bCumulative = true;
1941 double fBeta = GetDouble(); // scale
1942 double fAlpha = GetDouble(); // shape
1943 double fX = GetDouble(); // x
1944 if (fAlpha <= 0.0 || fBeta <= 0.0)
1945 PushIllegalArgument();
1946 else
1947 {
1948 if (bCumulative) // distribution
1949 PushDouble( GetGammaDist( fX, fAlpha, fBeta));
1950 else // density
1951 PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta));
1952 }
1953 }
1954
1955 void ScInterpreter::ScNormInv()
1956 {
1957 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScNormInv" );
1958 if ( MustHaveParamCount( GetByte(), 3 ) )
1959 {
1960 double sigma = GetDouble();
1961 double mue = GetDouble();
1962 double x = GetDouble();
1963 if (sigma <= 0.0 || x < 0.0 || x > 1.0)
1964 PushIllegalArgument();
1965 else if (x == 0.0 || x == 1.0)
1966 PushNoValue();
1967 else
1968 PushDouble(gaussinv(x)*sigma + mue);
1969 }
1970 }
1971
1972 void ScInterpreter::ScSNormInv()
1973 {
1974 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSNormInv" );
1975 double x = GetDouble();
1976 if (x < 0.0 || x > 1.0)
1977 PushIllegalArgument();
1978 else if (x == 0.0 || x == 1.0)
1979 PushNoValue();
1980 else
1981 PushDouble(gaussinv(x));
1982 }
1983
1984 void ScInterpreter::ScLogNormInv()
1985 {
1986 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScLogNormInv" );
1987 if ( MustHaveParamCount( GetByte(), 3 ) )
1988 {
1989 double sigma = GetDouble(); // Stdabw
1990 double mue = GetDouble(); // Mittelwert
1991 double y = GetDouble(); // y
1992 if (sigma <= 0.0 || y <= 0.0 || y >= 1.0)
1993 PushIllegalArgument();
1994 else
1995 PushDouble(exp(mue+sigma*gaussinv(y)));
1996 }
1997 }
1998
1999 class ScGammaDistFunction : public ScDistFunc
2000 {
2001 ScInterpreter& rInt;
2002 double fp, fAlpha, fBeta;
2003
2004 public:
2005 ScGammaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2006 rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2007
2008 double GetValue( double x ) const { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); }
2009 };
2010
2011 void ScInterpreter::ScGammaInv()
2012 {
2013 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGammaInv" );
2014 if ( !MustHaveParamCount( GetByte(), 3 ) )
2015 return;
2016 double fBeta = GetDouble();
2017 double fAlpha = GetDouble();
2018 double fP = GetDouble();
2019 if (fAlpha <= 0.0 || fBeta <= 0.0 || fP < 0.0 || fP >= 1.0 )
2020 {
2021 PushIllegalArgument();
2022 return;
2023 }
2024 if (fP == 0.0)
2025 PushInt(0);
2026 else
2027 {
2028 bool bConvError;
2029 ScGammaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2030 double fStart = fAlpha * fBeta;
2031 double fVal = lcl_IterateInverse( aFunc, fStart*0.5, fStart, bConvError );
2032 if (bConvError)
2033 SetError(errNoConvergence);
2034 PushDouble(fVal);
2035 }
2036 }
2037
2038 class ScBetaDistFunction : public ScDistFunc
2039 {
2040 ScInterpreter& rInt;
2041 double fp, fAlpha, fBeta;
2042
2043 public:
2044 ScBetaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2045 rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2046
2047 double GetValue( double x ) const { return fp - rInt.GetBetaDist(x, fAlpha, fBeta); }
2048 };
2049
2050 void ScInterpreter::ScBetaInv()
2051 {
2052 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScBetaInv" );
2053 BYTE nParamCount = GetByte();
2054 if ( !MustHaveParamCount( nParamCount, 3, 5 ) )
2055 return;
2056 double fP, fA, fB, fAlpha, fBeta;
2057 if (nParamCount == 5)
2058 fB = GetDouble();
2059 else
2060 fB = 1.0;
2061 if (nParamCount >= 4)
2062 fA = GetDouble();
2063 else
2064 fA = 0.0;
2065 fBeta = GetDouble();
2066 fAlpha = GetDouble();
2067 fP = GetDouble();
2068 if (fP < 0.0 || fP >= 1.0 || fA == fB || fAlpha <= 0.0 || fBeta <= 0.0)
2069 {
2070 PushIllegalArgument();
2071 return;
2072 }
2073 if (fP == 0.0)
2074 PushInt(0);
2075 else
2076 {
2077 bool bConvError;
2078 ScBetaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2079 // 0..1 as range for iteration so it isn't extended beyond the valid range
2080 double fVal = lcl_IterateInverse( aFunc, 0.0, 1.0, bConvError );
2081 if (bConvError)
2082 PushError( errNoConvergence);
2083 else
2084 PushDouble(fA + fVal*(fB-fA)); // scale to (A,B)
2085 }
2086 }
2087
2088 // Achtung: T, F und Chi
2089 // sind monoton fallend,
2090 // deshalb 1-Dist als Funktion
2091
2092 class ScTDistFunction : public ScDistFunc
2093 {
2094 ScInterpreter& rInt;
2095 double fp, fDF;
2096
2097 public:
2098 ScTDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2099 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2100
2101 double GetValue( double x ) const { return fp - 2 * rInt.GetTDist(x, fDF); }
2102 };
2103
2104 void ScInterpreter::ScTInv()
2105 {
2106 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTInv" );
2107 if ( !MustHaveParamCount( GetByte(), 2 ) )
2108 return;
2109 double fDF = ::rtl::math::approxFloor(GetDouble());
2110 double fP = GetDouble();
2111 if (fDF < 1.0 || fDF >= 1.0E5 || fP <= 0.0 || fP > 1.0 )
2112 {
2113 PushIllegalArgument();
2114 return;
2115 }
2116
2117 bool bConvError;
2118 ScTDistFunction aFunc( *this, fP, fDF );
2119 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2120 if (bConvError)
2121 SetError(errNoConvergence);
2122 PushDouble(fVal);
2123 }
2124
2125 class ScFDistFunction : public ScDistFunc
2126 {
2127 ScInterpreter& rInt;
2128 double fp, fF1, fF2;
2129
2130 public:
2131 ScFDistFunction( ScInterpreter& rI, double fpVal, double fF1Val, double fF2Val ) :
2132 rInt(rI), fp(fpVal), fF1(fF1Val), fF2(fF2Val) {}
2133
2134 double GetValue( double x ) const { return fp - rInt.GetFDist(x, fF1, fF2); }
2135 };
2136
2137 void ScInterpreter::ScFInv()
2138 {
2139 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFInv" );
2140 if ( !MustHaveParamCount( GetByte(), 3 ) )
2141 return;
2142 double fF2 = ::rtl::math::approxFloor(GetDouble());
2143 double fF1 = ::rtl::math::approxFloor(GetDouble());
2144 double fP = GetDouble();
2145 if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
2146 {
2147 PushIllegalArgument();
2148 return;
2149 }
2150
2151 bool bConvError;
2152 ScFDistFunction aFunc( *this, fP, fF1, fF2 );
2153 double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
2154 if (bConvError)
2155 SetError(errNoConvergence);
2156 PushDouble(fVal);
2157 }
2158
2159 class ScChiDistFunction : public ScDistFunc
2160 {
2161 ScInterpreter& rInt;
2162 double fp, fDF;
2163
2164 public:
2165 ScChiDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2166 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2167
2168 double GetValue( double x ) const { return fp - rInt.GetChiDist(x, fDF); }
2169 };
2170
2171 void ScInterpreter::ScChiInv()
2172 {
2173 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScChiInv" );
2174 if ( !MustHaveParamCount( GetByte(), 2 ) )
2175 return;
2176 double fDF = ::rtl::math::approxFloor(GetDouble());
2177 double fP = GetDouble();
2178 if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
2179 {
2180 PushIllegalArgument();
2181 return;
2182 }
2183
2184 bool bConvError;
2185 ScChiDistFunction aFunc( *this, fP, fDF );
2186 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2187 if (bConvError)
2188 SetError(errNoConvergence);
2189 PushDouble(fVal);
2190 }
2191
2192 /***********************************************/
2193 class ScChiSqDistFunction : public ScDistFunc
2194 {
2195 ScInterpreter& rInt;
2196 double fp, fDF;
2197
2198 public:
2199 ScChiSqDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2200 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2201
2202 double GetValue( double x ) const { return fp - rInt.GetChiSqDistCDF(x, fDF); }
2203 };
2204
2205
2206 void ScInterpreter::ScChiSqInv()
2207 {
2208 if ( !MustHaveParamCount( GetByte(), 2 ) )
2209 return;
2210 double fDF = ::rtl::math::approxFloor(GetDouble());
2211 double fP = GetDouble();
2212 if (fDF < 1.0 || fP < 0.0 || fP >= 1.0 )
2213 {
2214 PushIllegalArgument();
2215 return;
2216 }
2217
2218 bool bConvError;
2219 ScChiSqDistFunction aFunc( *this, fP, fDF );
2220 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2221 if (bConvError)
2222 SetError(errNoConvergence);
2223 PushDouble(fVal);
2224 }
2225
2226
2227 void ScInterpreter::ScConfidence()
2228 {
2229 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScConfidence" );
2230 if ( MustHaveParamCount( GetByte(), 3 ) )
2231 {
2232 double n = ::rtl::math::approxFloor(GetDouble());
2233 double sigma = GetDouble();
2234 double alpha = GetDouble();
2235 if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
2236 PushIllegalArgument();
2237 else
2238 PushDouble( gaussinv(1.0-alpha/2.0) * sigma/sqrt(n) );
2239 }
2240 }
2241
2242 void ScInterpreter::ScZTest()
2243 {
2244 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScZTest" );
2245 BYTE nParamCount = GetByte();
2246 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
2247 return;
2248 double sigma = 0.0, mue, x;
2249 if (nParamCount == 3)
2250 {
2251 sigma = GetDouble();
2252 if (sigma <= 0.0)
2253 {
2254 PushIllegalArgument();
2255 return;
2256 }
2257 }
2258 x = GetDouble();
2259
2260 double fSum = 0.0;
2261 double fSumSqr = 0.0;
2262 double fVal;
2263 double rValCount = 0.0;
2264 switch (GetStackType())
2265 {
2266 case formula::svDouble :
2267 {
2268 fVal = GetDouble();
2269 fSum += fVal;
2270 fSumSqr += fVal*fVal;
2271 rValCount++;
2272 }
2273 break;
2274 case svSingleRef :
2275 {
2276 ScAddress aAdr;
2277 PopSingleRef( aAdr );
2278 ScBaseCell* pCell = GetCell( aAdr );
2279 if (HasCellValueData(pCell))
2280 {
2281 fVal = GetCellValue( aAdr, pCell );
2282 fSum += fVal;
2283 fSumSqr += fVal*fVal;
2284 rValCount++;
2285 }
2286 }
2287 break;
2288 case svRefList :
2289 case formula::svDoubleRef :
2290 {
2291 short nParam = 1;
2292 size_t nRefInList = 0;
2293 while (nParam-- > 0)
2294 {
2295 ScRange aRange;
2296 USHORT nErr = 0;
2297 PopDoubleRef( aRange, nParam, nRefInList);
2298 ScValueIterator aValIter(pDok, aRange, glSubTotal);
2299 if (aValIter.GetFirst(fVal, nErr))
2300 {
2301 fSum += fVal;
2302 fSumSqr += fVal*fVal;
2303 rValCount++;
2304 while ((nErr == 0) && aValIter.GetNext(fVal, nErr))
2305 {
2306 fSum += fVal;
2307 fSumSqr += fVal*fVal;
2308 rValCount++;
2309 }
2310 SetError(nErr);
2311 }
2312 }
2313 }
2314 break;
2315 case svMatrix :
2316 {
2317 ScMatrixRef pMat = PopMatrix();
2318 if (pMat)
2319 {
2320 SCSIZE nCount = pMat->GetElementCount();
2321 if (pMat->IsNumeric())
2322 {
2323 for ( SCSIZE i = 0; i < nCount; i++ )
2324 {
2325 fVal= pMat->GetDouble(i);
2326 fSum += fVal;
2327 fSumSqr += fVal * fVal;
2328 rValCount++;
2329 }
2330 }
2331 else
2332 {
2333 for (SCSIZE i = 0; i < nCount; i++)
2334 if (!pMat->IsString(i))
2335 {
2336 fVal= pMat->GetDouble(i);
2337 fSum += fVal;
2338 fSumSqr += fVal * fVal;
2339 rValCount++;
2340 }
2341 }
2342 }
2343 }
2344 break;
2345 default : SetError(errIllegalParameter); break;
2346 }
2347 if (rValCount <= 1.0)
2348 PushError( errDivisionByZero);
2349 else
2350 {
2351 mue = fSum/rValCount;
2352 if (nParamCount != 3)
2353 sigma = (fSumSqr - fSum*fSum/rValCount)/(rValCount-1.0);
2354
2355 PushDouble(0.5 - gauss((mue-x)/sqrt(sigma/rValCount)));
2356 }
2357 }
2358 bool ScInterpreter::CalculateTest(BOOL _bTemplin
2359 ,const SCSIZE nC1, const SCSIZE nC2,const SCSIZE nR1,const SCSIZE nR2
2360 ,const ScMatrixRef& pMat1,const ScMatrixRef& pMat2
2361 ,double& fT,double& fF)
2362 {
2363 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateTest" );
2364 double fCount1 = 0.0;
2365 double fCount2 = 0.0;
2366 double fSum1 = 0.0;
2367 double fSumSqr1 = 0.0;
2368 double fSum2 = 0.0;
2369 double fSumSqr2 = 0.0;
2370 double fVal;
2371 SCSIZE i,j;
2372 for (i = 0; i < nC1; i++)
2373 for (j = 0; j < nR1; j++)
2374 {
2375 if (!pMat1->IsString(i,j))
2376 {
2377 fVal = pMat1->GetDouble(i,j);
2378 fSum1 += fVal;
2379 fSumSqr1 += fVal * fVal;
2380 fCount1++;
2381 }
2382 }
2383 for (i = 0; i < nC2; i++)
2384 for (j = 0; j < nR2; j++)
2385 {
2386 if (!pMat2->IsString(i,j))
2387 {
2388 fVal = pMat2->GetDouble(i,j);
2389 fSum2 += fVal;
2390 fSumSqr2 += fVal * fVal;
2391 fCount2++;
2392 }
2393 }
2394 if (fCount1 < 2.0 || fCount2 < 2.0)
2395 {
2396 PushNoValue();
2397 return false;
2398 } // if (fCount1 < 2.0 || fCount2 < 2.0)
2399 if ( _bTemplin )
2400 {
2401 double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0)/fCount1;
2402 double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0)/fCount2;
2403 if (fS1 + fS2 == 0.0)
2404 {
2405 PushNoValue();
2406 return false;
2407 }
2408 fT = fabs(fSum1/fCount1 - fSum2/fCount2)/sqrt(fS1+fS2);
2409 double c = fS1/(fS1+fS2);
2410 // s.u. fF = ::rtl::math::approxFloor(1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0)));
2411 // fF = ::rtl::math::approxFloor((fS1+fS2)*(fS1+fS2)/(fS1*fS1/(fCount1-1.0) + fS2*fS2/(fCount2-1.0)));
2412
2413 // GetTDist wird mit GetBetaDist berechnet und kommt auch mit nicht ganzzahligen
2414 // Freiheitsgraden klar. Dann stimmt das Ergebnis auch mit Excel ueberein (#52406#):
2415 fF = 1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0));
2416 }
2417 else
2418 {
2419 // laut Bronstein-Semendjajew
2420 double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1) / (fCount1 - 1.0); // Varianz
2421 double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2) / (fCount2 - 1.0);
2422 fT = fabs( fSum1/fCount1 - fSum2/fCount2 ) /
2423 sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
2424 sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
2425 fF = fCount1 + fCount2 - 2;
2426 }
2427 return true;
2428 }
2429 void ScInterpreter::ScTTest()
2430 {
2431 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTTest" );
2432 if ( !MustHaveParamCount( GetByte(), 4 ) )
2433 return;
2434 double fTyp = ::rtl::math::approxFloor(GetDouble());
2435 double fAnz = ::rtl::math::approxFloor(GetDouble());
2436 if (fAnz != 1.0 && fAnz != 2.0)
2437 {
2438 PushIllegalArgument();
2439 return;
2440 }
2441
2442 ScMatrixRef pMat2 = GetMatrix();
2443 ScMatrixRef pMat1 = GetMatrix();
2444 if (!pMat1 || !pMat2)
2445 {
2446 PushIllegalParameter();
2447 return;
2448 }
2449 double fT, fF;
2450 SCSIZE nC1, nC2;
2451 SCSIZE nR1, nR2;
2452 SCSIZE i, j;
2453 pMat1->GetDimensions(nC1, nR1);
2454 pMat2->GetDimensions(nC2, nR2);
2455 if (fTyp == 1.0)
2456 {
2457 if (nC1 != nC2 || nR1 != nR2)
2458 {
2459 PushIllegalArgument();
2460 return;
2461 }
2462 double fCount = 0.0;
2463 double fSum1 = 0.0;
2464 double fSum2 = 0.0;
2465 double fSumSqrD = 0.0;
2466 double fVal1, fVal2;
2467 for (i = 0; i < nC1; i++)
2468 for (j = 0; j < nR1; j++)
2469 {
2470 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
2471 {
2472 fVal1 = pMat1->GetDouble(i,j);
2473 fVal2 = pMat2->GetDouble(i,j);
2474 fSum1 += fVal1;
2475 fSum2 += fVal2;
2476 fSumSqrD += (fVal1 - fVal2)*(fVal1 - fVal2);
2477 fCount++;
2478 }
2479 }
2480 if (fCount < 1.0)
2481 {
2482 PushNoValue();
2483 return;
2484 }
2485 fT = sqrt(fCount-1.0) * fabs(fSum1 - fSum2) /
2486 sqrt(fCount * fSumSqrD - (fSum1-fSum2)*(fSum1-fSum2));
2487 fF = fCount - 1.0;
2488 }
2489 else if (fTyp == 2.0)
2490 {
2491 CalculateTest(FALSE,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF);
2492 }
2493 else if (fTyp == 3.0)
2494 {
2495 CalculateTest(TRUE,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF);
2496 }
2497
2498 else
2499 {
2500 PushIllegalArgument();
2501 return;
2502 }
2503 if (fAnz == 1.0)
2504 PushDouble(GetTDist(fT, fF));
2505 else
2506 PushDouble(2.0*GetTDist(fT, fF));
2507 }
2508
2509 void ScInterpreter::ScFTest()
2510 {
2511 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFTest" );
2512 if ( !MustHaveParamCount( GetByte(), 2 ) )
2513 return;
2514 ScMatrixRef pMat2 = GetMatrix();
2515 ScMatrixRef pMat1 = GetMatrix();
2516 if (!pMat1 || !pMat2)
2517 {
2518 PushIllegalParameter();
2519 return;
2520 }
2521 SCSIZE nC1, nC2;
2522 SCSIZE nR1, nR2;
2523 SCSIZE i, j;
2524 pMat1->GetDimensions(nC1, nR1);
2525 pMat2->GetDimensions(nC2, nR2);
2526 double fCount1 = 0.0;
2527 double fCount2 = 0.0;
2528 double fSum1 = 0.0;
2529 double fSumSqr1 = 0.0;
2530 double fSum2 = 0.0;
2531 double fSumSqr2 = 0.0;
2532 double fVal;
2533 for (i = 0; i < nC1; i++)
2534 for (j = 0; j < nR1; j++)
2535 {
2536 if (!pMat1->IsString(i,j))
2537 {
2538 fVal = pMat1->GetDouble(i,j);
2539 fSum1 += fVal;
2540 fSumSqr1 += fVal * fVal;
2541 fCount1++;
2542 }
2543 }
2544 for (i = 0; i < nC2; i++)
2545 for (j = 0; j < nR2; j++)
2546 {
2547 if (!pMat2->IsString(i,j))
2548 {
2549 fVal = pMat2->GetDouble(i,j);
2550 fSum2 += fVal;
2551 fSumSqr2 += fVal * fVal;
2552 fCount2++;
2553 }
2554 }
2555 if (fCount1 < 2.0 || fCount2 < 2.0)
2556 {
2557 PushNoValue();
2558 return;
2559 }
2560 double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0);
2561 double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0);
2562 if (fS1 == 0.0 || fS2 == 0.0)
2563 {
2564 PushNoValue();
2565 return;
2566 }
2567 double fF, fF1, fF2;
2568 if (fS1 > fS2)
2569 {
2570 fF = fS1/fS2;
2571 fF1 = fCount1-1.0;
2572 fF2 = fCount2-1.0;
2573 }
2574 else
2575 {
2576 fF = fS2/fS1;
2577 fF1 = fCount2-1.0;
2578 fF2 = fCount1-1.0;
2579 }
2580 PushDouble(2.0*GetFDist(fF, fF1, fF2));
2581 /*
2582 double Z = (pow(fF,1.0/3.0)*(1.0-2.0/(9.0*fF2)) - (1.0-2.0/(9.0*fF1))) /
2583 sqrt(2.0/(9.0*fF1) + pow(fF,2.0/3.0)*2.0/(9.0*fF2));
2584 PushDouble(1.0-2.0*gauss(Z));
2585 */
2586 }
2587
2588 void ScInterpreter::ScChiTest()
2589 {
2590 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScChiTest" );
2591 if ( !MustHaveParamCount( GetByte(), 2 ) )
2592 return;
2593 ScMatrixRef pMat2 = GetMatrix();
2594 ScMatrixRef pMat1 = GetMatrix();
2595 if (!pMat1 || !pMat2)
2596 {
2597 PushIllegalParameter();
2598 return;
2599 }
2600 SCSIZE nC1, nC2;
2601 SCSIZE nR1, nR2;
2602 pMat1->GetDimensions(nC1, nR1);
2603 pMat2->GetDimensions(nC2, nR2);
2604 if (nR1 != nR2 || nC1 != nC2)
2605 {
2606 PushIllegalArgument();
2607 return;
2608 }
2609 double fChi = 0.0;
2610 for (SCSIZE i = 0; i < nC1; i++)
2611 {
2612 for (SCSIZE j = 0; j < nR1; j++)
2613 {
2614 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
2615 {
2616 double fValX = pMat1->GetDouble(i,j);
2617 double fValE = pMat2->GetDouble(i,j);
2618 fChi += (fValX - fValE) * (fValX - fValE) / fValE;
2619 }
2620 else
2621 {
2622 PushIllegalArgument();
2623 return;
2624 }
2625 }
2626 }
2627 double fDF;
2628 if (nC1 == 1 || nR1 == 1)
2629 {
2630 fDF = (double)(nC1*nR1 - 1);
2631 if (fDF == 0.0)
2632 {
2633 PushNoValue();
2634 return;
2635 }
2636 }
2637 else
2638 fDF = (double)(nC1-1)*(double)(nR1-1);
2639 PushDouble(GetChiDist(fChi, fDF));
2640 /*
2641 double fX, fS, fT, fG;
2642 fX = 1.0;
2643 for (double fi = fDF; fi >= 2.0; fi -= 2.0)
2644 fX *= fChi/fi;
2645 fX *= exp(-fChi/2.0);
2646 if (fmod(fDF, 2.0) != 0.0)
2647 fX *= sqrt(2.0*fChi/F_PI);
2648 fS = 1.0;
2649 fT = 1.0;
2650 fG = fDF;
2651 while (fT >= 1.0E-7)
2652 {
2653 fG += 2.0;
2654 fT *= fChi/fG;
2655 fS += fT;
2656 }
2657 PushDouble(1.0 - fX*fS);
2658 */
2659 }
2660
2661 void ScInterpreter::ScKurt()
2662 {
2663 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScKurt" );
2664 double fSum,fCount,vSum;
2665 std::vector<double> values;
2666 if ( !CalculateSkew(fSum,fCount,vSum,values) )
2667 return;
2668
2669 if (fCount == 0.0)
2670 {
2671 PushError( errDivisionByZero);
2672 return;
2673 }
2674
2675 double fMean = fSum / fCount;
2676
2677 for (size_t i = 0; i < values.size(); i++)
2678 vSum += (values[i] - fMean) * (values[i] - fMean);
2679
2680 double fStdDev = sqrt(vSum / (fCount - 1.0));
2681 double dx = 0.0;
2682 double xpower4 = 0.0;
2683
2684 if (fStdDev == 0.0)
2685 {
2686 PushError( errDivisionByZero);
2687 return;
2688 }
2689
2690 for (size_t i = 0; i < values.size(); i++)
2691 {
2692 dx = (values[i] - fMean) / fStdDev;
2693 xpower4 = xpower4 + (dx * dx * dx * dx);
2694 }
2695
2696 double k_d = (fCount - 2.0) * (fCount - 3.0);
2697 double k_l = fCount * (fCount + 1.0) / ((fCount - 1.0) * k_d);
2698 double k_t = 3.0 * (fCount - 1.0) * (fCount - 1.0) / k_d;
2699
2700 PushDouble(xpower4 * k_l - k_t);
2701 }
2702
2703 void ScInterpreter::ScHarMean()
2704 {
2705 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScHarMean" );
2706 short nParamCount = GetByte();
2707 double nVal = 0.0;
2708 double nValCount = 0.0;
2709 ScAddress aAdr;
2710 ScRange aRange;
2711 size_t nRefInList = 0;
2712 while ((nGlobalError == 0) && (nParamCount-- > 0))
2713 {
2714 switch (GetStackType())
2715 {
2716 case formula::svDouble :
2717 {
2718 double x = GetDouble();
2719 if (x > 0.0)
2720 {
2721 nVal += 1.0/x;
2722 nValCount++;
2723 }
2724 else
2725 SetError( errIllegalArgument);
2726 break;
2727 }
2728 case svSingleRef :
2729 {
2730 PopSingleRef( aAdr );
2731 ScBaseCell* pCell = GetCell( aAdr );
2732 if (HasCellValueData(pCell))
2733 {
2734 double x = GetCellValue( aAdr, pCell );
2735 if (x > 0.0)
2736 {
2737 nVal += 1.0/x;
2738 nValCount++;
2739 }
2740 else
2741 SetError( errIllegalArgument);
2742 }
2743 break;
2744 }
2745 case formula::svDoubleRef :
2746 case svRefList :
2747 {
2748 USHORT nErr = 0;
2749 PopDoubleRef( aRange, nParamCount, nRefInList);
2750 double nCellVal;
2751 ScValueIterator aValIter(pDok, aRange, glSubTotal);
2752 if (aValIter.GetFirst(nCellVal, nErr))
2753 {
2754 if (nCellVal > 0.0)
2755 {
2756 nVal += 1.0/nCellVal;
2757 nValCount++;
2758 }
2759 else
2760 SetError( errIllegalArgument);
2761 SetError(nErr);
2762 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
2763 {
2764 if (nCellVal > 0.0)
2765 {
2766 nVal += 1.0/nCellVal;
2767 nValCount++;
2768 }
2769 else
2770 SetError( errIllegalArgument);
2771 }
2772 SetError(nErr);
2773 }
2774 }
2775 break;
2776 case svMatrix :
2777 {
2778 ScMatrixRef pMat = PopMatrix();
2779 if (pMat)
2780 {
2781 SCSIZE nCount = pMat->GetElementCount();
2782 if (pMat->IsNumeric())
2783 {
2784 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
2785 {
2786 double x = pMat->GetDouble(nElem);
2787 if (x > 0.0)
2788 {
2789 nVal += 1.0/x;
2790 nValCount++;
2791 }
2792 else
2793 SetError( errIllegalArgument);
2794 }
2795 }
2796 else
2797 {
2798 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
2799 if (!pMat->IsString(nElem))
2800 {
2801 double x = pMat->GetDouble(nElem);
2802 if (x > 0.0)
2803 {
2804 nVal += 1.0/x;
2805 nValCount++;
2806 }
2807 else
2808 SetError( errIllegalArgument);
2809 }
2810 }
2811 }
2812 }
2813 break;
2814 default : SetError(errIllegalParameter); break;
2815 }
2816 }
2817 if (nGlobalError == 0)
2818 PushDouble((double)nValCount/nVal);
2819 else
2820 PushError( nGlobalError);
2821 }
2822
2823 void ScInterpreter::ScGeoMean()
2824 {
2825 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGeoMean" );
2826 short nParamCount = GetByte();
2827 double nVal = 0.0;
2828 double nValCount = 0.0;
2829 ScAddress aAdr;
2830 ScRange aRange;
2831
2832 size_t nRefInList = 0;
2833 while ((nGlobalError == 0) && (nParamCount-- > 0))
2834 {
2835 switch (GetStackType())
2836 {
2837 case formula::svDouble :
2838 {
2839 double x = GetDouble();
2840 if (x > 0.0)
2841 {
2842 nVal += log(x);
2843 nValCount++;
2844 }
2845 else
2846 SetError( errIllegalArgument);
2847 break;
2848 }
2849 case svSingleRef :
2850 {
2851 PopSingleRef( aAdr );
2852 ScBaseCell* pCell = GetCell( aAdr );
2853 if (HasCellValueData(pCell))
2854 {
2855 double x = GetCellValue( aAdr, pCell );
2856 if (x > 0.0)
2857 {
2858 nVal += log(x);
2859 nValCount++;
2860 }
2861 else
2862 SetError( errIllegalArgument);
2863 }
2864 break;
2865 }
2866 case formula::svDoubleRef :
2867 case svRefList :
2868 {
2869 USHORT nErr = 0;
2870 PopDoubleRef( aRange, nParamCount, nRefInList);
2871 double nCellVal;
2872 ScValueIterator aValIter(pDok, aRange, glSubTotal);
2873 if (aValIter.GetFirst(nCellVal, nErr))
2874 {
2875 if (nCellVal > 0.0)
2876 {
2877 nVal += log(nCellVal);
2878 nValCount++;
2879 }
2880 else
2881 SetError( errIllegalArgument);
2882 SetError(nErr);
2883 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
2884 {
2885 if (nCellVal > 0.0)
2886 {
2887 nVal += log(nCellVal);
2888 nValCount++;
2889 }
2890 else
2891 SetError( errIllegalArgument);
2892 }
2893 SetError(nErr);
2894 }
2895 }
2896 break;
2897 case svMatrix :
2898 {
2899 ScMatrixRef pMat = PopMatrix();
2900 if (pMat)
2901 {
2902 SCSIZE nCount = pMat->GetElementCount();
2903 if (pMat->IsNumeric())
2904 {
2905 for (SCSIZE ui = 0; ui < nCount; ui++)
2906 {
2907 double x = pMat->GetDouble(ui);
2908 if (x > 0.0)
2909 {
2910 nVal += log(x);
2911 nValCount++;
2912 }
2913 else
2914 SetError( errIllegalArgument);
2915 }
2916 }
2917 else
2918 {
2919 for (SCSIZE ui = 0; ui < nCount; ui++)
2920 if (!pMat->IsString(ui))
2921 {
2922 double x = pMat->GetDouble(ui);
2923 if (x > 0.0)
2924 {
2925 nVal += log(x);
2926 nValCount++;
2927 }
2928 else
2929 SetError( errIllegalArgument);
2930 }
2931 }
2932 }
2933 }
2934 break;
2935 default : SetError(errIllegalParameter); break;
2936 }
2937 }
2938 if (nGlobalError == 0)
2939 PushDouble(exp(nVal / nValCount));
2940 else
2941 PushError( nGlobalError);
2942 }
2943
2944 void ScInterpreter::ScStandard()
2945 {
2946 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScStandard" );
2947 if ( MustHaveParamCount( GetByte(), 3 ) )
2948 {
2949 double sigma = GetDouble();
2950 double mue = GetDouble();
2951 double x = GetDouble();
2952 if (sigma < 0.0)
2953 PushError( errIllegalArgument);
2954 else if (sigma == 0.0)
2955 PushError( errDivisionByZero);
2956 else
2957 PushDouble((x-mue)/sigma);
2958 }
2959 }
2960 bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values)
2961 {
2962 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateSkew" );
2963 short nParamCount = GetByte();
2964 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
2965 return false;
2966
2967 fSum = 0.0;
2968 fCount = 0.0;
2969 vSum = 0.0;
2970 double fVal = 0.0;
2971 ScAddress aAdr;
2972 ScRange aRange;
2973 size_t nRefInList = 0;
2974 while (nParamCount-- > 0)
2975 {
2976 switch (GetStackType())
2977 {
2978 case formula::svDouble :
2979 {
2980 fVal = GetDouble();
2981 fSum += fVal;
2982 values.push_back(fVal);
2983 fCount++;
2984 }
2985 break;
2986 case svSingleRef :
2987 {
2988 PopSingleRef( aAdr );
2989 ScBaseCell* pCell = GetCell( aAdr );
2990 if (HasCellValueData(pCell))
2991 {
2992 fVal = GetCellValue( aAdr, pCell );
2993 fSum += fVal;
2994 values.push_back(fVal);
2995 fCount++;
2996 }
2997 }
2998 break;
2999 case formula::svDoubleRef :
3000 case svRefList :
3001 {
3002 PopDoubleRef( aRange, nParamCount, nRefInList);
3003 USHORT nErr = 0;
3004 ScValueIterator aValIter(pDok, aRange);
3005 if (aValIter.GetFirst(fVal, nErr))
3006 {
3007 fSum += fVal;
3008 values.push_back(fVal);
3009 fCount++;
3010 SetError(nErr);
3011 while ((nErr == 0) && aValIter.GetNext(fVal, nErr))
3012 {
3013 fSum += fVal;
3014 values.push_back(fVal);
3015 fCount++;
3016 }
3017 SetError(nErr);
3018 }
3019 }
3020 break;
3021 case svMatrix :
3022 {
3023 ScMatrixRef pMat = PopMatrix();
3024 if (pMat)
3025 {
3026 SCSIZE nCount = pMat->GetElementCount();
3027 if (pMat->IsNumeric())
3028 {
3029 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3030 {
3031 fVal = pMat->GetDouble(nElem);
3032 fSum += fVal;
3033 values.push_back(fVal);
3034 fCount++;
3035 }
3036 }
3037 else
3038 {
3039 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3040 if (!pMat->IsString(nElem))
3041 {
3042 fVal = pMat->GetDouble(nElem);
3043 fSum += fVal;
3044 values.push_back(fVal);
3045 fCount++;
3046 }
3047 }
3048 }
3049 }
3050 break;
3051 default :
3052 SetError(errIllegalParameter);
3053 break;
3054 }
3055 }
3056
3057 if (nGlobalError)
3058 {
3059 PushError( nGlobalError);
3060 return false;
3061 } // if (nGlobalError)
3062 return true;
3063 }
3064
3065 void ScInterpreter::ScSkew()
3066 {
3067 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSkew" );
3068 double fSum,fCount,vSum;
3069 std::vector<double> values;
3070 if ( !CalculateSkew(fSum,fCount,vSum,values) )
3071 return;
3072
3073 double fMean = fSum / fCount;
3074
3075 for (size_t i = 0; i < values.size(); i++)
3076 vSum += (values[i] - fMean) * (values[i] - fMean);
3077
3078 double fStdDev = sqrt(vSum / (fCount - 1.0));
3079 double dx = 0.0;
3080 double xcube = 0.0;
3081
3082 if (fStdDev == 0)
3083 {
3084 PushIllegalArgument();
3085 return;
3086 }
3087
3088 for (size_t i = 0; i < values.size(); i++)
3089 {
3090 dx = (values[i] - fMean) / fStdDev;
3091 xcube = xcube + (dx * dx * dx);
3092 }
3093
3094 PushDouble(((xcube * fCount) / (fCount - 1.0)) / (fCount - 2.0));
3095 }
3096
3097 double ScInterpreter::GetMedian( vector<double> & rArray )
3098 {
3099 size_t nSize = rArray.size();
3100 if (rArray.empty() || nSize == 0 || nGlobalError)
3101 {
3102 SetError( errNoValue);
3103 return 0.0;
3104 }
3105
3106 // Upper median.
3107 size_t nMid = nSize / 2;
3108 vector<double>::iterator iMid = rArray.begin() + nMid;
3109 ::std::nth_element( rArray.begin(), iMid, rArray.end());
3110 if (nSize & 1)
3111 return *iMid; // Lower and upper median are equal.
3112 else
3113 {
3114 double fUp = *iMid;
3115 // Lower median.
3116 iMid = rArray.begin() + nMid - 1;
3117 ::std::nth_element( rArray.begin(), iMid, rArray.end());
3118 return (fUp + *iMid) / 2;
3119 }
3120 }
3121
3122 void ScInterpreter::ScMedian()
3123 {
3124 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScMedian" );
3125 BYTE nParamCount = GetByte();
3126 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3127 return;
3128 vector<double> aArray;
3129 GetNumberSequenceArray( nParamCount, aArray);
3130 PushDouble( GetMedian( aArray));
3131 }
3132
3133 double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile )
3134 {
3135 size_t nSize = rArray.size();
3136 if (rArray.empty() || nSize == 0 || nGlobalError)
3137 {
3138 SetError( errNoValue);
3139 return 0.0;
3140 }
3141
3142 if (nSize == 1)
3143 return rArray[0];
3144 else
3145 {
3146 size_t nIndex = (size_t)::rtl::math::approxFloor( fPercentile * (nSize-1));
3147 double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1));
3148 DBG_ASSERT(nIndex < nSize, "GetPercentile: wrong index(1)");
3149 vector<double>::iterator iter = rArray.begin() + nIndex;
3150 ::std::nth_element( rArray.begin(), iter, rArray.end());
3151 if (fDiff == 0.0)
3152 return *iter;
3153 else
3154 {
3155 DBG_ASSERT(nIndex < nSize-1, "GetPercentile: wrong index(2)");
3156 double fVal = *iter;
3157 iter = rArray.begin() + nIndex+1;
3158 ::std::nth_element( rArray.begin(), iter, rArray.end());
3159 return fVal + fDiff * (*iter - fVal);
3160 }
3161 }
3162 }
3163
3164 void ScInterpreter::ScPercentile()
3165 {
3166 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPercentile" );
3167 if ( !MustHaveParamCount( GetByte(), 2 ) )
3168 return;
3169 double alpha = GetDouble();
3170 if (alpha < 0.0 || alpha > 1.0)
3171 {
3172 PushIllegalArgument();
3173 return;
3174 }
3175 vector<double> aArray;
3176 GetNumberSequenceArray( 1, aArray);
3177 PushDouble( GetPercentile( aArray, alpha));
3178 }
3179
3180 void ScInterpreter::ScQuartile()
3181 {
3182 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScQuartile" );
3183 if ( !MustHaveParamCount( GetByte(), 2 ) )
3184 return;
3185 double fFlag = ::rtl::math::approxFloor(GetDouble());
3186 if (fFlag < 0.0 || fFlag > 4.0)
3187 {
3188 PushIllegalArgument();
3189 return;
3190 }
3191 vector<double> aArray;
3192 GetNumberSequenceArray( 1, aArray);
3193 PushDouble( fFlag == 2.0 ? GetMedian( aArray) : GetPercentile( aArray, 0.25 * fFlag));
3194 }
3195
3196 void ScInterpreter::ScModalValue()
3197 {
3198 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScModalValue" );
3199 BYTE nParamCount = GetByte();
3200 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3201 return;
3202 vector<double> aSortArray;
3203 GetSortArray(nParamCount, aSortArray);
3204 SCSIZE nSize = aSortArray.size();
3205 if (aSortArray.empty() || nSize == 0 || nGlobalError)
3206 PushNoValue();
3207 else
3208 {
3209 SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
3210 double nOldVal = aSortArray[0];
3211 SCSIZE i;
3212
3213 for ( i = 1; i < nSize; i++)
3214 {
3215 if (aSortArray[i] == nOldVal)
3216 nCount++;
3217 else
3218 {
3219 nOldVal = aSortArray[i];
3220 if (nCount > nMax)
3221 {
3222 nMax = nCount;
3223 nMaxIndex = i-1;
3224 }
3225 nCount = 1;
3226 }
3227 }
3228 if (nCount > nMax)
3229 {
3230 nMax = nCount;
3231 nMaxIndex = i-1;
3232 }
3233 if (nMax == 1 && nCount == 1)
3234 PushNoValue();
3235 else if (nMax == 1)
3236 PushDouble(nOldVal);
3237 else
3238 PushDouble(aSortArray[nMaxIndex]);
3239 }
3240 }
3241
3242 void ScInterpreter::CalculateSmallLarge(BOOL bSmall)
3243 {
3244 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateSmallLarge" );
3245 if ( !MustHaveParamCount( GetByte(), 2 ) )
3246 return;
3247 double f = ::rtl::math::approxFloor(GetDouble());
3248 if (f < 1.0)
3249 {
3250 PushIllegalArgument();
3251 return;
3252 }
3253 SCSIZE k = static_cast<SCSIZE>(f);
3254 vector<double> aSortArray;
3255 /* TODO: using nth_element() is best for one single value, but LARGE/SMALL
3256 * actually are defined to return an array of values if an array of
3257 * positions was passed, in which case, depending on the number of values,
3258 * we may or will need a real sorted array again, see #i32345. */
3259 //GetSortArray(1, aSortArray);
3260 GetNumberSequenceArray(1, aSortArray);
3261 SCSIZE nSize = aSortArray.size();
3262 if (aSortArray.empty() || nSize == 0 || nGlobalError || nSize < k)
3263 PushNoValue();
3264 else
3265 {
3266 // TODO: the sorted case for array: PushDouble( aSortArray[ bSmall ? k-1 : nSize-k ] );
3267 vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
3268 ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
3269 PushDouble( *iPos);
3270 }
3271 }
3272
3273 void ScInterpreter::ScLarge()
3274 {
3275 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScLarge" );
3276 CalculateSmallLarge(FALSE);
3277 }
3278
3279 void ScInterpreter::ScSmall()
3280 {
3281 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSmall" );
3282 CalculateSmallLarge(TRUE);
3283 }
3284
3285 void ScInterpreter::ScPercentrank()
3286 {
3287 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPercentrank" );
3288 BYTE nParamCount = GetByte();
3289 if ( !MustHaveParamCount( nParamCount, 2 ) )
3290 return;
3291 #if 0
3292 /* wird nicht unterstuetzt
3293 double fPrec;
3294 if (nParamCount == 3)
3295 {
3296 fPrec = ::rtl::math::approxFloor(GetDouble());
3297 if (fPrec < 1.0)
3298 {
3299 PushIllegalArgument();
3300 return;
3301 }
3302 }
3303 else
3304 fPrec = 3.0;
3305 */
3306 #endif
3307 double fNum = GetDouble();
3308 vector<double> aSortArray;
3309 GetSortArray(1, aSortArray);
3310 SCSIZE nSize = aSortArray.size();
3311 if (aSortArray.empty() || nSize == 0 || nGlobalError)
3312 PushNoValue();
3313 else
3314 {
3315 if (fNum < aSortArray[0] || fNum > aSortArray[nSize-1])
3316 PushNoValue();
3317 else if ( nSize == 1 )
3318 PushDouble(1.0); // fNum == pSortArray[0], see test above
3319 else
3320 {
3321 double fRes;
3322 SCSIZE nOldCount = 0;
3323 double fOldVal = aSortArray[0];
3324 SCSIZE i;
3325 for (i = 1; i < nSize && aSortArray[i] < fNum; i++)
3326 {
3327 if (aSortArray[i] != fOldVal)
3328 {
3329 nOldCount = i;
3330 fOldVal = aSortArray[i];
3331 }
3332 }
3333 if (aSortArray[i] != fOldVal)
3334 nOldCount = i;
3335 if (fNum == aSortArray[i])
3336 fRes = (double)nOldCount/(double)(nSize-1);
3337 else
3338 {
3339 // #75312# nOldCount is the count of smaller entries
3340 // fNum is between pSortArray[nOldCount-1] and pSortArray[nOldCount]
3341 // use linear interpolation to find a position between the entries
3342
3343 if ( nOldCount == 0 )
3344 {
3345 DBG_ERROR("should not happen");
3346 fRes = 0.0;
3347 }
3348 else
3349 {
3350 double fFract = ( fNum - aSortArray[nOldCount-1] ) /
3351 ( aSortArray[nOldCount] - aSortArray[nOldCount-1] );
3352 fRes = ( (double)(nOldCount-1)+fFract )/(double)(nSize-1);
3353 }
3354 }
3355 PushDouble(fRes);
3356 }
3357 }
3358 }
3359
3360 void ScInterpreter::ScTrimMean()
3361 {
3362 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTrimMean" );
3363 if ( !MustHaveParamCount( GetByte(), 2 ) )
3364 return;
3365 double alpha = GetDouble();
3366 if (alpha < 0.0 || alpha >= 1.0)
3367 {
3368 PushIllegalArgument();
3369 return;
3370 }
3371 vector<double> aSortArray;
3372 GetSortArray(1, aSortArray);
3373 SCSIZE nSize = aSortArray.size();
3374 if (aSortArray.empty() || nSize == 0 || nGlobalError)
3375 PushNoValue();
3376 else
3377 {
3378 ULONG nIndex = (ULONG) ::rtl::math::approxFloor(alpha*(double)nSize);
3379 if (nIndex % 2 != 0)
3380 nIndex--;
3381 nIndex /= 2;
3382 DBG_ASSERT(nIndex < nSize, "ScTrimMean: falscher Index");
3383 double fSum = 0.0;
3384 for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
3385 fSum += aSortArray[i];
3386 PushDouble(fSum/(double)(nSize-2*nIndex));
3387 }
3388 }
3389
3390 void ScInterpreter::GetNumberSequenceArray( BYTE nParamCount, vector<double>& rArray )
3391 {
3392 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetSortArray" );
3393 ScAddress aAdr;
3394 ScRange aRange;
3395 short nParam = nParamCount;
3396 size_t nRefInList = 0;
3397 while (nParam-- > 0)
3398 {
3399 switch (GetStackType())
3400 {
3401 case formula::svDouble :
3402 rArray.push_back( PopDouble());
3403 break;
3404 case svSingleRef :
3405 {
3406 PopSingleRef( aAdr );
3407 ScBaseCell* pCell = GetCell( aAdr );
3408 if (HasCellValueData(pCell))
3409 rArray.push_back( GetCellValue( aAdr, pCell));
3410 }
3411 break;
3412 case formula::svDoubleRef :
3413 case svRefList :
3414 {
3415 PopDoubleRef( aRange, nParam, nRefInList);
3416 if (nGlobalError)
3417 break;
3418
3419 aRange.Justify();
3420 SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3421 nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3422 rArray.reserve( rArray.size() + nCellCount);
3423
3424 USHORT nErr = 0;
3425 double fCellVal;
3426 ScValueIterator aValIter(pDok, aRange);
3427 if (aValIter.GetFirst( fCellVal, nErr))
3428 {
3429 rArray.push_back( fCellVal);
3430 SetError(nErr);
3431 while ((nErr == 0) && aValIter.GetNext( fCellVal, nErr))
3432 rArray.push_back( fCellVal);
3433 SetError(nErr);
3434 }
3435 }
3436 break;
3437 case svMatrix :
3438 {
3439 ScMatrixRef pMat = PopMatrix();
3440 if (!pMat)
3441 break;
3442
3443 SCSIZE nCount = pMat->GetElementCount();
3444 rArray.reserve( rArray.size() + nCount);
3445 if (pMat->IsNumeric())
3446 {
3447 for (SCSIZE i = 0; i < nCount; ++i)
3448 rArray.push_back( pMat->GetDouble(i));
3449 }
3450 else
3451 {
3452 for (SCSIZE i = 0; i < nCount; ++i)
3453 if (!pMat->IsString(i))
3454 rArray.push_back( pMat->GetDouble(i));
3455 }
3456 }
3457 break;
3458 default :
3459 PopError();
3460 SetError( errIllegalParameter);
3461 break;
3462 }
3463 if (nGlobalError)
3464 break; // while
3465 }
3466 // nParam > 0 in case of error, clean stack environment and obtain earlier
3467 // error if there was one.
3468 while (nParam-- > 0)
3469 PopError();
3470 }
3471
3472 void ScInterpreter::GetSortArray( BYTE nParamCount, vector<double>& rSortArray, vector<long>* pIndexOrder )
3473 {
3474 GetNumberSequenceArray( nParamCount, rSortArray);
3475
3476 if (rSortArray.size() > MAX_ANZ_DOUBLE_FOR_SORT)
3477 SetError( errStackOverflow);
3478 else if (rSortArray.empty())
3479 SetError( errNoValue);
3480
3481 if (nGlobalError == 0)
3482 QuickSort( rSortArray, pIndexOrder);
3483 }
3484
3485 static void lcl_QuickSort( long nLo, long nHi, vector<double>& rSortArray, vector<long>* pIndexOrder )
3486 {
3487 // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().
3488
3489 using ::std::swap;
3490
3491 if (nHi - nLo == 1)
3492 {
3493 if (rSortArray[nLo] > rSortArray[nHi])
3494 {
3495 swap(rSortArray[nLo], rSortArray[nHi]);
3496 if (pIndexOrder)
3497 swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
3498 }
3499 return;
3500 }
3501
3502 long ni = nLo;
3503 long nj = nHi;
3504 do
3505 {
3506 double fLo = rSortArray[nLo];
3507 while (ni <= nHi && rSortArray[ni] < fLo) ni++;
3508 while (nj >= nLo && fLo < rSortArray[nj]) nj--;
3509 if (ni <= nj)
3510 {
3511 if (ni != nj)
3512 {
3513 swap(rSortArray[ni], rSortArray[nj]);
3514 if (pIndexOrder)
3515 swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
3516 }
3517
3518 ++ni;
3519 --nj;
3520 }
3521 }
3522 while (ni < nj);
3523
3524 if ((nj - nLo) < (nHi - ni))
3525 {
3526 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
3527 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
3528 }
3529 else
3530 {
3531 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
3532 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
3533 }
3534 }
3535
3536 void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<long>* pIndexOrder )
3537 {
3538 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::QuickSort" );
3539 long n = static_cast<long>(rSortArray.size());
3540
3541 if (pIndexOrder)
3542 {
3543 pIndexOrder->clear();
3544 pIndexOrder->reserve(n);
3545 for (long i = 0; i < n; ++i)
3546 pIndexOrder->push_back(i);
3547 }
3548
3549 if (n < 2)
3550 return;
3551
3552 size_t nValCount = rSortArray.size();
3553 for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
3554 {
3555 size_t nInd = rand() % (int) (nValCount-1);
3556 ::std::swap( rSortArray[i], rSortArray[nInd]);
3557 if (pIndexOrder)
3558 ::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
3559 }
3560
3561 lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
3562 }
3563
3564 void ScInterpreter::ScRank()
3565 {
3566 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScRank" );
3567 BYTE nParamCount = GetByte();
3568 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
3569 return;
3570 BOOL bDescending;
3571 if (nParamCount == 3)
3572 bDescending = GetBool();
3573 else
3574 bDescending = FALSE;
3575 double fCount = 1.0;
3576 BOOL bValid = FALSE;
3577 switch (GetStackType())
3578 {
3579 case formula::svDouble :
3580 {
3581 double x = GetDouble();
3582 double fVal = GetDouble();
3583 if (x == fVal)
3584 bValid = TRUE;
3585 break;
3586 }
3587 case svSingleRef :
3588 {
3589 ScAddress aAdr;
3590 PopSingleRef( aAdr );
3591 double fVal = GetDouble();
3592 ScBaseCell* pCell = GetCell( aAdr );
3593 if (HasCellValueData(pCell))
3594 {
3595 double x = GetCellValue( aAdr, pCell );
3596 if (x == fVal)
3597 bValid = TRUE;
3598 }
3599 break;
3600 }
3601 case formula::svDoubleRef :
3602 case svRefList :
3603 {
3604 ScRange aRange;
3605 short nParam = 1;
3606 size_t nRefInList = 0;
3607 while (nParam-- > 0)
3608 {
3609 USHORT nErr = 0;
3610 // Preserve stack until all RefList elements are done!
3611 USHORT nSaveSP = sp;
3612 PopDoubleRef( aRange, nParam, nRefInList);
3613 if (nParam)
3614 --sp; // simulate pop
3615 double fVal = GetDouble();
3616 if (nParam)
3617 sp = nSaveSP;
3618 double nCellVal;
3619 ScValueIterator aValIter(pDok, aRange, glSubTotal);
3620 if (aValIter.GetFirst(nCellVal, nErr))
3621 {
3622 if (nCellVal == fVal)
3623 bValid = TRUE;
3624 else if ((!bDescending && nCellVal > fVal) ||
3625 (bDescending && nCellVal < fVal) )
3626 fCount++;
3627 SetError(nErr);
3628 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
3629 {
3630 if (nCellVal == fVal)
3631 bValid = TRUE;
3632 else if ((!bDescending && nCellVal > fVal) ||
3633 (bDescending && nCellVal < fVal) )
3634 fCount++;
3635 }
3636 }
3637 SetError(nErr);
3638 }
3639 }
3640 break;
3641 case svMatrix :
3642 {
3643 ScMatrixRef pMat = PopMatrix();
3644 double fVal = GetDouble();
3645 if (pMat)
3646 {
3647 SCSIZE nCount = pMat->GetElementCount();
3648 if (pMat->IsNumeric())
3649 {
3650 for (SCSIZE i = 0; i < nCount; i++)
3651 {
3652 double x = pMat->GetDouble(i);
3653 if (x == fVal)
3654 bValid = TRUE;
3655 else if ((!bDescending && x > fVal) ||
3656 (bDescending && x < fVal) )
3657 fCount++;
3658 }
3659 }
3660 else
3661 {
3662 for (SCSIZE i = 0; i < nCount; i++)
3663 if (!pMat->IsString(i))
3664 {
3665 double x = pMat->GetDouble(i);
3666 if (x == fVal)
3667 bValid = TRUE;
3668 else if ((!bDescending && x > fVal) ||
3669 (bDescending && x < fVal) )
3670 fCount++;
3671 }
3672 }
3673 }
3674 }
3675 break;
3676 default : SetError(errIllegalParameter); break;
3677 }
3678 if (bValid)
3679 PushDouble(fCount);
3680 else
3681 PushNoValue();
3682 }
3683
3684 void ScInterpreter::ScAveDev()
3685 {
3686 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScAveDev" );
3687 BYTE nParamCount = GetByte();
3688 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3689 return;
3690 USHORT SaveSP = sp;
3691 double nMiddle = 0.0;
3692 double rVal = 0.0;
3693 double rValCount = 0.0;
3694 ScAddress aAdr;
3695 ScRange aRange;
3696 short nParam = nParamCount;
3697 size_t nRefInList = 0;
3698 while (nParam-- > 0)
3699 {
3700 switch (GetStackType())
3701 {
3702 case formula::svDouble :
3703 rVal += GetDouble();
3704 rValCount++;
3705 break;
3706 case svSingleRef :
3707 {
3708 PopSingleRef( aAdr );
3709 ScBaseCell* pCell = GetCell( aAdr );
3710 if (HasCellValueData(pCell))
3711 {
3712 rVal += GetCellValue( aAdr, pCell );
3713 rValCount++;
3714 }
3715 }
3716 break;
3717 case formula::svDoubleRef :
3718 case svRefList :
3719 {
3720 USHORT nErr = 0;
3721 double nCellVal;
3722 PopDoubleRef( aRange, nParam, nRefInList);
3723 ScValueIterator aValIter(pDok, aRange);
3724 if (aValIter.GetFirst(nCellVal, nErr))
3725 {
3726 rVal += nCellVal;
3727 rValCount++;
3728 SetError(nErr);
3729 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
3730 {
3731 rVal += nCellVal;
3732 rValCount++;
3733 }
3734 SetError(nErr);
3735 }
3736 }
3737 break;
3738 case svMatrix :
3739 {
3740 ScMatrixRef pMat = PopMatrix();
3741 if (pMat)
3742 {
3743 SCSIZE nCount = pMat->GetElementCount();
3744 if (pMat->IsNumeric())
3745 {
3746 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3747 {
3748 rVal += pMat->GetDouble(nElem);
3749 rValCount++;
3750 }
3751 }
3752 else
3753 {
3754 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3755 if (!pMat->IsString(nElem))
3756 {
3757 rVal += pMat->GetDouble(nElem);
3758 rValCount++;
3759 }
3760 }
3761 }
3762 }
3763 break;
3764 default :
3765 SetError(errIllegalParameter);
3766 break;
3767 }
3768 }
3769 if (nGlobalError)
3770 {
3771 PushError( nGlobalError);
3772 return;
3773 }
3774 nMiddle = rVal / rValCount;
3775 sp = SaveSP;
3776 rVal = 0.0;
3777 nParam = nParamCount;
3778 nRefInList = 0;
3779 while (nParam-- > 0)
3780 {
3781 switch (GetStackType())
3782 {
3783 case formula::svDouble :
3784 rVal += fabs(GetDouble() - nMiddle);
3785 break;
3786 case svSingleRef :
3787 {
3788 PopSingleRef( aAdr );
3789 ScBaseCell* pCell = GetCell( aAdr );
3790 if (HasCellValueData(pCell))
3791 rVal += fabs(GetCellValue( aAdr, pCell ) - nMiddle);
3792 }
3793 break;
3794 case formula::svDoubleRef :
3795 case svRefList :
3796 {
3797 USHORT nErr = 0;
3798 double nCellVal;
3799 PopDoubleRef( aRange, nParam, nRefInList);
3800 ScValueIterator aValIter(pDok, aRange);
3801 if (aValIter.GetFirst(nCellVal, nErr))
3802 {
3803 rVal += (fabs(nCellVal - nMiddle));
3804 while (aValIter.GetNext(nCellVal, nErr))
3805 rVal += fabs(nCellVal - nMiddle);
3806 }
3807 }
3808 break;
3809 case svMatrix :
3810 {
3811 ScMatrixRef pMat = PopMatrix();
3812 if (pMat)
3813 {
3814 SCSIZE nCount = pMat->GetElementCount();
3815 if (pMat->IsNumeric())
3816 {
3817 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3818 {
3819 rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
3820 }
3821 }
3822 else
3823 {
3824 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3825 {
3826 if (!pMat->IsString(nElem))
3827 rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
3828 }
3829 }
3830 }
3831 }
3832 break;
3833 default : SetError(errIllegalParameter); break;
3834 }
3835 }
3836 PushDouble(rVal / rValCount);
3837 }
3838
3839 void ScInterpreter::ScDevSq()
3840 {
3841 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScDevSq" );
3842 double nVal;
3843 double nValCount;
3844 GetStVarParams(nVal, nValCount);
3845 PushDouble(nVal);
3846 }
3847
3848 void ScInterpreter::ScProbability()
3849 {
3850 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScProbability" );
3851 BYTE nParamCount = GetByte();
3852 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
3853 return;
3854 double fUp, fLo;
3855 fUp = GetDouble();
3856 if (nParamCount == 4)
3857 fLo = GetDouble();
3858 else
3859 fLo = fUp;
3860 if (fLo > fUp)
3861 {
3862 double fTemp = fLo;
3863 fLo = fUp;
3864 fUp = fTemp;
3865 }
3866 ScMatrixRef pMatP = GetMatrix();
3867 ScMatrixRef pMatW = GetMatrix();
3868 if (!pMatP || !pMatW)
3869 PushIllegalParameter();
3870 else
3871 {
3872 SCSIZE nC1, nC2;
3873 SCSIZE nR1, nR2;
3874 pMatP->GetDimensions(nC1, nR1);
3875 pMatW->GetDimensions(nC2, nR2);
3876 if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
3877 nC2 == 0 || nR2 == 0)
3878 PushNA();
3879 else
3880 {
3881 double fSum = 0.0;
3882 double fRes = 0.0;
3883 BOOL bStop = FALSE;
3884 double fP, fW;
3885 SCSIZE nCount1 = nC1 * nR1;
3886 for ( SCSIZE i = 0; i < nCount1 && !bStop; i++ )
3887 {
3888 if (pMatP->IsValue(i) && pMatW->IsValue(i))
3889 {
3890 fP = pMatP->GetDouble(i);
3891 fW = pMatW->GetDouble(i);
3892 if (fP < 0.0 || fP > 1.0)
3893 bStop = TRUE;
3894 else
3895 {
3896 fSum += fP;
3897 if (fW >= fLo && fW <= fUp)
3898 fRes += fP;
3899 }
3900 }
3901 else
3902 SetError( errIllegalArgument);
3903 }
3904 if (bStop || fabs(fSum -1.0) > 1.0E-7)
3905 PushNoValue();
3906 else
3907 PushDouble(fRes);
3908 }
3909 }
3910 }
3911
3912 void ScInterpreter::ScCorrel()
3913 {
3914 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScCorrel" );
3915 // This is identical to ScPearson()
3916 ScPearson();
3917 }
3918
3919 void ScInterpreter::ScCovar()
3920 {
3921 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScCovar" );
3922 CalculatePearsonCovar(FALSE,FALSE);
3923 }
3924
3925 void ScInterpreter::ScPearson()
3926 {
3927 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPearson" );
3928 CalculatePearsonCovar(TRUE,FALSE);
3929 }
3930 void ScInterpreter::CalculatePearsonCovar(BOOL _bPearson,BOOL _bStexy)
3931 {
3932 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculatePearsonCovar" );
3933 if ( !MustHaveParamCount( GetByte(), 2 ) )
3934 return;
3935 ScMatrixRef pMat1 = GetMatrix();
3936 ScMatrixRef pMat2 = GetMatrix();
3937 if (!pMat1 || !pMat2)
3938 {
3939 PushIllegalParameter();
3940 return;
3941 }
3942 SCSIZE nC1, nC2;
3943 SCSIZE nR1, nR2;
3944 pMat1->GetDimensions(nC1, nR1);
3945 pMat2->GetDimensions(nC2, nR2);
3946 if (nR1 != nR2 || nC1 != nC2)
3947 {
3948 PushIllegalArgument();
3949 return;
3950 }
3951 /* #i78250#
3952 * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
3953 * but the latter produces wrong results if the absolute values are high,
3954 * for example above 10^8
3955 */
3956 double fCount = 0.0;
3957 double fSumX = 0.0;
3958 double fSumY = 0.0;
3959 double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
3960 double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
3961 double fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2
3962 for (SCSIZE i = 0; i < nC1; i++)
3963 {
3964 for (SCSIZE j = 0; j < nR1; j++)
3965 {
3966 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
3967 {
3968 double fValX = pMat1->GetDouble(i,j);
3969 double fValY = pMat2->GetDouble(i,j);
3970 fSumX += fValX;
3971 fSumY += fValY;
3972 fCount++;
3973 }
3974 }
3975 }
3976 if (fCount < (_bStexy ? 3.0 : 1.0)) // fCount==1 is handled by checking denominator later on
3977 PushNoValue();
3978 else
3979 {
3980 const double fMeanX = fSumX / fCount;
3981 const double fMeanY = fSumY / fCount;
3982 for (SCSIZE i = 0; i < nC1; i++)
3983 {
3984 for (SCSIZE j = 0; j < nR1; j++)
3985 {
3986 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
3987 {
3988 const double fValX = pMat1->GetDouble(i,j);
3989 const double fValY = pMat2->GetDouble(i,j);
3990 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
3991 if ( _bPearson )
3992 {
3993 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
3994 fSumSqrDeltaY += (fValY - fMeanY) * (fValY - fMeanY);
3995 }
3996 }
3997 }
3998 } // for (SCSIZE i = 0; i < nC1; i++)
3999 if ( _bPearson )
4000 {
4001 if (fSumSqrDeltaX == 0.0 || ( !_bStexy && fSumSqrDeltaY == 0.0) )
4002 PushError( errDivisionByZero);
4003 else if ( _bStexy )
4004 PushDouble( sqrt( (fSumSqrDeltaY - fSumDeltaXDeltaY *
4005 fSumDeltaXDeltaY / fSumSqrDeltaX) / (fCount-2)));
4006 else
4007 PushDouble( fSumDeltaXDeltaY / sqrt( fSumSqrDeltaX * fSumSqrDeltaY));
4008 } // if ( _bPearson )
4009 else
4010 {
4011 PushDouble( fSumDeltaXDeltaY / fCount);
4012 }
4013 }
4014 }
4015
4016 void ScInterpreter::ScRSQ()
4017 {
4018 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScRSQ" );
4019 // Same as ScPearson()*ScPearson()
4020 ScPearson();
4021 if (!nGlobalError)
4022 {
4023 switch (GetStackType())
4024 {
4025 case formula::svDouble:
4026 {
4027 double fVal = PopDouble();
4028 PushDouble( fVal * fVal);
4029 }
4030 break;
4031 default:
4032 PopError();
4033 PushNoValue();
4034 }
4035 }
4036 }
4037
4038 void ScInterpreter::ScSTEXY()
4039 {
4040 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSTEXY" );
4041 CalculatePearsonCovar(TRUE,TRUE);
4042 }
4043 void ScInterpreter::CalculateSlopeIntercept(BOOL bSlope)
4044 {
4045 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateSlopeIntercept" );
4046 if ( !MustHaveParamCount( GetByte(), 2 ) )
4047 return;
4048 ScMatrixRef pMat1 = GetMatrix();
4049 ScMatrixRef pMat2 = GetMatrix();
4050 if (!pMat1 || !pMat2)
4051 {
4052 PushIllegalParameter();
4053 return;
4054 }
4055 SCSIZE nC1, nC2;
4056 SCSIZE nR1, nR2;
4057 pMat1->GetDimensions(nC1, nR1);
4058 pMat2->GetDimensions(nC2, nR2);
4059 if (nR1 != nR2 || nC1 != nC2)
4060 {
4061 PushIllegalArgument();
4062 return;
4063 }
4064 // #i78250# numerical stability improved
4065 double fCount = 0.0;
4066 double fSumX = 0.0;
4067 double fSumY = 0.0;
4068 double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4069 double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4070 for (SCSIZE i = 0; i < nC1; i++)
4071 {
4072 for (SCSIZE j = 0; j < nR1; j++)
4073 {
4074 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4075 {
4076 double fValX = pMat1->GetDouble(i,j);
4077 double fValY = pMat2->GetDouble(i,j);
4078 fSumX += fValX;
4079 fSumY += fValY;
4080 fCount++;
4081 }
4082 }
4083 }
4084 if (fCount < 1.0)
4085 PushNoValue();
4086 else
4087 {
4088 double fMeanX = fSumX / fCount;
4089 double fMeanY = fSumY / fCount;
4090 for (SCSIZE i = 0; i < nC1; i++)
4091 {
4092 for (SCSIZE j = 0; j < nR1; j++)
4093 {
4094 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4095 {
4096 double fValX = pMat1->GetDouble(i,j);
4097 double fValY = pMat2->GetDouble(i,j);
4098 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4099 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4100 }
4101 }
4102 }
4103 if (fSumSqrDeltaX == 0.0)
4104 PushError( errDivisionByZero);
4105 else
4106 {
4107 if ( bSlope )
4108 PushDouble( fSumDeltaXDeltaY / fSumSqrDeltaX);
4109 else
4110 PushDouble( fMeanY - fSumDeltaXDeltaY / fSumSqrDeltaX * fMeanX);
4111 }
4112 }
4113 }
4114
4115 void ScInterpreter::ScSlope()
4116 {
4117 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSlope" );
4118 CalculateSlopeIntercept(TRUE);
4119 }
4120
4121 void ScInterpreter::ScIntercept()
4122 {
4123 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScIntercept" );
4124 CalculateSlopeIntercept(FALSE);
4125 }
4126
4127 void ScInterpreter::ScForecast()
4128 {
4129 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScForecast" );
4130 if ( !MustHaveParamCount( GetByte(), 3 ) )
4131 return;
4132 ScMatrixRef pMat1 = GetMatrix();
4133 ScMatrixRef pMat2 = GetMatrix();
4134 if (!pMat1 || !pMat2)
4135 {
4136 PushIllegalParameter();
4137 return;
4138 }
4139 SCSIZE nC1, nC2;
4140 SCSIZE nR1, nR2;
4141 pMat1->GetDimensions(nC1, nR1);
4142 pMat2->GetDimensions(nC2, nR2);
4143 if (nR1 != nR2 || nC1 != nC2)
4144 {
4145 PushIllegalArgument();
4146 return;
4147 }
4148 double fVal = GetDouble();
4149 // #i78250# numerical stability improved
4150 double fCount = 0.0;
4151 double fSumX = 0.0;
4152 double fSumY = 0.0;
4153 double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4154 double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4155 for (SCSIZE i = 0; i < nC1; i++)
4156 {
4157 for (SCSIZE j = 0; j < nR1; j++)
4158 {
4159 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4160 {
4161 double fValX = pMat1->GetDouble(i,j);
4162 double fValY = pMat2->GetDouble(i,j);
4163 fSumX += fValX;
4164 fSumY += fValY;
4165 fCount++;
4166 }
4167 }
4168 }
4169 if (fCount < 1.0)
4170 PushNoValue();
4171 else
4172 {
4173 double fMeanX = fSumX / fCount;
4174 double fMeanY = fSumY / fCount;
4175 for (SCSIZE i = 0; i < nC1; i++)
4176 {
4177 for (SCSIZE j = 0; j < nR1; j++)
4178 {
4179 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4180 {
4181 double fValX = pMat1->GetDouble(i,j);
4182 double fValY = pMat2->GetDouble(i,j);
4183 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4184 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4185 }
4186 }
4187 }
4188 if (fSumSqrDeltaX == 0.0)
4189 PushError( errDivisionByZero);
4190 else
4191 PushDouble( fMeanY + fSumDeltaXDeltaY / fSumSqrDeltaX * (fVal - fMeanX));
4192 }
4193 }
4194
vba interoperability for OOo