Attachment #80273 for issue #121561

View | Details | Raw Unified | Return to issue 121561
Collapse All | Expand All





#include <boost/math/special_functions/atanh.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/log1p.hpp>

using ::std::vector;

    return fSumNum/fSumDenom;
}

// The algorithm is based on tgamma in gamma.hpp
// in math library from http://www.boost.org
/** You must ensure fZ>0; fZ>171.624376956302 will overflow. */
double lcl_GetGammaHelper(double fZ)
{
    double fGamma = lcl_getLanczosSum(fZ);
    const double fg = 6.024680040776729583740234375;
    double fZgHelp = fZ + fg - 0.5;
    // avoid intermediate overflow
    double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25);
    fGamma *= fHalfpower;
    fGamma /= exp(fZgHelp);
    fGamma *= fHalfpower;
    if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ))
        fGamma = ::rtl::math::round(fGamma);
    return fGamma;
}

// The algorithm is based on tgamma in gamma.hpp
// in math library from http://www.boost.org
/** You must ensure fZ>0 */
double lcl_GetLogGammaHelper(double fZ)
{
    const double fg = 6.024680040776729583740234375;
    double fZgHelp = fZ + fg - 0.5;
    return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;
}

/** You must ensure non integer arguments for fZ<1 */
double ScInterpreter::GetGamma(double fZ)
{
    RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "bst", "ScInterpreter::GetGamma" );
    const double fLogPi = log(F_PI);
    const double fLogDblMax = log( ::std::numeric_limits<double>::max());

    if (fZ > fMaxGammaArgument)
    {
        SetError(errIllegalFPOperation);
        return HUGE_VAL;
    }
    return ::boost::math::tgamma(fZ);
    if (fZ >= 1.0)
        return lcl_GetGammaHelper(fZ);

    if (fZ >= 0.5)  // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
        return lcl_GetGammaHelper(fZ+1) / fZ;

    if (fZ >= -0.5) // shift to x>=1, might overflow
    {
        double fLogTest = lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log( fabs(fZ));
        if (fLogTest >= fLogDblMax)
        {
            SetError( errIllegalFPOperation);
            return HUGE_VAL;
        }
        return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
    }
    // fZ<-0.5
    // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) )
    double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( fabs( ::rtl::math::sin( F_PI*fZ)));
    if (fLogDivisor - fLogPi >= fLogDblMax)     // underflow
        return 0.0;

    if (fLogDivisor<0.0)
        if (fLogPi - fLogDivisor > fLogDblMax)  // overflow
        {
            SetError(errIllegalFPOperation);
            return HUGE_VAL;
        }

    return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( F_PI*fZ) < 0.0) ? -1.0 : 1.0);
}


/** You must ensure fZ>0 */
double ScInterpreter::GetLogGamma(double fZ)
{
    RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "bst", "ScInterpreter::GetLogGamma" );
    return ::boost::math::lgamma(fZ);
        return lcl_GetLogGammaHelper(fZ);
    if (fZ >= 1.0)
        return log(lcl_GetGammaHelper(fZ));
    if (fZ >= 0.5)
        return log( lcl_GetGammaHelper(fZ+1) / fZ);
    return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);
}

double ScInterpreter::GetFDist(double x, double fF1, double fF2)

Return to issue 121561

Lines 46-51 Link Here

(-)sc/source/core/tool/interpr3.cxx (-70 / +5 lines)
46		46
47	#include <boost/math/special_functions/atanh.hpp>	47	#include <boost/math/special_functions/atanh.hpp>
48	#include <boost/math/special_functions/expm1.hpp>	48	#include <boost/math/special_functions/expm1.hpp>
		49	#include <boost/math/special_functions/gamma.hpp>
49	#include <boost/math/special_functions/log1p.hpp>	50	#include <boost/math/special_functions/log1p.hpp>
50		51
51	using ::std::vector;	52	using ::std::vector;
Lines 577-667 Link Here
577	return fSumNum/fSumDenom;	578	return fSumNum/fSumDenom;
578	}	579	}
579		580
580	// The algorithm is based on tgamma in gamma.hpp
581	// in math library from http://www.boost.org
582	/** You must ensure fZ>0; fZ>171.624376956302 will overflow. */
583	double lcl_GetGammaHelper(double fZ)
584	{
585	double fGamma = lcl_getLanczosSum(fZ);
586	const double fg = 6.024680040776729583740234375;
587	double fZgHelp = fZ + fg - 0.5;
588	// avoid intermediate overflow
589	double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25);
590	fGamma *= fHalfpower;
591	fGamma /= exp(fZgHelp);
592	fGamma *= fHalfpower;
593	if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ))
594	fGamma = ::rtl::math::round(fGamma);
595	return fGamma;
596	}
597		581
598	// The algorithm is based on tgamma in gamma.hpp
599	// in math library from http://www.boost.org
600	/** You must ensure fZ>0 */
601	double lcl_GetLogGammaHelper(double fZ)
602	{
603	const double fg = 6.024680040776729583740234375;
604	double fZgHelp = fZ + fg - 0.5;
605	return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;
606	}
607
608	/** You must ensure non integer arguments for fZ<1 */	582	/** You must ensure non integer arguments for fZ<1 */
609	double ScInterpreter::GetGamma(double fZ)	583	double ScInterpreter::GetGamma(double fZ)
610	{	584	{
611	RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGamma" );	585	RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "bst", "ScInterpreter::GetGamma" );
612	const double fLogPi = log(F_PI);
613	const double fLogDblMax = log( ::std::numeric_limits<double>::max());
614
615	if (fZ > fMaxGammaArgument)	586	if (fZ > fMaxGammaArgument)
616	{	587	{
617	SetError(errIllegalFPOperation);	588	SetError(errIllegalFPOperation);
618	return HUGE_VAL;	589	return HUGE_VAL;
619	}	590	}
620		591	return ::boost::math::tgamma(fZ);
621	if (fZ >= 1.0)
622	return lcl_GetGammaHelper(fZ);
623
624	if (fZ >= 0.5) // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
625	return lcl_GetGammaHelper(fZ+1) / fZ;
626
627	if (fZ >= -0.5) // shift to x>=1, might overflow
628	{
629	double fLogTest = lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log( fabs(fZ));
630	if (fLogTest >= fLogDblMax)
631	{
632	SetError( errIllegalFPOperation);
633	return HUGE_VAL;
634	}
635	return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
636	}
637	// fZ<-0.5
638	// Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)sin(pix) )
639	double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( fabs( ::rtl::math::sin( F_PI*fZ)));
640	if (fLogDivisor - fLogPi >= fLogDblMax) // underflow
641	return 0.0;
642
643	if (fLogDivisor<0.0)
644	if (fLogPi - fLogDivisor > fLogDblMax) // overflow
645	{
646	SetError(errIllegalFPOperation);
647	return HUGE_VAL;
648	}
649
650	return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( F_PI*fZ) < 0.0) ? -1.0 : 1.0);
651	}	592	}
652		593
653		594
654	/** You must ensure fZ>0 */	595	/** You must ensure fZ>0 */
655	double ScInterpreter::GetLogGamma(double fZ)	596	double ScInterpreter::GetLogGamma(double fZ)
656	{	597	{
657	RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetLogGamma" );	598	RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "bst", "ScInterpreter::GetLogGamma" );
658	if (fZ >= fMaxGammaArgument)	599	return ::boost::math::lgamma(fZ);
659	return lcl_GetLogGammaHelper(fZ);
660	if (fZ >= 1.0)
661	return log(lcl_GetGammaHelper(fZ));
662	if (fZ >= 0.5)
663	return log( lcl_GetGammaHelper(fZ+1) / fZ);
664	return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);
665	}	600	}
666		601
667	double ScInterpreter::GetFDist(double x, double fF1, double fF2)	602	double ScInterpreter::GetFDist(double x, double fF1, double fF2)