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 <strstream> // can't use <sstream> because GCC 2.95.2 doesn't have it 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::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; // 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

Generated on Fri Aug 27 15:51:09 2004 for Crypto++ by doxygen 1.3.8