wmlpolynomialrootsdegn.inl

Wild Math Library数值计算库

// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2004.  All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.

//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::FindB (const Polynomial1<Real>& rkPoly,
    int iDigits)
{
    Real fBound = GetBound(rkPoly);
    return FindB(rkPoly,-fBound,fBound,iDigits);
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::FindN (const Polynomial1<Real>&, int)
{
    // TO DO:  Implement this.
    assert( false );
    return false;
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::FindE (const Polynomial1<Real>&, bool)
{
    // TO DO:  Implement this.
    assert( false );
    return false;
}
//----------------------------------------------------------------------------
template <class Real>
Real PolynomialRoots<Real>::GetBound (const Polynomial1<Real>& rkPoly)
{
    Polynomial1<Real> kCPoly = rkPoly;
    kCPoly.Compress(m_fEpsilon);
    int iDegree = kCPoly.GetDegree();
    if ( iDegree < 1 )
    {
        // polynomial is constant, return invalid bound
        return -(Real)1.0;
    }

    Real fInvCDeg = ((Real)1.0)/kCPoly[iDegree];
    Real fMax = (Real)0.0;
    for (int i = 0; i < iDegree; i++)
    {
        Real fTmp = Math<Real>::FAbs(kCPoly[i])*fInvCDeg;
        if ( fTmp > fMax )
            fMax = fTmp;
    }

    return (Real)1.0 + fMax;
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::FindB (const Polynomial1<Real>& rkPoly,
    Real fXMin, Real fXMax, int iDigits)
{
    // reallocate root array if necessary
    if ( rkPoly.GetDegree() > m_iMaxRoot )
    {
        m_iMaxRoot = rkPoly.GetDegree();
        delete[] m_afRoot;
        m_afRoot = new Real[m_iMaxRoot];
    }

    Real fRoot;
    if ( rkPoly.GetDegree() == 1 )
    {
        if ( Bisection(rkPoly,fXMin,fXMax,iDigits,fRoot) )
        {
            m_iCount = 1;
            m_afRoot[0] = fRoot;
            return true;
        }
        m_iCount = 0;
        return false;
    }

    // get roots of derivative polynomial
    Polynomial1<Real> kDeriv = rkPoly.GetDerivative();
    FindB(kDeriv,fXMin,fXMax,iDigits);

    int i, iNewCount = 0;
    Real* afNewRoot = new Real[m_iCount+1];

    if ( m_iCount > 0 )
    {
        // find root on [xmin,root[0]]
        if ( Bisection(rkPoly,fXMin,m_afRoot[0],iDigits,fRoot) )
            afNewRoot[iNewCount++] = fRoot;

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

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

    // copy to old buffer
    if ( iNewCount > 0 )
    {
        m_iCount = 1;
        m_afRoot[0] = afNewRoot[0];
        for (i = 1; i < iNewCount; i++)
        {
            Real fRootDiff = afNewRoot[i]-afNewRoot[i-1];
            if ( Math<Real>::FAbs(fRootDiff) > m_fEpsilon )
                m_afRoot[m_iCount++] = afNewRoot[i];
        }
    }
    else
    {
        m_iCount = 0;
    }

    delete[] afNewRoot;
    return m_iCount > 0;
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::FindN (const Polynomial1<Real>&, Real, Real, int)
{
    // TO DO:  Implement this.
    assert( false );
    return false;
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::AllRealPartsNegative (
    const Polynomial1<Real>& rkPoly)
{
    // make a copy of coefficients, later calls will change the copy
    int iDegree = rkPoly.GetDegree();
    const Real* afPolyCoeff = (const Real*)rkPoly;
    Real* afCoeff = new Real[iDegree+1];
    memcpy(afCoeff,afPolyCoeff,(iDegree+1)*sizeof(Real));

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

    return AllRealPartsNegative(iDegree,afCoeff);
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::AllRealPartsPositive (
    const Polynomial1<Real>& rkPoly)
{
    // make a copy of coefficients, later calls will change the copy
    int iDegree = rkPoly.GetDegree();
    const Real* afPolyCoeff = (const Real*)rkPoly;
    Real* afCoeff = new Real[iDegree+1];
    memcpy(afCoeff,afPolyCoeff,(iDegree+1)*sizeof(Real));

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

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

    return AllRealPartsNegative(iDegree,afCoeff);
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::AllRealPartsNegative (int iDegree, Real* afCoeff)
{
    // assert:  afCoeff[iDegree] = 1

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

    Real* afTmpCoeff = new Real[iDegree];
    afTmpCoeff[0] = ((Real)2.0)*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] *= (Real)2.0;
    }
    afTmpCoeff[iDegree-1] =
        ((Real)2.0)*afCoeff[iDegree-1]*afCoeff[iDegree-1];

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

    return AllRealPartsNegative(iNextDegree,afCoeff);
}
//----------------------------------------------------------------------------
template <class Real>
int PolynomialRoots<Real>::GetRootCount (const Polynomial1<Real>& rkPoly,
    Real fT0, Real fT1)
{
    int iDegree = rkPoly.GetDegree();
    const Real* afCoeff = (const Real*)rkPoly;

    if ( iDegree == 0 )
    {
        // polynomial is constant on the interval
        if ( afCoeff[0] != (Real)0.0 )
            return 0;
        else
            return -1;  // to indicate "infinitely many"
    }

    // generate the Sturm sequence
    std::vector<Polynomial1<Real>*> kSturm;
    Polynomial1<Real>* pkF0 = new Polynomial1<Real>(rkPoly);
    Polynomial1<Real>* pkF1 = new Polynomial1<Real>(pkF0->GetDerivative());
    kSturm.push_back(pkF0);
    kSturm.push_back(pkF1);

    while ( pkF1->GetDegree() > 0 )
    {
        Polynomial1<Real>* pkF2 = new Polynomial1<Real>(-1);
        Polynomial1<Real> kQuot;
        pkF0->Divide(*pkF1,kQuot,*pkF2,Math<Real>::EPSILON);
        *pkF2 = -(*pkF2);
        kSturm.push_back(pkF2);
        pkF0 = pkF1;
        pkF1 = pkF2;
    }

    int i;
    Real fValue0, fValue1;

    // count the sign changes at t0
    int iSignChanges0 = 0;
    if ( fT0 == -Math<Real>::MAX_REAL )
    {
        pkF0 = kSturm[0];
        if ( pkF0->GetDegree() & 1 )
            fValue0 = -(*pkF0)[pkF0->GetDegree()];
        else
            fValue0 = (*pkF0)[pkF0->GetDegree()];

        if ( Math<Real>::FAbs(fValue0) < m_fEpsilon )
            fValue0 = (Real)0.0;

        for (i = 1; i < (int)kSturm.size(); i++)
        {
            pkF1 = kSturm[i];

            if ( pkF1->GetDegree() & 1 )
                fValue1 = -(*pkF1)[pkF1->GetDegree()];
            else
                fValue1 = (*pkF1)[pkF1->GetDegree()];

            if ( Math<Real>::FAbs(fValue1) < m_fEpsilon )
                fValue1 = (Real)0.0;

            if ( fValue0*fValue1 < (Real)0.0 || fValue0 == (Real)0.0 )
                iSignChanges0++;