SSB_HighSpeed_Modem/hsmodem/fec/schifra_galois_field_polynomial.hpp
2021-01-15 22:15:58 +01:00

840 lines
26 KiB
C++
Executable File

/*
(**************************************************************************)
(* *)
(* Schifra *)
(* Reed-Solomon Error Correcting Code Library *)
(* *)
(* Release Version 0.0.1 *)
(* http://www.schifra.com *)
(* Copyright (c) 2000-2020 Arash Partow, All Rights Reserved. *)
(* *)
(* The Schifra Reed-Solomon error correcting code library and all its *)
(* components are supplied under the terms of the General Schifra License *)
(* agreement. The contents of the Schifra Reed-Solomon error correcting *)
(* code library and all its components may not be copied or disclosed *)
(* except in accordance with the terms of that agreement. *)
(* *)
(* URL: http://www.schifra.com/license.html *)
(* *)
(**************************************************************************)
*/
#ifndef INCLUDE_SCHIFRA_GALOIS_FIELD_POLYNOMIAL_HPP
#define INCLUDE_SCHIFRA_GALOIS_FIELD_POLYNOMIAL_HPP
#include <cassert>
#include <iostream>
#include <vector>
#include "schifra_galois_field.hpp"
#include "schifra_galois_field_element.hpp"
namespace schifra
{
namespace galois
{
class field_polynomial
{
public:
field_polynomial(const field& gfield);
field_polynomial(const field& gfield, const unsigned int& degree);
field_polynomial(const field& gfield, const unsigned int& degree, const field_element element[]);
field_polynomial(const field_polynomial& polynomial);
field_polynomial(const field_element& gfe);
~field_polynomial() {}
bool valid() const;
int deg() const;
const field& galois_field() const;
void set_degree(const unsigned int& x);
void simplify();
field_polynomial& operator = (const field_polynomial& polynomial);
field_polynomial& operator = (const field_element& element);
field_polynomial& operator += (const field_polynomial& element);
field_polynomial& operator += (const field_element& element);
field_polynomial& operator -= (const field_polynomial& element);
field_polynomial& operator -= (const field_element& element);
field_polynomial& operator *= (const field_polynomial& polynomial);
field_polynomial& operator *= (const field_element& element);
field_polynomial& operator /= (const field_polynomial& divisor);
field_polynomial& operator /= (const field_element& element);
field_polynomial& operator %= (const field_polynomial& divisor);
field_polynomial& operator %= (const unsigned int& power);
field_polynomial& operator ^= (const unsigned int& n);
field_polynomial& operator <<= (const unsigned int& n);
field_polynomial& operator >>= (const unsigned int& n);
field_element& operator[] (const std::size_t& term);
field_element operator() (const field_element& value);
field_element operator() (field_symbol value);
const field_element& operator[](const std::size_t& term) const;
const field_element operator()(const field_element& value) const;
const field_element operator()(field_symbol value) const;
bool operator==(const field_polynomial& polynomial) const;
bool operator!=(const field_polynomial& polynomial) const;
bool monic() const;
field_polynomial derivative() const;
friend std::ostream& operator << (std::ostream& os, const field_polynomial& polynomial);
private:
typedef std::vector<field_element>::iterator poly_iter;
typedef std::vector<field_element>::const_iterator const_poly_iter;
void simplify(field_polynomial& polynomial) const;
field& field_;
std::vector<field_element> poly_;
};
field_polynomial operator + (const field_polynomial& a, const field_polynomial& b);
field_polynomial operator + (const field_polynomial& a, const field_element& b);
field_polynomial operator + (const field_element& a, const field_polynomial& b);
field_polynomial operator + (const field_polynomial& a, const field_symbol& b);
field_polynomial operator + (const field_symbol& a, const field_polynomial& b);
field_polynomial operator - (const field_polynomial& a, const field_polynomial& b);
field_polynomial operator - (const field_polynomial& a, const field_element& b);
field_polynomial operator - (const field_element& a, const field_polynomial& b);
field_polynomial operator - (const field_polynomial& a, const field_symbol& b);
field_polynomial operator - (const field_symbol& a, const field_polynomial& b);
field_polynomial operator * (const field_polynomial& a, const field_polynomial& b);
field_polynomial operator * (const field_element& a, const field_polynomial& b);
field_polynomial operator * (const field_polynomial& a, const field_element& b);
field_polynomial operator / (const field_polynomial& a, const field_polynomial& b);
field_polynomial operator / (const field_polynomial& a, const field_element& b);
field_polynomial operator % (const field_polynomial& a, const field_polynomial& b);
field_polynomial operator % (const field_polynomial& a, const unsigned int& power);
field_polynomial operator ^ (const field_polynomial& a, const int& n);
field_polynomial operator <<(const field_polynomial& a, const unsigned int& n);
field_polynomial operator >>(const field_polynomial& a, const unsigned int& n);
field_polynomial gcd(const field_polynomial& a, const field_polynomial& b);
inline field_polynomial::field_polynomial(const field& gfield)
: field_(const_cast<field&>(gfield))
{
poly_.clear();
poly_.reserve(256);
}
inline field_polynomial::field_polynomial(const field& gfield, const unsigned int& degree)
: field_(const_cast<field&>(gfield))
{
poly_.reserve(256);
poly_.resize(degree + 1,field_element(field_,0));
}
inline field_polynomial::field_polynomial(const field& gfield, const unsigned int& degree, const field_element element[])
: field_(const_cast<field&>(gfield))
{
poly_.reserve(256);
if (element != NULL)
{
/*
It is assumed that element is an array of field elements
with size/element count of degree + 1.
*/
for (unsigned int i = 0; i <= degree; ++i)
{
poly_.push_back(element[i]);
}
}
else
poly_.resize(degree + 1, field_element(field_, 0));
}
inline field_polynomial::field_polynomial(const field_polynomial& polynomial)
: field_(const_cast<field&>(polynomial.field_)),
poly_ (polynomial.poly_)
{}
inline field_polynomial::field_polynomial(const field_element& element)
: field_(const_cast<field&>(element.galois_field()))
{
poly_.resize(1,element);
}
inline bool field_polynomial::valid() const
{
return (poly_.size() > 0);
}
inline int field_polynomial::deg() const
{
return static_cast<int>(poly_.size()) - 1;
}
inline const field& field_polynomial::galois_field() const
{
return field_;
}
inline void field_polynomial::set_degree(const unsigned int& x)
{
poly_.resize(x - 1,field_element(field_,0));
}
inline field_polynomial& field_polynomial::operator = (const field_polynomial& polynomial)
{
if ((this != &polynomial) && (&field_ == &(polynomial.field_)))
{
poly_ = polynomial.poly_;
}
return *this;
}
inline field_polynomial& field_polynomial::operator = (const field_element& element)
{
if (&field_ == &(element.galois_field()))
{
poly_.resize(1,element);
}
return *this;
}
inline field_polynomial& field_polynomial::operator += (const field_polynomial& polynomial)
{
if (&field_ == &(polynomial.field_))
{
if (poly_.size() < polynomial.poly_.size())
{
const_poly_iter it0 = polynomial.poly_.begin();
for (poly_iter it1 = poly_.begin(); it1 != poly_.end(); ++it0, ++it1)
{
(*it1) += (*it0);
}
while (it0 != polynomial.poly_.end())
{
poly_.push_back(*it0);
++it0;
}
}
else
{
poly_iter it0 = poly_.begin();
for (const_poly_iter it1 = polynomial.poly_.begin(); it1 != polynomial.poly_.end(); ++it0, ++it1)
{
(*it0) += (*it1);
}
}
simplify(*this);
}
return *this;
}
inline field_polynomial& field_polynomial::operator += (const field_element& element)
{
poly_[0] += element;
return *this;
}
inline field_polynomial& field_polynomial::operator -= (const field_polynomial& element)
{
return (*this += element);
}
inline field_polynomial& field_polynomial::operator -= (const field_element& element)
{
poly_[0] -= element;
return *this;
}
inline field_polynomial& field_polynomial::operator *= (const field_polynomial& polynomial)
{
if (&field_ == &(polynomial.field_))
{
field_polynomial product(field_,deg() + polynomial.deg() + 1);
poly_iter result_it = product.poly_.begin();
for (poly_iter it0 = poly_.begin(); it0 != poly_.end(); ++it0)
{
poly_iter current_result_it = result_it;
for (const_poly_iter it1 = polynomial.poly_.begin(); it1 != polynomial.poly_.end(); ++it1)
{
(*current_result_it) += (*it0) * (*it1);
++current_result_it;
}
++result_it;
}
simplify(product);
poly_ = product.poly_;
}
return *this;
}
inline field_polynomial& field_polynomial::operator *= (const field_element& element)
{
if (field_ == element.galois_field())
{
for (poly_iter it = poly_.begin(); it != poly_.end(); ++it)
{
(*it) *= element;
}
}
return *this;
}
inline field_polynomial& field_polynomial::operator /= (const field_polynomial& divisor)
{
if (
(&field_ == &divisor.field_) &&
(deg() >= divisor.deg()) &&
(divisor.deg() >= 0)
)
{
field_polynomial quotient (field_, deg() - divisor.deg() + 1);
field_polynomial remainder(field_, divisor.deg() - 1);
for (int i = static_cast<int>(deg()); i >= 0; i--)
{
if (i <= static_cast<int>(quotient.deg()))
{
quotient[i] = remainder[remainder.deg()] / divisor[divisor.deg()];
for (int j = static_cast<int>(remainder.deg()); j > 0; --j)
{
remainder[j] = remainder[j - 1] + (quotient[i] * divisor[j]);
}
remainder[0] = poly_[i] + (quotient[i] * divisor[0]);
}
else
{
for (int j = static_cast<int>(remainder.deg()); j > 0; --j)
{
remainder[j] = remainder[j - 1];
}
remainder[0] = poly_[i];
}
}
simplify(quotient);
poly_ = quotient.poly_;
}
return *this;
}
inline field_polynomial& field_polynomial::operator /= (const field_element& element)
{
if (field_ == element.galois_field())
{
for (poly_iter it = poly_.begin(); it != poly_.end(); ++it)
{
(*it) /= element;
}
}
return *this;
}
inline field_polynomial& field_polynomial::operator %= (const field_polynomial& divisor)
{
if (
(field_ == divisor.field_) &&
(deg() >= divisor.deg() ) &&
(divisor.deg() >= 0 )
)
{
field_polynomial quotient (field_, deg() - divisor.deg() + 1);
field_polynomial remainder(field_, divisor.deg() - 1);
for (int i = static_cast<int>(deg()); i >= 0; i--)
{
if (i <= static_cast<int>(quotient.deg()))
{
quotient[i] = remainder[remainder.deg()] / divisor[divisor.deg()];
for (int j = static_cast<int>(remainder.deg()); j > 0; --j)
{
remainder[j] = remainder[j - 1] + (quotient[i] * divisor[j]);
}
remainder[0] = poly_[i] + (quotient[i] * divisor[0]);
}
else
{
for (int j = static_cast<int>(remainder.deg()); j > 0; --j)
{
remainder[j] = remainder[j - 1];
}
remainder[0] = poly_[i];
}
}
poly_ = remainder.poly_;
}
return *this;
}
inline field_polynomial& field_polynomial::operator %= (const unsigned int& power)
{
if (poly_.size() >= power)
{
poly_.resize(power,field_element(field_,0));
simplify(*this);
}
return *this;
}
inline field_polynomial& field_polynomial::operator ^= (const unsigned int& n)
{
field_polynomial result = *this;
for (std::size_t i = 0; i < n; ++i)
{
result *= *this;
}
*this = result;
return *this;
}
inline field_polynomial& field_polynomial::operator <<= (const unsigned int& n)
{
if (poly_.size() > 0)
{
size_t initial_size = poly_.size();
poly_.resize(poly_.size() + n, field_element(field_,0));
for (size_t i = initial_size - 1; static_cast<int>(i) >= 0; --i)
{
poly_[i + n] = poly_[i];
}
for (unsigned int i = 0; i < n; ++i)
{
poly_[i] = 0;
}
}
return *this;
}
inline field_polynomial& field_polynomial::operator >>= (const unsigned int& n)
{
if (n <= poly_.size())
{
for (unsigned int i = 0; i <= deg() - n; ++i)
{
poly_[i] = poly_[i + n];
}
poly_.resize(poly_.size() - n,field_element(field_,0));
}
else if (static_cast<int>(n) >= (deg() + 1))
{
poly_.resize(0,field_element(field_,0));
}
return *this;
}
inline const field_element& field_polynomial::operator [] (const std::size_t& term) const
{
assert(term < poly_.size());
return poly_[term];
}
inline field_element& field_polynomial::operator [] (const std::size_t& term)
{
assert(term < poly_.size());
return poly_[term];
}
inline field_element field_polynomial::operator () (const field_element& value)
{
field_element result(field_,0);
if (!poly_.empty())
{
int i = 0;
field_symbol total_sum = 0 ;
field_symbol value_poly_form = value.poly();
for (poly_iter it = poly_.begin(); it != poly_.end(); ++it, ++i)
{
total_sum ^= field_.mul(field_.exp(value_poly_form,i), (*it).poly());
}
result = total_sum;
}
return result;
}
inline const field_element field_polynomial::operator () (const field_element& value) const
{
if (!poly_.empty())
{
int i = 0;
field_symbol total_sum = 0 ;
field_symbol value_poly_form = value.poly();
for (const_poly_iter it = poly_.begin(); it != poly_.end(); ++it, ++i)
{
total_sum ^= field_.mul(field_.exp(value_poly_form,i), (*it).poly());
}
return field_element(field_,total_sum);
}
return field_element(field_,0);
}
inline field_element field_polynomial::operator () (field_symbol value)
{
if (!poly_.empty())
{
int i = 0;
field_symbol total_sum = 0 ;
for (const_poly_iter it = poly_.begin(); it != poly_.end(); ++it, ++i)
{
total_sum ^= field_.mul(field_.exp(value,i), (*it).poly());
}
return field_element(field_,total_sum);
}
return field_element(field_,0);
}
inline const field_element field_polynomial::operator () (field_symbol value) const
{
if (!poly_.empty())
{
int i = 0;
field_symbol total_sum = 0 ;
for (const_poly_iter it = poly_.begin(); it != poly_.end(); ++it, ++i)
{
total_sum ^= field_.mul(field_.exp(value, i), (*it).poly());
}
return field_element(field_,total_sum);
}
return field_element(field_,0);
}
inline bool field_polynomial::operator == (const field_polynomial& polynomial) const
{
if (field_ == polynomial.field_)
{
if (poly_.size() != polynomial.poly_.size())
return false;
else
{
const_poly_iter it0 = polynomial.poly_.begin();
for (const_poly_iter it1 = poly_.begin(); it1 != poly_.end(); ++it0, ++it1)
{
if ((*it0) != (*it1))
return false;
}
return true;
}
}
else
return false;
}
inline bool field_polynomial::operator != (const field_polynomial& polynomial) const
{
return !(*this == polynomial);
}
inline field_polynomial field_polynomial::derivative() const
{
if ((*this).poly_.size() > 1)
{
field_polynomial deriv(field_,deg());
const std::size_t upper_bound = poly_.size() - 1;
for (std::size_t i = 0; i < upper_bound; i += 2)
{
deriv.poly_[i] = poly_[i + 1];
}
simplify(deriv);
return deriv;
}
return field_polynomial(field_,0);
}
inline bool field_polynomial::monic() const
{
return (poly_[poly_.size() - 1] == static_cast<galois::field_symbol>(1));
}
inline void field_polynomial::simplify()
{
simplify(*this);
}
inline void field_polynomial::simplify(field_polynomial& polynomial) const
{
std::size_t poly_size = polynomial.poly_.size();
if ((poly_size > 0) && (polynomial.poly_.back() == 0))
{
poly_iter it = polynomial.poly_.end ();
poly_iter begin = polynomial.poly_.begin();
std::size_t count = 0;
while ((begin != it) && (*(--it) == 0))
{
++count;
}
if (0 != count)
{
polynomial.poly_.resize(poly_size - count, field_element(field_,0));
}
}
}
inline field_polynomial operator + (const field_polynomial& a, const field_polynomial& b)
{
field_polynomial result = a;
result += b;
return result;
}
inline field_polynomial operator + (const field_polynomial& a, const field_element& b)
{
field_polynomial result = a;
result += b;
return result;
}
inline field_polynomial operator + (const field_element& a, const field_polynomial& b)
{
field_polynomial result = b;
result += a;
return result;
}
inline field_polynomial operator + (const field_polynomial& a, const field_symbol& b)
{
return a + field_element(a.galois_field(),b);
}
inline field_polynomial operator + (const field_symbol& a, const field_polynomial& b)
{
return b + field_element(b.galois_field(),a);
}
inline field_polynomial operator - (const field_polynomial& a, const field_polynomial& b)
{
field_polynomial result = a;
result -= b;
return result;
}
inline field_polynomial operator - (const field_polynomial& a, const field_element& b)
{
field_polynomial result = a;
result -= b;
return result;
}
inline field_polynomial operator - (const field_element& a, const field_polynomial& b)
{
field_polynomial result = b;
result -= a;
return result;
}
inline field_polynomial operator - (const field_polynomial& a, const field_symbol& b)
{
return a - field_element(a.galois_field(),b);
}
inline field_polynomial operator - (const field_symbol& a, const field_polynomial& b)
{
return b - field_element(b.galois_field(),a);
}
inline field_polynomial operator * (const field_polynomial& a, const field_polynomial& b)
{
field_polynomial result = a;
result *= b;
return result;
}
inline field_polynomial operator * (const field_element& a, const field_polynomial& b)
{
field_polynomial result = b;
result *= a;
return result;
}
inline field_polynomial operator * (const field_polynomial& a, const field_element& b)
{
field_polynomial result = a;
result *= b;
return result;
}
inline field_polynomial operator / (const field_polynomial& a, const field_polynomial& b)
{
field_polynomial result = a;
result /= b;
return result;
}
inline field_polynomial operator / (const field_polynomial& a, const field_element& b)
{
field_polynomial result = a;
result /= b;
return result;
}
inline field_polynomial operator % (const field_polynomial& a, const field_polynomial& b)
{
field_polynomial result = a;
result %= b;
return result;
}
inline field_polynomial operator % (const field_polynomial& a, const unsigned int& n)
{
field_polynomial result = a;
result %= n;
return result;
}
inline field_polynomial operator ^ (const field_polynomial& a, const int& n)
{
field_polynomial result = a;
result ^= n;
return result;
}
inline field_polynomial operator << (const field_polynomial& a, const unsigned int& n)
{
field_polynomial result = a;
result <<= n;
return result;
}
inline field_polynomial operator >> (const field_polynomial& a, const unsigned int& n)
{
field_polynomial result = a;
result >>= n;
return result;
}
inline field_polynomial gcd(const field_polynomial& a, const field_polynomial& b)
{
if (&a.galois_field() == &b.galois_field())
{
if ((!a.valid()) && (!b.valid()))
{
field_polynomial error_polynomial(a.galois_field());
return error_polynomial;
}
if (!a.valid()) return b;
if (!b.valid()) return a;
field_polynomial x = a % b;
field_polynomial y = b;
field_polynomial z = x;
while ((z = (y % x)).valid())
{
y = x;
x = z;
}
return x;
}
else
{
field_polynomial error_polynomial(a.galois_field());
return error_polynomial;
}
}
inline field_polynomial generate_X(const field& gfield)
{
const field_element xgfe[2] = {
galois::field_element(gfield, 0),
galois::field_element(gfield, 1)
};
field_polynomial X_(gfield,1,xgfe);
return X_;
}
inline std::ostream& operator << (std::ostream& os, const field_polynomial& polynomial)
{
if (polynomial.deg() >= 0)
{
/*
for (unsigned int i = 0; i < polynomial.poly_.size(); ++i)
{
os << polynomial.poly[i].index()
<< ((i != (polynomial.deg())) ? " " : "");
}
std::cout << " poly form: ";
*/
for (unsigned int i = 0; i < polynomial.poly_.size(); ++i)
{
os << polynomial.poly_[i].poly()
<< " "
<< "x^"
<< i
<< ((static_cast<int>(i) != (polynomial.deg())) ? " + " : "");
}
}
return os;
}
} // namespace galois
} // namespace schifra
#endif