1 // polynomi.h - originally written and placed in the public domain by Wei Dai
2 
3 /// \file polynomi.h
4 /// \brief Classes for polynomial basis and operations
5 
6 #ifndef CRYPTOPP_POLYNOMI_H
7 #define CRYPTOPP_POLYNOMI_H
8 
9 #include "cryptlib.h"
10 #include "secblock.h"
11 #include "algebra.h"
12 #include "misc.h"
13 
14 #include <iosfwd>
15 #include <vector>
16 
NAMESPACE_BEGIN(CryptoPP)17 NAMESPACE_BEGIN(CryptoPP)
18 
19 /// represents single-variable polynomials over arbitrary rings
20 /*!	\nosubgrouping */
21 template <class T> class PolynomialOver
22 {
23 public:
24 	/// \name ENUMS, EXCEPTIONS, and TYPEDEFS
25 	//@{
26 		/// division by zero exception
27 		class DivideByZero : public Exception
28 		{
29 		public:
30 			DivideByZero() : Exception(OTHER_ERROR, "PolynomialOver<T>: division by zero") {}
31 		};
32 
33 		/// specify the distribution for randomization functions
34 		class RandomizationParameter
35 		{
36 		public:
37 			RandomizationParameter(unsigned int coefficientCount, const typename T::RandomizationParameter &coefficientParameter )
38 				: m_coefficientCount(coefficientCount), m_coefficientParameter(coefficientParameter) {}
39 
40 		private:
41 			unsigned int m_coefficientCount;
42 			typename T::RandomizationParameter m_coefficientParameter;
43 			friend class PolynomialOver<T>;
44 		};
45 
46 		typedef T Ring;
47 		typedef typename T::Element CoefficientType;
48 	//@}
49 
50 	/// \name CREATORS
51 	//@{
52 		/// creates the zero polynomial
53 		PolynomialOver() {}
54 
55 		///
56 		PolynomialOver(const Ring &ring, unsigned int count)
57 			: m_coefficients((size_t)count, ring.Identity()) {}
58 
59 		/// copy constructor
60 		PolynomialOver(const PolynomialOver<Ring> &t)
61 			: m_coefficients(t.m_coefficients.size()) {*this = t;}
62 
63 		/// construct constant polynomial
64 		PolynomialOver(const CoefficientType &element)
65 			: m_coefficients(1, element) {}
66 
67 		/// construct polynomial with specified coefficients, starting from coefficient of x^0
68 		template <typename Iterator> PolynomialOver(Iterator begin, Iterator end)
69 			: m_coefficients(begin, end) {}
70 
71 		/// convert from string
72 		PolynomialOver(const char *str, const Ring &ring) {FromStr(str, ring);}
73 
74 		/// convert from big-endian byte array
75 		PolynomialOver(const byte *encodedPolynomialOver, unsigned int byteCount);
76 
77 		/// convert from Basic Encoding Rules encoded byte array
78 		explicit PolynomialOver(const byte *BEREncodedPolynomialOver);
79 
80 		/// convert from BER encoded byte array stored in a BufferedTransformation object
81 		explicit PolynomialOver(BufferedTransformation &bt);
82 
83 		/// create a random PolynomialOver<T>
84 		PolynomialOver(RandomNumberGenerator &rng, const RandomizationParameter &parameter, const Ring &ring)
85 			{Randomize(rng, parameter, ring);}
86 	//@}
87 
88 	/// \name ACCESSORS
89 	//@{
90 		/// the zero polynomial will return a degree of -1
91 		int Degree(const Ring &ring) const {return int(CoefficientCount(ring))-1;}
92 		///
93 		unsigned int CoefficientCount(const Ring &ring) const;
94 		/// return coefficient for x^i
95 		CoefficientType GetCoefficient(unsigned int i, const Ring &ring) const;
96 	//@}
97 
98 	/// \name MANIPULATORS
99 	//@{
100 		///
101 		PolynomialOver<Ring>&  operator=(const PolynomialOver<Ring>& t);
102 
103 		///
104 		void Randomize(RandomNumberGenerator &rng, const RandomizationParameter &parameter, const Ring &ring);
105 
106 		/// set the coefficient for x^i to value
107 		void SetCoefficient(unsigned int i, const CoefficientType &value, const Ring &ring);
108 
109 		///
110 		void Negate(const Ring &ring);
111 
112 		///
113 		void swap(PolynomialOver<Ring> &t);
114 	//@}
115 
116 
117 	/// \name BASIC ARITHMETIC ON POLYNOMIALS
118 	//@{
119 		bool Equals(const PolynomialOver<Ring> &t, const Ring &ring) const;
120 		bool IsZero(const Ring &ring) const {return CoefficientCount(ring)==0;}
121 
122 		PolynomialOver<Ring> Plus(const PolynomialOver<Ring>& t, const Ring &ring) const;
123 		PolynomialOver<Ring> Minus(const PolynomialOver<Ring>& t, const Ring &ring) const;
124 		PolynomialOver<Ring> Inverse(const Ring &ring) const;
125 
126 		PolynomialOver<Ring> Times(const PolynomialOver<Ring>& t, const Ring &ring) const;
127 		PolynomialOver<Ring> DividedBy(const PolynomialOver<Ring>& t, const Ring &ring) const;
128 		PolynomialOver<Ring> Modulo(const PolynomialOver<Ring>& t, const Ring &ring) const;
129 		PolynomialOver<Ring> MultiplicativeInverse(const Ring &ring) const;
130 		bool IsUnit(const Ring &ring) const;
131 
132 		PolynomialOver<Ring>& Accumulate(const PolynomialOver<Ring>& t, const Ring &ring);
133 		PolynomialOver<Ring>& Reduce(const PolynomialOver<Ring>& t, const Ring &ring);
134 
135 		///
136 		PolynomialOver<Ring> Doubled(const Ring &ring) const {return Plus(*this, ring);}
137 		///
138 		PolynomialOver<Ring> Squared(const Ring &ring) const {return Times(*this, ring);}
139 
140 		CoefficientType EvaluateAt(const CoefficientType &x, const Ring &ring) const;
141 
142 		PolynomialOver<Ring>& ShiftLeft(unsigned int n, const Ring &ring);
143 		PolynomialOver<Ring>& ShiftRight(unsigned int n, const Ring &ring);
144 
145 		/// calculate r and q such that (a == d*q + r) && (0 <= degree of r < degree of d)
146 		static void Divide(PolynomialOver<Ring> &r, PolynomialOver<Ring> &q, const PolynomialOver<Ring> &a, const PolynomialOver<Ring> &d, const Ring &ring);
147 	//@}
148 
149 	/// \name INPUT/OUTPUT
150 	//@{
151 		std::istream& Input(std::istream &in, const Ring &ring);
152 		std::ostream& Output(std::ostream &out, const Ring &ring) const;
153 	//@}
154 
155 private:
156 	void FromStr(const char *str, const Ring &ring);
157 
158 	std::vector<CoefficientType> m_coefficients;
159 };
160 
161 /// Polynomials over a fixed ring
162 /*! Having a fixed ring allows overloaded operators */
163 template <class T, int instance> class PolynomialOverFixedRing : private PolynomialOver<T>
164 {
165 	typedef PolynomialOver<T> B;
166 	typedef PolynomialOverFixedRing<T, instance> ThisType;
167 
168 public:
169 	typedef T Ring;
170 	typedef typename T::Element CoefficientType;
171 	typedef typename B::DivideByZero DivideByZero;
172 	typedef typename B::RandomizationParameter RandomizationParameter;
173 
174 	/// \name CREATORS
175 	//@{
176 		/// creates the zero polynomial
B(ms_fixedRing,count)177 		PolynomialOverFixedRing(unsigned int count = 0) : B(ms_fixedRing, count) {}
178 
179 		/// copy constructor
PolynomialOverFixedRing(const ThisType & t)180 		PolynomialOverFixedRing(const ThisType &t) : B(t) {}
181 
PolynomialOverFixedRing(const B & t)182 		explicit PolynomialOverFixedRing(const B &t) : B(t) {}
183 
184 		/// construct constant polynomial
PolynomialOverFixedRing(const CoefficientType & element)185 		PolynomialOverFixedRing(const CoefficientType &element) : B(element) {}
186 
187 		/// construct polynomial with specified coefficients, starting from coefficient of x^0
PolynomialOverFixedRing(Iterator first,Iterator last)188 		template <typename Iterator> PolynomialOverFixedRing(Iterator first, Iterator last)
189 			: B(first, last) {}
190 
191 		/// convert from string
PolynomialOverFixedRing(const char * str)192 		explicit PolynomialOverFixedRing(const char *str) : B(str, ms_fixedRing) {}
193 
194 		/// convert from big-endian byte array
PolynomialOverFixedRing(const byte * encodedPoly,unsigned int byteCount)195 		PolynomialOverFixedRing(const byte *encodedPoly, unsigned int byteCount) : B(encodedPoly, byteCount) {}
196 
197 		/// convert from Basic Encoding Rules encoded byte array
PolynomialOverFixedRing(const byte * BEREncodedPoly)198 		explicit PolynomialOverFixedRing(const byte *BEREncodedPoly) : B(BEREncodedPoly) {}
199 
200 		/// convert from BER encoded byte array stored in a BufferedTransformation object
PolynomialOverFixedRing(BufferedTransformation & bt)201 		explicit PolynomialOverFixedRing(BufferedTransformation &bt) : B(bt) {}
202 
203 		/// create a random PolynomialOverFixedRing
PolynomialOverFixedRing(RandomNumberGenerator & rng,const RandomizationParameter & parameter)204 		PolynomialOverFixedRing(RandomNumberGenerator &rng, const RandomizationParameter &parameter) : B(rng, parameter, ms_fixedRing) {}
205 
206 		static const ThisType &Zero();
207 		static const ThisType &One();
208 	//@}
209 
210 	/// \name ACCESSORS
211 	//@{
212 		/// the zero polynomial will return a degree of -1
Degree()213 		int Degree() const {return B::Degree(ms_fixedRing);}
214 		/// degree + 1
CoefficientCount()215 		unsigned int CoefficientCount() const {return B::CoefficientCount(ms_fixedRing);}
216 		/// return coefficient for x^i
GetCoefficient(unsigned int i)217 		CoefficientType GetCoefficient(unsigned int i) const {return B::GetCoefficient(i, ms_fixedRing);}
218 		/// return coefficient for x^i
219 		CoefficientType operator[](unsigned int i) const {return B::GetCoefficient(i, ms_fixedRing);}
220 	//@}
221 
222 	/// \name MANIPULATORS
223 	//@{
224 		///
225 		ThisType&  operator=(const ThisType& t) {B::operator=(t); return *this;}
226 		///
227 		ThisType&  operator+=(const ThisType& t) {Accumulate(t, ms_fixedRing); return *this;}
228 		///
229 		ThisType&  operator-=(const ThisType& t) {Reduce(t, ms_fixedRing); return *this;}
230 		///
231 		ThisType&  operator*=(const ThisType& t) {return *this = *this*t;}
232 		///
233 		ThisType&  operator/=(const ThisType& t) {return *this = *this/t;}
234 		///
235 		ThisType&  operator%=(const ThisType& t) {return *this = *this%t;}
236 
237 		///
238 		ThisType&  operator<<=(unsigned int n) {ShiftLeft(n, ms_fixedRing); return *this;}
239 		///
240 		ThisType&  operator>>=(unsigned int n) {ShiftRight(n, ms_fixedRing); return *this;}
241 
242 		/// set the coefficient for x^i to value
SetCoefficient(unsigned int i,const CoefficientType & value)243 		void SetCoefficient(unsigned int i, const CoefficientType &value) {B::SetCoefficient(i, value, ms_fixedRing);}
244 
245 		///
Randomize(RandomNumberGenerator & rng,const RandomizationParameter & parameter)246 		void Randomize(RandomNumberGenerator &rng, const RandomizationParameter &parameter) {B::Randomize(rng, parameter, ms_fixedRing);}
247 
248 		///
Negate()249 		void Negate() {B::Negate(ms_fixedRing);}
250 
swap(ThisType & t)251 		void swap(ThisType &t) {B::swap(t);}
252 	//@}
253 
254 	/// \name UNARY OPERATORS
255 	//@{
256 		///
257 		bool operator!() const {return CoefficientCount()==0;}
258 		///
259 		ThisType operator+() const {return *this;}
260 		///
261 		ThisType operator-() const {return ThisType(Inverse(ms_fixedRing));}
262 	//@}
263 
264 	/// \name BINARY OPERATORS
265 	//@{
266 		///
267 		friend ThisType operator>>(ThisType a, unsigned int n)	{return ThisType(a>>=n);}
268 		///
269 		friend ThisType operator<<(ThisType a, unsigned int n)	{return ThisType(a<<=n);}
270 	//@}
271 
272 	/// \name OTHER ARITHMETIC FUNCTIONS
273 	//@{
274 		///
MultiplicativeInverse()275 		ThisType MultiplicativeInverse() const {return ThisType(B::MultiplicativeInverse(ms_fixedRing));}
276 		///
IsUnit()277 		bool IsUnit() const {return B::IsUnit(ms_fixedRing);}
278 
279 		///
Doubled()280 		ThisType Doubled() const {return ThisType(B::Doubled(ms_fixedRing));}
281 		///
Squared()282 		ThisType Squared() const {return ThisType(B::Squared(ms_fixedRing));}
283 
EvaluateAt(const CoefficientType & x)284 		CoefficientType EvaluateAt(const CoefficientType &x) const {return B::EvaluateAt(x, ms_fixedRing);}
285 
286 		/// calculate r and q such that (a == d*q + r) && (0 <= r < abs(d))
Divide(ThisType & r,ThisType & q,const ThisType & a,const ThisType & d)287 		static void Divide(ThisType &r, ThisType &q, const ThisType &a, const ThisType &d)
288 			{B::Divide(r, q, a, d, ms_fixedRing);}
289 	//@}
290 
291 	/// \name INPUT/OUTPUT
292 	//@{
293 		///
294 		friend std::istream& operator>>(std::istream& in, ThisType &a)
295 			{return a.Input(in, ms_fixedRing);}
296 		///
297 		friend std::ostream& operator<<(std::ostream& out, const ThisType &a)
298 			{return a.Output(out, ms_fixedRing);}
299 	//@}
300 
301 private:
302 	struct NewOnePolynomial
303 	{
operatorNewOnePolynomial304 		ThisType * operator()() const
305 		{
306 			return new ThisType(ms_fixedRing.MultiplicativeIdentity());
307 		}
308 	};
309 
310 	static const Ring ms_fixedRing;
311 };
312 
313 /// Ring of polynomials over another ring
314 template <class T> class RingOfPolynomialsOver : public AbstractEuclideanDomain<PolynomialOver<T> >
315 {
316 public:
317 	typedef T CoefficientRing;
318 	typedef PolynomialOver<T> Element;
319 	typedef typename Element::CoefficientType CoefficientType;
320 	typedef typename Element::RandomizationParameter RandomizationParameter;
321 
RingOfPolynomialsOver(const CoefficientRing & ring)322 	RingOfPolynomialsOver(const CoefficientRing &ring) : m_ring(ring) {}
323 
RandomElement(RandomNumberGenerator & rng,const RandomizationParameter & parameter)324 	Element RandomElement(RandomNumberGenerator &rng, const RandomizationParameter &parameter)
325 		{return Element(rng, parameter, m_ring);}
326 
Equal(const Element & a,const Element & b)327 	bool Equal(const Element &a, const Element &b) const
328 		{return a.Equals(b, m_ring);}
329 
Identity()330 	const Element& Identity() const
331 		{return this->result = m_ring.Identity();}
332 
Add(const Element & a,const Element & b)333 	const Element& Add(const Element &a, const Element &b) const
334 		{return this->result = a.Plus(b, m_ring);}
335 
Accumulate(Element & a,const Element & b)336 	Element& Accumulate(Element &a, const Element &b) const
337 		{a.Accumulate(b, m_ring); return a;}
338 
Inverse(const Element & a)339 	const Element& Inverse(const Element &a) const
340 		{return this->result = a.Inverse(m_ring);}
341 
Subtract(const Element & a,const Element & b)342 	const Element& Subtract(const Element &a, const Element &b) const
343 		{return this->result = a.Minus(b, m_ring);}
344 
Reduce(Element & a,const Element & b)345 	Element& Reduce(Element &a, const Element &b) const
346 		{return a.Reduce(b, m_ring);}
347 
Double(const Element & a)348 	const Element& Double(const Element &a) const
349 		{return this->result = a.Doubled(m_ring);}
350 
MultiplicativeIdentity()351 	const Element& MultiplicativeIdentity() const
352 		{return this->result = m_ring.MultiplicativeIdentity();}
353 
Multiply(const Element & a,const Element & b)354 	const Element& Multiply(const Element &a, const Element &b) const
355 		{return this->result = a.Times(b, m_ring);}
356 
Square(const Element & a)357 	const Element& Square(const Element &a) const
358 		{return this->result = a.Squared(m_ring);}
359 
IsUnit(const Element & a)360 	bool IsUnit(const Element &a) const
361 		{return a.IsUnit(m_ring);}
362 
MultiplicativeInverse(const Element & a)363 	const Element& MultiplicativeInverse(const Element &a) const
364 		{return this->result = a.MultiplicativeInverse(m_ring);}
365 
Divide(const Element & a,const Element & b)366 	const Element& Divide(const Element &a, const Element &b) const
367 		{return this->result = a.DividedBy(b, m_ring);}
368 
Mod(const Element & a,const Element & b)369 	const Element& Mod(const Element &a, const Element &b) const
370 		{return this->result = a.Modulo(b, m_ring);}
371 
DivisionAlgorithm(Element & r,Element & q,const Element & a,const Element & d)372 	void DivisionAlgorithm(Element &r, Element &q, const Element &a, const Element &d) const
373 		{Element::Divide(r, q, a, d, m_ring);}
374 
375 	class InterpolationFailed : public Exception
376 	{
377 	public:
InterpolationFailed()378 		InterpolationFailed() : Exception(OTHER_ERROR, "RingOfPolynomialsOver<T>: interpolation failed") {}
379 	};
380 
381 	Element Interpolate(const CoefficientType x[], const CoefficientType y[], unsigned int n) const;
382 
383 	// a faster version of Interpolate(x, y, n).EvaluateAt(position)
384 	CoefficientType InterpolateAt(const CoefficientType &position, const CoefficientType x[], const CoefficientType y[], unsigned int n) const;
385 /*
386 	void PrepareBulkInterpolation(CoefficientType *w, const CoefficientType x[], unsigned int n) const;
387 	void PrepareBulkInterpolationAt(CoefficientType *v, const CoefficientType &position, const CoefficientType x[], const CoefficientType w[], unsigned int n) const;
388 	CoefficientType BulkInterpolateAt(const CoefficientType y[], const CoefficientType v[], unsigned int n) const;
389 */
390 protected:
391 	void CalculateAlpha(std::vector<CoefficientType> &alpha, const CoefficientType x[], const CoefficientType y[], unsigned int n) const;
392 
393 	CoefficientRing m_ring;
394 };
395 
396 template <class Ring, class Element>
397 void PrepareBulkPolynomialInterpolation(const Ring &ring, Element *w, const Element x[], unsigned int n);
398 template <class Ring, class Element>
399 void PrepareBulkPolynomialInterpolationAt(const Ring &ring, Element *v, const Element &position, const Element x[], const Element w[], unsigned int n);
400 template <class Ring, class Element>
401 Element BulkPolynomialInterpolateAt(const Ring &ring, const Element y[], const Element v[], unsigned int n);
402 
403 ///
404 template <class T, int instance>
405 inline bool operator==(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
406 	{return a.Equals(b, a.ms_fixedRing);}
407 ///
408 template <class T, int instance>
409 inline bool operator!=(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
410 	{return !(a==b);}
411 
412 ///
413 template <class T, int instance>
414 inline bool operator> (const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
415 	{return a.Degree() > b.Degree();}
416 ///
417 template <class T, int instance>
418 inline bool operator>=(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
419 	{return a.Degree() >= b.Degree();}
420 ///
421 template <class T, int instance>
422 inline bool operator< (const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
423 	{return a.Degree() < b.Degree();}
424 ///
425 template <class T, int instance>
426 inline bool operator<=(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
427 	{return a.Degree() <= b.Degree();}
428 
429 ///
430 template <class T, int instance>
431 inline CryptoPP::PolynomialOverFixedRing<T, instance> operator+(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
432 	{return CryptoPP::PolynomialOverFixedRing<T, instance>(a.Plus(b, a.ms_fixedRing));}
433 ///
434 template <class T, int instance>
435 inline CryptoPP::PolynomialOverFixedRing<T, instance> operator-(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
436 	{return CryptoPP::PolynomialOverFixedRing<T, instance>(a.Minus(b, a.ms_fixedRing));}
437 ///
438 template <class T, int instance>
439 inline CryptoPP::PolynomialOverFixedRing<T, instance> operator*(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
440 	{return CryptoPP::PolynomialOverFixedRing<T, instance>(a.Times(b, a.ms_fixedRing));}
441 ///
442 template <class T, int instance>
443 inline CryptoPP::PolynomialOverFixedRing<T, instance> operator/(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
444 	{return CryptoPP::PolynomialOverFixedRing<T, instance>(a.DividedBy(b, a.ms_fixedRing));}
445 ///
446 template <class T, int instance>
447 inline CryptoPP::PolynomialOverFixedRing<T, instance> operator%(const CryptoPP::PolynomialOverFixedRing<T, instance> &a, const CryptoPP::PolynomialOverFixedRing<T, instance> &b)
448 	{return CryptoPP::PolynomialOverFixedRing<T, instance>(a.Modulo(b, a.ms_fixedRing));}
449 
450 NAMESPACE_END
451 
NAMESPACE_BEGIN(std)452 NAMESPACE_BEGIN(std)
453 template<class T> inline void swap(CryptoPP::PolynomialOver<T> &a, CryptoPP::PolynomialOver<T> &b)
454 {
455 	a.swap(b);
456 }
swap(CryptoPP::PolynomialOverFixedRing<T,i> & a,CryptoPP::PolynomialOverFixedRing<T,i> & b)457 template<class T, int i> inline void swap(CryptoPP::PolynomialOverFixedRing<T,i> &a, CryptoPP::PolynomialOverFixedRing<T,i> &b)
458 {
459 	a.swap(b);
460 }
461 NAMESPACE_END
462 
463 #endif
464