xref: /aoo41x/main/sc/source/core/tool/interpr6.cxx (revision b3f79822)
1 /**************************************************************
2  *
3  * Licensed to the Apache Software Foundation (ASF) under one
4  * or more contributor license agreements.  See the NOTICE file
5  * distributed with this work for additional information
6  * regarding copyright ownership.  The ASF licenses this file
7  * to you under the Apache License, Version 2.0 (the
8  * "License"); you may not use this file except in compliance
9  * with the License.  You may obtain a copy of the License at
10  *
11  *   http://www.apache.org/licenses/LICENSE-2.0
12  *
13  * Unless required by applicable law or agreed to in writing,
14  * software distributed under the License is distributed on an
15  * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16  * KIND, either express or implied.  See the License for the
17  * specific language governing permissions and limitations
18  * under the License.
19  *
20  *************************************************************/
21 
22 
23 
24 // MARKER(update_precomp.py): autogen include statement, do not remove
25 #include "precompiled_sc.hxx"
26 
27 // #include <math.h>
28 
29 #include <tools/debug.hxx>
30 #include <rtl/logfile.hxx>
31 #include "interpre.hxx"
32 
33 double const fHalfMachEps = 0.5 * ::std::numeric_limits<double>::epsilon();
34 
35 // The idea how this group of gamma functions is calculated, is
36 // based on the Cephes library
37 // online http://www.moshier.net/#Cephes [called 2008-02]
38 
39 /** You must ensure fA>0.0 && fX>0.0
40     valid results only if fX > fA+1.0
41     uses continued fraction with odd items */
GetGammaContFraction(double fA,double fX)42 double ScInterpreter::GetGammaContFraction( double fA, double fX )
43 {
44     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaContFraction" );
45 
46     double const fBigInv = ::std::numeric_limits<double>::epsilon();
47     double const fBig = 1.0/fBigInv;
48     double fCount = 0.0;
49     double fNum = 0.0;  // dummy value
50     double fY = 1.0 - fA;
51     double fDenom = fX + 2.0-fA;
52     double fPk = 0.0;   // dummy value
53     double fPkm1 = fX + 1.0;
54     double fPkm2 = 1.0;
55     double fQk = 1.0;   // dummy value
56     double fQkm1 = fDenom * fX;
57     double fQkm2 = fX;
58     double fApprox = fPkm1/fQkm1;
59     bool bFinished = false;
60     double fR = 0.0;    // dummy value
61     do
62     {
63         fCount = fCount +1.0;
64         fY = fY+ 1.0;
65         fNum = fY * fCount;
66         fDenom = fDenom +2.0;
67         fPk = fPkm1 * fDenom  -  fPkm2 * fNum;
68         fQk = fQkm1 * fDenom  -  fQkm2 * fNum;
69         if (fQk != 0.0)
70         {
71             fR = fPk/fQk;
72             bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps);
73             fApprox = fR;
74         }
75         fPkm2 = fPkm1;
76         fPkm1 = fPk;
77         fQkm2 = fQkm1;
78         fQkm1 = fQk;
79         if (fabs(fPk) > fBig)
80         {
81             // reduce a fraction does not change the value
82             fPkm2 = fPkm2 * fBigInv;
83             fPkm1 = fPkm1 * fBigInv;
84             fQkm2 = fQkm2 * fBigInv;
85             fQkm1 = fQkm1 * fBigInv;
86         }
87     } while (!bFinished && fCount<10000);
88     // most iterations, if fX==fAlpha+1.0; approx sqrt(fAlpha) iterations then
89     if (!bFinished)
90     {
91         SetError(errNoConvergence);
92     }
93     return fApprox;
94 }
95 
96 /** You must ensure fA>0.0 && fX>0.0
97     valid results only if fX <= fA+1.0
98     uses power series */
GetGammaSeries(double fA,double fX)99 double ScInterpreter::GetGammaSeries( double fA, double fX )
100 {
101     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaSeries" );
102     double fDenomfactor = fA;
103     double fSummand = 1.0/fA;
104     double fSum = fSummand;
105     int nCount=1;
106     do
107     {
108         fDenomfactor = fDenomfactor + 1.0;
109         fSummand = fSummand * fX/fDenomfactor;
110         fSum = fSum + fSummand;
111         nCount = nCount+1;
112     } while ( fSummand/fSum > fHalfMachEps && nCount<=10000);
113     // large amount of iterations will be carried out for huge fAlpha, even
114     // if fX <= fAlpha+1.0
115     if (nCount>10000)
116     {
117         SetError(errNoConvergence);
118     }
119     return fSum;
120 }
121 
122 /** You must ensure fA>0.0 && fX>0.0) */
GetLowRegIGamma(double fA,double fX)123 double ScInterpreter::GetLowRegIGamma( double fA, double fX )
124 {
125     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetLowRegIGamma" );
126     double fLnFactor = fA * log(fX) - fX - GetLogGamma(fA);
127     double fFactor = exp(fLnFactor);    // Do we need more accuracy than exp(ln()) has?
128     if (fX>fA+1.0)  // includes fX>1.0; 1-GetUpRegIGamma, continued fraction
129         return 1.0 - fFactor * GetGammaContFraction(fA,fX);
130     else            // fX<=1.0 || fX<=fA+1.0, series
131         return fFactor * GetGammaSeries(fA,fX);
132 }
133 
134 /** You must ensure fA>0.0 && fX>0.0) */
GetUpRegIGamma(double fA,double fX)135 double ScInterpreter::GetUpRegIGamma( double fA, double fX )
136 {
137     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetUpRegIGamma" );
138 
139     double fLnFactor= fA*log(fX)-fX-GetLogGamma(fA);
140     double fFactor = exp(fLnFactor); //Do I need more accuracy than exp(ln()) has?;
141     if (fX>fA+1.0) // includes fX>1.0
142             return fFactor * GetGammaContFraction(fA,fX);
143     else //fX<=1 || fX<=fA+1, 1-GetLowRegIGamma, series
144             return 1.0 -fFactor * GetGammaSeries(fA,fX);
145 }
146 
147 /** Gamma distribution, probability density function.
148     fLambda is "scale" parameter
149     You must ensure fAlpha>0.0 and fLambda>0.0 */
GetGammaDistPDF(double fX,double fAlpha,double fLambda)150 double ScInterpreter::GetGammaDistPDF( double fX, double fAlpha, double fLambda )
151 {
152     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDistPDF" );
153     if (fX <= 0.0)
154         return 0.0;     // see ODFF
155     else
156     {
157         double fXr = fX / fLambda;
158         // use exp(ln()) only for large arguments because of less accuracy
159         if (fXr > 1.0)
160         {
161             const double fLogDblMax = log( ::std::numeric_limits<double>::max());
162             if (log(fXr) * (fAlpha-1.0) < fLogDblMax && fAlpha < fMaxGammaArgument)
163             {
164                 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
165             }
166             else
167             {
168                 return exp( (fAlpha-1.0) * log(fXr) - fXr - log(fLambda) - GetLogGamma(fAlpha));
169             }
170         }
171         else    // fXr near to zero
172         {
173             if (fAlpha<fMaxGammaArgument)
174             {
175                 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
176             }
177             else
178             {
179                 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / exp( GetLogGamma(fAlpha));
180             }
181         }
182     }
183 }
184 
185 /** Gamma distribution, cumulative distribution function.
186     fLambda is "scale" parameter
187     You must ensure fAlpha>0.0 and fLambda>0.0 */
GetGammaDist(double fX,double fAlpha,double fLambda)188 double ScInterpreter::GetGammaDist( double fX, double fAlpha, double fLambda )
189 {
190     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDist" );
191     if (fX <= 0.0)
192         return 0.0;
193     else
194         return GetLowRegIGamma( fAlpha, fX / fLambda);
195 }
196