Security Scol plugin
polynomi.cpp
1// polynomi.cpp - originally written and placed in the public domain by Wei Dai
2
3// Part of the code for polynomial evaluation and interpolation
4// originally came from Hal Finney's public domain secsplit.c.
5
6#include "pch.h"
7#include "polynomi.h"
8#include "secblock.h"
9
10#include <sstream>
11#include <iostream>
12
13NAMESPACE_BEGIN(CryptoPP)
14
15template <class T>
16void PolynomialOver<T>::Randomize(RandomNumberGenerator &rng, const RandomizationParameter &parameter, const Ring &ring)
17{
18 m_coefficients.resize(parameter.m_coefficientCount);
19 for (unsigned int i=0; i<m_coefficients.size(); ++i)
20 m_coefficients[i] = ring.RandomElement(rng, parameter.m_coefficientParameter);
21}
22
23template <class T>
24void PolynomialOver<T>::FromStr(const char *str, const Ring &ring)
25{
26 std::istringstream in((char *)str);
27 bool positive = true;
28 CoefficientType coef;
29 unsigned int power;
30
31 while (in)
32 {
33 std::ws(in);
34 if (in.peek() == 'x')
35 coef = ring.MultiplicativeIdentity();
36 else
37 in >> coef;
38
39 std::ws(in);
40 if (in.peek() == 'x')
41 {
42 in.get();
43 std::ws(in);
44 if (in.peek() == '^')
45 {
46 in.get();
47 in >> power;
48 }
49 else
50 power = 1;
51 }
52 else
53 power = 0;
54
55 if (!positive)
56 coef = ring.Inverse(coef);
57
58 SetCoefficient(power, coef, ring);
59
60 std::ws(in);
61 switch (in.get())
62 {
63 case '+':
64 positive = true;
65 break;
66 case '-':
67 positive = false;
68 break;
69 default:
70 return; // something's wrong with the input string
71 }
72 }
73}
74
75template <class T>
76unsigned int PolynomialOver<T>::CoefficientCount(const Ring &ring) const
77{
78 unsigned count = m_coefficients.size();
79 while (count && ring.Equal(m_coefficients[count-1], ring.Identity()))
80 count--;
81 const_cast<std::vector<CoefficientType> &>(m_coefficients).resize(count);
82 return count;
83}
84
85template <class T>
86typename PolynomialOver<T>::CoefficientType PolynomialOver<T>::GetCoefficient(unsigned int i, const Ring &ring) const
87{
88 return (i < m_coefficients.size()) ? m_coefficients[i] : ring.Identity();
89}
90
91template <class T>
92PolynomialOver<T>& PolynomialOver<T>::operator=(const PolynomialOver<T>& t)
93{
94 if (this != &t)
95 {
96 m_coefficients.resize(t.m_coefficients.size());
97 for (unsigned int i=0; i<m_coefficients.size(); i++)
98 m_coefficients[i] = t.m_coefficients[i];
99 }
100 return *this;
101}
102
103template <class T>
104PolynomialOver<T>& PolynomialOver<T>::Accumulate(const PolynomialOver<T>& t, const Ring &ring)
105{
106 unsigned int count = t.CoefficientCount(ring);
107
108 if (count > CoefficientCount(ring))
109 m_coefficients.resize(count, ring.Identity());
110
111 for (unsigned int i=0; i<count; i++)
112 ring.Accumulate(m_coefficients[i], t.GetCoefficient(i, ring));
113
114 return *this;
115}
116
117template <class T>
118PolynomialOver<T>& PolynomialOver<T>::Reduce(const PolynomialOver<T>& t, const Ring &ring)
119{
120 unsigned int count = t.CoefficientCount(ring);
121
122 if (count > CoefficientCount(ring))
123 m_coefficients.resize(count, ring.Identity());
124
125 for (unsigned int i=0; i<count; i++)
126 ring.Reduce(m_coefficients[i], t.GetCoefficient(i, ring));
127
128 return *this;
129}
130
131template <class T>
132typename PolynomialOver<T>::CoefficientType PolynomialOver<T>::EvaluateAt(const CoefficientType &x, const Ring &ring) const
133{
134 int degree = Degree(ring);
135
136 if (degree < 0)
137 return ring.Identity();
138
139 CoefficientType result = m_coefficients[degree];
140 for (int j=degree-1; j>=0; j--)
141 {
142 result = ring.Multiply(result, x);
143 ring.Accumulate(result, m_coefficients[j]);
144 }
145 return result;
146}
147
148template <class T>
149PolynomialOver<T>& PolynomialOver<T>::ShiftLeft(unsigned int n, const Ring &ring)
150{
151 unsigned int i = CoefficientCount(ring) + n;
152 m_coefficients.resize(i, ring.Identity());
153 while (i > n)
154 {
155 i--;
156 m_coefficients[i] = m_coefficients[i-n];
157 }
158 while (i)
159 {
160 i--;
161 m_coefficients[i] = ring.Identity();
162 }
163 return *this;
164}
165
166template <class T>
167PolynomialOver<T>& PolynomialOver<T>::ShiftRight(unsigned int n, const Ring &ring)
168{
169 unsigned int count = CoefficientCount(ring);
170 if (count > n)
171 {
172 for (unsigned int i=0; i<count-n; i++)
173 m_coefficients[i] = m_coefficients[i+n];
174 m_coefficients.resize(count-n, ring.Identity());
175 }
176 else
177 m_coefficients.resize(0, ring.Identity());
178 return *this;
179}
180
181template <class T>
182void PolynomialOver<T>::SetCoefficient(unsigned int i, const CoefficientType &value, const Ring &ring)
183{
184 if (i >= m_coefficients.size())
185 m_coefficients.resize(i+1, ring.Identity());
186 m_coefficients[i] = value;
187}
188
189template <class T>
190void PolynomialOver<T>::Negate(const Ring &ring)
191{
192 unsigned int count = CoefficientCount(ring);
193 for (unsigned int i=0; i<count; i++)
194 m_coefficients[i] = ring.Inverse(m_coefficients[i]);
195}
196
197template <class T>
198void PolynomialOver<T>::swap(PolynomialOver<T> &t)
199{
200 m_coefficients.swap(t.m_coefficients);
201}
202
203template <class T>
204bool PolynomialOver<T>::Equals(const PolynomialOver<T>& t, const Ring &ring) const
205{
206 unsigned int count = CoefficientCount(ring);
207
208 if (count != t.CoefficientCount(ring))
209 return false;
210
211 for (unsigned int i=0; i<count; i++)
212 if (!ring.Equal(m_coefficients[i], t.m_coefficients[i]))
213 return false;
214
215 return true;
216}
217
218template <class T>
219PolynomialOver<T> PolynomialOver<T>::Plus(const PolynomialOver<T>& t, const Ring &ring) const
220{
221 unsigned int i;
222 unsigned int count = CoefficientCount(ring);
223 unsigned int tCount = t.CoefficientCount(ring);
224
225 if (count > tCount)
226 {
227 PolynomialOver<T> result(ring, count);
228
229 for (i=0; i<tCount; i++)
230 result.m_coefficients[i] = ring.Add(m_coefficients[i], t.m_coefficients[i]);
231 for (; i<count; i++)
232 result.m_coefficients[i] = m_coefficients[i];
233
234 return result;
235 }
236 else
237 {
238 PolynomialOver<T> result(ring, tCount);
239
240 for (i=0; i<count; i++)
241 result.m_coefficients[i] = ring.Add(m_coefficients[i], t.m_coefficients[i]);
242 for (; i<tCount; i++)
243 result.m_coefficients[i] = t.m_coefficients[i];
244
245 return result;
246 }
247}
248
249template <class T>
250PolynomialOver<T> PolynomialOver<T>::Minus(const PolynomialOver<T>& t, const Ring &ring) const
251{
252 unsigned int i;
253 unsigned int count = CoefficientCount(ring);
254 unsigned int tCount = t.CoefficientCount(ring);
255
256 if (count > tCount)
257 {
258 PolynomialOver<T> result(ring, count);
259
260 for (i=0; i<tCount; i++)
261 result.m_coefficients[i] = ring.Subtract(m_coefficients[i], t.m_coefficients[i]);
262 for (; i<count; i++)
263 result.m_coefficients[i] = m_coefficients[i];
264
265 return result;
266 }
267 else
268 {
269 PolynomialOver<T> result(ring, tCount);
270
271 for (i=0; i<count; i++)
272 result.m_coefficients[i] = ring.Subtract(m_coefficients[i], t.m_coefficients[i]);
273 for (; i<tCount; i++)
274 result.m_coefficients[i] = ring.Inverse(t.m_coefficients[i]);
275
276 return result;
277 }
278}
279
280template <class T>
281PolynomialOver<T> PolynomialOver<T>::Inverse(const Ring &ring) const
282{
283 unsigned int count = CoefficientCount(ring);
284 PolynomialOver<T> result(ring, count);
285
286 for (unsigned int i=0; i<count; i++)
287 result.m_coefficients[i] = ring.Inverse(m_coefficients[i]);
288
289 return result;
290}
291
292template <class T>
293PolynomialOver<T> PolynomialOver<T>::Times(const PolynomialOver<T>& t, const Ring &ring) const
294{
295 if (IsZero(ring) || t.IsZero(ring))
296 return PolynomialOver<T>();
297
298 unsigned int count1 = CoefficientCount(ring), count2 = t.CoefficientCount(ring);
299 PolynomialOver<T> result(ring, count1 + count2 - 1);
300
301 for (unsigned int i=0; i<count1; i++)
302 for (unsigned int j=0; j<count2; j++)
303 ring.Accumulate(result.m_coefficients[i+j], ring.Multiply(m_coefficients[i], t.m_coefficients[j]));
304
305 return result;
306}
307
308template <class T>
309PolynomialOver<T> PolynomialOver<T>::DividedBy(const PolynomialOver<T>& t, const Ring &ring) const
310{
311 PolynomialOver<T> remainder, quotient;
312 Divide(remainder, quotient, *this, t, ring);
313 return quotient;
314}
315
316template <class T>
317PolynomialOver<T> PolynomialOver<T>::Modulo(const PolynomialOver<T>& t, const Ring &ring) const
318{
319 PolynomialOver<T> remainder, quotient;
320 Divide(remainder, quotient, *this, t, ring);
321 return remainder;
322}
323
324template <class T>
325PolynomialOver<T> PolynomialOver<T>::MultiplicativeInverse(const Ring &ring) const
326{
327 return Degree(ring)==0 ? ring.MultiplicativeInverse(m_coefficients[0]) : ring.Identity();
328}
329
330template <class T>
331bool PolynomialOver<T>::IsUnit(const Ring &ring) const
332{
333 return Degree(ring)==0 && ring.IsUnit(m_coefficients[0]);
334}
335
336template <class T>
337std::istream& PolynomialOver<T>::Input(std::istream &in, const Ring &ring)
338{
339 char c;
340 unsigned int length = 0;
341 SecBlock<char> str(length + 16);
342 bool paren = false;
343
344 std::ws(in);
345
346 if (in.peek() == '(')
347 {
348 paren = true;
349 in.get();
350 }
351
352 do
353 {
354 in.read(&c, 1);
355 str[length++] = c;
356 if (length >= str.size())
357 str.Grow(length + 16);
358 }
359 // if we started with a left paren, then read until we find a right paren,
360 // otherwise read until the end of the line
361 while (in && ((paren && c != ')') || (!paren && c != '\n')));
362
363 str[length-1] = '\0';
364 *this = PolynomialOver<T>(str, ring);
365
366 return in;
367}
368
369template <class T>
370std::ostream& PolynomialOver<T>::Output(std::ostream &out, const Ring &ring) const
371{
372 unsigned int i = CoefficientCount(ring);
373 if (i)
374 {
375 bool firstTerm = true;
376
377 while (i--)
378 {
379 if (m_coefficients[i] != ring.Identity())
380 {
381 if (firstTerm)
382 {
383 firstTerm = false;
384 if (!i || !ring.Equal(m_coefficients[i], ring.MultiplicativeIdentity()))
385 out << m_coefficients[i];
386 }
387 else
388 {
389 CoefficientType inverse = ring.Inverse(m_coefficients[i]);
390 std::ostringstream pstr, nstr;
391
392 pstr << m_coefficients[i];
393 nstr << inverse;
394
395 if (pstr.str().size() <= nstr.str().size())
396 {
397 out << " + ";
398 if (!i || !ring.Equal(m_coefficients[i], ring.MultiplicativeIdentity()))
399 out << m_coefficients[i];
400 }
401 else
402 {
403 out << " - ";
404 if (!i || !ring.Equal(inverse, ring.MultiplicativeIdentity()))
405 out << inverse;
406 }
407 }
408
409 switch (i)
410 {
411 case 0:
412 break;
413 case 1:
414 out << "x";
415 break;
416 default:
417 out << "x^" << i;
418 }
419 }
420 }
421 }
422 else
423 {
424 out << ring.Identity();
425 }
426 return out;
427}
428
429template <class T>
430void PolynomialOver<T>::Divide(PolynomialOver<T> &r, PolynomialOver<T> &q, const PolynomialOver<T> &a, const PolynomialOver<T> &d, const Ring &ring)
431{
432 unsigned int i = a.CoefficientCount(ring);
433 const int dDegree = d.Degree(ring);
434
435 if (dDegree < 0)
436 throw DivideByZero();
437
438 r = a;
439 q.m_coefficients.resize(STDMAX(0, int(i - dDegree)));
440
441 while (i > (unsigned int)dDegree)
442 {
443 --i;
444 q.m_coefficients[i-dDegree] = ring.Divide(r.m_coefficients[i], d.m_coefficients[dDegree]);
445 for (int j=0; j<=dDegree; j++)
446 ring.Reduce(r.m_coefficients[i-dDegree+j], ring.Multiply(q.m_coefficients[i-dDegree], d.m_coefficients[j]));
447 }
448
449 r.CoefficientCount(ring); // resize r.m_coefficients
450}
451
452// ********************************************************
453
454// helper function for Interpolate() and InterpolateAt()
455template <class T>
456void RingOfPolynomialsOver<T>::CalculateAlpha(std::vector<CoefficientType> &alpha, const CoefficientType x[], const CoefficientType y[], unsigned int n) const
457{
458 for (unsigned int j=0; j<n; ++j)
459 alpha[j] = y[j];
460
461 for (unsigned int k=1; k<n; ++k)
462 {
463 for (unsigned int j=n-1; j>=k; --j)
464 {
465 m_ring.Reduce(alpha[j], alpha[j-1]);
466
467 CoefficientType d = m_ring.Subtract(x[j], x[j-k]);
468 if (!m_ring.IsUnit(d))
469 throw InterpolationFailed();
470 alpha[j] = m_ring.Divide(alpha[j], d);
471 }
472 }
473}
474
475template <class T>
476typename RingOfPolynomialsOver<T>::Element RingOfPolynomialsOver<T>::Interpolate(const CoefficientType x[], const CoefficientType y[], unsigned int n) const
477{
478 CRYPTOPP_ASSERT(n > 0);
479
480 std::vector<CoefficientType> alpha(n);
481 CalculateAlpha(alpha, x, y, n);
482
483 std::vector<CoefficientType> coefficients((size_t)n, m_ring.Identity());
484 coefficients[0] = alpha[n-1];
485
486 for (int j=n-2; j>=0; --j)
487 {
488 for (unsigned int i=n-j-1; i>0; i--)
489 coefficients[i] = m_ring.Subtract(coefficients[i-1], m_ring.Multiply(coefficients[i], x[j]));
490
491 coefficients[0] = m_ring.Subtract(alpha[j], m_ring.Multiply(coefficients[0], x[j]));
492 }
493
494 return PolynomialOver<T>(coefficients.begin(), coefficients.end());
495}
496
497template <class T>
498typename RingOfPolynomialsOver<T>::CoefficientType RingOfPolynomialsOver<T>::InterpolateAt(const CoefficientType &position, const CoefficientType x[], const CoefficientType y[], unsigned int n) const
499{
500 CRYPTOPP_ASSERT(n > 0);
501
502 std::vector<CoefficientType> alpha(n);
503 CalculateAlpha(alpha, x, y, n);
504
505 CoefficientType result = alpha[n-1];
506 for (int j=n-2; j>=0; --j)
507 {
508 result = m_ring.Multiply(result, m_ring.Subtract(position, x[j]));
509 m_ring.Accumulate(result, alpha[j]);
510 }
511 return result;
512}
513
514template <class Ring, class Element>
515void PrepareBulkPolynomialInterpolation(const Ring &ring, Element *w, const Element x[], unsigned int n)
516{
517 for (unsigned int i=0; i<n; i++)
518 {
519 Element t = ring.MultiplicativeIdentity();
520 for (unsigned int j=0; j<n; j++)
521 if (i != j)
522 t = ring.Multiply(t, ring.Subtract(x[i], x[j]));
523 w[i] = ring.MultiplicativeInverse(t);
524 }
525}
526
527template <class Ring, class Element>
528void PrepareBulkPolynomialInterpolationAt(const Ring &ring, Element *v, const Element &position, const Element x[], const Element w[], unsigned int n)
529{
530 CRYPTOPP_ASSERT(n > 0);
531
532 std::vector<Element> a(2*n-1);
533 unsigned int i;
534
535 for (i=0; i<n; i++)
536 a[n-1+i] = ring.Subtract(position, x[i]);
537
538 for (i=n-1; i>1; i--)
539 a[i-1] = ring.Multiply(a[2*i], a[2*i-1]);
540
541 a[0] = ring.MultiplicativeIdentity();
542
543 for (i=0; i<n-1; i++)
544 {
545 std::swap(a[2*i+1], a[2*i+2]);
546 a[2*i+1] = ring.Multiply(a[i], a[2*i+1]);
547 a[2*i+2] = ring.Multiply(a[i], a[2*i+2]);
548 }
549
550 for (i=0; i<n; i++)
551 v[i] = ring.Multiply(a[n-1+i], w[i]);
552}
553
554template <class Ring, class Element>
555Element BulkPolynomialInterpolateAt(const Ring &ring, const Element y[], const Element v[], unsigned int n)
556{
557 Element result = ring.Identity();
558 for (unsigned int i=0; i<n; i++)
559 ring.Accumulate(result, ring.Multiply(y[i], v[i]));
560 return result;
561}
562
563// ********************************************************
564
565template <class T, int instance>
566const PolynomialOverFixedRing<T, instance> &PolynomialOverFixedRing<T, instance>::Zero()
567{
568 return Singleton<ThisType>().Ref();
569}
570
571template <class T, int instance>
572const PolynomialOverFixedRing<T, instance> &PolynomialOverFixedRing<T, instance>::One()
573{
574 return Singleton<ThisType, NewOnePolynomial>().Ref();
575}
576
577NAMESPACE_END
PolynomialOver::DivideByZero
division by zero exception
Definition polynomi.h:28
PolynomialOverFixedRing
Polynomials over a fixed ring.
Definition polynomi.h:164
PolynomialOver
represents single-variable polynomials over arbitrary rings
Definition polynomi.h:22
PolynomialOver::GetCoefficient
CoefficientType GetCoefficient(unsigned int i, const Ring &ring) const
return coefficient for x^i
Definition polynomi.cpp:86
PolynomialOver::Divide
static void Divide(PolynomialOver< Ring > &r, PolynomialOver< Ring > &q, const PolynomialOver< Ring > &a, const PolynomialOver< Ring > &d, const Ring &ring)
calculate r and q such that (a == d*q + r) && (0 <= degree of r < degree of d)
Definition polynomi.cpp:430
PolynomialOver::SetCoefficient
void SetCoefficient(unsigned int i, const CoefficientType &value, const Ring &ring)
set the coefficient for x^i to value
Definition polynomi.cpp:182
RandomNumberGenerator
Interface for random number generators.
Definition cryptlib.h:1435
RingOfPolynomialsOver
Ring of polynomials over another ring.
Definition polynomi.h:315
SecBlock
Secure memory block with allocator and cleanup.
Definition secblock.h:731
Singleton
Restricts the instantiation of a class to one static object without locks.
Definition misc.h:307
Singleton::Ref
CRYPTOPP_NOINLINE const T & Ref(CRYPTOPP_NOINLINE_DOTDOTDOT) const
Return a reference to the inner Singleton object.
Definition misc.h:327
STDMAX
const T & STDMAX(const T &a, const T &b)
Replacement function for std::max.
Definition misc.h:666
pch.h
Precompiled header file.
polynomi.h
Classes for polynomial basis and operations.
secblock.h
Classes and functions for secure memory allocations.
Generated by doxygen 1.9.8