Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

# polynomi.cpp

    00001 // polynomi.cpp - written and placed in the public domain by Wei Dai
    00002
    00003 // Part of the code for polynomial evaluation and interpolation
    00004 // originally came from Hal Finney's public domain secsplit.c.
    00005
    00006 #include "pch.h"
    00007 #include "polynomi.h"
    00008 #include "secblock.h"
    00009
    00010 #include <sstream>
    00011 #include <iostream>
    00012
    00013 NAMESPACE_BEGIN(CryptoPP)
    00014
    00015 template <class T>
    00016 void PolynomialOver<T>::Randomize(RandomNumberGenerator &rng, const RandomizationParameter &parameter, 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::istringstream 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;         // something's wrong with the input string
    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         // if we started with a left paren, then read until we find a right paren,
    00360         // otherwise read until the end of the line
    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);       // resize r.m_coefficients
    00450 }
    00451
    00452 // ********************************************************
    00453
    00454 // helper function for Interpolate() and InterpolateAt()
    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         return Singleton<ThisType>().Ref();
    00569 }
    00570
    00571 template <class T, int instance>
    00572 const PolynomialOverFixedRing<T, instance> &PolynomialOverFixedRing<T, instance>::One()
    00573 {
    00574         return Singleton<ThisType, NewOnePolynomial>().Ref();
    00575 }
    00576
    00577 NAMESPACE_END