polynomial.hpp

polynomial.hpp

#pragma once

// Copyright plasma-effect 2015
// Distributed under the Boost Software License, Version 1.0.
// (See http://www.boost.org/LICENSE_1_0.txt)

namespace plasma
{
	namespace polynomial_types
	{
		template<class T, class Derived>struct polynomial_concept
		{
			T operator()(T const& v)const
			{
				return static_cast<Derived const*>(this)->operator()(v);
			}
			operator Derived()const
			{
				return *static_cast<Derived const*>(this);
			}
			typedef T value_type;
		};

		template<class T>struct monomial :polynomial_concept<T, monomial<T>>
		{
			T operator()(T const& v)const
			{
				return v;
			}
			monomial() = default;
			monomial(monomial<T> const&) = default;
			monomial(monomial<T>&&) = default;
		};
		
		template<class T, class Poly>inline auto operator+(polynomial_concept<T, Poly>const& p)
		{
			return p;
		}
		template<class T, class Poly>struct polynomial_minus
			:polynomial_concept<T, polynomial_minus<T, Poly>>
		{
			Poly p_;

			polynomial_minus(Poly const& p):p_(p){}
			polynomial_minus(polynomial_minus<T, Poly>const&) = default;
			polynomial_minus(polynomial_minus<T, Poly>&&) = default;
			T operator()(T const&v)const
			{
				return -p_(v);
			}
		};
		template<class T, class Poly>inline auto operator-(polynomial_concept<T, Poly>const& p)
		{
			return polynomial_minus<T, Poly>(p);
		}

		template<class T, class Poly>struct polynomial_add_t
			:polynomial_concept<T, polynomial_add_t<T, Poly>>
		{
			T value_;
			Poly poly_;
			T operator()(T const& v)const
			{
				return value_ + poly_(v);
			}
			polynomial_add_t(T value, Poly poly) :value_(value), poly_(poly) {}
			polynomial_add_t(polynomial_add_t<T, Poly>const&) = default;
			polynomial_add_t(polynomial_add_t<T, Poly>&&) = default;
		};
		template<class Poly0, class Poly1,class T>struct polynomial_add2_t
			:polynomial_concept<T, polynomial_add2_t<Poly0, Poly1,T>>
		{
			Poly0 poly0_;
			Poly1 poly1_;
			T operator()(T const& v)const
			{
				return poly0_(v) + poly1_(v);
			}
			polynomial_add2_t(Poly0 poly0, Poly1 poly1) :poly0_(poly0), poly1_(poly1) {}
			polynomial_add2_t(polynomial_add2_t<Poly0, Poly1,T> const&) = default;
			polynomial_add2_t(polynomial_add2_t<Poly0, Poly1,T>&&) = default;
		};

		template<class T, class Poly>inline auto operator+(
			T const& v, polynomial_concept<T, Poly> const& p)
		{
			return polynomial_add_t<T, Poly>(v, p);
		}
		template<class T, class Poly>inline auto operator+(
			polynomial_concept<T, Poly> const& p, T const& v)
		{
			return polynomial_add_t<T, Poly>(v, p);
		}
		template<class Poly0, class Poly1, class T>inline auto operator+(
			polynomial_concept<T, Poly0> const& p0,
			polynomial_concept<T, Poly1> const& p1)
		{
			return polynomial_add2_t<Poly0, Poly1,T>(p0, p1);
		}

		template<class T, class Poly>struct polynomial_sub_t
			:polynomial_concept<T, polynomial_sub_t<T, Poly>>
		{
			T value_;
			Poly poly_;
			T operator()(T const& v)const
			{
				return value_ - poly_(v);
			}
			polynomial_sub_t(T value, Poly poly) :value_(value), poly_(poly) {}
			polynomial_sub_t(polynomial_sub_t<T, Poly>const&) = default;
			polynomial_sub_t(polynomial_sub_t<T, Poly>&&) = default;
		};
		template<class Poly0, class Poly1,class T>struct polynomial_sub2_t
			:polynomial_concept<T, polynomial_sub2_t<Poly0, Poly1,T>>
		{
			Poly0 poly0_;
			Poly1 poly1_;
			T operator()(T const& v)const
			{
				return poly0_(v) - poly1_(v);
			}
			polynomial_sub2_t(Poly0 poly0, Poly1 poly1) :poly0_(poly0), poly1_(poly1) {}
			polynomial_sub2_t(polynomial_sub2_t<Poly0, Poly1,T> const&) = default;
			polynomial_sub2_t(polynomial_sub2_t<Poly0, Poly1,T>&&) = default;
		};

		template<class T, class Poly>inline auto operator-(
			T const& v, polynomial_concept<T, Poly> const& p)
		{
			return polynomial_sub_t<T, Poly>(v, p);
		}
		template<class T, class Poly>inline auto operator-(
			polynomial_concept<T, Poly> const& p, T const& v)
		{
			return polynomial_sub_t<T, polynomial_minus<T,Poly>>(-v, -static_cast<Poly>(p));
		}
		template<class Poly0, class Poly1, class T>inline auto operator-(
			polynomial_concept<T, Poly0> const& p0,
			polynomial_concept<T, Poly1> const& p1)
		{
			return polynomial_sub2_t<Poly0, Poly1,T>(p0, p1);
		}



		template<class T, class Poly>struct polynomial_mul_t
			:polynomial_concept<T, polynomial_mul_t<T, Poly>>
		{
			T value_;
			Poly poly_;
			T operator()(T const& v)const
			{
				return value_ * poly_(v);
			}
			polynomial_mul_t(T value, Poly poly) :value_(value), poly_(poly) {}
			polynomial_mul_t(polynomial_mul_t<T, Poly>const&) = default;
			polynomial_mul_t(polynomial_mul_t<T, Poly>&&) = default;
		};
		template<class Poly0, class Poly1,class T>struct polynomial_mul2_t
			:polynomial_concept<T, polynomial_mul2_t<Poly0, Poly1,T>>
		{
			Poly0 poly0_;
			Poly1 poly1_;
			T operator()(T const& v)const
			{
				return poly0_(v) * poly1_(v);
			}
			polynomial_mul2_t(Poly0 poly0, Poly1 poly1) :poly0_(poly0), poly1_(poly1) {}
			polynomial_mul2_t(polynomial_mul2_t<Poly0, Poly1,T> const&) = default;
			polynomial_mul2_t(polynomial_mul2_t<Poly0, Poly1,T>&&) = default;
		};

		template<class T, class Poly>inline auto operator*(
			T const& v, polynomial_concept<T, Poly> const& p)
		{
			return polynomial_mul_t<T, Poly>(v, p);
		}
		template<class T, class Poly>inline auto operator*(
			polynomial_concept<T, Poly> const& p, T const& v)
		{
			return polynomial_mul_t<T, Poly>(v, p);
		}
		template<class Poly0, class Poly1, class T>inline auto operator*(
			polynomial_concept<T, Poly0> const& p0,
			polynomial_concept<T, Poly1> const& p1)
		{
			return polynomial_mul2_t<Poly0, Poly1,T>(p0, p1);
		}
	
		namespace detail
		{
			template<class T>constexpr T power(T const& v, unsigned int i)
			{
				return (i == 0 ? T(1) : (i == 2 ? v*v : (i % 2 == 0 ? power(power(v, i / 2), 2) : power(power(v, i / 2), 2) + v)));
			}
		}
		template<class Poly>struct polynomial_power_t
			:polynomial_concept<typename Poly::value_type, polynomial_power_t<Poly>>
		{
			Poly poly_;
			unsigned int I;
			polynomial_power_t(Poly const& poly, unsigned int i) :poly_(poly), I(i) {}
			polynomial_power_t(polynomial_power_t<Poly> const&) = default;
			polynomial_power_t(polynomial_power_t<Poly>&&) = default;

			typename Poly::value_type operator()(typename Poly::value_type const& v)const
			{
				return detail::power(poly_(v), I);
			}
		};
		template<class Poly, class T>inline auto operator^(polynomial_concept<T, Poly> const& p, unsigned int i)
		{
			return polynomial_power_t<Poly>(p, i);
		}

		template<class Poly1, class Poly2, class T>struct composition_t
			:polynomial_concept<T, composition_t<Poly1, Poly2, T>>
		{
			Poly1 p1_;
			Poly2 p2_;
			composition_t(Poly1 const& p1, Poly2 const& p2) :p1_(p1), p2_(p2) {}
			composition_t(composition_t<Poly1, Poly2, T> const&) = default;
			composition_t(composition_t<Poly1, Poly2, T>&&) = default;

			T operator()(T const& v)const
			{
				return p2_(p1_(v));
			}
		};
	
		template<class Poly, class T>struct self_composition_t
			:polynomial_concept<T, self_composition_t<Poly, T>>
		{
			Poly p_;
			unsigned int I;

			self_composition_t(Poly p, unsigned int i) :p_(p), I(i) {}
			self_composition_t(self_composition_t<Poly, T> const&) = default;
			self_composition_t(self_composition_t<Poly, T>&&) = default;

			T operator()(T const& v)const
			{
				return call(v, I);
			}
		private:
			T call(T const& v, unsigned int i)const
			{
				return (i == 0 ? v : p_(call(v, i - 1)));
			}
		};
	}
	template<class T>inline auto make_monomial()
	{
		return polynomial_types::monomial<T>();
	}
	template<class Poly1, class Poly2, class T>inline auto composition(
		polynomial_types::polynomial_concept<T, Poly1> const& first,
		polynomial_types::polynomial_concept<T, Poly2> const& second)
	{
		return polynomial_types::composition_t<Poly1, Poly2, T>(first, second);
	}
	template<class Poly, class T>inline auto self_composition(
		polynomial_types::polynomial_concept<T, Poly> const& p, unsigned int count)
	{
		return polynomial_types::self_composition_t<Poly,T>(p, count);
	}
}

test.cpp

#include<iostream>
#include"polynomial.hpp"

int main()
{
	const auto x = plasma::make_monomial<int>();
	
	const auto f = (x ^ 2) + 2 * x + 3;
	const auto g = plasma::self_composition((x ^ 2) - 1, 2);
	const auto h = plasma::composition(g, f);
	std::cout << f(4) << std::endl;
	std::cout << g(4) << std::endl;
	std::cout << h(4) << std::endl;
}