wmlpolynomial1.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>
Polynomial1<Real>::Polynomial1 (int iDegree)
{
    if ( iDegree >= 0 )
    {
        m_iDegree = iDegree;
        m_afCoeff = new Real[m_iDegree+1];
    }
    else
    {
        // default creation
        m_iDegree = -1;
        m_afCoeff = NULL;
    }
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>::Polynomial1 (const Polynomial1& rkPoly)
{
    m_iDegree = rkPoly.m_iDegree;
    m_afCoeff = new Real[m_iDegree+1];
    for (int i = 0; i <= m_iDegree; i++)
        m_afCoeff[i] = rkPoly.m_afCoeff[i];
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>::~Polynomial1 ()
{
    delete[] m_afCoeff;
}
//----------------------------------------------------------------------------
template <class Real>
void Polynomial1<Real>::SetDegree (int iDegree)
{
    m_iDegree = iDegree;
    delete[] m_afCoeff;

    if ( m_iDegree >= 0 )
        m_afCoeff = new Real[m_iDegree+1];
    else
        m_afCoeff = NULL;
}
//----------------------------------------------------------------------------
template <class Real>
int Polynomial1<Real>::GetDegree () const
{
    return m_iDegree;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>::operator const Real* () const
{
    return m_afCoeff;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>::operator Real* ()
{
    return m_afCoeff;
}
//----------------------------------------------------------------------------
template <class Real>
Real Polynomial1<Real>::operator[] (int i) const
{
    assert( 0 <= i && i <= m_iDegree );
    return m_afCoeff[i];
}
//----------------------------------------------------------------------------
template <class Real>
Real& Polynomial1<Real>::operator[] (int i)
{
    assert( 0 <= i && i <= m_iDegree );
    return m_afCoeff[i];
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator= (const Polynomial1& rkPoly)
{
    delete[] m_afCoeff;
    m_iDegree = rkPoly.m_iDegree;

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

    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Real Polynomial1<Real>::operator() (Real fT) const
{
    assert( m_iDegree >= 0 );

    Real fResult = m_afCoeff[m_iDegree];
    for (int i = m_iDegree-1; i >= 0; i--)
    {
        fResult *= fT;
        fResult += m_afCoeff[i];
    }
    return fResult;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator+ (const Polynomial1& rkPoly)
    const
{
    assert( m_iDegree >= 0 && rkPoly.m_iDegree >= 0 );

    Polynomial1 kSum;
    int i;

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

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

    return kSum;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator- (const Polynomial1& rkPoly)
    const
{
    assert( m_iDegree >= 0 && rkPoly.m_iDegree >= 0 );

    Polynomial1 kDiff;
    int i;

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

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

    return kDiff;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator* (const Polynomial1& rkPoly)
    const
{
    assert( m_iDegree >= 0 && rkPoly.m_iDegree >= 0 );

    Polynomial1 kProd(m_iDegree + rkPoly.m_iDegree);

    memset(kProd.m_afCoeff,0,(kProd.m_iDegree+1)*sizeof(Real));
    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;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator+ (Real fScalar) const
{
    assert( m_iDegree >= 0 );
    Polynomial1 kSum(m_iDegree);
    kSum.m_afCoeff[0] += fScalar;
    return kSum;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator- (Real fScalar) const
{
    assert( m_iDegree >= 0 );
    Polynomial1 kDiff(m_iDegree);
    kDiff.m_afCoeff[0] -= fScalar;
    return kDiff;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator* (Real fScalar) const
{
    assert( m_iDegree >= 0 );

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

    return kProd;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator/ (Real fScalar) const
{
    assert( m_iDegree >= 0 );

    Polynomial1 kProd(m_iDegree);
    int i;

    if ( fScalar != (Real)0.0 )
    {
        Real fInvScalar = ((Real)1.0)/fScalar;
        for (i = 0; i <= m_iDegree; i++)
            kProd.m_afCoeff[i] = fInvScalar*m_afCoeff[i];
    }
    else
    {
        for (i = 0; i <= m_iDegree; i++)
            kProd.m_afCoeff[i] = Math<Real>::MAX_REAL;
    }

    return kProd;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::operator- () const
{
    assert( m_iDegree >= 0 );

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

    return kNeg;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> operator* (Real fScalar,
    const Polynomial1<Real>& rkPoly)
{
    assert( rkPoly.GetDegree() >= 0 );

    Polynomial1<Real> kProd(rkPoly.GetDegree());
    for (int i = 0; i <= rkPoly.GetDegree(); i++)
        kProd[i] = fScalar*rkPoly[i];

    return kProd;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator += (const Polynomial1& rkPoly)
{
    assert( m_iDegree >= 0 );
    *this = *this + rkPoly;
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator -= (const Polynomial1& rkPoly)
{
    assert( m_iDegree >= 0 );
    *this = *this - rkPoly;
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator *= (const Polynomial1& rkPoly)
{
    assert( m_iDegree >= 0 );
    *this = (*this)*rkPoly;
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator += (Real fScalar)
{
    assert( m_iDegree >= 0 );
    m_afCoeff[0] += fScalar;
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator -= (Real fScalar)
{
    assert( m_iDegree >= 0 );
    m_afCoeff[0] -= fScalar;
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator *= (Real fScalar)
{
    assert( m_iDegree >= 0 );
    *this = (*this)*fScalar;
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real>& Polynomial1<Real>::operator /= (Real fScalar)
{
    assert( m_iDegree >= 0 );
    *this = (*this)/fScalar;
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::GetDerivative () const
{
    if ( m_iDegree > 0 )
    {
        Polynomial1 kDeriv(m_iDegree-1);
        for (int i0 = 0, i1 = 1; i0 < m_iDegree; i0++, i1++)
            kDeriv.m_afCoeff[i0] = i1*m_afCoeff[i1];
        return kDeriv;
    }
    else if ( m_iDegree == 0 )
    {
        Polynomial1 kConst(0);
        kConst.m_afCoeff[0] = (Real)0.0;
        return kConst;
    }
    return Polynomial1<Real>();  // invalid in, invalid out
}
//----------------------------------------------------------------------------
template <class Real>
Polynomial1<Real> Polynomial1<Real>::GetInversion () const
{
    Polynomial1 kInvPoly(m_iDegree);
    for (int i = 0; i <= m_iDegree; i++)
        kInvPoly.m_afCoeff[i] = m_afCoeff[m_iDegree-i];
    return kInvPoly;
}
//----------------------------------------------------------------------------
template <class Real>
void Polynomial1<Real>::Compress (Real fEpsilon)
{
    int i;
    for (i = m_iDegree; i >= 0; i--)
    {
        if ( Math<Real>::FAbs(m_afCoeff[i]) <= fEpsilon )
            m_iDegree--;
        else
            break;
    }

    if ( m_iDegree >= 0 )
    {
        Real fInvLeading = ((Real)1.0)/m_afCoeff[m_iDegree];
        m_afCoeff[m_iDegree] = (Real)1.0;
        for (i = 0; i < m_iDegree; i++)
            m_afCoeff[i] *= fInvLeading;
    }
}
//----------------------------------------------------------------------------
template <class Real>
void Polynomial1<Real>::Divide (const Polynomial1& rkDiv, Polynomial1& rkQuot,
    Polynomial1& rkRem, Real fEpsilon) const
{
    int iQuotDegree = m_iDegree - rkDiv.m_iDegree;
    if ( iQuotDegree >= 0 )
    {
        rkQuot.SetDegree(iQuotDegree);

        // temporary storage for the remainder
        Polynomial1 kTmp = *this;

        // do the division (Euclidean algorithm)
        Real fInv = ((Real)1.0)/rkDiv[rkDiv.m_iDegree];
        for (int iQ = iQuotDegree; iQ >= 0; iQ--)
        {
            int iR = rkDiv.m_iDegree + iQ;
            rkQuot[iQ] = fInv*kTmp[iR];
            for (iR--; iR >= iQ; iR--)
                kTmp[iR] -= rkQuot[iQ]*rkDiv[iR-iQ];
        }

        // calculate the correct degree for the remainder
        int iRemDeg = rkDiv.m_iDegree - 1;
        while ( iRemDeg > 0 && Math<Real>::FAbs(kTmp[iRemDeg]) < fEpsilon )
            iRemDeg--;

        if ( iRemDeg == 0 && Math<Real>::FAbs(kTmp[0]) < fEpsilon )
            kTmp[0] = (Real)0.0;

        rkRem.SetDegree(iRemDeg);
        memcpy(rkRem.m_afCoeff,kTmp.m_afCoeff,(iRemDeg+1)*sizeof(Real));
    }
    else
    {
        rkQuot.SetDegree(0);
        rkQuot[0] = (Real)0.0;
        rkRem = *this;
    }
}
//----------------------------------------------------------------------------