mgcspecialfunction.cpp

3D Game Engine Design Source Code非常棒

// Magic Software, Inc.
// http://www.magic-software.com
// Copyright (c) 2000, All Rights Reserved
//
// Source code from Magic Software is supplied under the terms of a license
// agreement and may not be copied or disclosed except in accordance with the
// terms of that agreement.  The various license agreements may be found at
// the Magic Software web site.  This file is subject to the license
//
// FREE SOURCE CODE
// http://www.magic-software.com/License/free.pdf

#include "MgcSpecialFunction.h"

//---------------------------------------------------------------------------
MgcReal MgcSpecialFunction::LogGamma (MgcReal fX)
{
    static const MgcReal s_afCoeff[6] =
    {
        76.18009173, -86.50532033, 24.01409822, -1.231739516,
        0.120858003e-2, -0.536382e-5
    };

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

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

    if ( fX == 0.0 )
        return 0.0;

    return MgcMath::INFINITY; // LogGamma not defined for x < 0
}
//---------------------------------------------------------------------------
MgcReal MgcSpecialFunction::IncompleteGammaCF (MgcReal fA, MgcReal fX)
{
    const int iMaxIterations = 100;
    const MgcReal fTolerance = 3e-07;

    MgcReal fA0 = 1.0, fA1 = fX;
    MgcReal fB0 = 0, fB1 = 1.0;
    MgcReal fGold = 0.0, fFac = 1.0;

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

    return MgcMath::INFINITY;  // numerical error if you get here
}
//---------------------------------------------------------------------------
MgcReal MgcSpecialFunction::IncompleteGamma (MgcReal fA, MgcReal fX)
{
    if ( fX < 1.0 + fA )
        return IncompleteGammaS(fA,fX);
    else
        return 1-IncompleteGammaCF(fA,fX);
}
//---------------------------------------------------------------------------
MgcReal MgcSpecialFunction::Erf (MgcReal fX)
{
    return 1-Erfc(fX);
}
//---------------------------------------------------------------------------
MgcReal MgcSpecialFunction::Erfc (MgcReal fX)
{
    static const MgcReal s_afCoeff[10] =
    {
        -1.26551223,  1.00002368, 0.37409196,  0.09678418, -0.18628806,
         0.27886807, -1.13520398, 1.48851587, -0.82215223,  0.17087277
    };

    MgcReal fZ = MgcMath::Abs(fX);
    MgcReal fT = 1.0/(1.0+0.5*fZ);
    MgcReal fSum = s_afCoeff[9];

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

    MgcReal fResult = fT*MgcMath::Exp(-fZ*fZ + fSum);

    return fX >= 0.0 ? fResult : 2.0 - fResult;
}
//---------------------------------------------------------------------------
MgcReal MgcSpecialFunction::ModBessel0 (MgcReal fX)
{
    if ( fX < 0.0 )  // function is even
        fX = -fX;

    MgcReal fT, fResult;
    int i;

    if ( fX <= 3.75 )
    {
        static const MgcReal s_afCoeff[7] =
        {
            +1.0000000, +3.5156229, +3.0899424, +1.2067492, +0.2659732,
            +0.0360768, +0.0045813
        };

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

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

    return fResult;
}
//---------------------------------------------------------------------------
MgcReal MgcSpecialFunction::ModBessel1 (MgcReal fX)
{
    int iSign;
    if ( fX > 0.0 )
    {
        iSign = 1;
    }
    else if ( fX < 0.0 )
    {
        fX = -fX;
        iSign = -1;
    }
    else
    {
        return 0.0;
    }

    MgcReal fT, fResult;
    int i;

    if ( fX <= 3.75 )
    {
        static const MgcReal s_afCoeff[7] =
        {
            +0.50000000, +0.87890549, +0.51498869, +0.15084934,
            +0.02658733, +0.00301532, +0.00032411
        };

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

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

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