- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# algebra.cpp

    00001 // algebra.cpp - written and placed in the public domain by Wei Dai
    00002
    00003 #include "pch.h"
    00004
    00005 #ifndef CRYPTOPP_ALGEBRA_CPP    // SunCC workaround: compiler could cause this file to be included twice
    00006 #define CRYPTOPP_ALGEBRA_CPP
    00007
    00008 #include "algebra.h"
    00009 #include "integer.h"
    00010
    00011 #include <vector>
    00012
    00013 NAMESPACE_BEGIN(CryptoPP)
    00014
    00015 template <class T> const T& AbstractGroup<T>::Double(const Element &a) const
    00016 {
    00017         return Add(a, a);
    00018 }
    00019
    00020 template <class T> const T& AbstractGroup<T>::Subtract(const Element &a, const Element &b) const
    00021 {
    00022         // make copy of a in case Inverse() overwrites it
    00023         Element a1(a);
    00024         return Add(a1, Inverse(b));
    00025 }
    00026
    00027 template <class T> T& AbstractGroup<T>::Accumulate(Element &a, const Element &b) const
    00028 {
    00029         return a = Add(a, b);
    00030 }
    00031
    00032 template <class T> T& AbstractGroup<T>::Reduce(Element &a, const Element &b) const
    00033 {
    00034         return a = Subtract(a, b);
    00035 }
    00036
    00037 template <class T> const T& AbstractRing<T>::Square(const Element &a) const
    00038 {
    00039         return Multiply(a, a);
    00040 }
    00041
    00042 template <class T> const T& AbstractRing<T>::Divide(const Element &a, const Element &b) const
    00043 {
    00044         // make copy of a in case MultiplicativeInverse() overwrites it
    00045         Element a1(a);
    00046         return Multiply(a1, MultiplicativeInverse(b));
    00047 }
    00048
    00049 template <class T> const T& AbstractEuclideanDomain<T>::Mod(const Element &a, const Element &b) const
    00050 {
    00051         Element q;
    00052         DivisionAlgorithm(result, q, a, b);
    00053         return result;
    00054 }
    00055
    00056 template <class T> const T& AbstractEuclideanDomain<T>::Gcd(const Element &a, const Element &b) const
    00057 {
    00058         Element g[3]={b, a};
    00059         unsigned int i0=0, i1=1, i2=2;
    00060
    00061         while (!Equal(g[i1], this->Identity()))
    00062         {
    00063                 g[i2] = Mod(g[i0], g[i1]);
    00064                 unsigned int t = i0; i0 = i1; i1 = i2; i2 = t;
    00065         }
    00066
    00067         return result = g[i0];
    00068 }
    00069
    00070 template <class T> const typename QuotientRing<T>::Element& QuotientRing<T>::MultiplicativeInverse(const Element &a) const
    00071 {
    00072         Element g[3]={m_modulus, a};
    00073         Element v[3]={m_domain.Identity(), m_domain.MultiplicativeIdentity()};
    00074         Element y;
    00075         unsigned int i0=0, i1=1, i2=2;
    00076
    00077         while (!Equal(g[i1], Identity()))
    00078         {
    00079                 // y = g[i0] / g[i1];
    00080                 // g[i2] = g[i0] % g[i1];
    00081                 m_domain.DivisionAlgorithm(g[i2], y, g[i0], g[i1]);
    00082                 // v[i2] = v[i0] - (v[i1] * y);
    00083                 v[i2] = m_domain.Subtract(v[i0], m_domain.Multiply(v[i1], y));
    00084                 unsigned int t = i0; i0 = i1; i1 = i2; i2 = t;
    00085         }
    00086
    00087         return m_domain.IsUnit(g[i0]) ? m_domain.Divide(v[i0], g[i0]) : m_domain.Identity();
    00088 }
    00089
    00090 template <class T> T AbstractGroup<T>::ScalarMultiply(const Element &base, const Integer &exponent) const
    00091 {
    00092         Element result;
    00093         SimultaneousMultiply(&result, base, &exponent, 1);
    00094         return result;
    00095 }
    00096
    00097 template <class T> T AbstractGroup<T>::CascadeScalarMultiply(const Element &x, const Integer &e1, const Element &y, const Integer &e2) const
    00098 {
    00099         const unsigned expLen = STDMAX(e1.BitCount(), e2.BitCount());
    00100         if (expLen==0)
    00101                 return Identity();
    00102
    00103         const unsigned w = (expLen <= 46 ? 1 : (expLen <= 260 ? 2 : 3));
    00104         const unsigned tableSize = 1<<w;
    00105         std::vector<Element> powerTable(tableSize << w);
    00106
    00107         powerTable[1] = x;
    00108         powerTable[tableSize] = y;
    00109         if (w==1)
    00110                 powerTable[3] = Add(x,y);
    00111         else
    00112         {
    00113                 powerTable[2] = Double(x);
    00114                 powerTable[2*tableSize] = Double(y);
    00115
    00116                 unsigned i, j;
    00117
    00118                 for (i=3; i<tableSize; i+=2)
    00119                         powerTable[i] = Add(powerTable[i-2], powerTable[2]);
    00120                 for (i=1; i<tableSize; i+=2)
    00121                         for (j=i+tableSize; j<(tableSize<<w); j+=tableSize)
    00122                                 powerTable[j] = Add(powerTable[j-tableSize], y);
    00123
    00124                 for (i=3*tableSize; i<(tableSize<<w); i+=2*tableSize)
    00125                         powerTable[i] = Add(powerTable[i-2*tableSize], powerTable[2*tableSize]);
    00126                 for (i=tableSize; i<(tableSize<<w); i+=2*tableSize)
    00127                         for (j=i+2; j<i+tableSize; j+=2)
    00128                                 powerTable[j] = Add(powerTable[j-1], x);
    00129         }
    00130
    00131         Element result;
    00132         unsigned power1 = 0, power2 = 0, prevPosition = expLen-1;
    00133         bool firstTime = true;
    00134
    00135         for (int i = expLen-1; i>=0; i--)
    00136         {
    00137                 power1 = 2*power1 + e1.GetBit(i);
    00138                 power2 = 2*power2 + e2.GetBit(i);
    00139
    00140                 if (i==0 || 2*power1 >= tableSize || 2*power2 >= tableSize)
    00141                 {
    00142                         unsigned squaresBefore = prevPosition-i;
    00143                         unsigned squaresAfter = 0;
    00144                         prevPosition = i;
    00145                         while ((power1 || power2) && power1%2 == 0 && power2%2==0)
    00146                         {
    00147                                 power1 /= 2;
    00148                                 power2 /= 2;
    00149                                 squaresBefore--;
    00150                                 squaresAfter++;
    00151                         }
    00152                         if (firstTime)
    00153                         {
    00154                                 result = powerTable[(power2<<w) + power1];
    00155                                 firstTime = false;
    00156                         }
    00157                         else
    00158                         {
    00159                                 while (squaresBefore--)
    00160                                         result = Double(result);
    00161                                 if (power1 || power2)
    00162                                         Accumulate(result, powerTable[(power2<<w) + power1]);
    00163                         }
    00164                         while (squaresAfter--)
    00165                                 result = Double(result);
    00166                         power1 = power2 = 0;
    00167                 }
    00168         }
    00169         return result;
    00170 }
    00171
    00172 template <class Element, class Iterator> Element GeneralCascadeMultiplication(const AbstractGroup<Element> &group, Iterator begin, Iterator end)
    00173 {
    00174         if (end-begin == 1)
    00175                 return group.ScalarMultiply(begin->base, begin->exponent);
    00176         else if (end-begin == 2)
    00177                 return group.CascadeScalarMultiply(begin->base, begin->exponent, (begin+1)->base, (begin+1)->exponent);
    00178         else
    00179         {
    00180                 Integer q, t;
    00181                 Iterator last = end;
    00182                 --last;
    00183
    00184                 std::make_heap(begin, end);
    00185                 std::pop_heap(begin, end);
    00186
    00187                 while (!!begin->exponent)
    00188                 {
    00189                         // last->exponent is largest exponent, begin->exponent is next largest
    00190                         t = last->exponent;
    00191                         Integer::Divide(last->exponent, q, t, begin->exponent);
    00192
    00193                         if (q == Integer::One())
    00194                                 group.Accumulate(begin->base, last->base);      // avoid overhead of ScalarMultiply()
    00195                         else
    00196                                 group.Accumulate(begin->base, group.ScalarMultiply(last->base, q));
    00197
    00198                         std::push_heap(begin, end);
    00199                         std::pop_heap(begin, end);
    00200                 }
    00201
    00202                 return group.ScalarMultiply(last->base, last->exponent);
    00203         }
    00204 }
    00205
    00206 struct WindowSlider
    00207 {
    00208         WindowSlider(const Integer &expIn, bool fastNegate, unsigned int windowSizeIn=0)
    00209                 : exp(expIn), windowModulus(Integer::One()), windowSize(windowSizeIn), windowBegin(0), fastNegate(fastNegate), firstTime(true), finished(false)
    00210         {
    00211                 if (windowSize == 0)
    00212                 {
    00213                         unsigned int expLen = exp.BitCount();
    00214                         windowSize = expLen <= 17 ? 1 : (expLen <= 24 ? 2 : (expLen <= 70 ? 3 : (expLen <= 197 ? 4 : (expLen <= 539 ? 5 : (expLen <= 1434 ? 6 : 7)))));
    00215                 }
    00216                 windowModulus <<= windowSize;
    00217         }
    00218
    00219         void FindNextWindow()
    00220         {
    00221                 unsigned int expLen = exp.WordCount() * WORD_BITS;
    00222                 unsigned int skipCount = firstTime ? 0 : windowSize;
    00223                 firstTime = false;
    00224                 while (!exp.GetBit(skipCount))
    00225                 {
    00226                         if (skipCount >= expLen)
    00227                         {
    00228                                 finished = true;
    00229                                 return;
    00230                         }
    00231                         skipCount++;
    00232                 }
    00233
    00234                 exp >>= skipCount;
    00235                 windowBegin += skipCount;
    00236                 expWindow = word32(exp % (word(1) << windowSize));
    00237
    00238                 if (fastNegate && exp.GetBit(windowSize))
    00239                 {
    00240                         negateNext = true;
    00241                         expWindow = (word32(1) << windowSize) - expWindow;
    00242                         exp += windowModulus;
    00243                 }
    00244                 else
    00245                         negateNext = false;
    00246         }
    00247
    00248         Integer exp, windowModulus;
    00249         unsigned int windowSize, windowBegin;
    00250         word32 expWindow;
    00251         bool fastNegate, negateNext, firstTime, finished;
    00252 };
    00253
    00254 template <class T>
    00255 void AbstractGroup<T>::SimultaneousMultiply(T *results, const T &base, const Integer *expBegin, unsigned int expCount) const
    00256 {
    00257         std::vector<std::vector<Element> > buckets(expCount);
    00258         std::vector<WindowSlider> exponents;
    00259         exponents.reserve(expCount);
    00260         unsigned int i;
    00261
    00262         for (i=0; i<expCount; i++)
    00263         {
    00264                 assert(expBegin->NotNegative());
    00265                 exponents.push_back(WindowSlider(*expBegin++, InversionIsFast(), 0));
    00266                 exponents[i].FindNextWindow();
    00267                 buckets[i].resize(1<<(exponents[i].windowSize-1), Identity());
    00268         }
    00269
    00270         unsigned int expBitPosition = 0;
    00271         Element g = base;
    00272         bool notDone = true;
    00273
    00274         while (notDone)
    00275         {
    00276                 notDone = false;
    00277                 for (i=0; i<expCount; i++)
    00278                 {
    00279                         if (!exponents[i].finished && expBitPosition == exponents[i].windowBegin)
    00280                         {
    00281                                 Element &bucket = buckets[i][exponents[i].expWindow/2];
    00282                                 if (exponents[i].negateNext)
    00283                                         Accumulate(bucket, Inverse(g));
    00284                                 else
    00285                                         Accumulate(bucket, g);
    00286                                 exponents[i].FindNextWindow();
    00287                         }
    00288                         notDone = notDone || !exponents[i].finished;
    00289                 }
    00290
    00291                 if (notDone)
    00292                 {
    00293                         g = Double(g);
    00294                         expBitPosition++;
    00295                 }
    00296         }
    00297
    00298         for (i=0; i<expCount; i++)
    00299         {
    00300                 Element &r = *results++;
    00301                 r = buckets[i][buckets[i].size()-1];
    00302                 if (buckets[i].size() > 1)
    00303                 {
    00304                         for (int j = (int)buckets[i].size()-2; j >= 1; j--)
    00305                         {
    00306                                 Accumulate(buckets[i][j], buckets[i][j+1]);
    00307                                 Accumulate(r, buckets[i][j]);
    00308                         }
    00309                         Accumulate(buckets[i][0], buckets[i][1]);
    00310                         r = Add(Double(r), buckets[i][0]);
    00311                 }
    00312         }
    00313 }
    00314
    00315 template <class T> T AbstractRing<T>::Exponentiate(const Element &base, const Integer &exponent) const
    00316 {
    00317         Element result;
    00318         SimultaneousExponentiate(&result, base, &exponent, 1);
    00319         return result;
    00320 }
    00321
    00322 template <class T> T AbstractRing<T>::CascadeExponentiate(const Element &x, const Integer &e1, const Element &y, const Integer &e2) const
    00323 {
    00324         return MultiplicativeGroup().AbstractGroup<T>::CascadeScalarMultiply(x, e1, y, e2);
    00325 }
    00326
    00327 template <class Element, class Iterator> Element GeneralCascadeExponentiation(const AbstractRing<Element> &ring, Iterator begin, Iterator end)
    00328 {
    00329         return GeneralCascadeMultiplication<Element>(ring.MultiplicativeGroup(), begin, end);
    00330 }
    00331
    00332 template <class T>
    00333 void AbstractRing<T>::SimultaneousExponentiate(T *results, const T &base, const Integer *exponents, unsigned int expCount) const
    00334 {
    00335         MultiplicativeGroup().AbstractGroup<T>::SimultaneousMultiply(results, base, exponents, expCount);
    00336 }
    00337
    00338 NAMESPACE_END
    00339
    00340 #endif