Interpolate functions using Lagrange interpolation method

lagrange_interpolator.cpp

#include <iostream>
#include <vector>
#include <iomanip>
#include <algorithm>
#include <iterator>
#include <string>
#include <sstream>
#include <iomanip>
#include <cmath>


template<class ValueType>
class rational
{
public:

	rational(): numerator(0), denominator(0){}

	rational(ValueType num, ValueType denom): numerator(num), denomirator(denom)
	{
		reduce();
	}

	rational<ValueType> operator+(rational<ValueType>& right)
	{
		if(right.denominator==this->denominator)
		{
			return rational<ValueType>(this->numerator+right.numerator, this->denominator);
		}
		else
		{
			rational<ValueType> result=*this;
			auto found_lcm=lcm(this->denominator, right.denominator);
			result.numerator*=lcm/this->denominator;
			result.numerator+=right.numerator*lcm/right->denominator;
			result.denominator=lcm;
			result.reduce();
			return result;
		}
	}

	inline rational<ValueType> operator-(rational<ValueType>& right)
	{
		rational<ValueType> temp=*this;
		temp.numerator=-temp.numerator;
		return temp+right;
	}

	inline rational<ValueType> operator*(rational<ValueType>& right)
	{
		return rational<ValueType>(this->numerator*right.numerator, this->denominator*right.denominator);
	}

	inline rational<ValueType> operator/(rational<ValueType>& right)
	{
		rational<ValueType> temp=*this;
		std::swap(temp.numerator, temp.denominator);
		return temp*right;
	}

	inline rational<ValueType> operator+(ValueType right)
	{
		rational<ValueType> temp(right, right);
		return temp+*this;
	}

	inline rational<ValueType> operator-(ValueType right)
	{
		rational<ValueType> temp(right, right);
		return temp-*this;
	}

	inline rational<ValueType> operator*(ValueType right)
	{
		return rational<ValueType>(this->numerator*right, this->denominator);
	}

	inline rational<ValueType> operator/(ValueType right)
	{
		return rational<ValueType>(this->numerator, this->denominator*right);
	}

	double to_double() const
	{
		return numerator/denomirator;
	}

private:

	inline void reduce()
	{
		auto found_gcd=binary_gcd(numerator, denominator);
		numerator/=found_gcd;
		denominator/=found_gcd;
	}

	inline ValueType lcm(ValueType left, ValueType right)
	{
		return std::abs(left*right)/binary_gcd(left, right);
	}

	ValueType binary_gcd(ValueType left, ValueType right)
	{
		if(left==0)
		{
			return right;
		}
		if(right==0)
		{
			return left;
		}
		if(left==right)
		{
			return left;
		}
		if(left==1 || right==1)
		{
			return 1;
		}
		if(left%2==0 && right%2==0)
		{
			return 2*binary_gcd(left/2, right/2);
		}
		if(left%2==0 && right%2!=0)
		{
			return binary_gcd(left/2, right);
		}
		if(right%2==0 && left%2!=0)
		{
			return binary_gcd(left, right/2);
		}
		if(left%2!=0 && right%2!=0 && right>left)
		{
			return binary_gcd((right-left)/2, left);
		}
		if(left%2!=0 && right%2!=0 && left>right)
		{
			return binary_gcd((left-right)/2, right);
		}
	}

	ValueType numerator;
	ValueType denominator;
};

template<>
class rational<float>{};

template<>
class rational<double>{};

template<>
class rational<long double>{};


template<class ValueType>
struct monomial
{
	ValueType coefficient;
	ValueType power;

	inline monomial operator*(monomial& right)
	{
		return std::move(monomial{this->coefficient*right.coefficient, this->power+right.power});
	}

	inline monomial operator*(ValueType right)
	{
		return std::move(monomial{this->coefficient*right, this->power});
	}

	monomial operator/(monomial& right)
	{
		monomial result=*this;
		result.coefficient=result.coefficient/right.coefficient;
		result.power=result.power-right.power;
		return std::move(result);
	}

	std::string to_string(unsigned short prec=2, char expr_variable='x') const
	{
		std::stringstream target;
		target<<std::fixed<<std::setprecision(prec);
		if(coefficient==-1 && power==0)
		{
			target<<(ValueType)-1;
		}
		else if(coefficient==1 && power==0)
		{
			target<<(ValueType)1;
		}
		else if(coefficient!=1)
		{
			target<<coefficient;
		}
		if(power==1)
		{
			target<<expr_variable;
		}
		else if(power!=0)
		{
			target<<expr_variable<<'^'<<power;
		}
		return target.str();
	}

	std::string to_string_detailed(unsigned short prec=2, char expr_variable='x') const
	{
		std::stringstream target;
		target<<std::fixed<<std::setprecision(prec);
		target<<coefficient<<expr_variable<<'^'<<power;
		return target.str();
	}
};


template<class ValueType>
class polynomial
{
public:

	struct division_result
	{
		polynomial<ValueType> result;
		polynomial<ValueType> remainder;
	};

	polynomial()
	{
		parts.push_back(monomial<ValueType>{0, 0});
	}

	polynomial(std::vector<monomial<ValueType>> mp): parts(mp)
	{
		this->compose();
		smart_remove_dead_parts(this->parts);
	}

	polynomial<ValueType> operator+(polynomial<ValueType>& right)
	{
		auto result=this->addition_impl(right);
		smart_remove_dead_parts(result.parts);
		return std::move(result);
	}

	polynomial<ValueType> operator-(polynomial<ValueType>& right)
	{
		polynomial<ValueType> result=right;
		result.inverse();
		result=*this+result;
		//smart_remove_dead_parts(result.parts);
		return std::move(result);
	}

	polynomial<ValueType> operator*(polynomial<ValueType>& right)
	{
		polynomial<ValueType> result;
		this->compose();
		right.compose();
		for(auto& it: right.parts)
		{
			for(auto& st: this->parts)
			{
				result.parts.push_back(it*st);
			}
		}
		result.compose();
		smart_remove_dead_parts(result.parts);
		return std::move(result);
	}

	polynomial<ValueType> operator*(ValueType right)
	{
		auto result=*this;
		for(auto& it: result.parts)
		{
			it=it*right;
		}
		return std::move(result);
	}

	division_result operator/(polynomial<ValueType>& right)
	{
		if(this->highest_power()<right.highest_power())
		{
			return std::move(division_result
			{
				std::move(polynomial<ValueType>()),
				right
			});
		}
		auto reconstructed=reconstruct_parts(*this, this->highest_power());
		auto reconstructed_right=reconstruct_parts(right, this->highest_power());
		reconstructed.do_sort();
		reconstructed_right.do_sort();
		polynomial<ValueType> result;
		for(unsigned int it=1; it<=this->highest_power(); it++)
		{
			result.parts.push_back(monomial<ValueType>{0, (double)it});
		}
		polynomial<ValueType> remainder;
		for(unsigned int it=1; it<=this->highest_power(); it++)
		{
			remainder.parts.push_back(monomial<ValueType>{0, (double)it});
		}
		polynomial<ValueType> temp_right(right.parts);
		unsigned int pos=0;
		while(true)
		{
			auto current_left=reconstructed.parts[pos];
			auto current_right=reconstructed_right.parts[pos];
			auto current_result=current_left/current_right;
			//std::cout<<"    current_right: "<<current_left.to_string_detailed(0)<<std::endl;
			//std::cout<<"    current_left: "<<current_right.to_string_detailed(0)<<std::endl;
			//std::cout<<"    current_result: "<<current_result.to_string_detailed(0)<<std::endl;
			result=result+polynomial<ValueType>({current_result});
			//std::cout<<"    result: "<<result.to_string(0)<<std::endl;
			temp_right=polynomial<ValueType>({current_result})*reconstructed_right;
			//std::cout<<"    temp_right: "<<temp_right.to_string(0)<<std::endl;
			//std::cout<<"	reconstructed: "<<reconstructed.to_string(0)<<std::endl;
			reconstructed=reconstructed-temp_right;
			//std::cout<<"    reconstructed: "<<reconstructed.to_string(0)<<std::endl;
			reconstructed=std::move(reconstruct_parts(reconstructed, reconstructed.highest_power()));
			//std::cout<<"    reconstructed: "<<reconstructed.to_string(0)<<std::endl;
			//std::cout<<std::endl;
			if(current_result.power==reconstructed_right.highest_power())
			{
				break;
			}
		}
		remainder=std::move(reconstructed);
		return std::move(division_result
		{
			std::move(result), 
			std::move(remainder)
		});
	}

	inline double highest_power()
	{
		do_sort();
		return (*(parts.end()-1)).power;
	}

	void compose()
	{
		do_sort();
		unsigned int pos=1;
		while(pos!=parts.size())
		{
			if(parts[pos-1].power==parts[pos].power)
			{
				parts[pos-1].coefficient+=parts[pos].coefficient;
				parts.erase(parts.begin()+pos);
			}
			else
			{
				pos++;
			}
		}
	}

	inline void inverse()
	{
		for(auto& it: parts)
		{
			it.coefficient=-it.coefficient;
		}
	}

	inline void do_reverse_sort()
	{
		std::sort(parts.begin(), parts.end(), 
			[](monomial<ValueType>& left, monomial<ValueType>& right)
		{
			return (left.power<right.power);
		});
	}

	inline void do_sort()
	{
		std::sort(parts.begin(), parts.end(), 
			[](monomial<ValueType>& left, monomial<ValueType>& right)
		{
			return (left.power>right.power);
		});
	}

	std::string to_string(unsigned short prec=2, char expr_variable='x')
	{
		do_sort();
		std::string result;
		for(unsigned int it=0; it<parts.size(); it++)
		{
			if(it!=0 && parts[it].coefficient>=0)
			{
				result+='+';
			}
			result+=parts[it].to_string(prec, expr_variable);
		}
		return result;
	}

	inline const std::vector<monomial<ValueType>>& get_parts() const
	{
		return parts;
	}

private:

	polynomial<ValueType> addition_impl(polynomial<ValueType>& right)
	{
		polynomial<ValueType> result=*this;
		for(auto& it: right.parts)
		{
			result.parts.push_back(it);
		}
		result.compose();
		return std::move(result);
	}

	static polynomial<ValueType> reconstruct_parts(polynomial<ValueType>& target, 
		unsigned int max_power)
	{
		polynomial<ValueType> result;
		target.do_sort();
		for(unsigned int it=1; it<=max_power; it++)
		{
			result.parts.push_back(monomial<ValueType>{0, (ValueType)it});
		}
		return std::move(result.addition_impl(target));
	}

	static void smart_remove_dead_parts(std::vector<monomial<ValueType>>& target_parts)
	{
		unsigned int pos=0;
		while(pos!=target_parts.size())
		{
			if(target_parts[pos].coefficient==0)
			{
				target_parts.erase(target_parts.begin()+pos);
			}
			else
			{
				pos++;
			}
		}
		if(target_parts.size()==0)
		{
			target_parts.push_back(monomial<ValueType>{0, 0});
		}
	}

	std::vector<monomial<ValueType>> parts;
};


using namespace std;

template<class ValueType>
class polynomial_function
{
public:

	polynomial_function(){}

	polynomial_function(polynomial<ValueType> target): target_polynomial(target){}

	ValueType operator()(ValueType variable)
	{
		ValueType result=0;
		for(auto& it: target_polynomial.get_parts())
		{
			result+=it.coefficient*std::pow(variable, it.power);
		}
		return result;
	}

	inline auto& get_polynomial()
	{
		return target_polynomial;
	}

private:

	polynomial<ValueType> target_polynomial;
};


template<class ValueType>
class lagrange_polynomial_interpolator
{
public:

	struct data_point
	{
		ValueType variable;
		ValueType result;
	};

	lagrange_polynomial_interpolator(){}

	inline void add_data_point(data_point point)
	{
		data.push_back(point);
	}

	void calc_basis_polynomials()
	{
		for(unsigned int it=0; it<data.size(); it++)
		{
			polynomial<ValueType> temp({{1, 0}});
			for(unsigned int st=0; st<data.size(); st++)
			{
				if(st!=it)
				{
					temp=temp*(polynomial<ValueType>({{1, 1}, {-(data[st].variable), 0}})*(1/(data[it].variable-data[st].variable)));
				}
			}
			basis_polynomials.push_back(temp);
		}
	}

	void calc_result_function()
	{
		polynomial<ValueType> temp;
		for(unsigned int it=0; it<data.size(); it++)
		{
			auto ctmp=basis_polynomials[it]*data[it].result;
			temp=temp+ctmp;
		}
		target_fn=polynomial_function<ValueType>(temp);
	}

	inline auto& get_basis_polynomials()
	{
		return basis_polynomials;
	}

	inline auto& get_function()
	{
		return target_fn;
	}

private:
	
	std::vector<data_point> data;
	polynomial_function<ValueType> target_fn;
	std::vector<polynomial<ValueType>> basis_polynomials;
};

using namespace std;

void test_case()
{
	polynomial_function<double> result_fn(polynomial<double>({{1, 2}}));
	cout<<"Target function: f(x)="<<result_fn.get_polynomial().to_string(2)<<endl;
	lagrange_polynomial_interpolator<double> interpolator;
	const unsigned int iterations=4;
	for(unsigned int it=0; it<iterations; it++)
	{
		cout<<"f("<<it<<")="<<result_fn(it)<<endl;
		interpolator.add_data_point({(double)it, result_fn((double)it)});
	}
	interpolator.calc_basis_polynomials();
	cout<<endl;
	cout<<"Calculated basis polynomials"<<endl;
	for(auto& it: interpolator.get_basis_polynomials())
	{
		cout<<it.to_string(2)<<endl;
	}

	/*
	cout<<endl;
	for(auto& it: interpolator.get_basis_polynomials())
	{
		cout<<polynomial_function<double>(it)(0)<<endl;
	}
	*/
	
	interpolator.calc_result_function();
	cout<<endl;
	cout<<"Calculated function: L(x)="<<interpolator.get_function().get_polynomial().to_string(2)<<endl;
	for(unsigned int it=0; it<iterations; it++)
	{
		cout<<"L("<<it<<")="<<interpolator.get_function()((double)it)<<endl;
	}
}

//Test case for long division
/*
void test_case()
{
	using polynomial_t=polynomial<double>;
	polynomial_t left({{3, 3}, {-2, 2}, {4, 1}, {-3, 0}});
	polynomial_t right({{1, 2}, {3, 1}, {3, 0}});
	auto div_res=left/right;
	std::cout<<left.to_string(0)<<std::endl;
	std::cout<<right.to_string(0)<<std::endl;
	std::cout<<std::endl;
	std::cout<<"result: "<<div_res.result.to_string(0)<<std::endl;
	std::cout<<"remainder: "<<div_res.remainder.to_string(0)<<std::endl;
}
*/

int main(int argc, char** argv)
{
	test_case();
	system("pause");
	return 0;
}