Fabcoin Core  0.16.2
P2P Digital Currency
polynomi.cpp
Go to the documentation of this file.
1 // polynomi.cpp - 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 
14 
15 template <class T>
16 void 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 
23 template <class T>
24 void 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 
75 template <class T>
76 unsigned 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 
85 template <class T>
87 {
88  return (i < m_coefficients.size()) ? m_coefficients[i] : ring.Identity();
89 }
90 
91 template <class 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 
103 template <class T>
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 
117 template <class T>
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 
131 template <class T>
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 
148 template <class T>
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 
166 template <class T>
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 
181 template <class T>
182 void 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 
189 template <class T>
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 
197 template <class T>
199 {
200  m_coefficients.swap(t.m_coefficients);
201 }
202 
203 template <class T>
204 bool 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 
218 template <class T>
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 
249 template <class T>
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 
280 template <class T>
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 
292 template <class T>
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 
308 template <class T>
310 {
311  PolynomialOver<T> remainder, quotient;
312  Divide(remainder, quotient, *this, t, ring);
313  return quotient;
314 }
315 
316 template <class T>
318 {
319  PolynomialOver<T> remainder, quotient;
320  Divide(remainder, quotient, *this, t, ring);
321  return remainder;
322 }
323 
324 template <class T>
326 {
327  return Degree(ring)==0 ? ring.MultiplicativeInverse(m_coefficients[0]) : ring.Identity();
328 }
329 
330 template <class T>
331 bool PolynomialOver<T>::IsUnit(const Ring &ring) const
332 {
333  return Degree(ring)==0 && ring.IsUnit(m_coefficients[0]);
334 }
335 
336 template <class T>
337 std::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 
369 template <class T>
370 std::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 
429 template <class T>
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()
455 template <class T>
456 void 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 
475 template <class T>
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 
497 template <class T>
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 
514 template <class Ring, class Element>
515 void 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 
527 template <class Ring, class Element>
528 void 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 
554 template <class Ring, class Element>
555 Element 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 
565 template <class T, int instance>
567 {
568  return Singleton<ThisType>().Ref();
569 }
570 
571 template <class T, int instance>
573 {
575 }
576 
PolynomialOver< Ring > & operator=(const PolynomialOver< Ring > &t)
Definition: polynomi.cpp:92
specify the distribution for randomization functions
Definition: polynomi.h:38
Element::CoefficientType CoefficientType
Definition: polynomi.h:323
void SetCoefficient(unsigned int i, const CoefficientType &value, const Ring &ring)
set the coefficient for x^i to value
Definition: polynomi.cpp:182
Restricts the instantiation of a class to one static object without locks.
Definition: misc.h:274
void swap(dev::eth::Watch &_a, dev::eth::Watch &_b)
Definition: Interface.h:284
PolynomialOver< Ring > Minus(const PolynomialOver< Ring > &t, const Ring &ring) const
Definition: polynomi.cpp:250
#define T(i, x)
#define NAMESPACE_BEGIN(x)
Definition: config.h:200
PolynomialOver< Ring > & ShiftRight(unsigned int n, const Ring &ring)
Definition: polynomi.cpp:167
size_t count
Definition: ExecStats.cpp:37
Secure memory block with allocator and cleanup.
Definition: secblock.h:437
void swap(PolynomialOver< Ring > &t)
Definition: polynomi.cpp:198
CoefficientType GetCoefficient(unsigned int i, const Ring &ring) const
return coefficient for x^i
Definition: polynomi.cpp:86
size_type size() const
Provides the count of elements in the SecBlock.
Definition: secblock.h:524
#define c(i)
Interface for random number generators.
Definition: cryptlib.h:1188
bool Equals(const PolynomialOver< Ring > &t, const Ring &ring) const
Definition: polynomi.cpp:204
PolynomialOver< Ring > Inverse(const Ring &ring) const
Definition: polynomi.cpp:281
unsigned int CoefficientCount(const Ring &ring) const
Definition: polynomi.cpp:76
represents single-variable polynomials over arbitrary rings
Definition: polynomi.h:25
T::Element CoefficientType
Definition: polynomi.h:51
bool IsUnit(const Ring &ring) const
Definition: polynomi.cpp:331
Classes and functions for secure memory allocations.
CoefficientType EvaluateAt(const CoefficientType &x, const Ring &ring) const
Definition: polynomi.cpp:132
bool IsZero(const Ring &ring) const
Definition: polynomi.h:124
#define a(i)
PolynomialOver< Ring > Plus(const PolynomialOver< Ring > &t, const Ring &ring) const
Definition: polynomi.cpp:219
Element Interpolate(const CoefficientType x[], const CoefficientType y[], unsigned int n) const
Definition: polynomi.cpp:476
#define x(i)
void FromStr(const char *str, const Ring &ring)
Definition: polynomi.cpp:24
CoefficientType InterpolateAt(const CoefficientType &position, const CoefficientType x[], const CoefficientType y[], unsigned int n) const
Definition: polynomi.cpp:498
void CalculateAlpha(std::vector< CoefficientType > &alpha, const CoefficientType x[], const CoefficientType y[], unsigned int n) const
Definition: polynomi.cpp:456
PolynomialOver< Ring > & Accumulate(const PolynomialOver< Ring > &t, const Ring &ring)
Definition: polynomi.cpp:104
PolynomialOver< Ring > Modulo(const PolynomialOver< Ring > &t, const Ring &ring) const
Definition: polynomi.cpp:317
void Divide(word *R, word *Q, word *T, const word *A, size_t NA, const word *B, size_t NB)
Definition: integer.cpp:2692
Polynomials over a fixed ring.
Definition: polynomi.h:167
static const ThisType & Zero()
Definition: polynomi.cpp:566
PolynomialOver< Ring > & ShiftLeft(unsigned int n, const Ring &ring)
Definition: polynomi.cpp:149
PolynomialOver< Ring > DividedBy(const PolynomialOver< Ring > &t, const Ring &ring) const
Definition: polynomi.cpp:309
#define CRYPTOPP_ASSERT(exp)
Definition: trap.h:92
Classes for polynomial basis and operations.
void PrepareBulkPolynomialInterpolationAt(const Ring &ring, Element *v, const Element &position, const Element x[], const Element w[], unsigned int n)
Definition: polynomi.cpp:528
static const ThisType & One()
Definition: polynomi.cpp:572
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< Ring > & Reduce(const PolynomialOver< Ring > &t, const Ring &ring)
Definition: polynomi.cpp:118
std::ostream & Output(std::ostream &out, const Ring &ring) const
Definition: polynomi.cpp:370
const T & STDMAX(const T &a, const T &b)
Replacement function for std::max.
Definition: misc.h:487
#define NAMESPACE_END
Definition: config.h:201
void Grow(size_type newSize)
Change size and preserve contents.
Definition: secblock.h:673
std::vector< CoefficientType > m_coefficients
Definition: polynomi.h:162
#define d(i)
Definition: sha.cpp:732
division by zero exception
Definition: polynomi.h:31
void PrepareBulkPolynomialInterpolation(const Ring &ring, Element *w, const Element x[], unsigned int n)
Definition: polynomi.cpp:515
PolynomialOver< Ring > Times(const PolynomialOver< Ring > &t, const Ring &ring) const
Definition: polynomi.cpp:293
void Negate(const Ring &ring)
Definition: polynomi.cpp:190
Element BulkPolynomialInterpolateAt(const Ring &ring, const Element y[], const Element v[], unsigned int n)
Definition: polynomi.cpp:555
int Degree(const Ring &ring) const
the zero polynomial will return a degree of -1
Definition: polynomi.h:95
PolynomialOver< Ring > MultiplicativeInverse(const Ring &ring) const
Definition: polynomi.cpp:325
std::istream & Input(std::istream &in, const Ring &ring)
Definition: polynomi.cpp:337