Polynomial.h

// This file is a part of the OpenSurgSim project.
// Copyright 2013-2016, SimQuest Solutions Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef SURGSIM_MATH_POLYNOMIAL_H
#define SURGSIM_MATH_POLYNOMIAL_H

#include <iostream>

namespace SurgSim
{
namespace Math
{

namespace
{
double polynomial_epsilon = 1.0e-09;
}

template <typename T>
bool isNearZero(const T& value, const T& epsilon = static_cast<T>(polynomial_epsilon));

template <typename T, int N> class Polynomial
{
    static_assert(N >= 0, "Polynomials must have degree >= 0.");
    static_assert(N <= 3, "Polynomials of degree > 3 are not yet supported.");
};

template <typename T>
class Polynomial<T, 0>
{
public:
    Polynomial();

    explicit Polynomial(const T& a0);

    T evaluate(const T& x) const;

    T operator()(const T& x) const;

    T& operator[](const size_t i);
    const T& operator[](const size_t i) const;
    Polynomial<T, 0> operator- () const;
    Polynomial<T, 0> operator+ (const Polynomial<T, 0>& rhs) const;
    Polynomial<T, 0>& operator+= (const Polynomial<T, 0>& rhs);
    Polynomial<T, 0> operator- (const Polynomial<T, 0>& rhs) const;
    Polynomial<T, 0>& operator-= (const Polynomial<T, 0>& rhs);

    bool isNearZero(const T& epsilon = static_cast<T>(polynomial_epsilon)) const;

    bool isApprox(const Polynomial<T, 0>& p, const T& epsilon) const;

    T getCoefficient(size_t i) const;

    void setCoefficient(size_t i, const T& value);

private:
    T m_a0;
};

template <typename T>
class Polynomial<T, 1>
{
public:
    Polynomial();

    Polynomial(const T& a0, const T& a1);

    T evaluate(const T& x) const;

    T operator()(const T& x) const;

    T& operator[](const size_t i);
    const T& operator[](const size_t i) const;
    Polynomial<T, 1> operator- () const;
    Polynomial<T, 1> operator+ (const Polynomial<T, 1>& rhs) const;
    Polynomial<T, 1>& operator+= (const Polynomial<T, 1>& rhs);
    Polynomial<T, 1> operator- (const Polynomial<T, 1>& rhs) const;
    Polynomial<T, 1>& operator-= (const Polynomial<T, 1>& rhs);

    Polynomial<T, 0> derivative() const;

    bool isNearZero(const T& epsilon = static_cast<T>(polynomial_epsilon)) const;

    bool isApprox(const Polynomial<T, 1>& p, const T& epsilon) const;

    T getCoefficient(size_t i) const;

    void setCoefficient(size_t i, const T& value);

private:
    T m_a0;
    T m_a1;
};

template <typename T>
class Polynomial<T, 2>
{
public:
    Polynomial();

    Polynomial(const T& a0, const T& a1, const T& a2);

    T discriminant() const;

    T evaluate(const T& x) const;

    T operator()(const T& x) const;

    T& operator[](const size_t i);
    const T& operator[](const size_t i) const;
    Polynomial<T, 2> operator- () const;
    Polynomial<T, 2> operator+ (const Polynomial<T, 2>& rhs) const;
    Polynomial<T, 2>& operator+= (const Polynomial<T, 2>& rhs);
    Polynomial<T, 2> operator- (const Polynomial<T, 2>& rhs) const;
    Polynomial<T, 2>& operator-= (const Polynomial<T, 2>& rhs);

    Polynomial<T, 1> derivative() const;

    bool isNearZero(const T& epsilon = static_cast<T>(polynomial_epsilon)) const;

    bool isApprox(const Polynomial<T, 2>& p, const T& epsilon) const;

    T getCoefficient(size_t i) const;

    void setCoefficient(size_t i, const T& value);

private:
    T m_a0;
    T m_a1;
    T m_a2;
};

template <typename T>
class Polynomial<T, 3>
{
public:
    Polynomial();

    Polynomial(const T& a0, const T& a1, const T& a2, const T& a3);

    T evaluate(const T& x) const;

    T operator()(const T& x) const;

    T& operator[](const size_t i);
    const T& operator[](const size_t i) const;
    Polynomial<T, 3> operator- () const;
    Polynomial<T, 3> operator+ (const Polynomial<T, 3>& rhs) const;
    Polynomial<T, 3>& operator+= (const Polynomial<T, 3>& rhs);
    Polynomial<T, 3> operator- (const Polynomial<T, 3>& rhs) const;
    Polynomial<T, 3>& operator-= (const Polynomial<T, 3>& rhs);

    Polynomial<T, 2> derivative() const;

    bool isNearZero(const T& epsilon = static_cast<T>(polynomial_epsilon)) const;

    bool isApprox(const Polynomial<T, 3>& p, const T& epsilon) const;

    T getCoefficient(size_t i) const;

    void setCoefficient(size_t i, const T& value);

private:
    T m_a0;
    T m_a1;
    T m_a2;
    T m_a3;
};

// Operators

template <typename T, int N, int M>
Polynomial < T, N + M > operator*(const Polynomial<T, N>& p, const Polynomial<T, M>& q);

template <typename T>
Polynomial<T, 2> operator*(const Polynomial<T, 1>& p, const Polynomial<T, 1>& q);

template <typename T>
Polynomial<T, 3> operator*(const Polynomial<T, 2>& p, const Polynomial<T, 1>& q);

template <typename T>
Polynomial<T, 3> operator*(const Polynomial<T, 1>& p, const Polynomial<T, 2>& q);

template <typename T>
Polynomial<T, 0> square(const Polynomial<T, 0>& p);

template <typename T>
Polynomial<T, 2> square(const Polynomial<T, 1>& p);

template <typename T, int N>
std::ostream& operator<<(std::ostream& stream, const Polynomial<T, N>& p);

}; // Math
}; // SurgSim

#include "SurgSim/Math/Polynomial-inl.h"

#endif // SURGSIM_MATH_POLYNOMIAL_H