00001
00002
00003
00004
00005
00006
#include "pch.h"
00007
#include "polynomi.h"
00008
#include "secblock.h"
00009
00010
#include <strstream>
00011
#include <iostream>
00012
00013 NAMESPACE_BEGIN(CryptoPP)
00014
00015 template <class T>
00016
void PolynomialOver<T>::Randomize(
RandomNumberGenerator &rng, const RandomizationParameter ¶meter, const Ring &ring)
00017 {
00018 m_coefficients.resize(parameter.m_coefficientCount);
00019
for (
unsigned int i=0; i<m_coefficients.size(); ++i)
00020 m_coefficients[i] = ring.RandomElement(rng, parameter.m_coefficientParameter);
00021 }
00022
00023
template <
class T>
00024
void PolynomialOver<T>::FromStr(
const char *str,
const Ring &ring)
00025 {
00026 std::istrstream in((
char *)str);
00027
bool positive =
true;
00028 CoefficientType coef;
00029
unsigned int power;
00030
00031
while (in)
00032 {
00033 std::ws(in);
00034
if (in.peek() ==
'x')
00035 coef = ring.MultiplicativeIdentity();
00036
else
00037 in >> coef;
00038
00039 std::ws(in);
00040
if (in.peek() ==
'x')
00041 {
00042 in.get();
00043 std::ws(in);
00044
if (in.peek() ==
'^')
00045 {
00046 in.get();
00047 in >> power;
00048 }
00049
else
00050 power = 1;
00051 }
00052
else
00053 power = 0;
00054
00055
if (!positive)
00056 coef = ring.Inverse(coef);
00057
00058 SetCoefficient(power, coef, ring);
00059
00060 std::ws(in);
00061
switch (in.get())
00062 {
00063
case '+':
00064 positive =
true;
00065
break;
00066
case '-':
00067 positive =
false;
00068
break;
00069
default:
00070
return;
00071 }
00072 }
00073 }
00074
00075
template <
class T>
00076
unsigned int PolynomialOver<T>::CoefficientCount(
const Ring &ring)
const
00077
{
00078
unsigned count = m_coefficients.size();
00079
while (count && ring.Equal(m_coefficients[count-1], ring.Identity()))
00080 count--;
00081 const_cast<std::vector<CoefficientType> &>(m_coefficients).resize(count);
00082
return count;
00083 }
00084
00085
template <
class T>
00086 typename PolynomialOver<T>::CoefficientType
PolynomialOver<T>::GetCoefficient(
unsigned int i,
const Ring &ring)
const
00087
{
00088
return (i < m_coefficients.size()) ? m_coefficients[i] : ring.Identity();
00089 }
00090
00091
template <
class T>
00092
PolynomialOver<T>&
PolynomialOver<T>::operator=(
const PolynomialOver<T>& t)
00093 {
00094
if (
this != &t)
00095 {
00096 m_coefficients.resize(t.m_coefficients.size());
00097
for (
unsigned int i=0; i<m_coefficients.size(); i++)
00098 m_coefficients[i] = t.m_coefficients[i];
00099 }
00100
return *
this;
00101 }
00102
00103
template <
class T>
00104
PolynomialOver<T>&
PolynomialOver<T>::Accumulate(
const PolynomialOver<T>& t,
const Ring &ring)
00105 {
00106
unsigned int count = t.CoefficientCount(ring);
00107
00108
if (count > CoefficientCount(ring))
00109 m_coefficients.resize(count, ring.Identity());
00110
00111
for (
unsigned int i=0; i<count; i++)
00112 ring.Accumulate(m_coefficients[i], t.GetCoefficient(i, ring));
00113
00114
return *
this;
00115 }
00116
00117
template <
class T>
00118
PolynomialOver<T>&
PolynomialOver<T>::Reduce(
const PolynomialOver<T>& t,
const Ring &ring)
00119 {
00120
unsigned int count = t.CoefficientCount(ring);
00121
00122
if (count > CoefficientCount(ring))
00123 m_coefficients.resize(count, ring.Identity());
00124
00125
for (
unsigned int i=0; i<count; i++)
00126 ring.Reduce(m_coefficients[i], t.GetCoefficient(i, ring));
00127
00128
return *
this;
00129 }
00130
00131
template <
class T>
00132
typename PolynomialOver<T>::CoefficientType
PolynomialOver<T>::EvaluateAt(
const CoefficientType &x,
const Ring &ring)
const
00133
{
00134
int degree =
Degree(ring);
00135
00136
if (degree < 0)
00137
return ring.Identity();
00138
00139 CoefficientType result = m_coefficients[degree];
00140
for (
int j=degree-1; j>=0; j--)
00141 {
00142 result = ring.Multiply(result, x);
00143 ring.Accumulate(result, m_coefficients[j]);
00144 }
00145
return result;
00146 }
00147
00148
template <
class T>
00149
PolynomialOver<T>&
PolynomialOver<T>::ShiftLeft(
unsigned int n,
const Ring &ring)
00150 {
00151
unsigned int i = CoefficientCount(ring) + n;
00152 m_coefficients.resize(i, ring.Identity());
00153
while (i > n)
00154 {
00155 i--;
00156 m_coefficients[i] = m_coefficients[i-n];
00157 }
00158
while (i)
00159 {
00160 i--;
00161 m_coefficients[i] = ring.Identity();
00162 }
00163
return *
this;
00164 }
00165
00166
template <
class T>
00167
PolynomialOver<T>&
PolynomialOver<T>::ShiftRight(
unsigned int n,
const Ring &ring)
00168 {
00169
unsigned int count = CoefficientCount(ring);
00170
if (count > n)
00171 {
00172
for (
unsigned int i=0; i<count-n; i++)
00173 m_coefficients[i] = m_coefficients[i+n];
00174 m_coefficients.resize(count-n, ring.Identity());
00175 }
00176
else
00177 m_coefficients.resize(0, ring.Identity());
00178
return *
this;
00179 }
00180
00181
template <
class T>
00182 void PolynomialOver<T>::SetCoefficient(
unsigned int i,
const CoefficientType &value,
const Ring &ring)
00183 {
00184
if (i >= m_coefficients.size())
00185 m_coefficients.resize(i+1, ring.Identity());
00186 m_coefficients[i] = value;
00187 }
00188
00189
template <
class T>
00190
void PolynomialOver<T>::Negate(
const Ring &ring)
00191 {
00192
unsigned int count = CoefficientCount(ring);
00193
for (
unsigned int i=0; i<count; i++)
00194 m_coefficients[i] = ring.Inverse(m_coefficients[i]);
00195 }
00196
00197
template <
class T>
00198
void PolynomialOver<T>::swap(
PolynomialOver<T> &t)
00199 {
00200 m_coefficients.swap(t.m_coefficients);
00201 }
00202
00203
template <
class T>
00204
bool PolynomialOver<T>::Equals(
const PolynomialOver<T>& t,
const Ring &ring)
const
00205
{
00206
unsigned int count = CoefficientCount(ring);
00207
00208
if (count != t.CoefficientCount(ring))
00209
return false;
00210
00211
for (
unsigned int i=0; i<count; i++)
00212
if (!ring.Equal(m_coefficients[i], t.m_coefficients[i]))
00213
return false;
00214
00215
return true;
00216 }
00217
00218
template <
class T>
00219
PolynomialOver<T> PolynomialOver<T>::Plus(
const PolynomialOver<T>& t,
const Ring &ring)
const
00220
{
00221
unsigned int i;
00222
unsigned int count = CoefficientCount(ring);
00223
unsigned int tCount = t.CoefficientCount(ring);
00224
00225
if (count > tCount)
00226 {
00227
PolynomialOver<T> result(ring, count);
00228
00229
for (i=0; i<tCount; i++)
00230 result.
m_coefficients[i] = ring.Add(m_coefficients[i], t.m_coefficients[i]);
00231
for (; i<count; i++)
00232 result.
m_coefficients[i] = m_coefficients[i];
00233
00234
return result;
00235 }
00236
else
00237 {
00238
PolynomialOver<T> result(ring, tCount);
00239
00240
for (i=0; i<count; i++)
00241 result.
m_coefficients[i] = ring.Add(m_coefficients[i], t.m_coefficients[i]);
00242
for (; i<tCount; i++)
00243 result.
m_coefficients[i] = t.m_coefficients[i];
00244
00245
return result;
00246 }
00247 }
00248
00249
template <
class T>
00250
PolynomialOver<T> PolynomialOver<T>::Minus(
const PolynomialOver<T>& t,
const Ring &ring)
const
00251
{
00252
unsigned int i;
00253
unsigned int count = CoefficientCount(ring);
00254
unsigned int tCount = t.CoefficientCount(ring);
00255
00256
if (count > tCount)
00257 {
00258
PolynomialOver<T> result(ring, count);
00259
00260
for (i=0; i<tCount; i++)
00261 result.
m_coefficients[i] = ring.Subtract(m_coefficients[i], t.m_coefficients[i]);
00262
for (; i<count; i++)
00263 result.
m_coefficients[i] = m_coefficients[i];
00264
00265
return result;
00266 }
00267
else
00268 {
00269
PolynomialOver<T> result(ring, tCount);
00270
00271
for (i=0; i<count; i++)
00272 result.
m_coefficients[i] = ring.Subtract(m_coefficients[i], t.m_coefficients[i]);
00273
for (; i<tCount; i++)
00274 result.
m_coefficients[i] = ring.Inverse(t.m_coefficients[i]);
00275
00276
return result;
00277 }
00278 }
00279
00280
template <
class T>
00281
PolynomialOver<T> PolynomialOver<T>::Inverse(
const Ring &ring)
const
00282
{
00283
unsigned int count = CoefficientCount(ring);
00284
PolynomialOver<T> result(ring, count);
00285
00286
for (
unsigned int i=0; i<count; i++)
00287 result.
m_coefficients[i] = ring.Inverse(m_coefficients[i]);
00288
00289
return result;
00290 }
00291
00292
template <
class T>
00293
PolynomialOver<T> PolynomialOver<T>::Times(
const PolynomialOver<T>& t,
const Ring &ring)
const
00294
{
00295
if (IsZero(ring) || t.IsZero(ring))
00296
return PolynomialOver<T>();
00297
00298
unsigned int count1 = CoefficientCount(ring), count2 = t.CoefficientCount(ring);
00299
PolynomialOver<T> result(ring, count1 + count2 - 1);
00300
00301
for (
unsigned int i=0; i<count1; i++)
00302
for (
unsigned int j=0; j<count2; j++)
00303 ring.Accumulate(result.
m_coefficients[i+j], ring.Multiply(m_coefficients[i], t.m_coefficients[j]));
00304
00305
return result;
00306 }
00307
00308
template <
class T>
00309
PolynomialOver<T> PolynomialOver<T>::DividedBy(
const PolynomialOver<T>& t,
const Ring &ring)
const
00310
{
00311
PolynomialOver<T> remainder, quotient;
00312
Divide(remainder, quotient, *
this, t, ring);
00313
return quotient;
00314 }
00315
00316
template <
class T>
00317
PolynomialOver<T> PolynomialOver<T>::Modulo(
const PolynomialOver<T>& t,
const Ring &ring)
const
00318
{
00319
PolynomialOver<T> remainder, quotient;
00320
Divide(remainder, quotient, *
this, t, ring);
00321
return remainder;
00322 }
00323
00324
template <
class T>
00325
PolynomialOver<T> PolynomialOver<T>::MultiplicativeInverse(
const Ring &ring)
const
00326
{
00327
return Degree(ring)==0 ? ring.MultiplicativeInverse(m_coefficients[0]) : ring.Identity();
00328 }
00329
00330
template <
class T>
00331
bool PolynomialOver<T>::IsUnit(
const Ring &ring)
const
00332
{
00333
return Degree(ring)==0 && ring.IsUnit(m_coefficients[0]);
00334 }
00335
00336
template <
class T>
00337 std::istream&
PolynomialOver<T>::Input(std::istream &in,
const Ring &ring)
00338 {
00339
char c;
00340
unsigned int length = 0;
00341
SecBlock<char> str(length + 16);
00342
bool paren =
false;
00343
00344 std::ws(in);
00345
00346
if (in.peek() ==
'(')
00347 {
00348 paren =
true;
00349 in.get();
00350 }
00351
00352
do
00353 {
00354 in.read(&c, 1);
00355 str[length++] = c;
00356
if (length >= str.
size())
00357 str.
Grow(length + 16);
00358 }
00359
00360
00361
while (in && ((paren && c !=
')') || (!paren && c !=
'\n')));
00362
00363 str[length-1] =
'\0';
00364 *
this =
PolynomialOver<T>(str, ring);
00365
00366
return in;
00367 }
00368
00369
template <
class T>
00370 std::ostream&
PolynomialOver<T>::Output(std::ostream &out,
const Ring &ring)
const
00371
{
00372
unsigned int i = CoefficientCount(ring);
00373
if (i)
00374 {
00375
bool firstTerm =
true;
00376
00377
while (i--)
00378 {
00379
if (m_coefficients[i] != ring.Identity())
00380 {
00381
if (firstTerm)
00382 {
00383 firstTerm =
false;
00384
if (!i || !ring.Equal(m_coefficients[i], ring.MultiplicativeIdentity()))
00385 out << m_coefficients[i];
00386 }
00387
else
00388 {
00389 CoefficientType inverse = ring.Inverse(m_coefficients[i]);
00390 std::ostrstream pstr, nstr;
00391
00392 pstr << m_coefficients[i];
00393 nstr << inverse;
00394
00395
if (pstr.pcount() <= nstr.pcount())
00396 {
00397 out <<
" + ";
00398
if (!i || !ring.Equal(m_coefficients[i], ring.MultiplicativeIdentity()))
00399 out << m_coefficients[i];
00400 }
00401
else
00402 {
00403 out <<
" - ";
00404
if (!i || !ring.Equal(inverse, ring.MultiplicativeIdentity()))
00405 out << inverse;
00406 }
00407 }
00408
00409
switch (i)
00410 {
00411
case 0:
00412
break;
00413
case 1:
00414 out <<
"x";
00415
break;
00416
default:
00417 out <<
"x^" << i;
00418 }
00419 }
00420 }
00421 }
00422
else
00423 {
00424 out << ring.Identity();
00425 }
00426
return out;
00427 }
00428
00429
template <
class T>
00430
void PolynomialOver<T>::Divide(
PolynomialOver<T> &r,
PolynomialOver<T> &q,
const PolynomialOver<T> &a,
const PolynomialOver<T> &d,
const Ring &ring)
00431 {
00432
unsigned int i = a.CoefficientCount(ring);
00433
const int dDegree = d.Degree(ring);
00434
00435
if (dDegree < 0)
00436
throw DivideByZero();
00437
00438 r = a;
00439 q.
m_coefficients.resize(STDMAX(0,
int(i - dDegree)));
00440
00441
while (i > (
unsigned int)dDegree)
00442 {
00443 --i;
00444 q.
m_coefficients[i-dDegree] = ring.Divide(r.
m_coefficients[i], d.m_coefficients[dDegree]);
00445
for (
int j=0; j<=dDegree; j++)
00446 ring.Reduce(r.
m_coefficients[i-dDegree+j], ring.Multiply(q.
m_coefficients[i-dDegree], d.m_coefficients[j]));
00447 }
00448
00449 r.
CoefficientCount(ring);
00450 }
00451
00452
00453
00454
00455
template <
class T>
00456
void RingOfPolynomialsOver<T>::CalculateAlpha(std::vector<CoefficientType> &alpha,
const CoefficientType x[],
const CoefficientType y[],
unsigned int n)
const
00457
{
00458
for (
unsigned int j=0; j<n; ++j)
00459 alpha[j] = y[j];
00460
00461
for (
unsigned int k=1; k<n; ++k)
00462 {
00463
for (
unsigned int j=n-1; j>=k; --j)
00464 {
00465 m_ring.Reduce(alpha[j], alpha[j-1]);
00466
00467 CoefficientType d = m_ring.Subtract(x[j], x[j-k]);
00468
if (!m_ring.IsUnit(d))
00469
throw InterpolationFailed();
00470 alpha[j] = m_ring.Divide(alpha[j], d);
00471 }
00472 }
00473 }
00474
00475
template <
class T>
00476
typename RingOfPolynomialsOver<T>::Element
RingOfPolynomialsOver<T>::Interpolate(
const CoefficientType x[],
const CoefficientType y[],
unsigned int n)
const
00477
{
00478 assert(n > 0);
00479
00480 std::vector<CoefficientType> alpha(n);
00481 CalculateAlpha(alpha, x, y, n);
00482
00483 std::vector<CoefficientType> coefficients((size_t)n, m_ring.Identity());
00484 coefficients[0] = alpha[n-1];
00485
00486
for (
int j=n-2; j>=0; --j)
00487 {
00488
for (
unsigned int i=n-j-1; i>0; i--)
00489 coefficients[i] = m_ring.Subtract(coefficients[i-1], m_ring.Multiply(coefficients[i], x[j]));
00490
00491 coefficients[0] = m_ring.Subtract(alpha[j], m_ring.Multiply(coefficients[0], x[j]));
00492 }
00493
00494
return PolynomialOver<T>(coefficients.begin(), coefficients.end());
00495 }
00496
00497
template <
class T>
00498
typename RingOfPolynomialsOver<T>::CoefficientType
RingOfPolynomialsOver<T>::InterpolateAt(
const CoefficientType &position,
const CoefficientType x[],
const CoefficientType y[],
unsigned int n)
const
00499
{
00500 assert(n > 0);
00501
00502 std::vector<CoefficientType> alpha(n);
00503 CalculateAlpha(alpha, x, y, n);
00504
00505 CoefficientType result = alpha[n-1];
00506
for (
int j=n-2; j>=0; --j)
00507 {
00508 result = m_ring.Multiply(result, m_ring.Subtract(position, x[j]));
00509 m_ring.Accumulate(result, alpha[j]);
00510 }
00511
return result;
00512 }
00513
00514
template <
class Ring,
class Element>
00515
void PrepareBulkPolynomialInterpolation(
const Ring &ring, Element *w,
const Element x[],
unsigned int n)
00516 {
00517
for (
unsigned int i=0; i<n; i++)
00518 {
00519 Element t = ring.MultiplicativeIdentity();
00520
for (
unsigned int j=0; j<n; j++)
00521
if (i != j)
00522 t = ring.Multiply(t, ring.Subtract(x[i], x[j]));
00523 w[i] = ring.MultiplicativeInverse(t);
00524 }
00525 }
00526
00527
template <
class Ring,
class Element>
00528
void PrepareBulkPolynomialInterpolationAt(
const Ring &ring, Element *v,
const Element &position,
const Element x[],
const Element w[],
unsigned int n)
00529 {
00530 assert(n > 0);
00531
00532 std::vector<Element> a(2*n-1);
00533
unsigned int i;
00534
00535
for (i=0; i<n; i++)
00536 a[n-1+i] = ring.Subtract(position, x[i]);
00537
00538
for (i=n-1; i>1; i--)
00539 a[i-1] = ring.Multiply(a[2*i], a[2*i-1]);
00540
00541 a[0] = ring.MultiplicativeIdentity();
00542
00543
for (i=0; i<n-1; i++)
00544 {
00545 std::swap(a[2*i+1], a[2*i+2]);
00546 a[2*i+1] = ring.Multiply(a[i], a[2*i+1]);
00547 a[2*i+2] = ring.Multiply(a[i], a[2*i+2]);
00548 }
00549
00550
for (i=0; i<n; i++)
00551 v[i] = ring.Multiply(a[n-1+i], w[i]);
00552 }
00553
00554
template <
class Ring,
class Element>
00555 Element BulkPolynomialInterpolateAt(
const Ring &ring,
const Element y[],
const Element v[],
unsigned int n)
00556 {
00557 Element result = ring.Identity();
00558
for (
unsigned int i=0; i<n; i++)
00559 ring.Accumulate(result, ring.Multiply(y[i], v[i]));
00560
return result;
00561 }
00562
00563
00564
00565
template <
class T,
int instance>
00566
const PolynomialOverFixedRing<T, instance> &
PolynomialOverFixedRing<T, instance>::Zero()
00567 {
00568
static const PolynomialOverFixedRing<T, instance> zero;
00569
return zero;
00570 }
00571
00572
template <
class T,
int instance>
00573
const PolynomialOverFixedRing<T, instance> &
PolynomialOverFixedRing<T, instance>::One()
00574 {
00575
static const PolynomialOverFixedRing<T, instance> one = fixedRing.MultiplicativeIdentity();
00576
return one;
00577 }
00578
00579 NAMESPACE_END