sc/source/core/tool/interpr6.cxx

   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: interpr6.cxx,v $
  10  * $Revision: 1.8 $
  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 <math.h>
  35
  36 #include <tools/debug.hxx>
  37 #include <rtl/logfile.hxx>
  38 #include "interpre.hxx"
  39
  40 double const fHalfMachEps = 0.5 * ::std::numeric_limits<double>::epsilon();
  41
  42 // The idea how this group of gamma functions is calculated, is
  43 // based on the Cephes library
  44 // online http://www.moshier.net/#Cephes [called 2008-02]
  45
  46 /** You must ensure fA>0.0 && fX>0.0
  47     valid results only if fX > fA+1.0
  48     uses continued fraction with odd items */
  49 double ScInterpreter::GetGammaContFraction( double fA, double fX )
  50 {
  51     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaContFraction" );
  52
  53     double const fBigInv = ::std::numeric_limits<double>::epsilon();
  54     double const fBig = 1.0/fBigInv;
  55     double fCount = 0.0;
  56     double fNum = 0.0;  // dummy value
  57     double fY = 1.0 - fA;
  58     double fDenom = fX + 2.0-fA;
  59     double fPk = 0.0;   // dummy value
  60     double fPkm1 = fX + 1.0;
  61     double fPkm2 = 1.0;
  62     double fQk = 1.0;   // dummy value
  63     double fQkm1 = fDenom * fX;
  64     double fQkm2 = fX;
  65     double fApprox = fPkm1/fQkm1;
  66     bool bFinished = false;
  67     double fR = 0.0;    // dummy value
  68     do
  69     {
  70         fCount = fCount +1.0;
  71         fY = fY+ 1.0;
  72         fNum = fY * fCount;
  73         fDenom = fDenom +2.0;
  74         fPk = fPkm1 * fDenom  -  fPkm2 * fNum;
  75         fQk = fQkm1 * fDenom  -  fQkm2 * fNum;
  76         if (fQk != 0.0)
  77         {
  78             fR = fPk/fQk;
  79             bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps);
  80             fApprox = fR;
  81         }
  82         fPkm2 = fPkm1;
  83         fPkm1 = fPk;
  84         fQkm2 = fQkm1;
  85         fQkm1 = fQk;
  86         if (fabs(fPk) > fBig)
  87         {
  88             // reduce a fraction does not change the value
  89             fPkm2 = fPkm2 * fBigInv;
  90             fPkm1 = fPkm1 * fBigInv;
  91             fQkm2 = fQkm2 * fBigInv;
  92             fQkm1 = fQkm1 * fBigInv;
  93         }
  94     } while (!bFinished && fCount<10000);
  95     // most iterations, if fX==fAlpha+1.0; approx sqrt(fAlpha) iterations then
  96     if (!bFinished)
  97     {
  98         SetError(errNoConvergence);
  99     }
 100     return fApprox;
 101 }
 102
 103 /** You must ensure fA>0.0 && fX>0.0
 104     valid results only if fX <= fA+1.0
 105     uses power series */
 106 double ScInterpreter::GetGammaSeries( double fA, double fX )
 107 {
 108     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaSeries" );
 109     double fDenomfactor = fA;
 110     double fSummand = 1.0/fA;
 111     double fSum = fSummand;
 112     int nCount=1;
 113     do
 114     {
 115         fDenomfactor = fDenomfactor + 1.0;
 116         fSummand = fSummand * fX/fDenomfactor;
 117         fSum = fSum + fSummand;
 118         nCount = nCount+1;
 119     } while ( fSummand/fSum > fHalfMachEps && nCount<=10000);
 120     // large amount of iterations will be carried out for huge fAlpha, even
 121     // if fX <= fAlpha+1.0
 122     if (nCount>10000)
 123     {
 124         SetError(errNoConvergence);
 125     }
 126     return fSum;
 127 }
 128
 129 /** You must ensure fA>0.0 && fX>0.0) */
 130 double ScInterpreter::GetLowRegIGamma( double fA, double fX )
 131 {
 132     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetLowRegIGamma" );
 133     double fLnFactor = fA * log(fX) - fX - GetLogGamma(fA);
 134     double fFactor = exp(fLnFactor);    // Do we need more accuracy than exp(ln()) has?
 135     if (fX>fA+1.0)  // includes fX>1.0; 1-GetUpRegIGamma, continued fraction
 136         return 1.0 - fFactor * GetGammaContFraction(fA,fX);
 137     else            // fX<=1.0 || fX<=fA+1.0, series
 138         return fFactor * GetGammaSeries(fA,fX);
 139 }
 140
 141 /** You must ensure fA>0.0 && fX>0.0) */
 142 double ScInterpreter::GetUpRegIGamma( double fA, double fX )
 143 {
 144     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetUpRegIGamma" );
 145
 146     double fLnFactor= fA*log(fX)-fX-GetLogGamma(fA);
 147     double fFactor = exp(fLnFactor); //Do I need more accuracy than exp(ln()) has?;
 148     if (fX>fA+1.0) // includes fX>1.0
 149             return fFactor * GetGammaContFraction(fA,fX);
 150     else //fX<=1 || fX<=fA+1, 1-GetLowRegIGamma, series
 151             return 1.0 -fFactor * GetGammaSeries(fA,fX);
 152 }
 153
 154 /** Gamma distribution, probability density function.
 155     fLambda is "scale" parameter
 156     You must ensure fAlpha>0.0 and fLambda>0.0 */
 157 double ScInterpreter::GetGammaDistPDF( double fX, double fAlpha, double fLambda )
 158 {
 159     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDistPDF" );
 160     if (fX <= 0.0)
 161         return 0.0;     // see ODFF
 162     else
 163     {
 164         double fXr = fX / fLambda;
 165         // use exp(ln()) only for large arguments because of less accuracy
 166         if (fXr > 1.0)
 167         {
 168             const double fLogDblMax = log( ::std::numeric_limits<double>::max());
 169             if (log(fXr) * (fAlpha-1.0) < fLogDblMax && fAlpha < fMaxGammaArgument)
 170             {
 171                 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
 172             }
 173             else
 174             {
 175                 return exp( (fAlpha-1.0) * log(fXr) - fXr - log(fLambda) - GetLogGamma(fAlpha));
 176             }
 177         }
 178         else    // fXr near to zero
 179         {
 180             if (fAlpha<fMaxGammaArgument)
 181             {
 182                 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
 183             }
 184             else
 185             {
 186                 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / exp( GetLogGamma(fAlpha));
 187             }
 188         }
 189     }
 190 }
 191
 192 /** Gamma distribution, cumulative distribution function.
 193     fLambda is "scale" parameter
 194     You must ensure fAlpha>0.0 and fLambda>0.0 */
 195 double ScInterpreter::GetGammaDist( double fX, double fAlpha, double fLambda )
 196 {
 197     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDist" );
 198     if (fX <= 0.0)
 199         return 0.0;
 200     else
 201         return GetLowRegIGamma( fAlpha, fX / fLambda);
 202 }