wmlmath.inl

Wild Math Library数值计算库

    fResult += (Real)0.999866;
    fResult *= fValue;
    return fResult;
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::FastInvTan1 (Real fValue)
{
    Real fVSqr = fValue*fValue;
    Real fResult = (Real)0.0028662257;
    fResult *= fVSqr;
    fResult -= (Real)0.0161657367;
    fResult *= fVSqr;
    fResult += (Real)0.0429096138;
    fResult *= fVSqr;
    fResult -= (Real)0.0752896400;
    fResult *= fVSqr;
    fResult += (Real)0.1065626393;
    fResult *= fVSqr;
    fResult -= (Real)0.1420889944;
    fResult *= fVSqr;
    fResult += (Real)0.1999355085;
    fResult *= fVSqr;
    fResult -= (Real)0.3333314528;
    fResult *= fVSqr;
    fResult += (Real)1.0;
    fResult *= fValue;
    return fResult;
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::FastInvSqrt (Real fValue)
{
    // TO DO.  This routine was designed for 'float'.  Come up with an
    // equivalent one for 'double' and specialize the templates for 'float'
    // and 'double'.
    float fFValue = (float)fValue;
    float fHalf = 0.5f*fFValue;
    int i  = *(int*)&fFValue;
    i = 0x5f3759df - (i >> 1);
    fFValue = *(float*)&i;
    fFValue = fFValue*(1.5f - fHalf*fFValue*fFValue);
    return (Real)fFValue;
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::LogGamma (Real fX)
{
    const Real afCoeff[6] =
    {
        +(Real)76.18009173,
        -(Real)86.50532033,
        +(Real)24.01409822,
        -(Real)1.231739516,
        +(Real)0.120858003e-02,
        -(Real)0.536382000e-05
    };

    fX -= (Real)1.0;
    Real fTmp = fX + (Real)5.5;
    fTmp -= (fX+((Real)0.5))*Math::Log(fTmp);
    Real fSeries = (Real)1.0;
    for (int j = 0; j <= 5; j++)
    {
        fX += (Real)1.0;
        fSeries += afCoeff[j]/fX;
    }
    return -fTmp + Math::Log(((Real)2.50662827465)*fSeries);
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::Gamma (Real fX)
{
    return Math::Exp(LogGamma(fX));
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::IncompleteGammaS (Real fA, Real fX)
{
    const int iMaxIterations = 100;
    const Real fTolerance = (Real)3e-07;

    if ( fX > (Real)0.0 )
    {
        Real fAp = fA;
        Real fSum = ((Real)1.0)/fA, fDel = fSum;
        for (int i = 1; i <= iMaxIterations; i++)
        {
            fAp += (Real)1.0;
            fDel *= fX/fAp;
            fSum += fDel;
            if ( Math::FAbs(fDel) < Math::FAbs(fSum)*fTolerance )
            {
                Real fArg = -fX+fA*Math::Log(fX)-LogGamma(fA);
                return fSum*Math::Exp(fArg);
            }
        }
    }

    if ( fX == (Real)0.0 )
        return (Real)0.0;

    return Math::MAX_REAL; // LogGamma not defined for x < 0
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::IncompleteGammaCF (Real fA, Real fX)
{
    const int iMaxIterations = 100;
    const Real fTolerance = (Real)3e-07;

    Real fA0 = (Real)1.0, fA1 = fX;
    Real fB0 = 0, fB1 = (Real)1.0;
    Real fGold = (Real)0.0, fFac = (Real)1.0;

    for (int i = 1; i <= iMaxIterations; i++)
    {
        Real fI = (Real) i;
        Real fImA = fI - fA;
        fA0 = (fA1 + fA0*fImA)*fFac;
        fB0 = (fB1 + fB0*fImA)*fFac;
        Real fItF = fI*fFac;
        fA1 = fX*fA0 + fItF*fA1;
        fB1 = fX*fB0 + fItF*fB1;
        if ( fA1 != (Real)0.0 )
        {
            fFac = ((Real)1.0)/fA1;
            Real fG = fB1*fFac;
            if ( Math::FAbs((fG-fGold)/fG) < fTolerance)
            {
                Real fArg = -fX + fA*Math::Log(fX) - LogGamma(fA);
                return fG*Math::Exp(fArg);
            }
            fGold = fG;
        }
    }

    return Math::MAX_REAL;  // numerical error if you get here
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::IncompleteGamma (Real fA, Real fX)
{
    if ( fX < (Real)1.0 + fA )
        return IncompleteGammaS(fA,fX);
    else
        return (Real)1.0-IncompleteGammaCF(fA,fX);
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::Erf (Real fX)
{
    return (Real)1.0-Erfc(fX);
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::Erfc (Real fX)
{
    const Real afCoeff[10] =
    {
        -(Real)1.26551223,
        +(Real)1.00002368,
        +(Real)0.37409196,
        +(Real)0.09678418,
        -(Real)0.18628806,
        +(Real)0.27886807,
        -(Real)1.13520398,
        +(Real)1.48851587,
        -(Real)0.82215223,
        +(Real)0.17087277
    };

    Real fZ = Math::FAbs(fX);
    Real fT = ((Real)1.0)/((Real)1.0+((Real)0.5)*fZ);
    Real fSum = afCoeff[9];

    for (int i = 9; i >= 0; i--)
        fSum = fT*fSum + afCoeff[i];

    Real fResult = fT*Math::Exp(-fZ*fZ + fSum);

    return ( fX >= (Real)0.0 ? fResult : (Real)2.0 - fResult );
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::ModBessel0 (Real fX)
{
    if ( fX < (Real)0.0 )  // function is even
        fX = -fX;

    Real fT, fResult;
    int i;

    if ( fX <= (Real)3.75 )
    {
        const Real afCoeff[7] =
        {
            (Real)1.0000000,
            (Real)3.5156229,
            (Real)3.0899424,
            (Real)1.2067492,
            (Real)0.2659732,
            (Real)0.0360768,
            (Real)0.0045813
        };

        fT = fX/(Real)3.75;
        Real fT2 = fT*fT;
        fResult = afCoeff[6];
        for (i = 5; i >= 0; i--)
        {
            fResult *= fT2;
            fResult += afCoeff[i];
        }
        // |error| < 1.6e-07
    }
    else
    {
        const Real afCoeff[9] =
        {
            +(Real)0.39894228,
            +(Real)0.01328592,
            +(Real)0.00225319,
            -(Real)0.00157565,
            +(Real)0.00916281,
            -(Real)0.02057706,
            +(Real)0.02635537,
            -(Real)0.01647633,
            +(Real)0.00392377
        };

        fT = fX/(Real)3.75;
        Real fInvT = ((Real)1.0)/fT;
        fResult = afCoeff[8];
        for (i = 7; i >= 0; i--)
        {
            fResult *= fInvT;
            fResult += afCoeff[i];
        }
        fResult *= Math::Exp(fX);
        fResult /= Math::Sqrt(fX);
        // |error| < 1.9e-07
    }

    return fResult;
}
//----------------------------------------------------------------------------
template <class Real>
Real Math<Real>::ModBessel1 (Real fX)
{
    int iSign;
    if ( fX > (Real)0.0 )
    {
        iSign = 1;
    }
    else if ( fX < (Real)0.0 )
    {
        fX = -fX;
        iSign = -1;
    }
    else
    {
        return (Real)0.0;
    }

    Real fT, fResult;
    int i;

    if ( fX <= (Real)3.75 )
    {
        const Real afCoeff[7] =
        {
            (Real)0.50000000,
            (Real)0.87890549,
            (Real)0.51498869,
            (Real)0.15084934,
            (Real)0.02658733,
            (Real)0.00301532,
            (Real)0.00032411
        };

        fT = fX/(Real)3.75;
        Real fT2 = fT*fT;
        fResult = afCoeff[6];
        for (i = 5; i >= 0; i--)
        {
            fResult *= fT2;
            fResult += afCoeff[i];
        }
        fResult *= fX;
        // |error| < 8e-09
    }
    else
    {
        const Real afCoeff[9] =
        {
            +(Real)0.39894228,
            -(Real)0.03988024,
            -(Real)0.00362018,
            +(Real)0.00163801,
            -(Real)0.01031555,
            +(Real)0.02282967,
            -(Real)0.02895312,
            +(Real)0.01787654,
            -(Real)0.00420059
        };

        fT = fX/(Real)3.75;
        Real fInvT = ((Real)1.0)/fT;
        fResult = afCoeff[8];
        for (i = 7; i >= 0; i--)
        {
            fResult *= fInvT;
            fResult += afCoeff[i];
        }
        fResult *= Math::Exp(fX);
        fResult /= Math::Sqrt(fX);
        // |error| < 2.2e-07
    }

    fResult *= iSign;
    return fResult;
}
//----------------------------------------------------------------------------