You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
577 lines
14 KiB
577 lines
14 KiB
// polynomi.cpp - written and placed in the public domain by Wei Dai |
|
|
|
// Part of the code for polynomial evaluation and interpolation |
|
// originally came from Hal Finney's public domain secsplit.c. |
|
|
|
#include "pch.h" |
|
#include "polynomi.h" |
|
#include "secblock.h" |
|
|
|
#include <sstream> |
|
#include <iostream> |
|
|
|
NAMESPACE_BEGIN(CryptoPP) |
|
|
|
template <class T> |
|
void PolynomialOver<T>::Randomize(RandomNumberGenerator &rng, const RandomizationParameter ¶meter, const Ring &ring) |
|
{ |
|
m_coefficients.resize(parameter.m_coefficientCount); |
|
for (unsigned int i=0; i<m_coefficients.size(); ++i) |
|
m_coefficients[i] = ring.RandomElement(rng, parameter.m_coefficientParameter); |
|
} |
|
|
|
template <class T> |
|
void PolynomialOver<T>::FromStr(const char *str, const Ring &ring) |
|
{ |
|
std::istringstream in((char *)str); |
|
bool positive = true; |
|
CoefficientType coef; |
|
unsigned int power; |
|
|
|
while (in) |
|
{ |
|
std::ws(in); |
|
if (in.peek() == 'x') |
|
coef = ring.MultiplicativeIdentity(); |
|
else |
|
in >> coef; |
|
|
|
std::ws(in); |
|
if (in.peek() == 'x') |
|
{ |
|
in.get(); |
|
std::ws(in); |
|
if (in.peek() == '^') |
|
{ |
|
in.get(); |
|
in >> power; |
|
} |
|
else |
|
power = 1; |
|
} |
|
else |
|
power = 0; |
|
|
|
if (!positive) |
|
coef = ring.Inverse(coef); |
|
|
|
SetCoefficient(power, coef, ring); |
|
|
|
std::ws(in); |
|
switch (in.get()) |
|
{ |
|
case '+': |
|
positive = true; |
|
break; |
|
case '-': |
|
positive = false; |
|
break; |
|
default: |
|
return; // something's wrong with the input string |
|
} |
|
} |
|
} |
|
|
|
template <class T> |
|
unsigned int PolynomialOver<T>::CoefficientCount(const Ring &ring) const |
|
{ |
|
unsigned count = m_coefficients.size(); |
|
while (count && ring.Equal(m_coefficients[count-1], ring.Identity())) |
|
count--; |
|
const_cast<std::vector<CoefficientType> &>(m_coefficients).resize(count); |
|
return count; |
|
} |
|
|
|
template <class T> |
|
typename PolynomialOver<T>::CoefficientType PolynomialOver<T>::GetCoefficient(unsigned int i, const Ring &ring) const |
|
{ |
|
return (i < m_coefficients.size()) ? m_coefficients[i] : ring.Identity(); |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T>& PolynomialOver<T>::operator=(const PolynomialOver<T>& t) |
|
{ |
|
if (this != &t) |
|
{ |
|
m_coefficients.resize(t.m_coefficients.size()); |
|
for (unsigned int i=0; i<m_coefficients.size(); i++) |
|
m_coefficients[i] = t.m_coefficients[i]; |
|
} |
|
return *this; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T>& PolynomialOver<T>::Accumulate(const PolynomialOver<T>& t, const Ring &ring) |
|
{ |
|
unsigned int count = t.CoefficientCount(ring); |
|
|
|
if (count > CoefficientCount(ring)) |
|
m_coefficients.resize(count, ring.Identity()); |
|
|
|
for (unsigned int i=0; i<count; i++) |
|
ring.Accumulate(m_coefficients[i], t.GetCoefficient(i, ring)); |
|
|
|
return *this; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T>& PolynomialOver<T>::Reduce(const PolynomialOver<T>& t, const Ring &ring) |
|
{ |
|
unsigned int count = t.CoefficientCount(ring); |
|
|
|
if (count > CoefficientCount(ring)) |
|
m_coefficients.resize(count, ring.Identity()); |
|
|
|
for (unsigned int i=0; i<count; i++) |
|
ring.Reduce(m_coefficients[i], t.GetCoefficient(i, ring)); |
|
|
|
return *this; |
|
} |
|
|
|
template <class T> |
|
typename PolynomialOver<T>::CoefficientType PolynomialOver<T>::EvaluateAt(const CoefficientType &x, const Ring &ring) const |
|
{ |
|
int degree = Degree(ring); |
|
|
|
if (degree < 0) |
|
return ring.Identity(); |
|
|
|
CoefficientType result = m_coefficients[degree]; |
|
for (int j=degree-1; j>=0; j--) |
|
{ |
|
result = ring.Multiply(result, x); |
|
ring.Accumulate(result, m_coefficients[j]); |
|
} |
|
return result; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T>& PolynomialOver<T>::ShiftLeft(unsigned int n, const Ring &ring) |
|
{ |
|
unsigned int i = CoefficientCount(ring) + n; |
|
m_coefficients.resize(i, ring.Identity()); |
|
while (i > n) |
|
{ |
|
i--; |
|
m_coefficients[i] = m_coefficients[i-n]; |
|
} |
|
while (i) |
|
{ |
|
i--; |
|
m_coefficients[i] = ring.Identity(); |
|
} |
|
return *this; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T>& PolynomialOver<T>::ShiftRight(unsigned int n, const Ring &ring) |
|
{ |
|
unsigned int count = CoefficientCount(ring); |
|
if (count > n) |
|
{ |
|
for (unsigned int i=0; i<count-n; i++) |
|
m_coefficients[i] = m_coefficients[i+n]; |
|
m_coefficients.resize(count-n, ring.Identity()); |
|
} |
|
else |
|
m_coefficients.resize(0, ring.Identity()); |
|
return *this; |
|
} |
|
|
|
template <class T> |
|
void PolynomialOver<T>::SetCoefficient(unsigned int i, const CoefficientType &value, const Ring &ring) |
|
{ |
|
if (i >= m_coefficients.size()) |
|
m_coefficients.resize(i+1, ring.Identity()); |
|
m_coefficients[i] = value; |
|
} |
|
|
|
template <class T> |
|
void PolynomialOver<T>::Negate(const Ring &ring) |
|
{ |
|
unsigned int count = CoefficientCount(ring); |
|
for (unsigned int i=0; i<count; i++) |
|
m_coefficients[i] = ring.Inverse(m_coefficients[i]); |
|
} |
|
|
|
template <class T> |
|
void PolynomialOver<T>::swap(PolynomialOver<T> &t) |
|
{ |
|
m_coefficients.swap(t.m_coefficients); |
|
} |
|
|
|
template <class T> |
|
bool PolynomialOver<T>::Equals(const PolynomialOver<T>& t, const Ring &ring) const |
|
{ |
|
unsigned int count = CoefficientCount(ring); |
|
|
|
if (count != t.CoefficientCount(ring)) |
|
return false; |
|
|
|
for (unsigned int i=0; i<count; i++) |
|
if (!ring.Equal(m_coefficients[i], t.m_coefficients[i])) |
|
return false; |
|
|
|
return true; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T> PolynomialOver<T>::Plus(const PolynomialOver<T>& t, const Ring &ring) const |
|
{ |
|
unsigned int i; |
|
unsigned int count = CoefficientCount(ring); |
|
unsigned int tCount = t.CoefficientCount(ring); |
|
|
|
if (count > tCount) |
|
{ |
|
PolynomialOver<T> result(ring, count); |
|
|
|
for (i=0; i<tCount; i++) |
|
result.m_coefficients[i] = ring.Add(m_coefficients[i], t.m_coefficients[i]); |
|
for (; i<count; i++) |
|
result.m_coefficients[i] = m_coefficients[i]; |
|
|
|
return result; |
|
} |
|
else |
|
{ |
|
PolynomialOver<T> result(ring, tCount); |
|
|
|
for (i=0; i<count; i++) |
|
result.m_coefficients[i] = ring.Add(m_coefficients[i], t.m_coefficients[i]); |
|
for (; i<tCount; i++) |
|
result.m_coefficients[i] = t.m_coefficients[i]; |
|
|
|
return result; |
|
} |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T> PolynomialOver<T>::Minus(const PolynomialOver<T>& t, const Ring &ring) const |
|
{ |
|
unsigned int i; |
|
unsigned int count = CoefficientCount(ring); |
|
unsigned int tCount = t.CoefficientCount(ring); |
|
|
|
if (count > tCount) |
|
{ |
|
PolynomialOver<T> result(ring, count); |
|
|
|
for (i=0; i<tCount; i++) |
|
result.m_coefficients[i] = ring.Subtract(m_coefficients[i], t.m_coefficients[i]); |
|
for (; i<count; i++) |
|
result.m_coefficients[i] = m_coefficients[i]; |
|
|
|
return result; |
|
} |
|
else |
|
{ |
|
PolynomialOver<T> result(ring, tCount); |
|
|
|
for (i=0; i<count; i++) |
|
result.m_coefficients[i] = ring.Subtract(m_coefficients[i], t.m_coefficients[i]); |
|
for (; i<tCount; i++) |
|
result.m_coefficients[i] = ring.Inverse(t.m_coefficients[i]); |
|
|
|
return result; |
|
} |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T> PolynomialOver<T>::Inverse(const Ring &ring) const |
|
{ |
|
unsigned int count = CoefficientCount(ring); |
|
PolynomialOver<T> result(ring, count); |
|
|
|
for (unsigned int i=0; i<count; i++) |
|
result.m_coefficients[i] = ring.Inverse(m_coefficients[i]); |
|
|
|
return result; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T> PolynomialOver<T>::Times(const PolynomialOver<T>& t, const Ring &ring) const |
|
{ |
|
if (IsZero(ring) || t.IsZero(ring)) |
|
return PolynomialOver<T>(); |
|
|
|
unsigned int count1 = CoefficientCount(ring), count2 = t.CoefficientCount(ring); |
|
PolynomialOver<T> result(ring, count1 + count2 - 1); |
|
|
|
for (unsigned int i=0; i<count1; i++) |
|
for (unsigned int j=0; j<count2; j++) |
|
ring.Accumulate(result.m_coefficients[i+j], ring.Multiply(m_coefficients[i], t.m_coefficients[j])); |
|
|
|
return result; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T> PolynomialOver<T>::DividedBy(const PolynomialOver<T>& t, const Ring &ring) const |
|
{ |
|
PolynomialOver<T> remainder, quotient; |
|
Divide(remainder, quotient, *this, t, ring); |
|
return quotient; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T> PolynomialOver<T>::Modulo(const PolynomialOver<T>& t, const Ring &ring) const |
|
{ |
|
PolynomialOver<T> remainder, quotient; |
|
Divide(remainder, quotient, *this, t, ring); |
|
return remainder; |
|
} |
|
|
|
template <class T> |
|
PolynomialOver<T> PolynomialOver<T>::MultiplicativeInverse(const Ring &ring) const |
|
{ |
|
return Degree(ring)==0 ? ring.MultiplicativeInverse(m_coefficients[0]) : ring.Identity(); |
|
} |
|
|
|
template <class T> |
|
bool PolynomialOver<T>::IsUnit(const Ring &ring) const |
|
{ |
|
return Degree(ring)==0 && ring.IsUnit(m_coefficients[0]); |
|
} |
|
|
|
template <class T> |
|
std::istream& PolynomialOver<T>::Input(std::istream &in, const Ring &ring) |
|
{ |
|
char c; |
|
unsigned int length = 0; |
|
SecBlock<char> str(length + 16); |
|
bool paren = false; |
|
|
|
std::ws(in); |
|
|
|
if (in.peek() == '(') |
|
{ |
|
paren = true; |
|
in.get(); |
|
} |
|
|
|
do |
|
{ |
|
in.read(&c, 1); |
|
str[length++] = c; |
|
if (length >= str.size()) |
|
str.Grow(length + 16); |
|
} |
|
// if we started with a left paren, then read until we find a right paren, |
|
// otherwise read until the end of the line |
|
while (in && ((paren && c != ')') || (!paren && c != '\n'))); |
|
|
|
str[length-1] = '\0'; |
|
*this = PolynomialOver<T>(str, ring); |
|
|
|
return in; |
|
} |
|
|
|
template <class T> |
|
std::ostream& PolynomialOver<T>::Output(std::ostream &out, const Ring &ring) const |
|
{ |
|
unsigned int i = CoefficientCount(ring); |
|
if (i) |
|
{ |
|
bool firstTerm = true; |
|
|
|
while (i--) |
|
{ |
|
if (m_coefficients[i] != ring.Identity()) |
|
{ |
|
if (firstTerm) |
|
{ |
|
firstTerm = false; |
|
if (!i || !ring.Equal(m_coefficients[i], ring.MultiplicativeIdentity())) |
|
out << m_coefficients[i]; |
|
} |
|
else |
|
{ |
|
CoefficientType inverse = ring.Inverse(m_coefficients[i]); |
|
std::ostringstream pstr, nstr; |
|
|
|
pstr << m_coefficients[i]; |
|
nstr << inverse; |
|
|
|
if (pstr.str().size() <= nstr.str().size()) |
|
{ |
|
out << " + "; |
|
if (!i || !ring.Equal(m_coefficients[i], ring.MultiplicativeIdentity())) |
|
out << m_coefficients[i]; |
|
} |
|
else |
|
{ |
|
out << " - "; |
|
if (!i || !ring.Equal(inverse, ring.MultiplicativeIdentity())) |
|
out << inverse; |
|
} |
|
} |
|
|
|
switch (i) |
|
{ |
|
case 0: |
|
break; |
|
case 1: |
|
out << "x"; |
|
break; |
|
default: |
|
out << "x^" << i; |
|
} |
|
} |
|
} |
|
} |
|
else |
|
{ |
|
out << ring.Identity(); |
|
} |
|
return out; |
|
} |
|
|
|
template <class T> |
|
void PolynomialOver<T>::Divide(PolynomialOver<T> &r, PolynomialOver<T> &q, const PolynomialOver<T> &a, const PolynomialOver<T> &d, const Ring &ring) |
|
{ |
|
unsigned int i = a.CoefficientCount(ring); |
|
const int dDegree = d.Degree(ring); |
|
|
|
if (dDegree < 0) |
|
throw DivideByZero(); |
|
|
|
r = a; |
|
q.m_coefficients.resize(STDMAX(0, int(i - dDegree))); |
|
|
|
while (i > (unsigned int)dDegree) |
|
{ |
|
--i; |
|
q.m_coefficients[i-dDegree] = ring.Divide(r.m_coefficients[i], d.m_coefficients[dDegree]); |
|
for (int j=0; j<=dDegree; j++) |
|
ring.Reduce(r.m_coefficients[i-dDegree+j], ring.Multiply(q.m_coefficients[i-dDegree], d.m_coefficients[j])); |
|
} |
|
|
|
r.CoefficientCount(ring); // resize r.m_coefficients |
|
} |
|
|
|
// ******************************************************** |
|
|
|
// helper function for Interpolate() and InterpolateAt() |
|
template <class T> |
|
void RingOfPolynomialsOver<T>::CalculateAlpha(std::vector<CoefficientType> &alpha, const CoefficientType x[], const CoefficientType y[], unsigned int n) const |
|
{ |
|
for (unsigned int j=0; j<n; ++j) |
|
alpha[j] = y[j]; |
|
|
|
for (unsigned int k=1; k<n; ++k) |
|
{ |
|
for (unsigned int j=n-1; j>=k; --j) |
|
{ |
|
m_ring.Reduce(alpha[j], alpha[j-1]); |
|
|
|
CoefficientType d = m_ring.Subtract(x[j], x[j-k]); |
|
if (!m_ring.IsUnit(d)) |
|
throw InterpolationFailed(); |
|
alpha[j] = m_ring.Divide(alpha[j], d); |
|
} |
|
} |
|
} |
|
|
|
template <class T> |
|
typename RingOfPolynomialsOver<T>::Element RingOfPolynomialsOver<T>::Interpolate(const CoefficientType x[], const CoefficientType y[], unsigned int n) const |
|
{ |
|
assert(n > 0); |
|
|
|
std::vector<CoefficientType> alpha(n); |
|
CalculateAlpha(alpha, x, y, n); |
|
|
|
std::vector<CoefficientType> coefficients((size_t)n, m_ring.Identity()); |
|
coefficients[0] = alpha[n-1]; |
|
|
|
for (int j=n-2; j>=0; --j) |
|
{ |
|
for (unsigned int i=n-j-1; i>0; i--) |
|
coefficients[i] = m_ring.Subtract(coefficients[i-1], m_ring.Multiply(coefficients[i], x[j])); |
|
|
|
coefficients[0] = m_ring.Subtract(alpha[j], m_ring.Multiply(coefficients[0], x[j])); |
|
} |
|
|
|
return PolynomialOver<T>(coefficients.begin(), coefficients.end()); |
|
} |
|
|
|
template <class T> |
|
typename RingOfPolynomialsOver<T>::CoefficientType RingOfPolynomialsOver<T>::InterpolateAt(const CoefficientType &position, const CoefficientType x[], const CoefficientType y[], unsigned int n) const |
|
{ |
|
assert(n > 0); |
|
|
|
std::vector<CoefficientType> alpha(n); |
|
CalculateAlpha(alpha, x, y, n); |
|
|
|
CoefficientType result = alpha[n-1]; |
|
for (int j=n-2; j>=0; --j) |
|
{ |
|
result = m_ring.Multiply(result, m_ring.Subtract(position, x[j])); |
|
m_ring.Accumulate(result, alpha[j]); |
|
} |
|
return result; |
|
} |
|
|
|
template <class Ring, class Element> |
|
void PrepareBulkPolynomialInterpolation(const Ring &ring, Element *w, const Element x[], unsigned int n) |
|
{ |
|
for (unsigned int i=0; i<n; i++) |
|
{ |
|
Element t = ring.MultiplicativeIdentity(); |
|
for (unsigned int j=0; j<n; j++) |
|
if (i != j) |
|
t = ring.Multiply(t, ring.Subtract(x[i], x[j])); |
|
w[i] = ring.MultiplicativeInverse(t); |
|
} |
|
} |
|
|
|
template <class Ring, class Element> |
|
void PrepareBulkPolynomialInterpolationAt(const Ring &ring, Element *v, const Element &position, const Element x[], const Element w[], unsigned int n) |
|
{ |
|
assert(n > 0); |
|
|
|
std::vector<Element> a(2*n-1); |
|
unsigned int i; |
|
|
|
for (i=0; i<n; i++) |
|
a[n-1+i] = ring.Subtract(position, x[i]); |
|
|
|
for (i=n-1; i>1; i--) |
|
a[i-1] = ring.Multiply(a[2*i], a[2*i-1]); |
|
|
|
a[0] = ring.MultiplicativeIdentity(); |
|
|
|
for (i=0; i<n-1; i++) |
|
{ |
|
std::swap(a[2*i+1], a[2*i+2]); |
|
a[2*i+1] = ring.Multiply(a[i], a[2*i+1]); |
|
a[2*i+2] = ring.Multiply(a[i], a[2*i+2]); |
|
} |
|
|
|
for (i=0; i<n; i++) |
|
v[i] = ring.Multiply(a[n-1+i], w[i]); |
|
} |
|
|
|
template <class Ring, class Element> |
|
Element BulkPolynomialInterpolateAt(const Ring &ring, const Element y[], const Element v[], unsigned int n) |
|
{ |
|
Element result = ring.Identity(); |
|
for (unsigned int i=0; i<n; i++) |
|
ring.Accumulate(result, ring.Multiply(y[i], v[i])); |
|
return result; |
|
} |
|
|
|
// ******************************************************** |
|
|
|
template <class T, int instance> |
|
const PolynomialOverFixedRing<T, instance> &PolynomialOverFixedRing<T, instance>::Zero() |
|
{ |
|
return Singleton<ThisType>().Ref(); |
|
} |
|
|
|
template <class T, int instance> |
|
const PolynomialOverFixedRing<T, instance> &PolynomialOverFixedRing<T, instance>::One() |
|
{ |
|
return Singleton<ThisType, NewOnePolynomial>().Ref(); |
|
} |
|
|
|
NAMESPACE_END
|
|
|