?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
update dev300-m57
[ooovba.git] / sc / source / core / tool / interpr3.cxx
blobb86e266151560622ca67be13471efc38acd2e153
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 {
936 if (fB < 1.0 && fX == 1.0)
937 {
938 SetError(errIllegalArgument);
939 return HUGE_VAL;
940 }
941 else
942 return 0.0;
943 }
944
945 // normal cases; result x^(a-1)*(1-x)^(b-1)/Beta(a,b)
946 const double fLogDblMax = log( ::std::numeric_limits<double>::max());
947 const double fLogDblMin = log( ::std::numeric_limits<double>::min());
948 double fLogY = (fX < 0.1) ? ::rtl::math::log1p(-fX) : log(0.5-fX+0.5);
949 double fLogX = log(fX);
950 double fAm1 = fA-1.0;
951 double fBm1 = fB-1.0;
952 double fLogBeta = GetLogBeta(fA,fB);
953 // check whether parts over- or underflow
954 if ( fAm1 * fLogX < fLogDblMax && fAm1 * fLogX > fLogDblMin
955 && fBm1 * fLogY < fLogDblMax && fBm1* fLogY > fLogDblMin
956 && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin )
957 return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);
958 else // need logarithm;
959 // might overflow as a whole, but seldom, not worth to pre-detect it
960 return exp((fA-1.0)*fLogX + (fB-1.0)* fLogY - fLogBeta);
961 }
962
963
964 /*
965 x^a * (1-x)^b
966 I_x(a,b) = ---------------- * result of ContFrac
967 a * Beta(a,b)
968 */
969 double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)
970 { // like old version
971 double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, cfnew, cf;
972 a1 = 1.0; b1 = 1.0;
973 b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;
974 if (b2 == 0.0)
975 {
976 a2 = 0.0;
977 fnorm = 1.0;
978 cf = 1.0;
979 }
980 else
981 {
982 a2 = 1.0;
983 fnorm = 1.0/b2;
984 cf = a2*fnorm;
985 }
986 cfnew = 1.0;
987 double rm = 1.0;
988
989 const double fMaxIter = 50000.0;
990 // loop security, normal cases converge in less than 100 iterations.
991 // FIXME: You will get so much iteratons for fX near mean,
992 // I do not know a better algorithm.
993 bool bfinished = false;
994 do
995 {
996 apl2m = fA + 2.0*rm;
997 d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m);
998 d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0));
999 a1 = (a2+d2m*a1)*fnorm;
1000 b1 = (b2+d2m*b1)*fnorm;
1001 a2 = a1 + d2m1*a2*fnorm;
1002 b2 = b1 + d2m1*b2*fnorm;
1003 if (b2 != 0.0)
1004 {
1005 fnorm = 1.0/b2;
1006 cfnew = a2*fnorm;
1007 bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);
1008 }
1009 cf = cfnew;
1010 rm += 1.0;
1011 }
1012 while (rm < fMaxIter && !bfinished);
1013 return cf;
1014 }
1015
1016 // cumulative distribution function, normalized
1017 double ScInterpreter::GetBetaDist(double fXin, double fAlpha, double fBeta)
1018 {
1019 // special cases
1020 if (fXin <= 0.0) // values are valid, see spec
1021 return 0.0;
1022 if (fXin >= 1.0) // values are valid, see spec
1023 return 1.0;
1024 if (fBeta == 1.0)
1025 return pow(fXin, fAlpha);
1026 if (fAlpha == 1.0)
1027 // 1.0 - pow(1.0-fX,fBeta) is not accurate enough
1028 return -::rtl::math::expm1(fBeta * ::rtl::math::log1p(-fXin));
1029 //FIXME: need special algorithm for fX near fP for large fA,fB
1030 double fResult;
1031 // I use always continued fraction, power series are neither
1032 // faster nor more accurate.
1033 double fY = (0.5-fXin)+0.5;
1034 double flnY = ::rtl::math::log1p(-fXin);
1035 double fX = fXin;
1036 double flnX = log(fXin);
1037 double fA = fAlpha;
1038 double fB = fBeta;
1039 bool bReflect = fXin > fAlpha/(fAlpha+fBeta);
1040 if (bReflect)
1041 {
1042 fA = fBeta;
1043 fB = fAlpha;
1044 fX = fY;
1045 fY = fXin;
1046 flnX = flnY;
1047 flnY = log(fXin);
1048 }
1049 fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);
1050 fResult = fResult/fA;
1051 double fP = fA/(fA+fB);
1052 double fQ = fB/(fA+fB);
1053 double fTemp;
1054 if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental
1055 fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;
1056 else
1057 fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));
1058 fResult *= fTemp;
1059 if (bReflect)
1060 fResult = 0.5 - fResult + 0.5;
1061 if (fResult > 1.0) // ensure valid range
1062 fResult = 1.0;
1063 if (fResult < 0.0)
1064 fResult = 0.0;
1065 return fResult;
1066 }
1067
1068 void ScInterpreter::ScBetaDist()
1069 {
1070 BYTE nParamCount = GetByte();
1071 if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547#
1072 return;
1073 double fLowerBound, fUpperBound;
1074 double alpha, beta, x;
1075 bool bIsCumulative;
1076 if (nParamCount == 6)
1077 bIsCumulative = GetBool();
1078 else
1079 bIsCumulative = true;
1080 if (nParamCount >= 5)
1081 fUpperBound = GetDouble();
1082 else
1083 fUpperBound = 1.0;
1084 if (nParamCount >= 4)
1085 fLowerBound = GetDouble();
1086 else
1087 fLowerBound = 0.0;
1088 beta = GetDouble();
1089 alpha = GetDouble();
1090 x = GetDouble();
1091 double fScale = fUpperBound - fLowerBound;
1092 if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0)
1093 {
1094 PushIllegalArgument();
1095 return;
1096 }
1097 if (bIsCumulative) // cumulative distribution function
1098 {
1099 // special cases
1100 if (x < fLowerBound)
1101 {
1102 PushDouble(0.0); return; //see spec
1103 }
1104 if (x > fUpperBound)
1105 {
1106 PushDouble(1.0); return; //see spec
1107 }
1108 // normal cases
1109 x = (x-fLowerBound)/fScale; // convert to standard form
1110 PushDouble(GetBetaDist(x, alpha, beta));
1111 return;
1112 }
1113 else // probability density function
1114 {
1115 if (x < fLowerBound || x > fUpperBound)
1116 {
1117 PushDouble(0.0);
1118 return;
1119 }
1120 x = (x-fLowerBound)/fScale;
1121 PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
1122 return;
1123 }
1124 }
1125
1126 void ScInterpreter::ScPhi()
1127 {
1128 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPhi" );
1129 PushDouble(phi(GetDouble()));
1130 }
1131
1132 void ScInterpreter::ScGauss()
1133 {
1134 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGauss" );
1135 PushDouble(gauss(GetDouble()));
1136 }
1137
1138 void ScInterpreter::ScFisher()
1139 {
1140 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFisher" );
1141 double fVal = GetDouble();
1142 if (fabs(fVal) >= 1.0)
1143 PushIllegalArgument();
1144 else
1145 PushDouble( ::rtl::math::atanh( fVal));
1146 }
1147
1148 void ScInterpreter::ScFisherInv()
1149 {
1150 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFisherInv" );
1151 PushDouble( tanh( GetDouble()));
1152 }
1153
1154 void ScInterpreter::ScFact()
1155 {
1156 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFact" );
1157 double nVal = GetDouble();
1158 if (nVal < 0.0)
1159 PushIllegalArgument();
1160 else
1161 PushDouble(Fakultaet(nVal));
1162 }
1163
1164 void ScInterpreter::ScKombin()
1165 {
1166 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScKombin" );
1167 if ( MustHaveParamCount( GetByte(), 2 ) )
1168 {
1169 double k = ::rtl::math::approxFloor(GetDouble());
1170 double n = ::rtl::math::approxFloor(GetDouble());
1171 if (k < 0.0 || n < 0.0 || k > n)
1172 PushIllegalArgument();
1173 else
1174 PushDouble(BinomKoeff(n, k));
1175 }
1176 }
1177
1178 void ScInterpreter::ScKombin2()
1179 {
1180 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScKombin2" );
1181 if ( MustHaveParamCount( GetByte(), 2 ) )
1182 {
1183 double k = ::rtl::math::approxFloor(GetDouble());
1184 double n = ::rtl::math::approxFloor(GetDouble());
1185 if (k < 0.0 || n < 0.0 || k > n)
1186 PushIllegalArgument();
1187 else
1188 PushDouble(BinomKoeff(n + k - 1, k));
1189 }
1190 }
1191
1192 void ScInterpreter::ScVariationen()
1193 {
1194 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScVariationen" );
1195 if ( MustHaveParamCount( GetByte(), 2 ) )
1196 {
1197 double k = ::rtl::math::approxFloor(GetDouble());
1198 double n = ::rtl::math::approxFloor(GetDouble());
1199 if (n < 0.0 || k < 0.0 || k > n)
1200 PushIllegalArgument();
1201 else if (k == 0.0)
1202 PushInt(1); // (n! / (n - 0)!) == 1
1203 else
1204 {
1205 double nVal = n;
1206 for (ULONG i = (ULONG)k-1; i >= 1; i--)
1207 nVal *= n-(double)i;
1208 PushDouble(nVal);
1209 }
1210 }
1211 }
1212
1213 void ScInterpreter::ScVariationen2()
1214 {
1215 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScVariationen2" );
1216 if ( MustHaveParamCount( GetByte(), 2 ) )
1217 {
1218 double k = ::rtl::math::approxFloor(GetDouble());
1219 double n = ::rtl::math::approxFloor(GetDouble());
1220 if (n < 0.0 || k < 0.0 || k > n)
1221 PushIllegalArgument();
1222 else
1223 PushDouble(pow(n,k));
1224 }
1225 }
1226
1227 void ScInterpreter::ScB()
1228 {
1229 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScB" );
1230 BYTE nParamCount = GetByte();
1231 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1232 return ;
1233 if (nParamCount == 3)
1234 {
1235 double x = ::rtl::math::approxFloor(GetDouble());
1236 double p = GetDouble();
1237 double n = ::rtl::math::approxFloor(GetDouble());
1238 if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1239 PushIllegalArgument();
1240 else
1241 {
1242 double q = 1.0 - p;
1243 double fFactor = pow(q, n);
1244 if (fFactor == 0.0)
1245 {
1246 fFactor = pow(p, n);
1247 if (fFactor == 0.0)
1248 PushNoValue();
1249 else
1250 {
1251 ULONG max = (ULONG) (n - x);
1252 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1253 fFactor *= (n-i)/(i+1)*q/p;
1254 PushDouble(fFactor);
1255 }
1256 }
1257 else
1258 {
1259 ULONG max = (ULONG) x;
1260 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1261 fFactor *= (n-i)/(i+1)*p/q;
1262 PushDouble(fFactor);
1263 }
1264 }
1265 }
1266 else if (nParamCount == 4)
1267 {
1268 double xe = GetDouble();
1269 double xs = GetDouble();
1270 double p = GetDouble();
1271 double n = GetDouble();
1272 // alter Stand 300-SC
1273 // if ((xs < n) && (xe < n) && (p < 1.0))
1274 // {
1275 // double Varianz = sqrt(n * p * (1.0 - p));
1276 // xs = fabs(xs - (n * p /* / 2.0 STE */ ));
1277 // xe = fabs(xe - (n * p /* / 2.0 STE */ ));
1278 //// STE double nVal = gauss((xs + 0.5) / Varianz) + gauss((xe + 0.5) / Varianz);
1279 // double nVal = fabs(gauss(xs / Varianz) - gauss(xe / Varianz));
1280 // PushDouble(nVal);
1281 // }
1282 if (xe <= n && xs <= xe &&
1283 p < 1.0 && p > 0.0 && n >= 0.0 && xs >= 0.0 )
1284 {
1285 double q = 1.0 - p;
1286 double fFactor = pow(q, n);
1287 if (fFactor == 0.0)
1288 {
1289 fFactor = pow(p, n);
1290 if (fFactor == 0.0)
1291 PushNoValue();
1292 else
1293 {
1294 double fSum = 0.0;
1295 ULONG max;
1296 if (xe < (ULONG) n)
1297 max = (ULONG) (n-xe)-1;
1298 else
1299 max = 0;
1300 ULONG i;
1301 for (i = 0; i < max && fFactor > 0.0; i++)
1302 fFactor *= (n-i)/(i+1)*q/p;
1303 if (xs < (ULONG) n)
1304 max = (ULONG) (n-xs);
1305 else
1306 fSum = fFactor;
1307 for (; i < max && fFactor > 0.0; i++)
1308 {
1309 fFactor *= (n-i)/(i+1)*q/p;
1310 fSum += fFactor;
1311 }
1312 PushDouble(fSum);
1313 }
1314 }
1315 else
1316 {
1317 ULONG max;
1318 double fSum;
1319 if ( (ULONG) xs == 0)
1320 {
1321 fSum = fFactor;
1322 max = 0;
1323 }
1324 else
1325 {
1326 max = (ULONG) xs-1;
1327 fSum = 0.0;
1328 }
1329 ULONG i;
1330 for (i = 0; i < max && fFactor > 0.0; i++)
1331 fFactor *= (n-i)/(i+1)*p/q;
1332 if ((ULONG)xe == 0) // beide 0
1333 fSum = fFactor;
1334 else
1335 max = (ULONG) xe;
1336 for (; i < max && fFactor > 0.0; i++)
1337 {
1338 fFactor *= (n-i)/(i+1)*p/q;
1339 fSum += fFactor;
1340 }
1341 PushDouble(fSum);
1342 }
1343 }
1344 else
1345 PushIllegalArgument();
1346 }
1347 }
1348
1349 void ScInterpreter::ScBinomDist()
1350 {
1351 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScBinomDist" );
1352 if ( MustHaveParamCount( GetByte(), 4 ) )
1353 {
1354 double kum = GetDouble(); // 0 oder 1
1355 double p = GetDouble(); // p
1356 double n = ::rtl::math::approxFloor(GetDouble()); // n
1357 double x = ::rtl::math::approxFloor(GetDouble()); // x
1358 double fFactor, q, fSum;
1359 if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
1360 PushIllegalArgument();
1361 else if (kum == 0.0) // Dichte
1362 {
1363 q = 1.0 - p;
1364 fFactor = pow(q, n);
1365 if (fFactor == 0.0)
1366 {
1367 fFactor = pow(p, n);
1368 if (fFactor == 0.0)
1369 PushNoValue();
1370 else
1371 {
1372 ULONG max = (ULONG) (n - x);
1373 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1374 fFactor *= (n-i)/(i+1)*q/p;
1375 PushDouble(fFactor);
1376 }
1377 }
1378 else
1379 {
1380 ULONG max = (ULONG) x;
1381 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1382 fFactor *= (n-i)/(i+1)*p/q;
1383 PushDouble(fFactor);
1384 }
1385 }
1386 else // Verteilung
1387 {
1388 if (n == x)
1389 PushDouble(1.0);
1390 else
1391 {
1392 q = 1.0 - p;
1393 fFactor = pow(q, n);
1394 if (fFactor == 0.0)
1395 {
1396 fFactor = pow(p, n);
1397 if (fFactor == 0.0)
1398 PushNoValue();
1399 else
1400 {
1401 fSum = 1.0 - fFactor;
1402 ULONG max = (ULONG) (n - x) - 1;
1403 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1404 {
1405 fFactor *= (n-i)/(i+1)*q/p;
1406 fSum -= fFactor;
1407 }
1408 if (fSum < 0.0)
1409 PushDouble(0.0);
1410 else
1411 PushDouble(fSum);
1412 }
1413 }
1414 else
1415 {
1416 fSum = fFactor;
1417 ULONG max = (ULONG) x;
1418 for (ULONG i = 0; i < max && fFactor > 0.0; i++)
1419 {
1420 fFactor *= (n-i)/(i+1)*p/q;
1421 fSum += fFactor;
1422 }
1423 PushDouble(fSum);
1424 }
1425 }
1426 }
1427 }
1428 }
1429
1430 void ScInterpreter::ScCritBinom()
1431 {
1432 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScCritBinom" );
1433 if ( MustHaveParamCount( GetByte(), 3 ) )
1434 {
1435 double alpha = GetDouble(); // alpha
1436 double p = GetDouble(); // p
1437 double n = ::rtl::math::approxFloor(GetDouble());
1438 if (n < 0.0 || alpha <= 0.0 || alpha >= 1.0 || p < 0.0 || p > 1.0)
1439 PushIllegalArgument();
1440 else
1441 {
1442 double q = 1.0 - p;
1443 double fFactor = pow(q,n);
1444 if (fFactor == 0.0)
1445 {
1446 fFactor = pow(p, n);
1447 if (fFactor == 0.0)
1448 PushNoValue();
1449 else
1450 {
1451 double fSum = 1.0 - fFactor; ULONG max = (ULONG) n;
1452 ULONG i;
1453
1454 for ( i = 0; i < max && fSum >= alpha; i++)
1455 {
1456 fFactor *= (n-i)/(i+1)*q/p;
1457 fSum -= fFactor;
1458 }
1459 PushDouble(n-i);
1460 }
1461 }
1462 else
1463 {
1464 double fSum = fFactor; ULONG max = (ULONG) n;
1465 ULONG i;
1466
1467 for ( i = 0; i < max && fSum < alpha; i++)
1468 {
1469 fFactor *= (n-i)/(i+1)*p/q;
1470 fSum += fFactor;
1471 }
1472 PushDouble(i);
1473 }
1474 }
1475 }
1476 }
1477
1478 void ScInterpreter::ScNegBinomDist()
1479 {
1480 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScNegBinomDist" );
1481 if ( MustHaveParamCount( GetByte(), 3 ) )
1482 {
1483 double p = GetDouble(); // p
1484 double r = GetDouble(); // r
1485 double x = GetDouble(); // x
1486 if (r < 0.0 || x < 0.0 || p < 0.0 || p > 1.0)
1487 PushIllegalArgument();
1488 else
1489 {
1490 double q = 1.0 - p;
1491 double fFactor = pow(p,r);
1492 for (double i = 0.0; i < x; i++)
1493 fFactor *= (i+r)/(i+1.0)*q;
1494 PushDouble(fFactor);
1495 }
1496 }
1497 }
1498
1499 void ScInterpreter::ScNormDist()
1500 {
1501 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScNormDist" );
1502 if ( MustHaveParamCount( GetByte(), 4 ) )
1503 {
1504 double kum = GetDouble(); // 0 oder 1
1505 double sigma = GetDouble(); // Stdabw
1506 double mue = GetDouble(); // Mittelwert
1507 double x = GetDouble(); // x
1508 if (sigma < 0.0)
1509 PushError( errIllegalArgument);
1510 else if (sigma == 0.0)
1511 PushError( errDivisionByZero);
1512 else if (kum == 0.0) // Dichte
1513 PushDouble(phi((x-mue)/sigma)/sigma);
1514 else // Verteilung
1515 PushDouble(0.5 + gauss((x-mue)/sigma));
1516 }
1517 }
1518
1519 void ScInterpreter::ScLogNormDist()
1520 {
1521 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScLogNormDist" );
1522 if ( MustHaveParamCount( GetByte(), 3 ) )
1523 {
1524 double sigma = GetDouble(); // Stdabw
1525 double mue = GetDouble(); // Mittelwert
1526 double x = GetDouble(); // x
1527 if (sigma < 0.0)
1528 PushError( errIllegalArgument);
1529 else if (sigma == 0.0)
1530 PushError( errDivisionByZero);
1531 else if (x <= 0.0)
1532 PushIllegalArgument();
1533 else
1534 PushDouble(0.5 + gauss((log(x)-mue)/sigma));
1535 }
1536 }
1537
1538 void ScInterpreter::ScStdNormDist()
1539 {
1540 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScStdNormDist" );
1541 PushDouble(0.5 + gauss(GetDouble()));
1542 }
1543
1544 void ScInterpreter::ScExpDist()
1545 {
1546 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScExpDist" );
1547 if ( MustHaveParamCount( GetByte(), 3 ) )
1548 {
1549 double kum = GetDouble(); // 0 oder 1
1550 double lambda = GetDouble(); // lambda
1551 double x = GetDouble(); // x
1552 if (lambda <= 0.0)
1553 PushIllegalArgument();
1554 else if (kum == 0.0) // Dichte
1555 {
1556 if (x >= 0.0)
1557 PushDouble(lambda * exp(-lambda*x));
1558 else
1559 PushInt(0);
1560 }
1561 else // Verteilung
1562 {
1563 if (x > 0.0)
1564 PushDouble(1.0 - exp(-lambda*x));
1565 else
1566 PushInt(0);
1567 }
1568 }
1569 }
1570
1571 void ScInterpreter::ScTDist()
1572 {
1573 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTDist" );
1574 if ( !MustHaveParamCount( GetByte(), 3 ) )
1575 return;
1576 double fFlag = ::rtl::math::approxFloor(GetDouble());
1577 double fDF = ::rtl::math::approxFloor(GetDouble());
1578 double T = GetDouble();
1579 if (fDF < 1.0 || T < 0.0 || (fFlag != 1.0 && fFlag != 2.0) )
1580 {
1581 PushIllegalArgument();
1582 return;
1583 }
1584 double R = GetTDist(T, fDF);
1585 if (fFlag == 1.0)
1586 PushDouble(R);
1587 else
1588 PushDouble(2.0*R);
1589 }
1590
1591 void ScInterpreter::ScFDist()
1592 {
1593 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFDist" );
1594 if ( !MustHaveParamCount( GetByte(), 3 ) )
1595 return;
1596 double fF2 = ::rtl::math::approxFloor(GetDouble());
1597 double fF1 = ::rtl::math::approxFloor(GetDouble());
1598 double fF = GetDouble();
1599 if (fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10)
1600 {
1601 PushIllegalArgument();
1602 return;
1603 }
1604 PushDouble(GetFDist(fF, fF1, fF2));
1605 }
1606
1607 void ScInterpreter::ScChiDist()
1608 {
1609 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScChiDist" );
1610 double fResult;
1611 if ( !MustHaveParamCount( GetByte(), 2 ) )
1612 return;
1613 double fDF = ::rtl::math::approxFloor(GetDouble());
1614 double fChi = GetDouble();
1615 if (fDF < 1.0) // x<=0 returns 1, see ODFF 6.17.10
1616 {
1617 PushIllegalArgument();
1618 return;
1619 }
1620 fResult = GetChiDist( fChi, fDF);
1621 if (nGlobalError)
1622 {
1623 PushError( nGlobalError);
1624 return;
1625 }
1626 PushDouble(fResult);
1627 }
1628
1629 void ScInterpreter::ScWeibull()
1630 {
1631 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScWeibull" );
1632 if ( MustHaveParamCount( GetByte(), 4 ) )
1633 {
1634 double kum = GetDouble(); // 0 oder 1
1635 double beta = GetDouble(); // beta
1636 double alpha = GetDouble(); // alpha
1637 double x = GetDouble(); // x
1638 if (alpha <= 0.0 || beta <= 0.0 || x < 0.0)
1639 PushIllegalArgument();
1640 else if (kum == 0.0) // Dichte
1641 PushDouble(alpha/pow(beta,alpha)*pow(x,alpha-1.0)*
1642 exp(-pow(x/beta,alpha)));
1643 else // Verteilung
1644 PushDouble(1.0 - exp(-pow(x/beta,alpha)));
1645 }
1646 }
1647
1648 void ScInterpreter::ScPoissonDist()
1649 {
1650 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPoissonDist" );
1651 BYTE nParamCount = GetByte();
1652 if ( MustHaveParamCount( nParamCount, 2, 3 ) )
1653 {
1654 double kum = (nParamCount == 3 ? GetDouble() : 1); // 0 oder 1
1655 double lambda = GetDouble(); // Mittelwert
1656 double x = ::rtl::math::approxFloor(GetDouble()); // x
1657 if (lambda < 0.0 || x < 0.0)
1658 PushIllegalArgument();
1659 else if (kum == 0.0) // Dichte
1660 {
1661 if (lambda == 0.0)
1662 PushInt(0);
1663 else
1664 {
1665 double fPoissonVar = 1.0;
1666 for ( double f = 0.0; f < x; ++f )
1667 fPoissonVar *= lambda / ( f + 1.0 );
1668 PushDouble( fPoissonVar*exp( -lambda ) );
1669 }
1670 }
1671 else // Verteilung
1672 {
1673 if (lambda == 0.0)
1674 PushInt(1);
1675 else
1676 {
1677 double sum = 1.0;
1678 double fFak = 1.0;
1679 ULONG nEnd = (ULONG) x;
1680 for (ULONG i = 1; i <= nEnd; i++)
1681 {
1682 fFak *= (double)i;
1683 sum += pow( lambda, (double)i ) / fFak;
1684 }
1685 sum *= exp(-lambda);
1686 PushDouble(sum);
1687 }
1688 }
1689 }
1690 }
1691
1692 /** Local function used in the calculation of the hypergeometric distribution.
1693 */
1694 void lcl_PutFactorialElements( ::std::vector< double >& cn, double fLower, double fUpper, double fBase )
1695 {
1696 for ( double i = fLower; i <= fUpper; ++i )
1697 {
1698 double fVal = fBase - i;
1699 if ( fVal > 1.0 )
1700 cn.push_back( fVal );
1701 }
1702 }
1703
1704 /** Calculates a value of the hypergeometric distribution.
1705
1706 The algorithm is designed to avoid unnecessary multiplications and division
1707 by expanding all factorial elements (9 of them). It is done by excluding
1708 those ranges that overlap in the numerator and the denominator. This allows
1709 for a fast calculation for large values which would otherwise cause an overflow
1710 in the intermediate values.
1711
1712 @author Kohei Yoshida <kohei@openoffice.org>
1713
1714 @see #i47296#
1715
1716 */
1717 void ScInterpreter::ScHypGeomDist()
1718 {
1719 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScHypGeomDist" );
1720 const size_t nMaxArraySize = 500000; // arbitrary max array size
1721
1722 if ( !MustHaveParamCount( GetByte(), 4 ) )
1723 return;
1724
1725 double N = ::rtl::math::approxFloor(GetDouble());
1726 double M = ::rtl::math::approxFloor(GetDouble());
1727 double n = ::rtl::math::approxFloor(GetDouble());
1728 double x = ::rtl::math::approxFloor(GetDouble());
1729
1730 if( (x < 0.0) || (n < x) || (M < x) || (N < n) || (N < M) || (x < n - N + M) )
1731 {
1732 PushIllegalArgument();
1733 return;
1734 }
1735
1736 typedef ::std::vector< double > HypContainer;
1737 HypContainer cnNumer, cnDenom;
1738
1739 size_t nEstContainerSize = static_cast<size_t>( x + ::std::min( n, M ) );
1740 size_t nMaxSize = ::std::min( cnNumer.max_size(), nMaxArraySize );
1741 if ( nEstContainerSize > nMaxSize )
1742 {
1743 PushNoValue();
1744 return;
1745 }
1746 cnNumer.reserve( nEstContainerSize + 10 );
1747 cnDenom.reserve( nEstContainerSize + 10 );
1748
1749 // Trim coefficient C first
1750 double fCNumVarUpper = N - n - M + x - 1.0;
1751 double fCDenomVarLower = 1.0;
1752 if ( N - n - M + x >= M - x + 1.0 )
1753 {
1754 fCNumVarUpper = M - x - 1.0;
1755 fCDenomVarLower = N - n - 2.0*(M - x) + 1.0;
1756 }
1757
1758 #ifdef DBG_UTIL
1759 double fCNumLower = N - n - fCNumVarUpper;
1760 #endif
1761 double fCDenomUpper = N - n - M + x + 1.0 - fCDenomVarLower;
1762
1763 double fDNumVarLower = n - M;
1764
1765 if ( n >= M + 1.0 )
1766 {
1767 if ( N - M < n + 1.0 )
1768 {
1769 // Case 1
1770
1771 if ( N - n < n + 1.0 )
1772 {
1773 // no overlap
1774 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1775 lcl_PutFactorialElements( cnDenom, 0.0, N - n - 1.0, N );
1776 }
1777 else
1778 {
1779 // overlap
1780 DBG_ASSERT( fCNumLower < n + 1.0, "ScHypGeomDist: wrong assertion" );
1781 lcl_PutFactorialElements( cnNumer, N - 2.0*n, fCNumVarUpper, N - n );
1782 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1783 }
1784
1785 DBG_ASSERT( fCDenomUpper <= N - M, "ScHypGeomDist: wrong assertion" );
1786
1787 if ( fCDenomUpper < n - x + 1.0 )
1788 // no overlap
1789 lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
1790 else
1791 {
1792 // overlap
1793 lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1794
1795 fCDenomUpper = n - x;
1796 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1797 }
1798 }
1799 else
1800 {
1801 // Case 2
1802
1803 if ( n > M - 1.0 )
1804 {
1805 // no overlap
1806 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1807 lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
1808 }
1809 else
1810 {
1811 lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
1812 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1813 }
1814
1815 DBG_ASSERT( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
1816
1817 if ( fCDenomUpper < n - x + 1.0 )
1818 // no overlap
1819 lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - n + x, N - M + 1.0 );
1820 else
1821 {
1822 lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1823 fCDenomUpper = n - x;
1824 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1825 }
1826 }
1827
1828 DBG_ASSERT( fCDenomUpper <= M, "ScHypGeomDist: wrong assertion" );
1829 }
1830 else
1831 {
1832 if ( N - M < M + 1.0 )
1833 {
1834 // Case 3
1835
1836 if ( N - n < M + 1.0 )
1837 {
1838 // No overlap
1839 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1840 lcl_PutFactorialElements( cnDenom, 0.0, N - M - 1.0, N );
1841 }
1842 else
1843 {
1844 lcl_PutFactorialElements( cnNumer, N - n - M, fCNumVarUpper, N - n );
1845 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1846 }
1847
1848 if ( n - x + 1.0 > fCDenomUpper )
1849 // No overlap
1850 lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
1851 else
1852 {
1853 // Overlap
1854 lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1855
1856 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1857 fCDenomUpper = n - x;
1858 }
1859 }
1860 else
1861 {
1862 // Case 4
1863
1864 DBG_ASSERT( M >= n - x, "ScHypGeomDist: wrong assertion" );
1865 DBG_ASSERT( M - x <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
1866
1867 if ( N - n < N - M + 1.0 )
1868 {
1869 // No overlap
1870 lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
1871 lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
1872 }
1873 else
1874 {
1875 // Overlap
1876 DBG_ASSERT( fCNumLower <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
1877
1878 lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
1879 lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
1880 }
1881
1882 if ( n - x + 1.0 > fCDenomUpper )
1883 // No overlap
1884 lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - n + x, N - M + 1.0 );
1885 else if ( M >= fCDenomUpper )
1886 {
1887 lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
1888
1889 fCDenomUpper = n - x;
1890 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1891 }
1892 else
1893 {
1894 DBG_ASSERT( M <= fCDenomUpper, "ScHypGeomDist: wrong assertion" );
1895 lcl_PutFactorialElements( cnDenom, fCDenomVarLower, N - n - 2.0*M + x,
1896 N - n - M + x + 1.0 );
1897
1898 fCDenomUpper = n - x;
1899 fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
1900 }
1901 }
1902
1903 DBG_ASSERT( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
1904
1905 fDNumVarLower = 0.0;
1906 }
1907
1908 double nDNumVarUpper = fCDenomUpper < x + 1.0 ? n - x - 1.0 : n - fCDenomUpper - 1.0;
1909 double nDDenomVarLower = fCDenomUpper < x + 1.0 ? fCDenomVarLower : N - n - M + 1.0;
1910 lcl_PutFactorialElements( cnNumer, fDNumVarLower, nDNumVarUpper, n );
1911 lcl_PutFactorialElements( cnDenom, nDDenomVarLower, N - n - M + x, N - n - M + x + 1.0 );
1912
1913 ::std::sort( cnNumer.begin(), cnNumer.end() );
1914 ::std::sort( cnDenom.begin(), cnDenom.end() );
1915 HypContainer::reverse_iterator it1 = cnNumer.rbegin(), it1End = cnNumer.rend();
1916 HypContainer::reverse_iterator it2 = cnDenom.rbegin(), it2End = cnDenom.rend();
1917
1918 double fFactor = 1.0;
1919 for ( ; it1 != it1End || it2 != it2End; )
1920 {
1921 double fEnum = 1.0, fDenom = 1.0;
1922 if ( it1 != it1End )
1923 fEnum = *it1++;
1924 if ( it2 != it2End )
1925 fDenom = *it2++;
1926 fFactor *= fEnum / fDenom;
1927 }
1928
1929 PushDouble(fFactor);
1930 }
1931
1932 void ScInterpreter::ScGammaDist()
1933 {
1934 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGammaDist" );
1935 BYTE nParamCount = GetByte();
1936 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
1937 return;
1938 double bCumulative;
1939 if (nParamCount == 4)
1940 bCumulative = GetBool();
1941 else
1942 bCumulative = true;
1943 double fBeta = GetDouble(); // scale
1944 double fAlpha = GetDouble(); // shape
1945 double fX = GetDouble(); // x
1946 if (fAlpha <= 0.0 || fBeta <= 0.0)
1947 PushIllegalArgument();
1948 else
1949 {
1950 if (bCumulative) // distribution
1951 PushDouble( GetGammaDist( fX, fAlpha, fBeta));
1952 else // density
1953 PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta));
1954 }
1955 }
1956
1957 void ScInterpreter::ScNormInv()
1958 {
1959 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScNormInv" );
1960 if ( MustHaveParamCount( GetByte(), 3 ) )
1961 {
1962 double sigma = GetDouble();
1963 double mue = GetDouble();
1964 double x = GetDouble();
1965 if (sigma <= 0.0 || x < 0.0 || x > 1.0)
1966 PushIllegalArgument();
1967 else if (x == 0.0 || x == 1.0)
1968 PushNoValue();
1969 else
1970 PushDouble(gaussinv(x)*sigma + mue);
1971 }
1972 }
1973
1974 void ScInterpreter::ScSNormInv()
1975 {
1976 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSNormInv" );
1977 double x = GetDouble();
1978 if (x < 0.0 || x > 1.0)
1979 PushIllegalArgument();
1980 else if (x == 0.0 || x == 1.0)
1981 PushNoValue();
1982 else
1983 PushDouble(gaussinv(x));
1984 }
1985
1986 void ScInterpreter::ScLogNormInv()
1987 {
1988 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScLogNormInv" );
1989 if ( MustHaveParamCount( GetByte(), 3 ) )
1990 {
1991 double sigma = GetDouble(); // Stdabw
1992 double mue = GetDouble(); // Mittelwert
1993 double y = GetDouble(); // y
1994 if (sigma <= 0.0 || y <= 0.0 || y >= 1.0)
1995 PushIllegalArgument();
1996 else
1997 PushDouble(exp(mue+sigma*gaussinv(y)));
1998 }
1999 }
2000
2001 class ScGammaDistFunction : public ScDistFunc
2002 {
2003 ScInterpreter& rInt;
2004 double fp, fAlpha, fBeta;
2005
2006 public:
2007 ScGammaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2008 rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2009
2010 double GetValue( double x ) const { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); }
2011 };
2012
2013 void ScInterpreter::ScGammaInv()
2014 {
2015 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGammaInv" );
2016 if ( !MustHaveParamCount( GetByte(), 3 ) )
2017 return;
2018 double fBeta = GetDouble();
2019 double fAlpha = GetDouble();
2020 double fP = GetDouble();
2021 if (fAlpha <= 0.0 || fBeta <= 0.0 || fP < 0.0 || fP >= 1.0 )
2022 {
2023 PushIllegalArgument();
2024 return;
2025 }
2026 if (fP == 0.0)
2027 PushInt(0);
2028 else
2029 {
2030 bool bConvError;
2031 ScGammaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2032 double fStart = fAlpha * fBeta;
2033 double fVal = lcl_IterateInverse( aFunc, fStart*0.5, fStart, bConvError );
2034 if (bConvError)
2035 SetError(errNoConvergence);
2036 PushDouble(fVal);
2037 }
2038 }
2039
2040 class ScBetaDistFunction : public ScDistFunc
2041 {
2042 ScInterpreter& rInt;
2043 double fp, fAlpha, fBeta;
2044
2045 public:
2046 ScBetaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
2047 rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
2048
2049 double GetValue( double x ) const { return fp - rInt.GetBetaDist(x, fAlpha, fBeta); }
2050 };
2051
2052 void ScInterpreter::ScBetaInv()
2053 {
2054 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScBetaInv" );
2055 BYTE nParamCount = GetByte();
2056 if ( !MustHaveParamCount( nParamCount, 3, 5 ) )
2057 return;
2058 double fP, fA, fB, fAlpha, fBeta;
2059 if (nParamCount == 5)
2060 fB = GetDouble();
2061 else
2062 fB = 1.0;
2063 if (nParamCount >= 4)
2064 fA = GetDouble();
2065 else
2066 fA = 0.0;
2067 fBeta = GetDouble();
2068 fAlpha = GetDouble();
2069 fP = GetDouble();
2070 if (fP < 0.0 || fP >= 1.0 || fA == fB || fAlpha <= 0.0 || fBeta <= 0.0)
2071 {
2072 PushIllegalArgument();
2073 return;
2074 }
2075 if (fP == 0.0)
2076 PushInt(0);
2077 else
2078 {
2079 bool bConvError;
2080 ScBetaDistFunction aFunc( *this, fP, fAlpha, fBeta );
2081 // 0..1 as range for iteration so it isn't extended beyond the valid range
2082 double fVal = lcl_IterateInverse( aFunc, 0.0, 1.0, bConvError );
2083 if (bConvError)
2084 PushError( errNoConvergence);
2085 else
2086 PushDouble(fA + fVal*(fB-fA)); // scale to (A,B)
2087 }
2088 }
2089
2090 // Achtung: T, F und Chi
2091 // sind monoton fallend,
2092 // deshalb 1-Dist als Funktion
2093
2094 class ScTDistFunction : public ScDistFunc
2095 {
2096 ScInterpreter& rInt;
2097 double fp, fDF;
2098
2099 public:
2100 ScTDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2101 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2102
2103 double GetValue( double x ) const { return fp - 2 * rInt.GetTDist(x, fDF); }
2104 };
2105
2106 void ScInterpreter::ScTInv()
2107 {
2108 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTInv" );
2109 if ( !MustHaveParamCount( GetByte(), 2 ) )
2110 return;
2111 double fDF = ::rtl::math::approxFloor(GetDouble());
2112 double fP = GetDouble();
2113 if (fDF < 1.0 || fDF >= 1.0E5 || fP <= 0.0 || fP > 1.0 )
2114 {
2115 PushIllegalArgument();
2116 return;
2117 }
2118
2119 bool bConvError;
2120 ScTDistFunction aFunc( *this, fP, fDF );
2121 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2122 if (bConvError)
2123 SetError(errNoConvergence);
2124 PushDouble(fVal);
2125 }
2126
2127 class ScFDistFunction : public ScDistFunc
2128 {
2129 ScInterpreter& rInt;
2130 double fp, fF1, fF2;
2131
2132 public:
2133 ScFDistFunction( ScInterpreter& rI, double fpVal, double fF1Val, double fF2Val ) :
2134 rInt(rI), fp(fpVal), fF1(fF1Val), fF2(fF2Val) {}
2135
2136 double GetValue( double x ) const { return fp - rInt.GetFDist(x, fF1, fF2); }
2137 };
2138
2139 void ScInterpreter::ScFInv()
2140 {
2141 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFInv" );
2142 if ( !MustHaveParamCount( GetByte(), 3 ) )
2143 return;
2144 double fF2 = ::rtl::math::approxFloor(GetDouble());
2145 double fF1 = ::rtl::math::approxFloor(GetDouble());
2146 double fP = GetDouble();
2147 if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
2148 {
2149 PushIllegalArgument();
2150 return;
2151 }
2152
2153 bool bConvError;
2154 ScFDistFunction aFunc( *this, fP, fF1, fF2 );
2155 double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
2156 if (bConvError)
2157 SetError(errNoConvergence);
2158 PushDouble(fVal);
2159 }
2160
2161 class ScChiDistFunction : public ScDistFunc
2162 {
2163 ScInterpreter& rInt;
2164 double fp, fDF;
2165
2166 public:
2167 ScChiDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2168 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2169
2170 double GetValue( double x ) const { return fp - rInt.GetChiDist(x, fDF); }
2171 };
2172
2173 void ScInterpreter::ScChiInv()
2174 {
2175 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScChiInv" );
2176 if ( !MustHaveParamCount( GetByte(), 2 ) )
2177 return;
2178 double fDF = ::rtl::math::approxFloor(GetDouble());
2179 double fP = GetDouble();
2180 if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
2181 {
2182 PushIllegalArgument();
2183 return;
2184 }
2185
2186 bool bConvError;
2187 ScChiDistFunction aFunc( *this, fP, fDF );
2188 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2189 if (bConvError)
2190 SetError(errNoConvergence);
2191 PushDouble(fVal);
2192 }
2193
2194 /***********************************************/
2195 class ScChiSqDistFunction : public ScDistFunc
2196 {
2197 ScInterpreter& rInt;
2198 double fp, fDF;
2199
2200 public:
2201 ScChiSqDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
2202 rInt(rI), fp(fpVal), fDF(fDFVal) {}
2203
2204 double GetValue( double x ) const { return fp - rInt.GetChiSqDistCDF(x, fDF); }
2205 };
2206
2207
2208 void ScInterpreter::ScChiSqInv()
2209 {
2210 if ( !MustHaveParamCount( GetByte(), 2 ) )
2211 return;
2212 double fDF = ::rtl::math::approxFloor(GetDouble());
2213 double fP = GetDouble();
2214 if (fDF < 1.0 || fP < 0.0 || fP >= 1.0 )
2215 {
2216 PushIllegalArgument();
2217 return;
2218 }
2219
2220 bool bConvError;
2221 ScChiSqDistFunction aFunc( *this, fP, fDF );
2222 double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
2223 if (bConvError)
2224 SetError(errNoConvergence);
2225 PushDouble(fVal);
2226 }
2227
2228
2229 void ScInterpreter::ScConfidence()
2230 {
2231 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScConfidence" );
2232 if ( MustHaveParamCount( GetByte(), 3 ) )
2233 {
2234 double n = ::rtl::math::approxFloor(GetDouble());
2235 double sigma = GetDouble();
2236 double alpha = GetDouble();
2237 if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
2238 PushIllegalArgument();
2239 else
2240 PushDouble( gaussinv(1.0-alpha/2.0) * sigma/sqrt(n) );
2241 }
2242 }
2243
2244 void ScInterpreter::ScZTest()
2245 {
2246 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScZTest" );
2247 BYTE nParamCount = GetByte();
2248 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
2249 return;
2250 double sigma = 0.0, mue, x;
2251 if (nParamCount == 3)
2252 {
2253 sigma = GetDouble();
2254 if (sigma <= 0.0)
2255 {
2256 PushIllegalArgument();
2257 return;
2258 }
2259 }
2260 x = GetDouble();
2261
2262 double fSum = 0.0;
2263 double fSumSqr = 0.0;
2264 double fVal;
2265 double rValCount = 0.0;
2266 switch (GetStackType())
2267 {
2268 case formula::svDouble :
2269 {
2270 fVal = GetDouble();
2271 fSum += fVal;
2272 fSumSqr += fVal*fVal;
2273 rValCount++;
2274 }
2275 break;
2276 case svSingleRef :
2277 {
2278 ScAddress aAdr;
2279 PopSingleRef( aAdr );
2280 ScBaseCell* pCell = GetCell( aAdr );
2281 if (HasCellValueData(pCell))
2282 {
2283 fVal = GetCellValue( aAdr, pCell );
2284 fSum += fVal;
2285 fSumSqr += fVal*fVal;
2286 rValCount++;
2287 }
2288 }
2289 break;
2290 case svRefList :
2291 case formula::svDoubleRef :
2292 {
2293 short nParam = 1;
2294 size_t nRefInList = 0;
2295 while (nParam-- > 0)
2296 {
2297 ScRange aRange;
2298 USHORT nErr = 0;
2299 PopDoubleRef( aRange, nParam, nRefInList);
2300 ScValueIterator aValIter(pDok, aRange, glSubTotal);
2301 if (aValIter.GetFirst(fVal, nErr))
2302 {
2303 fSum += fVal;
2304 fSumSqr += fVal*fVal;
2305 rValCount++;
2306 while ((nErr == 0) && aValIter.GetNext(fVal, nErr))
2307 {
2308 fSum += fVal;
2309 fSumSqr += fVal*fVal;
2310 rValCount++;
2311 }
2312 SetError(nErr);
2313 }
2314 }
2315 }
2316 break;
2317 case svMatrix :
2318 {
2319 ScMatrixRef pMat = PopMatrix();
2320 if (pMat)
2321 {
2322 SCSIZE nCount = pMat->GetElementCount();
2323 if (pMat->IsNumeric())
2324 {
2325 for ( SCSIZE i = 0; i < nCount; i++ )
2326 {
2327 fVal= pMat->GetDouble(i);
2328 fSum += fVal;
2329 fSumSqr += fVal * fVal;
2330 rValCount++;
2331 }
2332 }
2333 else
2334 {
2335 for (SCSIZE i = 0; i < nCount; i++)
2336 if (!pMat->IsString(i))
2337 {
2338 fVal= pMat->GetDouble(i);
2339 fSum += fVal;
2340 fSumSqr += fVal * fVal;
2341 rValCount++;
2342 }
2343 }
2344 }
2345 }
2346 break;
2347 default : SetError(errIllegalParameter); break;
2348 }
2349 if (rValCount <= 1.0)
2350 PushError( errDivisionByZero);
2351 else
2352 {
2353 mue = fSum/rValCount;
2354 if (nParamCount != 3)
2355 sigma = (fSumSqr - fSum*fSum/rValCount)/(rValCount-1.0);
2356
2357 PushDouble(0.5 - gauss((mue-x)/sqrt(sigma/rValCount)));
2358 }
2359 }
2360 bool ScInterpreter::CalculateTest(BOOL _bTemplin
2361 ,const SCSIZE nC1, const SCSIZE nC2,const SCSIZE nR1,const SCSIZE nR2
2362 ,const ScMatrixRef& pMat1,const ScMatrixRef& pMat2
2363 ,double& fT,double& fF)
2364 {
2365 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateTest" );
2366 double fCount1 = 0.0;
2367 double fCount2 = 0.0;
2368 double fSum1 = 0.0;
2369 double fSumSqr1 = 0.0;
2370 double fSum2 = 0.0;
2371 double fSumSqr2 = 0.0;
2372 double fVal;
2373 SCSIZE i,j;
2374 for (i = 0; i < nC1; i++)
2375 for (j = 0; j < nR1; j++)
2376 {
2377 if (!pMat1->IsString(i,j))
2378 {
2379 fVal = pMat1->GetDouble(i,j);
2380 fSum1 += fVal;
2381 fSumSqr1 += fVal * fVal;
2382 fCount1++;
2383 }
2384 }
2385 for (i = 0; i < nC2; i++)
2386 for (j = 0; j < nR2; j++)
2387 {
2388 if (!pMat2->IsString(i,j))
2389 {
2390 fVal = pMat2->GetDouble(i,j);
2391 fSum2 += fVal;
2392 fSumSqr2 += fVal * fVal;
2393 fCount2++;
2394 }
2395 }
2396 if (fCount1 < 2.0 || fCount2 < 2.0)
2397 {
2398 PushNoValue();
2399 return false;
2400 } // if (fCount1 < 2.0 || fCount2 < 2.0)
2401 if ( _bTemplin )
2402 {
2403 double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0)/fCount1;
2404 double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0)/fCount2;
2405 if (fS1 + fS2 == 0.0)
2406 {
2407 PushNoValue();
2408 return false;
2409 }
2410 fT = fabs(fSum1/fCount1 - fSum2/fCount2)/sqrt(fS1+fS2);
2411 double c = fS1/(fS1+fS2);
2412 // s.u. fF = ::rtl::math::approxFloor(1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0)));
2413 // fF = ::rtl::math::approxFloor((fS1+fS2)*(fS1+fS2)/(fS1*fS1/(fCount1-1.0) + fS2*fS2/(fCount2-1.0)));
2414
2415 // GetTDist wird mit GetBetaDist berechnet und kommt auch mit nicht ganzzahligen
2416 // Freiheitsgraden klar. Dann stimmt das Ergebnis auch mit Excel ueberein (#52406#):
2417 fF = 1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0));
2418 }
2419 else
2420 {
2421 // laut Bronstein-Semendjajew
2422 double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1) / (fCount1 - 1.0); // Varianz
2423 double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2) / (fCount2 - 1.0);
2424 fT = fabs( fSum1/fCount1 - fSum2/fCount2 ) /
2425 sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
2426 sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
2427 fF = fCount1 + fCount2 - 2;
2428 }
2429 return true;
2430 }
2431 void ScInterpreter::ScTTest()
2432 {
2433 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTTest" );
2434 if ( !MustHaveParamCount( GetByte(), 4 ) )
2435 return;
2436 double fTyp = ::rtl::math::approxFloor(GetDouble());
2437 double fAnz = ::rtl::math::approxFloor(GetDouble());
2438 if (fAnz != 1.0 && fAnz != 2.0)
2439 {
2440 PushIllegalArgument();
2441 return;
2442 }
2443
2444 ScMatrixRef pMat2 = GetMatrix();
2445 ScMatrixRef pMat1 = GetMatrix();
2446 if (!pMat1 || !pMat2)
2447 {
2448 PushIllegalParameter();
2449 return;
2450 }
2451 double fT, fF;
2452 SCSIZE nC1, nC2;
2453 SCSIZE nR1, nR2;
2454 SCSIZE i, j;
2455 pMat1->GetDimensions(nC1, nR1);
2456 pMat2->GetDimensions(nC2, nR2);
2457 if (fTyp == 1.0)
2458 {
2459 if (nC1 != nC2 || nR1 != nR2)
2460 {
2461 PushIllegalArgument();
2462 return;
2463 }
2464 double fCount = 0.0;
2465 double fSum1 = 0.0;
2466 double fSum2 = 0.0;
2467 double fSumSqrD = 0.0;
2468 double fVal1, fVal2;
2469 for (i = 0; i < nC1; i++)
2470 for (j = 0; j < nR1; j++)
2471 {
2472 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
2473 {
2474 fVal1 = pMat1->GetDouble(i,j);
2475 fVal2 = pMat2->GetDouble(i,j);
2476 fSum1 += fVal1;
2477 fSum2 += fVal2;
2478 fSumSqrD += (fVal1 - fVal2)*(fVal1 - fVal2);
2479 fCount++;
2480 }
2481 }
2482 if (fCount < 1.0)
2483 {
2484 PushNoValue();
2485 return;
2486 }
2487 fT = sqrt(fCount-1.0) * fabs(fSum1 - fSum2) /
2488 sqrt(fCount * fSumSqrD - (fSum1-fSum2)*(fSum1-fSum2));
2489 fF = fCount - 1.0;
2490 }
2491 else if (fTyp == 2.0)
2492 {
2493 CalculateTest(FALSE,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF);
2494 }
2495 else if (fTyp == 3.0)
2496 {
2497 CalculateTest(TRUE,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF);
2498 }
2499
2500 else
2501 {
2502 PushIllegalArgument();
2503 return;
2504 }
2505 if (fAnz == 1.0)
2506 PushDouble(GetTDist(fT, fF));
2507 else
2508 PushDouble(2.0*GetTDist(fT, fF));
2509 }
2510
2511 void ScInterpreter::ScFTest()
2512 {
2513 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScFTest" );
2514 if ( !MustHaveParamCount( GetByte(), 2 ) )
2515 return;
2516 ScMatrixRef pMat2 = GetMatrix();
2517 ScMatrixRef pMat1 = GetMatrix();
2518 if (!pMat1 || !pMat2)
2519 {
2520 PushIllegalParameter();
2521 return;
2522 }
2523 SCSIZE nC1, nC2;
2524 SCSIZE nR1, nR2;
2525 SCSIZE i, j;
2526 pMat1->GetDimensions(nC1, nR1);
2527 pMat2->GetDimensions(nC2, nR2);
2528 double fCount1 = 0.0;
2529 double fCount2 = 0.0;
2530 double fSum1 = 0.0;
2531 double fSumSqr1 = 0.0;
2532 double fSum2 = 0.0;
2533 double fSumSqr2 = 0.0;
2534 double fVal;
2535 for (i = 0; i < nC1; i++)
2536 for (j = 0; j < nR1; j++)
2537 {
2538 if (!pMat1->IsString(i,j))
2539 {
2540 fVal = pMat1->GetDouble(i,j);
2541 fSum1 += fVal;
2542 fSumSqr1 += fVal * fVal;
2543 fCount1++;
2544 }
2545 }
2546 for (i = 0; i < nC2; i++)
2547 for (j = 0; j < nR2; j++)
2548 {
2549 if (!pMat2->IsString(i,j))
2550 {
2551 fVal = pMat2->GetDouble(i,j);
2552 fSum2 += fVal;
2553 fSumSqr2 += fVal * fVal;
2554 fCount2++;
2555 }
2556 }
2557 if (fCount1 < 2.0 || fCount2 < 2.0)
2558 {
2559 PushNoValue();
2560 return;
2561 }
2562 double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0);
2563 double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0);
2564 if (fS1 == 0.0 || fS2 == 0.0)
2565 {
2566 PushNoValue();
2567 return;
2568 }
2569 double fF, fF1, fF2;
2570 if (fS1 > fS2)
2571 {
2572 fF = fS1/fS2;
2573 fF1 = fCount1-1.0;
2574 fF2 = fCount2-1.0;
2575 }
2576 else
2577 {
2578 fF = fS2/fS1;
2579 fF1 = fCount2-1.0;
2580 fF2 = fCount1-1.0;
2581 }
2582 PushDouble(2.0*GetFDist(fF, fF1, fF2));
2583 /*
2584 double Z = (pow(fF,1.0/3.0)*(1.0-2.0/(9.0*fF2)) - (1.0-2.0/(9.0*fF1))) /
2585 sqrt(2.0/(9.0*fF1) + pow(fF,2.0/3.0)*2.0/(9.0*fF2));
2586 PushDouble(1.0-2.0*gauss(Z));
2587 */
2588 }
2589
2590 void ScInterpreter::ScChiTest()
2591 {
2592 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScChiTest" );
2593 if ( !MustHaveParamCount( GetByte(), 2 ) )
2594 return;
2595 ScMatrixRef pMat2 = GetMatrix();
2596 ScMatrixRef pMat1 = GetMatrix();
2597 if (!pMat1 || !pMat2)
2598 {
2599 PushIllegalParameter();
2600 return;
2601 }
2602 SCSIZE nC1, nC2;
2603 SCSIZE nR1, nR2;
2604 pMat1->GetDimensions(nC1, nR1);
2605 pMat2->GetDimensions(nC2, nR2);
2606 if (nR1 != nR2 || nC1 != nC2)
2607 {
2608 PushIllegalArgument();
2609 return;
2610 }
2611 double fChi = 0.0;
2612 for (SCSIZE i = 0; i < nC1; i++)
2613 {
2614 for (SCSIZE j = 0; j < nR1; j++)
2615 {
2616 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
2617 {
2618 double fValX = pMat1->GetDouble(i,j);
2619 double fValE = pMat2->GetDouble(i,j);
2620 fChi += (fValX - fValE) * (fValX - fValE) / fValE;
2621 }
2622 else
2623 {
2624 PushIllegalArgument();
2625 return;
2626 }
2627 }
2628 }
2629 double fDF;
2630 if (nC1 == 1 || nR1 == 1)
2631 {
2632 fDF = (double)(nC1*nR1 - 1);
2633 if (fDF == 0.0)
2634 {
2635 PushNoValue();
2636 return;
2637 }
2638 }
2639 else
2640 fDF = (double)(nC1-1)*(double)(nR1-1);
2641 PushDouble(GetChiDist(fChi, fDF));
2642 /*
2643 double fX, fS, fT, fG;
2644 fX = 1.0;
2645 for (double fi = fDF; fi >= 2.0; fi -= 2.0)
2646 fX *= fChi/fi;
2647 fX *= exp(-fChi/2.0);
2648 if (fmod(fDF, 2.0) != 0.0)
2649 fX *= sqrt(2.0*fChi/F_PI);
2650 fS = 1.0;
2651 fT = 1.0;
2652 fG = fDF;
2653 while (fT >= 1.0E-7)
2654 {
2655 fG += 2.0;
2656 fT *= fChi/fG;
2657 fS += fT;
2658 }
2659 PushDouble(1.0 - fX*fS);
2660 */
2661 }
2662
2663 void ScInterpreter::ScKurt()
2664 {
2665 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScKurt" );
2666 double fSum,fCount,vSum;
2667 std::vector<double> values;
2668 if ( !CalculateSkew(fSum,fCount,vSum,values) )
2669 return;
2670
2671 if (fCount == 0.0)
2672 {
2673 PushError( errDivisionByZero);
2674 return;
2675 }
2676
2677 double fMean = fSum / fCount;
2678
2679 for (size_t i = 0; i < values.size(); i++)
2680 vSum += (values[i] - fMean) * (values[i] - fMean);
2681
2682 double fStdDev = sqrt(vSum / (fCount - 1.0));
2683 double dx = 0.0;
2684 double xpower4 = 0.0;
2685
2686 if (fStdDev == 0.0)
2687 {
2688 PushError( errDivisionByZero);
2689 return;
2690 }
2691
2692 for (size_t i = 0; i < values.size(); i++)
2693 {
2694 dx = (values[i] - fMean) / fStdDev;
2695 xpower4 = xpower4 + (dx * dx * dx * dx);
2696 }
2697
2698 double k_d = (fCount - 2.0) * (fCount - 3.0);
2699 double k_l = fCount * (fCount + 1.0) / ((fCount - 1.0) * k_d);
2700 double k_t = 3.0 * (fCount - 1.0) * (fCount - 1.0) / k_d;
2701
2702 PushDouble(xpower4 * k_l - k_t);
2703 }
2704
2705 void ScInterpreter::ScHarMean()
2706 {
2707 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScHarMean" );
2708 short nParamCount = GetByte();
2709 double nVal = 0.0;
2710 double nValCount = 0.0;
2711 ScAddress aAdr;
2712 ScRange aRange;
2713 size_t nRefInList = 0;
2714 while ((nGlobalError == 0) && (nParamCount-- > 0))
2715 {
2716 switch (GetStackType())
2717 {
2718 case formula::svDouble :
2719 {
2720 double x = GetDouble();
2721 if (x > 0.0)
2722 {
2723 nVal += 1.0/x;
2724 nValCount++;
2725 }
2726 else
2727 SetError( errIllegalArgument);
2728 break;
2729 }
2730 case svSingleRef :
2731 {
2732 PopSingleRef( aAdr );
2733 ScBaseCell* pCell = GetCell( aAdr );
2734 if (HasCellValueData(pCell))
2735 {
2736 double x = GetCellValue( aAdr, pCell );
2737 if (x > 0.0)
2738 {
2739 nVal += 1.0/x;
2740 nValCount++;
2741 }
2742 else
2743 SetError( errIllegalArgument);
2744 }
2745 break;
2746 }
2747 case formula::svDoubleRef :
2748 case svRefList :
2749 {
2750 USHORT nErr = 0;
2751 PopDoubleRef( aRange, nParamCount, nRefInList);
2752 double nCellVal;
2753 ScValueIterator aValIter(pDok, aRange, glSubTotal);
2754 if (aValIter.GetFirst(nCellVal, nErr))
2755 {
2756 if (nCellVal > 0.0)
2757 {
2758 nVal += 1.0/nCellVal;
2759 nValCount++;
2760 }
2761 else
2762 SetError( errIllegalArgument);
2763 SetError(nErr);
2764 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
2765 {
2766 if (nCellVal > 0.0)
2767 {
2768 nVal += 1.0/nCellVal;
2769 nValCount++;
2770 }
2771 else
2772 SetError( errIllegalArgument);
2773 }
2774 SetError(nErr);
2775 }
2776 }
2777 break;
2778 case svMatrix :
2779 {
2780 ScMatrixRef pMat = PopMatrix();
2781 if (pMat)
2782 {
2783 SCSIZE nCount = pMat->GetElementCount();
2784 if (pMat->IsNumeric())
2785 {
2786 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
2787 {
2788 double x = pMat->GetDouble(nElem);
2789 if (x > 0.0)
2790 {
2791 nVal += 1.0/x;
2792 nValCount++;
2793 }
2794 else
2795 SetError( errIllegalArgument);
2796 }
2797 }
2798 else
2799 {
2800 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
2801 if (!pMat->IsString(nElem))
2802 {
2803 double x = pMat->GetDouble(nElem);
2804 if (x > 0.0)
2805 {
2806 nVal += 1.0/x;
2807 nValCount++;
2808 }
2809 else
2810 SetError( errIllegalArgument);
2811 }
2812 }
2813 }
2814 }
2815 break;
2816 default : SetError(errIllegalParameter); break;
2817 }
2818 }
2819 if (nGlobalError == 0)
2820 PushDouble((double)nValCount/nVal);
2821 else
2822 PushError( nGlobalError);
2823 }
2824
2825 void ScInterpreter::ScGeoMean()
2826 {
2827 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScGeoMean" );
2828 short nParamCount = GetByte();
2829 double nVal = 0.0;
2830 double nValCount = 0.0;
2831 ScAddress aAdr;
2832 ScRange aRange;
2833
2834 size_t nRefInList = 0;
2835 while ((nGlobalError == 0) && (nParamCount-- > 0))
2836 {
2837 switch (GetStackType())
2838 {
2839 case formula::svDouble :
2840 {
2841 double x = GetDouble();
2842 if (x > 0.0)
2843 {
2844 nVal += log(x);
2845 nValCount++;
2846 }
2847 else
2848 SetError( errIllegalArgument);
2849 break;
2850 }
2851 case svSingleRef :
2852 {
2853 PopSingleRef( aAdr );
2854 ScBaseCell* pCell = GetCell( aAdr );
2855 if (HasCellValueData(pCell))
2856 {
2857 double x = GetCellValue( aAdr, pCell );
2858 if (x > 0.0)
2859 {
2860 nVal += log(x);
2861 nValCount++;
2862 }
2863 else
2864 SetError( errIllegalArgument);
2865 }
2866 break;
2867 }
2868 case formula::svDoubleRef :
2869 case svRefList :
2870 {
2871 USHORT nErr = 0;
2872 PopDoubleRef( aRange, nParamCount, nRefInList);
2873 double nCellVal;
2874 ScValueIterator aValIter(pDok, aRange, glSubTotal);
2875 if (aValIter.GetFirst(nCellVal, nErr))
2876 {
2877 if (nCellVal > 0.0)
2878 {
2879 nVal += log(nCellVal);
2880 nValCount++;
2881 }
2882 else
2883 SetError( errIllegalArgument);
2884 SetError(nErr);
2885 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
2886 {
2887 if (nCellVal > 0.0)
2888 {
2889 nVal += log(nCellVal);
2890 nValCount++;
2891 }
2892 else
2893 SetError( errIllegalArgument);
2894 }
2895 SetError(nErr);
2896 }
2897 }
2898 break;
2899 case svMatrix :
2900 {
2901 ScMatrixRef pMat = PopMatrix();
2902 if (pMat)
2903 {
2904 SCSIZE nCount = pMat->GetElementCount();
2905 if (pMat->IsNumeric())
2906 {
2907 for (SCSIZE ui = 0; ui < nCount; ui++)
2908 {
2909 double x = pMat->GetDouble(ui);
2910 if (x > 0.0)
2911 {
2912 nVal += log(x);
2913 nValCount++;
2914 }
2915 else
2916 SetError( errIllegalArgument);
2917 }
2918 }
2919 else
2920 {
2921 for (SCSIZE ui = 0; ui < nCount; ui++)
2922 if (!pMat->IsString(ui))
2923 {
2924 double x = pMat->GetDouble(ui);
2925 if (x > 0.0)
2926 {
2927 nVal += log(x);
2928 nValCount++;
2929 }
2930 else
2931 SetError( errIllegalArgument);
2932 }
2933 }
2934 }
2935 }
2936 break;
2937 default : SetError(errIllegalParameter); break;
2938 }
2939 }
2940 if (nGlobalError == 0)
2941 PushDouble(exp(nVal / nValCount));
2942 else
2943 PushError( nGlobalError);
2944 }
2945
2946 void ScInterpreter::ScStandard()
2947 {
2948 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScStandard" );
2949 if ( MustHaveParamCount( GetByte(), 3 ) )
2950 {
2951 double sigma = GetDouble();
2952 double mue = GetDouble();
2953 double x = GetDouble();
2954 if (sigma < 0.0)
2955 PushError( errIllegalArgument);
2956 else if (sigma == 0.0)
2957 PushError( errDivisionByZero);
2958 else
2959 PushDouble((x-mue)/sigma);
2960 }
2961 }
2962 bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values)
2963 {
2964 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateSkew" );
2965 short nParamCount = GetByte();
2966 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
2967 return false;
2968
2969 fSum = 0.0;
2970 fCount = 0.0;
2971 vSum = 0.0;
2972 double fVal = 0.0;
2973 ScAddress aAdr;
2974 ScRange aRange;
2975 size_t nRefInList = 0;
2976 while (nParamCount-- > 0)
2977 {
2978 switch (GetStackType())
2979 {
2980 case formula::svDouble :
2981 {
2982 fVal = GetDouble();
2983 fSum += fVal;
2984 values.push_back(fVal);
2985 fCount++;
2986 }
2987 break;
2988 case svSingleRef :
2989 {
2990 PopSingleRef( aAdr );
2991 ScBaseCell* pCell = GetCell( aAdr );
2992 if (HasCellValueData(pCell))
2993 {
2994 fVal = GetCellValue( aAdr, pCell );
2995 fSum += fVal;
2996 values.push_back(fVal);
2997 fCount++;
2998 }
2999 }
3000 break;
3001 case formula::svDoubleRef :
3002 case svRefList :
3003 {
3004 PopDoubleRef( aRange, nParamCount, nRefInList);
3005 USHORT nErr = 0;
3006 ScValueIterator aValIter(pDok, aRange);
3007 if (aValIter.GetFirst(fVal, nErr))
3008 {
3009 fSum += fVal;
3010 values.push_back(fVal);
3011 fCount++;
3012 SetError(nErr);
3013 while ((nErr == 0) && aValIter.GetNext(fVal, nErr))
3014 {
3015 fSum += fVal;
3016 values.push_back(fVal);
3017 fCount++;
3018 }
3019 SetError(nErr);
3020 }
3021 }
3022 break;
3023 case svMatrix :
3024 {
3025 ScMatrixRef pMat = PopMatrix();
3026 if (pMat)
3027 {
3028 SCSIZE nCount = pMat->GetElementCount();
3029 if (pMat->IsNumeric())
3030 {
3031 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3032 {
3033 fVal = pMat->GetDouble(nElem);
3034 fSum += fVal;
3035 values.push_back(fVal);
3036 fCount++;
3037 }
3038 }
3039 else
3040 {
3041 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3042 if (!pMat->IsString(nElem))
3043 {
3044 fVal = pMat->GetDouble(nElem);
3045 fSum += fVal;
3046 values.push_back(fVal);
3047 fCount++;
3048 }
3049 }
3050 }
3051 }
3052 break;
3053 default :
3054 SetError(errIllegalParameter);
3055 break;
3056 }
3057 }
3058
3059 if (nGlobalError)
3060 {
3061 PushError( nGlobalError);
3062 return false;
3063 } // if (nGlobalError)
3064 return true;
3065 }
3066
3067 void ScInterpreter::ScSkew()
3068 {
3069 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSkew" );
3070 double fSum,fCount,vSum;
3071 std::vector<double> values;
3072 if ( !CalculateSkew(fSum,fCount,vSum,values) )
3073 return;
3074
3075 double fMean = fSum / fCount;
3076
3077 for (size_t i = 0; i < values.size(); i++)
3078 vSum += (values[i] - fMean) * (values[i] - fMean);
3079
3080 double fStdDev = sqrt(vSum / (fCount - 1.0));
3081 double dx = 0.0;
3082 double xcube = 0.0;
3083
3084 if (fStdDev == 0)
3085 {
3086 PushIllegalArgument();
3087 return;
3088 }
3089
3090 for (size_t i = 0; i < values.size(); i++)
3091 {
3092 dx = (values[i] - fMean) / fStdDev;
3093 xcube = xcube + (dx * dx * dx);
3094 }
3095
3096 PushDouble(((xcube * fCount) / (fCount - 1.0)) / (fCount - 2.0));
3097 }
3098
3099 double ScInterpreter::GetMedian( vector<double> & rArray )
3100 {
3101 size_t nSize = rArray.size();
3102 if (rArray.empty() || nSize == 0 || nGlobalError)
3103 {
3104 SetError( errNoValue);
3105 return 0.0;
3106 }
3107
3108 // Upper median.
3109 size_t nMid = nSize / 2;
3110 vector<double>::iterator iMid = rArray.begin() + nMid;
3111 ::std::nth_element( rArray.begin(), iMid, rArray.end());
3112 if (nSize & 1)
3113 return *iMid; // Lower and upper median are equal.
3114 else
3115 {
3116 double fUp = *iMid;
3117 // Lower median.
3118 iMid = rArray.begin() + nMid - 1;
3119 ::std::nth_element( rArray.begin(), iMid, rArray.end());
3120 return (fUp + *iMid) / 2;
3121 }
3122 }
3123
3124 void ScInterpreter::ScMedian()
3125 {
3126 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScMedian" );
3127 BYTE nParamCount = GetByte();
3128 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3129 return;
3130 vector<double> aArray;
3131 GetNumberSequenceArray( nParamCount, aArray);
3132 PushDouble( GetMedian( aArray));
3133 }
3134
3135 double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile )
3136 {
3137 size_t nSize = rArray.size();
3138 if (rArray.empty() || nSize == 0 || nGlobalError)
3139 {
3140 SetError( errNoValue);
3141 return 0.0;
3142 }
3143
3144 if (nSize == 1)
3145 return rArray[0];
3146 else
3147 {
3148 size_t nIndex = (size_t)::rtl::math::approxFloor( fPercentile * (nSize-1));
3149 double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1));
3150 DBG_ASSERT(nIndex < nSize, "GetPercentile: wrong index(1)");
3151 vector<double>::iterator iter = rArray.begin() + nIndex;
3152 ::std::nth_element( rArray.begin(), iter, rArray.end());
3153 if (fDiff == 0.0)
3154 return *iter;
3155 else
3156 {
3157 DBG_ASSERT(nIndex < nSize-1, "GetPercentile: wrong index(2)");
3158 double fVal = *iter;
3159 iter = rArray.begin() + nIndex+1;
3160 ::std::nth_element( rArray.begin(), iter, rArray.end());
3161 return fVal + fDiff * (*iter - fVal);
3162 }
3163 }
3164 }
3165
3166 void ScInterpreter::ScPercentile()
3167 {
3168 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPercentile" );
3169 if ( !MustHaveParamCount( GetByte(), 2 ) )
3170 return;
3171 double alpha = GetDouble();
3172 if (alpha < 0.0 || alpha > 1.0)
3173 {
3174 PushIllegalArgument();
3175 return;
3176 }
3177 vector<double> aArray;
3178 GetNumberSequenceArray( 1, aArray);
3179 PushDouble( GetPercentile( aArray, alpha));
3180 }
3181
3182 void ScInterpreter::ScQuartile()
3183 {
3184 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScQuartile" );
3185 if ( !MustHaveParamCount( GetByte(), 2 ) )
3186 return;
3187 double fFlag = ::rtl::math::approxFloor(GetDouble());
3188 if (fFlag < 0.0 || fFlag > 4.0)
3189 {
3190 PushIllegalArgument();
3191 return;
3192 }
3193 vector<double> aArray;
3194 GetNumberSequenceArray( 1, aArray);
3195 PushDouble( fFlag == 2.0 ? GetMedian( aArray) : GetPercentile( aArray, 0.25 * fFlag));
3196 }
3197
3198 void ScInterpreter::ScModalValue()
3199 {
3200 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScModalValue" );
3201 BYTE nParamCount = GetByte();
3202 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3203 return;
3204 vector<double> aSortArray;
3205 GetSortArray(nParamCount, aSortArray);
3206 SCSIZE nSize = aSortArray.size();
3207 if (aSortArray.empty() || nSize == 0 || nGlobalError)
3208 PushNoValue();
3209 else
3210 {
3211 SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
3212 double nOldVal = aSortArray[0];
3213 SCSIZE i;
3214
3215 for ( i = 1; i < nSize; i++)
3216 {
3217 if (aSortArray[i] == nOldVal)
3218 nCount++;
3219 else
3220 {
3221 nOldVal = aSortArray[i];
3222 if (nCount > nMax)
3223 {
3224 nMax = nCount;
3225 nMaxIndex = i-1;
3226 }
3227 nCount = 1;
3228 }
3229 }
3230 if (nCount > nMax)
3231 {
3232 nMax = nCount;
3233 nMaxIndex = i-1;
3234 }
3235 if (nMax == 1 && nCount == 1)
3236 PushNoValue();
3237 else if (nMax == 1)
3238 PushDouble(nOldVal);
3239 else
3240 PushDouble(aSortArray[nMaxIndex]);
3241 }
3242 }
3243
3244 void ScInterpreter::CalculateSmallLarge(BOOL bSmall)
3245 {
3246 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateSmallLarge" );
3247 if ( !MustHaveParamCount( GetByte(), 2 ) )
3248 return;
3249 double f = ::rtl::math::approxFloor(GetDouble());
3250 if (f < 1.0)
3251 {
3252 PushIllegalArgument();
3253 return;
3254 }
3255 SCSIZE k = static_cast<SCSIZE>(f);
3256 vector<double> aSortArray;
3257 /* TODO: using nth_element() is best for one single value, but LARGE/SMALL
3258 * actually are defined to return an array of values if an array of
3259 * positions was passed, in which case, depending on the number of values,
3260 * we may or will need a real sorted array again, see #i32345. */
3261 //GetSortArray(1, aSortArray);
3262 GetNumberSequenceArray(1, aSortArray);
3263 SCSIZE nSize = aSortArray.size();
3264 if (aSortArray.empty() || nSize == 0 || nGlobalError || nSize < k)
3265 PushNoValue();
3266 else
3267 {
3268 // TODO: the sorted case for array: PushDouble( aSortArray[ bSmall ? k-1 : nSize-k ] );
3269 vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
3270 ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
3271 PushDouble( *iPos);
3272 }
3273 }
3274
3275 void ScInterpreter::ScLarge()
3276 {
3277 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScLarge" );
3278 CalculateSmallLarge(FALSE);
3279 }
3280
3281 void ScInterpreter::ScSmall()
3282 {
3283 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSmall" );
3284 CalculateSmallLarge(TRUE);
3285 }
3286
3287 void ScInterpreter::ScPercentrank()
3288 {
3289 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPercentrank" );
3290 BYTE nParamCount = GetByte();
3291 if ( !MustHaveParamCount( nParamCount, 2 ) )
3292 return;
3293 #if 0
3294 /* wird nicht unterstuetzt
3295 double fPrec;
3296 if (nParamCount == 3)
3297 {
3298 fPrec = ::rtl::math::approxFloor(GetDouble());
3299 if (fPrec < 1.0)
3300 {
3301 PushIllegalArgument();
3302 return;
3303 }
3304 }
3305 else
3306 fPrec = 3.0;
3307 */
3308 #endif
3309 double fNum = GetDouble();
3310 vector<double> aSortArray;
3311 GetSortArray(1, aSortArray);
3312 SCSIZE nSize = aSortArray.size();
3313 if (aSortArray.empty() || nSize == 0 || nGlobalError)
3314 PushNoValue();
3315 else
3316 {
3317 if (fNum < aSortArray[0] || fNum > aSortArray[nSize-1])
3318 PushNoValue();
3319 else if ( nSize == 1 )
3320 PushDouble(1.0); // fNum == pSortArray[0], see test above
3321 else
3322 {
3323 double fRes;
3324 SCSIZE nOldCount = 0;
3325 double fOldVal = aSortArray[0];
3326 SCSIZE i;
3327 for (i = 1; i < nSize && aSortArray[i] < fNum; i++)
3328 {
3329 if (aSortArray[i] != fOldVal)
3330 {
3331 nOldCount = i;
3332 fOldVal = aSortArray[i];
3333 }
3334 }
3335 if (aSortArray[i] != fOldVal)
3336 nOldCount = i;
3337 if (fNum == aSortArray[i])
3338 fRes = (double)nOldCount/(double)(nSize-1);
3339 else
3340 {
3341 // #75312# nOldCount is the count of smaller entries
3342 // fNum is between pSortArray[nOldCount-1] and pSortArray[nOldCount]
3343 // use linear interpolation to find a position between the entries
3344
3345 if ( nOldCount == 0 )
3346 {
3347 DBG_ERROR("should not happen");
3348 fRes = 0.0;
3349 }
3350 else
3351 {
3352 double fFract = ( fNum - aSortArray[nOldCount-1] ) /
3353 ( aSortArray[nOldCount] - aSortArray[nOldCount-1] );
3354 fRes = ( (double)(nOldCount-1)+fFract )/(double)(nSize-1);
3355 }
3356 }
3357 PushDouble(fRes);
3358 }
3359 }
3360 }
3361
3362 void ScInterpreter::ScTrimMean()
3363 {
3364 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScTrimMean" );
3365 if ( !MustHaveParamCount( GetByte(), 2 ) )
3366 return;
3367 double alpha = GetDouble();
3368 if (alpha < 0.0 || alpha >= 1.0)
3369 {
3370 PushIllegalArgument();
3371 return;
3372 }
3373 vector<double> aSortArray;
3374 GetSortArray(1, aSortArray);
3375 SCSIZE nSize = aSortArray.size();
3376 if (aSortArray.empty() || nSize == 0 || nGlobalError)
3377 PushNoValue();
3378 else
3379 {
3380 ULONG nIndex = (ULONG) ::rtl::math::approxFloor(alpha*(double)nSize);
3381 if (nIndex % 2 != 0)
3382 nIndex--;
3383 nIndex /= 2;
3384 DBG_ASSERT(nIndex < nSize, "ScTrimMean: falscher Index");
3385 double fSum = 0.0;
3386 for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
3387 fSum += aSortArray[i];
3388 PushDouble(fSum/(double)(nSize-2*nIndex));
3389 }
3390 }
3391
3392 void ScInterpreter::GetNumberSequenceArray( BYTE nParamCount, vector<double>& rArray )
3393 {
3394 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::GetSortArray" );
3395 ScAddress aAdr;
3396 ScRange aRange;
3397 short nParam = nParamCount;
3398 size_t nRefInList = 0;
3399 while (nParam-- > 0)
3400 {
3401 switch (GetStackType())
3402 {
3403 case formula::svDouble :
3404 rArray.push_back( PopDouble());
3405 break;
3406 case svSingleRef :
3407 {
3408 PopSingleRef( aAdr );
3409 ScBaseCell* pCell = GetCell( aAdr );
3410 if (HasCellValueData(pCell))
3411 rArray.push_back( GetCellValue( aAdr, pCell));
3412 }
3413 break;
3414 case formula::svDoubleRef :
3415 case svRefList :
3416 {
3417 PopDoubleRef( aRange, nParam, nRefInList);
3418 if (nGlobalError)
3419 break;
3420
3421 aRange.Justify();
3422 SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
3423 nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
3424 rArray.reserve( rArray.size() + nCellCount);
3425
3426 USHORT nErr = 0;
3427 double fCellVal;
3428 ScValueIterator aValIter(pDok, aRange);
3429 if (aValIter.GetFirst( fCellVal, nErr))
3430 {
3431 rArray.push_back( fCellVal);
3432 SetError(nErr);
3433 while ((nErr == 0) && aValIter.GetNext( fCellVal, nErr))
3434 rArray.push_back( fCellVal);
3435 SetError(nErr);
3436 }
3437 }
3438 break;
3439 case svMatrix :
3440 {
3441 ScMatrixRef pMat = PopMatrix();
3442 if (!pMat)
3443 break;
3444
3445 SCSIZE nCount = pMat->GetElementCount();
3446 rArray.reserve( rArray.size() + nCount);
3447 if (pMat->IsNumeric())
3448 {
3449 for (SCSIZE i = 0; i < nCount; ++i)
3450 rArray.push_back( pMat->GetDouble(i));
3451 }
3452 else
3453 {
3454 for (SCSIZE i = 0; i < nCount; ++i)
3455 if (!pMat->IsString(i))
3456 rArray.push_back( pMat->GetDouble(i));
3457 }
3458 }
3459 break;
3460 default :
3461 PopError();
3462 SetError( errIllegalParameter);
3463 break;
3464 }
3465 if (nGlobalError)
3466 break; // while
3467 }
3468 // nParam > 0 in case of error, clean stack environment and obtain earlier
3469 // error if there was one.
3470 while (nParam-- > 0)
3471 PopError();
3472 }
3473
3474 void ScInterpreter::GetSortArray( BYTE nParamCount, vector<double>& rSortArray, vector<long>* pIndexOrder )
3475 {
3476 GetNumberSequenceArray( nParamCount, rSortArray);
3477
3478 if (rSortArray.size() > MAX_ANZ_DOUBLE_FOR_SORT)
3479 SetError( errStackOverflow);
3480 else if (rSortArray.empty())
3481 SetError( errNoValue);
3482
3483 if (nGlobalError == 0)
3484 QuickSort( rSortArray, pIndexOrder);
3485 }
3486
3487 static void lcl_QuickSort( long nLo, long nHi, vector<double>& rSortArray, vector<long>* pIndexOrder )
3488 {
3489 // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().
3490
3491 using ::std::swap;
3492
3493 if (nHi - nLo == 1)
3494 {
3495 if (rSortArray[nLo] > rSortArray[nHi])
3496 {
3497 swap(rSortArray[nLo], rSortArray[nHi]);
3498 if (pIndexOrder)
3499 swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
3500 }
3501 return;
3502 }
3503
3504 long ni = nLo;
3505 long nj = nHi;
3506 do
3507 {
3508 double fLo = rSortArray[nLo];
3509 while (ni <= nHi && rSortArray[ni] < fLo) ni++;
3510 while (nj >= nLo && fLo < rSortArray[nj]) nj--;
3511 if (ni <= nj)
3512 {
3513 if (ni != nj)
3514 {
3515 swap(rSortArray[ni], rSortArray[nj]);
3516 if (pIndexOrder)
3517 swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
3518 }
3519
3520 ++ni;
3521 --nj;
3522 }
3523 }
3524 while (ni < nj);
3525
3526 if ((nj - nLo) < (nHi - ni))
3527 {
3528 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
3529 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
3530 }
3531 else
3532 {
3533 if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
3534 if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
3535 }
3536 }
3537
3538 void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<long>* pIndexOrder )
3539 {
3540 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::QuickSort" );
3541 long n = static_cast<long>(rSortArray.size());
3542
3543 if (pIndexOrder)
3544 {
3545 pIndexOrder->clear();
3546 pIndexOrder->reserve(n);
3547 for (long i = 0; i < n; ++i)
3548 pIndexOrder->push_back(i);
3549 }
3550
3551 if (n < 2)
3552 return;
3553
3554 size_t nValCount = rSortArray.size();
3555 for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
3556 {
3557 size_t nInd = rand() % (int) (nValCount-1);
3558 ::std::swap( rSortArray[i], rSortArray[nInd]);
3559 if (pIndexOrder)
3560 ::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
3561 }
3562
3563 lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
3564 }
3565
3566 void ScInterpreter::ScRank()
3567 {
3568 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScRank" );
3569 BYTE nParamCount = GetByte();
3570 if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
3571 return;
3572 BOOL bDescending;
3573 if (nParamCount == 3)
3574 bDescending = GetBool();
3575 else
3576 bDescending = FALSE;
3577 double fCount = 1.0;
3578 BOOL bValid = FALSE;
3579 switch (GetStackType())
3580 {
3581 case formula::svDouble :
3582 {
3583 double x = GetDouble();
3584 double fVal = GetDouble();
3585 if (x == fVal)
3586 bValid = TRUE;
3587 break;
3588 }
3589 case svSingleRef :
3590 {
3591 ScAddress aAdr;
3592 PopSingleRef( aAdr );
3593 double fVal = GetDouble();
3594 ScBaseCell* pCell = GetCell( aAdr );
3595 if (HasCellValueData(pCell))
3596 {
3597 double x = GetCellValue( aAdr, pCell );
3598 if (x == fVal)
3599 bValid = TRUE;
3600 }
3601 break;
3602 }
3603 case formula::svDoubleRef :
3604 case svRefList :
3605 {
3606 ScRange aRange;
3607 short nParam = 1;
3608 size_t nRefInList = 0;
3609 while (nParam-- > 0)
3610 {
3611 USHORT nErr = 0;
3612 // Preserve stack until all RefList elements are done!
3613 USHORT nSaveSP = sp;
3614 PopDoubleRef( aRange, nParam, nRefInList);
3615 if (nParam)
3616 --sp; // simulate pop
3617 double fVal = GetDouble();
3618 if (nParam)
3619 sp = nSaveSP;
3620 double nCellVal;
3621 ScValueIterator aValIter(pDok, aRange, glSubTotal);
3622 if (aValIter.GetFirst(nCellVal, nErr))
3623 {
3624 if (nCellVal == fVal)
3625 bValid = TRUE;
3626 else if ((!bDescending && nCellVal > fVal) ||
3627 (bDescending && nCellVal < fVal) )
3628 fCount++;
3629 SetError(nErr);
3630 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
3631 {
3632 if (nCellVal == fVal)
3633 bValid = TRUE;
3634 else if ((!bDescending && nCellVal > fVal) ||
3635 (bDescending && nCellVal < fVal) )
3636 fCount++;
3637 }
3638 }
3639 SetError(nErr);
3640 }
3641 }
3642 break;
3643 case svMatrix :
3644 {
3645 ScMatrixRef pMat = PopMatrix();
3646 double fVal = GetDouble();
3647 if (pMat)
3648 {
3649 SCSIZE nCount = pMat->GetElementCount();
3650 if (pMat->IsNumeric())
3651 {
3652 for (SCSIZE i = 0; i < nCount; i++)
3653 {
3654 double x = pMat->GetDouble(i);
3655 if (x == fVal)
3656 bValid = TRUE;
3657 else if ((!bDescending && x > fVal) ||
3658 (bDescending && x < fVal) )
3659 fCount++;
3660 }
3661 }
3662 else
3663 {
3664 for (SCSIZE i = 0; i < nCount; i++)
3665 if (!pMat->IsString(i))
3666 {
3667 double x = pMat->GetDouble(i);
3668 if (x == fVal)
3669 bValid = TRUE;
3670 else if ((!bDescending && x > fVal) ||
3671 (bDescending && x < fVal) )
3672 fCount++;
3673 }
3674 }
3675 }
3676 }
3677 break;
3678 default : SetError(errIllegalParameter); break;
3679 }
3680 if (bValid)
3681 PushDouble(fCount);
3682 else
3683 PushNoValue();
3684 }
3685
3686 void ScInterpreter::ScAveDev()
3687 {
3688 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScAveDev" );
3689 BYTE nParamCount = GetByte();
3690 if ( !MustHaveParamCountMin( nParamCount, 1 ) )
3691 return;
3692 USHORT SaveSP = sp;
3693 double nMiddle = 0.0;
3694 double rVal = 0.0;
3695 double rValCount = 0.0;
3696 ScAddress aAdr;
3697 ScRange aRange;
3698 short nParam = nParamCount;
3699 size_t nRefInList = 0;
3700 while (nParam-- > 0)
3701 {
3702 switch (GetStackType())
3703 {
3704 case formula::svDouble :
3705 rVal += GetDouble();
3706 rValCount++;
3707 break;
3708 case svSingleRef :
3709 {
3710 PopSingleRef( aAdr );
3711 ScBaseCell* pCell = GetCell( aAdr );
3712 if (HasCellValueData(pCell))
3713 {
3714 rVal += GetCellValue( aAdr, pCell );
3715 rValCount++;
3716 }
3717 }
3718 break;
3719 case formula::svDoubleRef :
3720 case svRefList :
3721 {
3722 USHORT nErr = 0;
3723 double nCellVal;
3724 PopDoubleRef( aRange, nParam, nRefInList);
3725 ScValueIterator aValIter(pDok, aRange);
3726 if (aValIter.GetFirst(nCellVal, nErr))
3727 {
3728 rVal += nCellVal;
3729 rValCount++;
3730 SetError(nErr);
3731 while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr))
3732 {
3733 rVal += nCellVal;
3734 rValCount++;
3735 }
3736 SetError(nErr);
3737 }
3738 }
3739 break;
3740 case svMatrix :
3741 {
3742 ScMatrixRef pMat = PopMatrix();
3743 if (pMat)
3744 {
3745 SCSIZE nCount = pMat->GetElementCount();
3746 if (pMat->IsNumeric())
3747 {
3748 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3749 {
3750 rVal += pMat->GetDouble(nElem);
3751 rValCount++;
3752 }
3753 }
3754 else
3755 {
3756 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3757 if (!pMat->IsString(nElem))
3758 {
3759 rVal += pMat->GetDouble(nElem);
3760 rValCount++;
3761 }
3762 }
3763 }
3764 }
3765 break;
3766 default :
3767 SetError(errIllegalParameter);
3768 break;
3769 }
3770 }
3771 if (nGlobalError)
3772 {
3773 PushError( nGlobalError);
3774 return;
3775 }
3776 nMiddle = rVal / rValCount;
3777 sp = SaveSP;
3778 rVal = 0.0;
3779 nParam = nParamCount;
3780 nRefInList = 0;
3781 while (nParam-- > 0)
3782 {
3783 switch (GetStackType())
3784 {
3785 case formula::svDouble :
3786 rVal += fabs(GetDouble() - nMiddle);
3787 break;
3788 case svSingleRef :
3789 {
3790 PopSingleRef( aAdr );
3791 ScBaseCell* pCell = GetCell( aAdr );
3792 if (HasCellValueData(pCell))
3793 rVal += fabs(GetCellValue( aAdr, pCell ) - nMiddle);
3794 }
3795 break;
3796 case formula::svDoubleRef :
3797 case svRefList :
3798 {
3799 USHORT nErr = 0;
3800 double nCellVal;
3801 PopDoubleRef( aRange, nParam, nRefInList);
3802 ScValueIterator aValIter(pDok, aRange);
3803 if (aValIter.GetFirst(nCellVal, nErr))
3804 {
3805 rVal += (fabs(nCellVal - nMiddle));
3806 while (aValIter.GetNext(nCellVal, nErr))
3807 rVal += fabs(nCellVal - nMiddle);
3808 }
3809 }
3810 break;
3811 case svMatrix :
3812 {
3813 ScMatrixRef pMat = PopMatrix();
3814 if (pMat)
3815 {
3816 SCSIZE nCount = pMat->GetElementCount();
3817 if (pMat->IsNumeric())
3818 {
3819 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3820 {
3821 rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
3822 }
3823 }
3824 else
3825 {
3826 for (SCSIZE nElem = 0; nElem < nCount; nElem++)
3827 {
3828 if (!pMat->IsString(nElem))
3829 rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
3830 }
3831 }
3832 }
3833 }
3834 break;
3835 default : SetError(errIllegalParameter); break;
3836 }
3837 }
3838 PushDouble(rVal / rValCount);
3839 }
3840
3841 void ScInterpreter::ScDevSq()
3842 {
3843 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScDevSq" );
3844 double nVal;
3845 double nValCount;
3846 GetStVarParams(nVal, nValCount);
3847 PushDouble(nVal);
3848 }
3849
3850 void ScInterpreter::ScProbability()
3851 {
3852 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScProbability" );
3853 BYTE nParamCount = GetByte();
3854 if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
3855 return;
3856 double fUp, fLo;
3857 fUp = GetDouble();
3858 if (nParamCount == 4)
3859 fLo = GetDouble();
3860 else
3861 fLo = fUp;
3862 if (fLo > fUp)
3863 {
3864 double fTemp = fLo;
3865 fLo = fUp;
3866 fUp = fTemp;
3867 }
3868 ScMatrixRef pMatP = GetMatrix();
3869 ScMatrixRef pMatW = GetMatrix();
3870 if (!pMatP || !pMatW)
3871 PushIllegalParameter();
3872 else
3873 {
3874 SCSIZE nC1, nC2;
3875 SCSIZE nR1, nR2;
3876 pMatP->GetDimensions(nC1, nR1);
3877 pMatW->GetDimensions(nC2, nR2);
3878 if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
3879 nC2 == 0 || nR2 == 0)
3880 PushNA();
3881 else
3882 {
3883 double fSum = 0.0;
3884 double fRes = 0.0;
3885 BOOL bStop = FALSE;
3886 double fP, fW;
3887 SCSIZE nCount1 = nC1 * nR1;
3888 for ( SCSIZE i = 0; i < nCount1 && !bStop; i++ )
3889 {
3890 if (pMatP->IsValue(i) && pMatW->IsValue(i))
3891 {
3892 fP = pMatP->GetDouble(i);
3893 fW = pMatW->GetDouble(i);
3894 if (fP < 0.0 || fP > 1.0)
3895 bStop = TRUE;
3896 else
3897 {
3898 fSum += fP;
3899 if (fW >= fLo && fW <= fUp)
3900 fRes += fP;
3901 }
3902 }
3903 else
3904 SetError( errIllegalArgument);
3905 }
3906 if (bStop || fabs(fSum -1.0) > 1.0E-7)
3907 PushNoValue();
3908 else
3909 PushDouble(fRes);
3910 }
3911 }
3912 }
3913
3914 void ScInterpreter::ScCorrel()
3915 {
3916 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScCorrel" );
3917 // This is identical to ScPearson()
3918 ScPearson();
3919 }
3920
3921 void ScInterpreter::ScCovar()
3922 {
3923 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScCovar" );
3924 CalculatePearsonCovar(FALSE,FALSE);
3925 }
3926
3927 void ScInterpreter::ScPearson()
3928 {
3929 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScPearson" );
3930 CalculatePearsonCovar(TRUE,FALSE);
3931 }
3932 void ScInterpreter::CalculatePearsonCovar(BOOL _bPearson,BOOL _bStexy)
3933 {
3934 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculatePearsonCovar" );
3935 if ( !MustHaveParamCount( GetByte(), 2 ) )
3936 return;
3937 ScMatrixRef pMat1 = GetMatrix();
3938 ScMatrixRef pMat2 = GetMatrix();
3939 if (!pMat1 || !pMat2)
3940 {
3941 PushIllegalParameter();
3942 return;
3943 }
3944 SCSIZE nC1, nC2;
3945 SCSIZE nR1, nR2;
3946 pMat1->GetDimensions(nC1, nR1);
3947 pMat2->GetDimensions(nC2, nR2);
3948 if (nR1 != nR2 || nC1 != nC2)
3949 {
3950 PushIllegalArgument();
3951 return;
3952 }
3953 /* #i78250#
3954 * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
3955 * but the latter produces wrong results if the absolute values are high,
3956 * for example above 10^8
3957 */
3958 double fCount = 0.0;
3959 double fSumX = 0.0;
3960 double fSumY = 0.0;
3961 double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
3962 double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
3963 double fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2
3964 for (SCSIZE i = 0; i < nC1; i++)
3965 {
3966 for (SCSIZE j = 0; j < nR1; j++)
3967 {
3968 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
3969 {
3970 double fValX = pMat1->GetDouble(i,j);
3971 double fValY = pMat2->GetDouble(i,j);
3972 fSumX += fValX;
3973 fSumY += fValY;
3974 fCount++;
3975 }
3976 }
3977 }
3978 if (fCount < (_bStexy ? 3.0 : 1.0)) // fCount==1 is handled by checking denominator later on
3979 PushNoValue();
3980 else
3981 {
3982 const double fMeanX = fSumX / fCount;
3983 const double fMeanY = fSumY / fCount;
3984 for (SCSIZE i = 0; i < nC1; i++)
3985 {
3986 for (SCSIZE j = 0; j < nR1; j++)
3987 {
3988 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
3989 {
3990 const double fValX = pMat1->GetDouble(i,j);
3991 const double fValY = pMat2->GetDouble(i,j);
3992 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
3993 if ( _bPearson )
3994 {
3995 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
3996 fSumSqrDeltaY += (fValY - fMeanY) * (fValY - fMeanY);
3997 }
3998 }
3999 }
4000 } // for (SCSIZE i = 0; i < nC1; i++)
4001 if ( _bPearson )
4002 {
4003 if (fSumSqrDeltaX == 0.0 || ( !_bStexy && fSumSqrDeltaY == 0.0) )
4004 PushError( errDivisionByZero);
4005 else if ( _bStexy )
4006 PushDouble( sqrt( (fSumSqrDeltaY - fSumDeltaXDeltaY *
4007 fSumDeltaXDeltaY / fSumSqrDeltaX) / (fCount-2)));
4008 else
4009 PushDouble( fSumDeltaXDeltaY / sqrt( fSumSqrDeltaX * fSumSqrDeltaY));
4010 } // if ( _bPearson )
4011 else
4012 {
4013 PushDouble( fSumDeltaXDeltaY / fCount);
4014 }
4015 }
4016 }
4017
4018 void ScInterpreter::ScRSQ()
4019 {
4020 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScRSQ" );
4021 // Same as ScPearson()*ScPearson()
4022 ScPearson();
4023 if (!nGlobalError)
4024 {
4025 switch (GetStackType())
4026 {
4027 case formula::svDouble:
4028 {
4029 double fVal = PopDouble();
4030 PushDouble( fVal * fVal);
4031 }
4032 break;
4033 default:
4034 PopError();
4035 PushNoValue();
4036 }
4037 }
4038 }
4039
4040 void ScInterpreter::ScSTEXY()
4041 {
4042 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSTEXY" );
4043 CalculatePearsonCovar(TRUE,TRUE);
4044 }
4045 void ScInterpreter::CalculateSlopeIntercept(BOOL bSlope)
4046 {
4047 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::CalculateSlopeIntercept" );
4048 if ( !MustHaveParamCount( GetByte(), 2 ) )
4049 return;
4050 ScMatrixRef pMat1 = GetMatrix();
4051 ScMatrixRef pMat2 = GetMatrix();
4052 if (!pMat1 || !pMat2)
4053 {
4054 PushIllegalParameter();
4055 return;
4056 }
4057 SCSIZE nC1, nC2;
4058 SCSIZE nR1, nR2;
4059 pMat1->GetDimensions(nC1, nR1);
4060 pMat2->GetDimensions(nC2, nR2);
4061 if (nR1 != nR2 || nC1 != nC2)
4062 {
4063 PushIllegalArgument();
4064 return;
4065 }
4066 // #i78250# numerical stability improved
4067 double fCount = 0.0;
4068 double fSumX = 0.0;
4069 double fSumY = 0.0;
4070 double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4071 double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4072 for (SCSIZE i = 0; i < nC1; i++)
4073 {
4074 for (SCSIZE j = 0; j < nR1; j++)
4075 {
4076 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4077 {
4078 double fValX = pMat1->GetDouble(i,j);
4079 double fValY = pMat2->GetDouble(i,j);
4080 fSumX += fValX;
4081 fSumY += fValY;
4082 fCount++;
4083 }
4084 }
4085 }
4086 if (fCount < 1.0)
4087 PushNoValue();
4088 else
4089 {
4090 double fMeanX = fSumX / fCount;
4091 double fMeanY = fSumY / fCount;
4092 for (SCSIZE i = 0; i < nC1; i++)
4093 {
4094 for (SCSIZE j = 0; j < nR1; j++)
4095 {
4096 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4097 {
4098 double fValX = pMat1->GetDouble(i,j);
4099 double fValY = pMat2->GetDouble(i,j);
4100 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4101 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4102 }
4103 }
4104 }
4105 if (fSumSqrDeltaX == 0.0)
4106 PushError( errDivisionByZero);
4107 else
4108 {
4109 if ( bSlope )
4110 PushDouble( fSumDeltaXDeltaY / fSumSqrDeltaX);
4111 else
4112 PushDouble( fMeanY - fSumDeltaXDeltaY / fSumSqrDeltaX * fMeanX);
4113 }
4114 }
4115 }
4116
4117 void ScInterpreter::ScSlope()
4118 {
4119 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScSlope" );
4120 CalculateSlopeIntercept(TRUE);
4121 }
4122
4123 void ScInterpreter::ScIntercept()
4124 {
4125 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScIntercept" );
4126 CalculateSlopeIntercept(FALSE);
4127 }
4128
4129 void ScInterpreter::ScForecast()
4130 {
4131 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "Eike.Rathke@sun.com", "ScInterpreter::ScForecast" );
4132 if ( !MustHaveParamCount( GetByte(), 3 ) )
4133 return;
4134 ScMatrixRef pMat1 = GetMatrix();
4135 ScMatrixRef pMat2 = GetMatrix();
4136 if (!pMat1 || !pMat2)
4137 {
4138 PushIllegalParameter();
4139 return;
4140 }
4141 SCSIZE nC1, nC2;
4142 SCSIZE nR1, nR2;
4143 pMat1->GetDimensions(nC1, nR1);
4144 pMat2->GetDimensions(nC2, nR2);
4145 if (nR1 != nR2 || nC1 != nC2)
4146 {
4147 PushIllegalArgument();
4148 return;
4149 }
4150 double fVal = GetDouble();
4151 // #i78250# numerical stability improved
4152 double fCount = 0.0;
4153 double fSumX = 0.0;
4154 double fSumY = 0.0;
4155 double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
4156 double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
4157 for (SCSIZE i = 0; i < nC1; i++)
4158 {
4159 for (SCSIZE j = 0; j < nR1; j++)
4160 {
4161 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4162 {
4163 double fValX = pMat1->GetDouble(i,j);
4164 double fValY = pMat2->GetDouble(i,j);
4165 fSumX += fValX;
4166 fSumY += fValY;
4167 fCount++;
4168 }
4169 }
4170 }
4171 if (fCount < 1.0)
4172 PushNoValue();
4173 else
4174 {
4175 double fMeanX = fSumX / fCount;
4176 double fMeanY = fSumY / fCount;
4177 for (SCSIZE i = 0; i < nC1; i++)
4178 {
4179 for (SCSIZE j = 0; j < nR1; j++)
4180 {
4181 if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
4182 {
4183 double fValX = pMat1->GetDouble(i,j);
4184 double fValY = pMat2->GetDouble(i,j);
4185 fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
4186 fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
4187 }
4188 }
4189 }
4190 if (fSumSqrDeltaX == 0.0)
4191 PushError( errDivisionByZero);
4192 else
4193 PushDouble( fMeanY + fSumDeltaXDeltaY / fSumSqrDeltaX * (fVal - fMeanX));
4194 }
4195 }
4196
vba interoperability for OOo