// 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 "SM_MathPCH.h"
#include "MgcRTLib.h"
#include "MgcPolynomial.h"

int MgcPolynomial::DIGITS_ACCURACY = 6;
float MgcPolynomial::ZERO_TOLERANCE = 1e-06f;
const float MgcPolynomial::ms_fInvLog2 = 1.0f/logf(2.0f);
const float MgcPolynomial::ms_fLog10 = logf(10.0f);
const float MgcPolynomial::ms_fThird = 1.0f/3.0f;
const float MgcPolynomial::ms_fSqrt3 = sqrtf(3.0f);
const float MgcPolynomial::ms_fTwentySeventh = 1.0f/27.0f;

//----------------------------------------------------------------------------
MgcPolynomial::MgcPolynomial (int iDegree)
{
    if ( iDegree >= 0 )
    {
        m_iDegree = iDegree;
        m_afCoeff = new float[m_iDegree+1];
    }
    else
    {
        // default creation
        m_iDegree = -1;
        m_afCoeff = 0;
    }
}
//----------------------------------------------------------------------------
MgcPolynomial::MgcPolynomial (const MgcPolynomial& rkPoly)
{
    m_iDegree = rkPoly.m_iDegree;
    m_afCoeff = new float[m_iDegree+1];
    for (int i = 0; i <= m_iDegree; i++)
        m_afCoeff[i] = rkPoly.m_afCoeff[i];
}
//----------------------------------------------------------------------------
MgcPolynomial::~MgcPolynomial ()
{
    delete[] m_afCoeff;
}
//----------------------------------------------------------------------------
float& MgcPolynomial::operator[] (int i) const
{
    assert( 0 <= i && i <= m_iDegree );
    return ((float*)m_afCoeff)[i];
}
//----------------------------------------------------------------------------
MgcPolynomial& MgcPolynomial::operator= (const MgcPolynomial& rkPoly)
{
    delete[] m_afCoeff;
    m_iDegree = rkPoly.m_iDegree;

    if ( m_iDegree >= 0 )
    {
        m_afCoeff = new float[m_iDegree+1];
        for (int i = 0; i <= m_iDegree; i++)
            m_afCoeff[i] = rkPoly.m_afCoeff[i];
    }

    return *this;
}
//----------------------------------------------------------------------------
void MgcPolynomial::SetDegree (int iDegree)
{
    m_iDegree = iDegree;
    delete[] m_afCoeff;

    if ( m_iDegree >= 0 )
        m_afCoeff = new float[m_iDegree+1];
    else
        m_afCoeff = 0;
}
//----------------------------------------------------------------------------
int MgcPolynomial::GetDegree () const
{
    return m_iDegree;
}
//----------------------------------------------------------------------------
float MgcPolynomial::operator() (float fT) const
{
    assert( m_iDegree >= 0 );

    float fResult = m_afCoeff[m_iDegree];
    for (int i = m_iDegree-1; i >= 0; i--)
    {
        fResult *= fT;
        fResult += m_afCoeff[i];
    }
    return fResult;
}
//----------------------------------------------------------------------------
MgcPolynomial MgcPolynomial::GetDerivative () const
{
    if ( m_iDegree > 0 )
    {
        MgcPolynomial kDeriv(m_iDegree-1);
        for (int i0 = 0, i1 = 1; i0 < m_iDegree; i0++, i1++)
            kDeriv[i0] = i1*m_afCoeff[i1];
        return kDeriv;
    }
    else
    {
        return MgcPolynomial(-1);
    }
}
//----------------------------------------------------------------------------
MgcPolynomial MgcPolynomial::GetInversion () const
{
    MgcPolynomial kInvPoly(m_iDegree);
    for (int i = 0; i <= m_iDegree; i++)
        kInvPoly[i] = m_afCoeff[m_iDegree-i];
    return kInvPoly;
}
//----------------------------------------------------------------------------
MgcPolynomial MgcPolynomial::operator+ (const MgcPolynomial& rkPoly) const
{
    assert( m_iDegree >= 0 && rkPoly.m_iDegree >= 0 );

    MgcPolynomial kSum(-1);
    int i;

    if ( m_iDegree > rkPoly.m_iDegree )
    {
        kSum.SetDegree(m_iDegree);
        for (i = 0; i <= rkPoly.m_iDegree; i++)
            kSum[i] = m_afCoeff[i] + rkPoly.m_afCoeff[i];
        for (i = rkPoly.m_iDegree+1; i <= m_iDegree; i++)
            kSum[i] = m_afCoeff[i];

    }
    else
    {
        kSum.SetDegree(rkPoly.m_iDegree);
        for (i = 0; i <= m_iDegree; i++)
            kSum[i] = m_afCoeff[i] + rkPoly.m_afCoeff[i];
        for (i = m_iDegree+1; i <= m_iDegree; i++)
            kSum[i] = rkPoly.m_afCoeff[i];
    }

    return kSum;
}
//----------------------------------------------------------------------------
MgcPolynomial MgcPolynomial::operator- (const MgcPolynomial& rkPoly) const
{
    assert( m_iDegree >= 0 && rkPoly.m_iDegree >= 0 );

    MgcPolynomial kDiff(-1);
    int i;

    if ( m_iDegree > rkPoly.m_iDegree )
    {
        kDiff.SetDegree(m_iDegree);
        for (i = 0; i <= rkPoly.m_iDegree; i++)
            kDiff[i] = m_afCoeff[i] - rkPoly.m_afCoeff[i];
        for (i = rkPoly.m_iDegree+1; i <= m_iDegree; i++)
            kDiff[i] = m_afCoeff[i];

    }
    else
    {
        kDiff.SetDegree(rkPoly.m_iDegree);
        for (i = 0; i <= m_iDegree; i++)
            kDiff[i] = m_afCoeff[i] - rkPoly.m_afCoeff[i];
        for (i = m_iDegree+1; i <= m_iDegree; i++)
            kDiff[i] = -rkPoly.m_afCoeff[i];
    }

    return kDiff;
}
//----------------------------------------------------------------------------
MgcPolynomial MgcPolynomial::operator* (const MgcPolynomial& rkPoly) const
{
    assert( m_iDegree >= 0 && rkPoly.m_iDegree >= 0 );

    MgcPolynomial kProd(m_iDegree + rkPoly.m_iDegree);

    memset(kProd.m_afCoeff,0,(kProd.m_iDegree+1)*sizeof(float));
    for (int i0 = 0; i0 <= m_iDegree; i0++)
    {
        for (int i1 = 0; i1 <= rkPoly.m_iDegree; i1++)
        {
            int i2 = i0 + i1;
            kProd.m_afCoeff[i2] += m_afCoeff[i0]*rkPoly.m_afCoeff[i1];
        }
    }

    return kProd;
}
//----------------------------------------------------------------------------
MgcPolynomial MgcPolynomial::operator* (float fScalar) const
{
    assert( m_iDegree >= 0 );

    MgcPolynomial kProd(m_iDegree);
    for (int i = 0; i <= m_iDegree; i++)
        kProd[i] = fScalar*m_afCoeff[i];

    return kProd;
}
//----------------------------------------------------------------------------
MgcPolynomial MgcPolynomial::operator- () const
{
    assert( m_iDegree >= 0 );

    MgcPolynomial kNeg(m_iDegree);
    for (int i = 0; i <= m_iDegree; i++)
        kNeg[i] = -m_afCoeff[i];

    return kNeg;
}
//----------------------------------------------------------------------------
MgcPolynomial operator* (float fScalar, const MgcPolynomial& rkPoly)
{
    assert( rkPoly.m_iDegree >= 0 );

    MgcPolynomial kProd(rkPoly.m_iDegree);
    for (int i = 0; i <= rkPoly.m_iDegree; i++)
        kProd[i] = fScalar*rkPoly.m_afCoeff[i];

    return kProd;
}
//----------------------------------------------------------------------------
bool MgcPolynomial::Bisection (float fXMin, float fXMax,
    float& rfRoot) const
{
    float fP0 = (*this)(fXMin);
    float fP1 = (*this)(fXMax);

    // check for endpoint zeros
    if ( fabsf(fP0) <= ZERO_TOLERANCE )
    {
        rfRoot = fXMin;
        return true;
    }
    if ( fabsf(fP1) <= ZERO_TOLERANCE )
    {
        rfRoot = fXMax;
        return true;
    }

    if ( fP0*fP1 > 0.0f )
        return false;

    // determine number of iterations to get 'digits' accuracy.
    float fTmp0 = logf(fXMax-fXMin);
    float fTmp1 = DIGITS_ACCURACY*ms_fLog10;
    float fArg = (fTmp0 + fTmp1)*ms_fInvLog2;
    int iMaxIter = (int)ceilf(fArg);
    
    for (int i = 0; i < iMaxIter; i++)
    {
        rfRoot = 0.5f*(fXMin + fXMax);
        float fP = (*this)(rfRoot);
        if ( fabsf(fP) <= ZERO_TOLERANCE )
            return true;
        
        if ( fP*fP0 < 0.0f )
        {
            fXMax = rfRoot;
            fP1 = fP;
        }
        else
        {
            fXMin = rfRoot;
            fP0 = fP;
        }
    }

    return true;
}
//----------------------------------------------------------------------------
void MgcPolynomial::GetRootsOnInterval (float fXMin, float fXMax,
    int& riCount, float* afRoot) const
{
    float fRoot;

    if ( m_iDegree == 1 )
    {
        if ( Bisection(fXMin,fXMax,fRoot) )
        {
            riCount = 1;
            afRoot[0] = fRoot;
        }
        else
        {
            riCount = 0;
        }
        return;
    }

    // get roots of derivative polynomial
    MgcPolynomial kDeriv = GetDerivative();
    kDeriv.GetRootsOnInterval(fXMin,fXMax,riCount,afRoot);

    int i, iNewCount = 0;
    float* afNewRoot = new float[riCount+1];

    if ( riCount > 0 )
    {
        // find root on [xmin,root[0]]
        if ( Bisection(fXMin,afRoot[0],fRoot) )
            afNewRoot[iNewCount++] = fRoot;

        // find root on [root[i],root[i+1]] for 0 <= i <= count-2
        for (i = 0; i <= riCount-2; i++)
        {
            if ( Bisection(afRoot[i],afRoot[i+1],fRoot) )
                afNewRoot[iNewCount++] = fRoot;
        }

        // find root on [root[count-1],xmax]
        if ( Bisection(afRoot[riCount-1],fXMax,fRoot) )
            afNewRoot[iNewCount++] = fRoot;
    }
    else
    {
        // polynomial is monotone on [xmin,xmax], has at most one root
        if ( Bisection(fXMin,fXMax,fRoot) )
            afNewRoot[iNewCount++] = fRoot;
    }

    // copy to old buffer
    if ( iNewCount > 0 )
    {
        riCount = 1;
        afRoot[0] = afNewRoot[0];
        for (i = 1; i < iNewCount; i++)
        {
            if ( fabsf(afNewRoot[i]-afNewRoot[i-1]) > ZERO_TOLERANCE )
                afRoot[riCount++] = afNewRoot[i];
        }
    }
    else
    {
        riCount = 0;
    }

    delete[] afNewRoot;
}
//----------------------------------------------------------------------------
void MgcPolynomial::GetAllRoots (int& riCount, float* afRoot) const
{
    float fBound, fTmp;
    int i;

    float fAbs0 = fabsf(m_afCoeff[0]);
    float fAbsD = fabsf(m_afCoeff[m_iDegree]);

    if ( fAbsD >= fAbs0 )
    {
        fBound = fAbs0;
        for (i = 0; i < m_iDegree; i++)
        {
            fTmp = fabsf(m_afCoeff[i]) + fAbsD;
            if ( fBound < fTmp )
                fBound = fTmp;
        }
        fBound /= fAbsD;

        GetRootsOnInterval(-fBound,fBound,riCount,afRoot);
    }
    else
    {
        fBound = fAbsD;
        for (i = 0; i < m_iDegree; i++)
        {
            fTmp = fabsf(m_afCoeff[i]) + fAbs0;
            if ( fBound < fTmp )
                fBound = fTmp;
        }
        fBound /= fAbs0;

        // use inverted polynomial to find inverse roots
        float* afInvRoot = new float[m_iDegree];
        MgcPolynomial kInvPoly = GetInversion();
        kInvPoly.GetRootsOnInterval(-fBound,fBound,riCount,afInvRoot);
        for (i = 0; i < riCount; i++)
            afRoot[i] = 1.0f/afInvRoot[riCount-1-i];
        delete[] afInvRoot;
    }
}
//----------------------------------------------------------------------------
float MgcPolynomial::GetRootBound () const
{
    float fBound, fTmp;
    int i;

    float fAbs0 = fabsf(m_afCoeff[0]);
    float fAbsD = fabsf(m_afCoeff[m_iDegree]);

    if ( fAbsD >= fAbs0 )
    {
        fBound = fAbs0;
        for (i = 0; i < m_iDegree; i++)
        {
            fTmp = fabsf(m_afCoeff[i]) + fAbsD;
            if ( fBound < fTmp )
                fBound = fTmp;
        }
        fBound /= fAbsD;
    }
    else
    {
        fBound = fAbsD;
        for (i = 0; i < m_iDegree; i++)
        {
            fTmp = fabsf(m_afCoeff[i]) + fAbs0;
            if ( fBound < fTmp )
                fBound = fTmp;
        }
        fBound /= fAbs0;
    }

    return fBound;
}
//----------------------------------------------------------------------------
bool MgcPolynomial::AllRealPartsNegative (int iDegree, float* afCoeff)
{
    // assert:  afCoeff[iDegree] = 1

    if ( afCoeff[iDegree-1] <= 0.0f )
        return false;
    if ( iDegree == 1 )
        return true;

    float* afTmpCoeff = new float[iDegree];
    afTmpCoeff[0] = 2.0f*afCoeff[0]*afCoeff[iDegree-1];
    int i;
    for (i = 1; i <= iDegree-2; i++) 
    {
        afTmpCoeff[i] = afCoeff[iDegree-1]*afCoeff[i];
        if ( ((iDegree-i) % 2) == 0 )
            afTmpCoeff[i] -= afCoeff[i-1];
        afTmpCoeff[i] *= 2.0f;
    }
    afTmpCoeff[iDegree-1] = 2.0f*afCoeff[iDegree-1]*afCoeff[iDegree-1];

    int iNextDegree;
    for (iNextDegree = iDegree-1; iNextDegree >= 0; iNextDegree--)
    {
        if ( afTmpCoeff[iNextDegree] != 0.0f )
            break;
    }
    for (i = 0; i <= iNextDegree-1; i++)
        afCoeff[i] = afTmpCoeff[i]/afTmpCoeff[iNextDegree];
    delete[] afTmpCoeff;

    return AllRealPartsNegative(iNextDegree,afCoeff);
}
//----------------------------------------------------------------------------
bool MgcPolynomial::AllRealPartsNegative () const
{
    // make a copy of coefficients, later calls will change the copy
    float* afCoeff = new float[m_iDegree+1];
    memcpy(afCoeff,m_afCoeff,(m_iDegree+1)*sizeof(float));

    // make polynomial monic
    if ( afCoeff[m_iDegree] != 1.0f )
    {
        float fInv = 1.0f/afCoeff[m_iDegree];
        for (int i = 0; i < m_iDegree; i++)
            afCoeff[i] *= fInv;
        afCoeff[m_iDegree] = 1.0f;
    }

    return AllRealPartsNegative(m_iDegree,afCoeff);
}
//----------------------------------------------------------------------------
bool MgcPolynomial::AllRealPartsPositive () const
{
    // make a copy of coefficients, later calls will change the copy
    float* afCoeff = new float[m_iDegree+1];
    memcpy(afCoeff,m_afCoeff,(m_iDegree+1)*sizeof(float));

    // make polynomial monic
    int i;
    if ( afCoeff[m_iDegree] != 1.0f )
    {
        float fInv = 1.0f/afCoeff[m_iDegree];
        for (i = 0; i < m_iDegree; i++)
            afCoeff[i] *= fInv;
        afCoeff[m_iDegree] = 1.0f;
    }

    // reflect z -> -z
    int iSign = -1;
    for (i = m_iDegree-1; i >= 0; i--, iSign = -iSign)
        afCoeff[i] *= iSign;

    return AllRealPartsNegative(m_iDegree,afCoeff);
}
//----------------------------------------------------------------------------
bool MgcPolynomial::RootsDegree2 (int& riCount, float afRoot[2]) const
{
    // compute real roots to x^2+c[1]*x+c[0] = 0
    if ( m_iDegree != 2 )
        return false;

    // make polynomial monic
    float afCoeff[2] = { m_afCoeff[0], m_afCoeff[1] };
    if ( m_afCoeff[2] != 1.0f )
    {
        float fInv = 1.0f/m_afCoeff[2];
        afCoeff[0] *= fInv;
        afCoeff[1] *= fInv;
    }

    float fDiscr = afCoeff[1]*afCoeff[1]-4.0f*afCoeff[0];
    if ( fabsf(fDiscr) <= ZERO_TOLERANCE )
        fDiscr = 0.0f;

    if ( fDiscr >= 0.0f )
    {
        fDiscr = sqrtf(fDiscr);
        afRoot[0] = 0.5f*(-afCoeff[1]-fDiscr);
        afRoot[1] = 0.5f*(-afCoeff[1]+fDiscr);
        riCount = 2;
    }
    else
    {
        riCount = 0;
    }

    return true;
}
//----------------------------------------------------------------------------
bool MgcPolynomial::RootsDegree3 (int& riCount, float afRoot[3]) const
{
    // compute real roots to x^3+c[2]*x^2+c[1]*x+c[0] = 0
    if ( m_iDegree != 3 )
        return false;

    // make polynomial monic
    float afCoeff[3] = { m_afCoeff[0], m_afCoeff[1], m_afCoeff[2] };
    if ( m_afCoeff[3] != 1.0f )
    {
        float fInv = 1.0f/m_afCoeff[3];
        afCoeff[0] *= fInv;
        afCoeff[1] *= fInv;
        afCoeff[2] *= fInv;
    }

    // convert to y^3+a*y+b = 0 by x = y-c[2]/3 and
    float fA = ms_fThird*(3.0f*afCoeff[1]-afCoeff[2]*afCoeff[2]);
    float fB = ms_fTwentySeventh*(2.0f*afCoeff[2]*afCoeff[2]*afCoeff[2] -
        9.0f*afCoeff[1]*afCoeff[2]+27.0f*afCoeff[0]);
    float fOffset = ms_fThird*afCoeff[2];

    float fDiscr = 0.25f*fB*fB + ms_fTwentySeventh*fA*fA*fA;
    if ( fabsf(fDiscr) <= ZERO_TOLERANCE )
        fDiscr = 0.0f;

    float fHalfB = 0.5f*fB;
    if ( fDiscr > 0.0f )  // 1 real, 2 complex roots
    {
        fDiscr = sqrtf(fDiscr);
        float fTemp = -fHalfB + fDiscr;
        if ( fTemp >= 0.0f )
            afRoot[0] = powf(fTemp,ms_fThird);
        else
            afRoot[0] = -powf(-fTemp,ms_fThird);
        fTemp = -fHalfB - fDiscr;
        if ( fTemp >= 0.0f )
            afRoot[0] += powf(fTemp,ms_fThird);
        else
            afRoot[0] -= powf(-fTemp,ms_fThird);
        afRoot[0] -= fOffset;
        riCount = 1;
    }
    else if ( fDiscr < 0.0f ) 
    {
        float fDist = sqrtf(-ms_fThird*fA);
        float fAngle = ms_fThird*atan2f(sqrtf(-fDiscr),
            -fHalfB);
        float fCos = cosf(fAngle);
        float fSin = sinf(fAngle);
        afRoot[0] = 2.0f*fDist*fCos-fOffset;
        afRoot[1] = -fDist*(fCos+ms_fSqrt3*fSin)-fOffset;
        afRoot[2] = -fDist*(fCos-ms_fSqrt3*fSin)-fOffset;
        riCount = 3;
    }
    else 
    {
        float fTemp;
        if ( fHalfB >= 0.0f )
            fTemp = -powf(fHalfB,ms_fThird);
        else
            fTemp = powf(-fHalfB,ms_fThird);
        afRoot[0] = 2.0f*fTemp-fOffset;
        afRoot[1] = -fTemp-fOffset;
        afRoot[2] = afRoot[1];
        riCount = 3;
    }

    return true;
}
//----------------------------------------------------------------------------
bool MgcPolynomial::RootsDegree4 (int& riCount, float afRoot[4]) const
{
    // compute real roots to x^4+c[3]*x^3+c[2]*x^2+c[1]*x+c[0] = 0
    if ( m_iDegree != 4 )
        return false;

    // make polynomial monic
    float afCoeff[4] = { m_afCoeff[0], m_afCoeff[1], m_afCoeff[2],
        m_afCoeff[3] };
    if ( m_afCoeff[4] != 1.0f )
    {
        float fInv = 1.0f/m_afCoeff[4];
        afCoeff[0] *= fInv;
        afCoeff[1] *= fInv;
        afCoeff[2] *= fInv;
        afCoeff[3] *= fInv;
    }

    // reduction to resolvent cubic polynomial
    MgcPolynomial kResolve(3);
    kResolve[3] = 1.0f;
    kResolve[2] = -afCoeff[2];
    kResolve[1] = afCoeff[3]*afCoeff[1]-4.0f*afCoeff[0];
    kResolve[0] = -afCoeff[3]*afCoeff[3]*afCoeff[0] +
        4.0f*afCoeff[2]*afCoeff[0]-afCoeff[1]*afCoeff[1];
    int iResolveCount;
    float afResolveRoot[3];
    kResolve.RootsDegree3(iResolveCount,afResolveRoot);
    float fY = afResolveRoot[0];

    riCount = 0;
    float fDiscr = 0.25f*afCoeff[3]*afCoeff[3]-afCoeff[2]+fY;
    if ( fabsf(fDiscr) <= ZERO_TOLERANCE )
        fDiscr = 0.0f;

    if ( fDiscr > 0.0f ) 
    {
        float fR = sqrtf(fDiscr);
        float fT1 = 0.75f*afCoeff[3]*afCoeff[3]-fR*fR-2.0f*afCoeff[2];
        float fT2 = (4.0f*afCoeff[3]*afCoeff[2]-8.0f*afCoeff[1]-
            afCoeff[3]*afCoeff[3]*afCoeff[3])/(4.0f*fR);

        float fTplus = fT1+fT2;
        float fTminus = fT1-fT2;
        if ( fabsf(fTplus) <= ZERO_TOLERANCE ) 
            fTplus = 0.0f;
        if ( fabsf(fTminus) <= ZERO_TOLERANCE ) 
            fTminus = 0.0f;

        if ( fTplus >= 0.0f )
        {
            float fD = sqrtf(fTplus);
            afRoot[0] = -0.25f*afCoeff[3]+0.5f*(fR+fD);
            afRoot[1] = -0.25f*afCoeff[3]+0.5f*(fR-fD);
            riCount += 2;
        }
        if ( fTminus >= 0.0f )
        {
            float fE = sqrtf(fTminus);
            afRoot[riCount++] = -0.25f*afCoeff[3]+0.5f*(fE-fR);
            afRoot[riCount++] = -0.25f*afCoeff[3]-0.5f*(fE+fR);
        }
    }
    else if ( fDiscr < 0.0f )
    {
        riCount = 0;
    }
    else
    {
        float fT2 = fY*fY-4.0f*afCoeff[0];
        if ( fT2 >= -ZERO_TOLERANCE ) 
        {
            if ( fT2 < 0.0f ) // round to zero
                fT2 = 0.0f;
            fT2 = 2.0f*sqrtf(fT2);
            float fT1 = 0.75f*afCoeff[3]*afCoeff[3]-2.0f*afCoeff[2];
            if ( fT1+fT2 >= ZERO_TOLERANCE ) 
            {
                float fD = sqrtf(fT1+fT2);
                afRoot[0] = -0.25f*afCoeff[3]+0.5f*fD;
                afRoot[1] = -0.25f*afCoeff[3]-0.5f*fD;
                riCount += 2;
            }
            if ( fT1-fT2 >= ZERO_TOLERANCE ) 
            {
                float fE = sqrtf(fT1-fT2);
                afRoot[riCount++] = -0.25f*afCoeff[3]+0.5f*fE;
                afRoot[riCount++] = -0.25f*afCoeff[3]-0.5f*fE;
            }
        }
    }

    return true;
}
//----------------------------------------------------------------------------
float MgcPolynomial::SpecialCubeRoot (float fA, float fB, float fC)
{
    // Solve A*r^3 + B*r = C where A > 0 and B > 0.
    //
    // Let r = D*sinh(u) where D = sqrt(4*B/(3*A)).  Then
    // sinh(3*u) = 4*[sinh(u)]^3+3*sinh(u) = E where E = 4*C/(A*D^3).
    // sinh(3*u) = E has solution u = (1/3)*log(E+sqrt(E^2+1)).  This
    // leads to sinh(u) = ((E+sqrt(E^2+1))^{1/3}-(E+sqrt(E^2+1))^{-1/3})/2.
    // Therefore,
    //
    //     r = D*((E+sqrt(E^2+1))^{1/3}-(E+sqrt(E^2+1))^{-1/3})/2.

    float fD = sqrtf(4.0f*ms_fThird*fB/fA);
    float fE = 4.0f*fC/(fA*fD*fD*fD);
    float fF = powf(fE+sqrtf(fE*fE+1.0f),ms_fThird);
    float fRoot = 0.5f*fD*(fF-1.0f/fF);

    return fRoot;
}
//----------------------------------------------------------------------------
