source: trunk/src-cryptopp/algebra.cpp

Last change on this file was e230cb0, checked in by David Stainton <dstainton415@…>, at 2016-10-12T13:27:29Z

Add cryptopp from tag CRYPTOPP_5_6_5
  • Property mode set to 100644
File size: 9.4 KB
Line 
1// algebra.cpp - written and placed in the public domain by Wei Dai
2
3#include "pch.h"
4
5#ifndef CRYPTOPP_ALGEBRA_CPP    // SunCC workaround: compiler could cause this file to be included twice
6#define CRYPTOPP_ALGEBRA_CPP
7
8#include "algebra.h"
9#include "integer.h"
10
11#include <vector>
12
13NAMESPACE_BEGIN(CryptoPP)
14
15template <class T> const T& AbstractGroup<T>::Double(const Element &a) const
16{
17        return this->Add(a, a);
18}
19
20template <class T> const T& AbstractGroup<T>::Subtract(const Element &a, const Element &b) const
21{
22        // make copy of a in case Inverse() overwrites it
23        Element a1(a);
24        return this->Add(a1, Inverse(b));
25}
26
27template <class T> T& AbstractGroup<T>::Accumulate(Element &a, const Element &b) const
28{
29        return a = this->Add(a, b);
30}
31
32template <class T> T& AbstractGroup<T>::Reduce(Element &a, const Element &b) const
33{
34        return a = this->Subtract(a, b);
35}
36
37template <class T> const T& AbstractRing<T>::Square(const Element &a) const
38{
39        return this->Multiply(a, a);
40}
41
42template <class T> const T& AbstractRing<T>::Divide(const Element &a, const Element &b) const
43{
44        // make copy of a in case MultiplicativeInverse() overwrites it
45        Element a1(a);
46        return this->Multiply(a1, this->MultiplicativeInverse(b));
47}
48
49template <class T> const T& AbstractEuclideanDomain<T>::Mod(const Element &a, const Element &b) const
50{
51        Element q;
52        this->DivisionAlgorithm(result, q, a, b);
53        return result;
54}
55
56template <class T> const T& AbstractEuclideanDomain<T>::Gcd(const Element &a, const Element &b) const
57{
58        Element g[3]={b, a};
59        unsigned int i0=0, i1=1, i2=2;
60
61        while (!this->Equal(g[i1], this->Identity()))
62        {
63                g[i2] = this->Mod(g[i0], g[i1]);
64                unsigned int t = i0; i0 = i1; i1 = i2; i2 = t;
65        }
66
67        return result = g[i0];
68}
69
70template <class T> const typename QuotientRing<T>::Element& QuotientRing<T>::MultiplicativeInverse(const Element &a) const
71{
72        Element g[3]={m_modulus, a};
73        Element v[3]={m_domain.Identity(), m_domain.MultiplicativeIdentity()};
74        Element y;
75        unsigned int i0=0, i1=1, i2=2;
76
77        while (!this->Equal(g[i1], this->Identity()))
78        {
79                // y = g[i0] / g[i1];
80                // g[i2] = g[i0] % g[i1];
81                m_domain.DivisionAlgorithm(g[i2], y, g[i0], g[i1]);
82                // v[i2] = v[i0] - (v[i1] * y);
83                v[i2] = m_domain.Subtract(v[i0], m_domain.Multiply(v[i1], y));
84                unsigned int t = i0; i0 = i1; i1 = i2; i2 = t;
85        }
86
87        return m_domain.IsUnit(g[i0]) ? m_domain.Divide(v[i0], g[i0]) : m_domain.Identity();
88}
89
90template <class T> T AbstractGroup<T>::ScalarMultiply(const Element &base, const Integer &exponent) const
91{
92        Element result;
93        this->SimultaneousMultiply(&result, base, &exponent, 1);
94        return result;
95}
96
97template <class T> T AbstractGroup<T>::CascadeScalarMultiply(const Element &x, const Integer &e1, const Element &y, const Integer &e2) const
98{
99        const unsigned expLen = STDMAX(e1.BitCount(), e2.BitCount());
100        if (expLen==0)
101                return this->Identity();
102
103        const unsigned w = (expLen <= 46 ? 1 : (expLen <= 260 ? 2 : 3));
104        const unsigned tableSize = 1<<w;
105        std::vector<Element> powerTable(tableSize << w);
106
107        powerTable[1] = x;
108        powerTable[tableSize] = y;
109        if (w==1)
110                powerTable[3] = this->Add(x,y);
111        else
112        {
113                powerTable[2] = this->Double(x);
114                powerTable[2*tableSize] = this->Double(y);
115
116                unsigned i, j;
117
118                for (i=3; i<tableSize; i+=2)
119                        powerTable[i] = Add(powerTable[i-2], powerTable[2]);
120                for (i=1; i<tableSize; i+=2)
121                        for (j=i+tableSize; j<(tableSize<<w); j+=tableSize)
122                                powerTable[j] = Add(powerTable[j-tableSize], y);
123
124                for (i=3*tableSize; i<(tableSize<<w); i+=2*tableSize)
125                        powerTable[i] = Add(powerTable[i-2*tableSize], powerTable[2*tableSize]);
126                for (i=tableSize; i<(tableSize<<w); i+=2*tableSize)
127                        for (j=i+2; j<i+tableSize; j+=2)
128                                powerTable[j] = Add(powerTable[j-1], x);
129        }
130
131        Element result;
132        unsigned power1 = 0, power2 = 0, prevPosition = expLen-1;
133        bool firstTime = true;
134
135        for (int i = expLen-1; i>=0; i--)
136        {
137                power1 = 2*power1 + e1.GetBit(i);
138                power2 = 2*power2 + e2.GetBit(i);
139
140                if (i==0 || 2*power1 >= tableSize || 2*power2 >= tableSize)
141                {
142                        unsigned squaresBefore = prevPosition-i;
143                        unsigned squaresAfter = 0;
144                        prevPosition = i;
145                        while ((power1 || power2) && power1%2 == 0 && power2%2==0)
146                        {
147                                power1 /= 2;
148                                power2 /= 2;
149                                squaresBefore--;
150                                squaresAfter++;
151                        }
152                        if (firstTime)
153                        {
154                                result = powerTable[(power2<<w) + power1];
155                                firstTime = false;
156                        }
157                        else
158                        {
159                                while (squaresBefore--)
160                                        result = this->Double(result);
161                                if (power1 || power2)
162                                        Accumulate(result, powerTable[(power2<<w) + power1]);
163                        }
164                        while (squaresAfter--)
165                                result = this->Double(result);
166                        power1 = power2 = 0;
167                }
168        }
169        return result;
170}
171
172template <class Element, class Iterator> Element GeneralCascadeMultiplication(const AbstractGroup<Element> &group, Iterator begin, Iterator end)
173{
174        if (end-begin == 1)
175                return group.ScalarMultiply(begin->base, begin->exponent);
176        else if (end-begin == 2)
177                return group.CascadeScalarMultiply(begin->base, begin->exponent, (begin+1)->base, (begin+1)->exponent);
178        else
179        {
180                Integer q, t;
181                Iterator last = end;
182                --last;
183
184                std::make_heap(begin, end);
185                std::pop_heap(begin, end);
186
187                while (!!begin->exponent)
188                {
189                        // last->exponent is largest exponent, begin->exponent is next largest
190                        t = last->exponent;
191                        Integer::Divide(last->exponent, q, t, begin->exponent);
192
193                        if (q == Integer::One())
194                                group.Accumulate(begin->base, last->base);      // avoid overhead of ScalarMultiply()
195                        else
196                                group.Accumulate(begin->base, group.ScalarMultiply(last->base, q));
197
198                        std::push_heap(begin, end);
199                        std::pop_heap(begin, end);
200                }
201
202                return group.ScalarMultiply(last->base, last->exponent);
203        }
204}
205
206struct WindowSlider
207{
208        WindowSlider(const Integer &expIn, bool fastNegate, unsigned int windowSizeIn=0)
209                : exp(expIn), windowModulus(Integer::One()), windowSize(windowSizeIn), windowBegin(0), expWindow(0)
210                , fastNegate(fastNegate), negateNext(false), firstTime(true), finished(false)
211        {
212                if (windowSize == 0)
213                {
214                        unsigned int expLen = exp.BitCount();
215                        windowSize = expLen <= 17 ? 1 : (expLen <= 24 ? 2 : (expLen <= 70 ? 3 : (expLen <= 197 ? 4 : (expLen <= 539 ? 5 : (expLen <= 1434 ? 6 : 7)))));
216                }
217                windowModulus <<= windowSize;
218        }
219
220        void FindNextWindow()
221        {
222                unsigned int expLen = exp.WordCount() * WORD_BITS;
223                unsigned int skipCount = firstTime ? 0 : windowSize;
224                firstTime = false;
225                while (!exp.GetBit(skipCount))
226                {
227                        if (skipCount >= expLen)
228                        {
229                                finished = true;
230                                return;
231                        }
232                        skipCount++;
233                }
234
235                exp >>= skipCount;
236                windowBegin += skipCount;
237                expWindow = word32(exp % (word(1) << windowSize));
238
239                if (fastNegate && exp.GetBit(windowSize))
240                {
241                        negateNext = true;
242                        expWindow = (word32(1) << windowSize) - expWindow;
243                        exp += windowModulus;
244                }
245                else
246                        negateNext = false;
247        }
248
249        Integer exp, windowModulus;
250        unsigned int windowSize, windowBegin;
251        word32 expWindow;
252        bool fastNegate, negateNext, firstTime, finished;
253};
254
255template <class T>
256void AbstractGroup<T>::SimultaneousMultiply(T *results, const T &base, const Integer *expBegin, unsigned int expCount) const
257{
258        std::vector<std::vector<Element> > buckets(expCount);
259        std::vector<WindowSlider> exponents;
260        exponents.reserve(expCount);
261        unsigned int i;
262
263        for (i=0; i<expCount; i++)
264        {
265                CRYPTOPP_ASSERT(expBegin->NotNegative());
266                exponents.push_back(WindowSlider(*expBegin++, InversionIsFast(), 0));
267                exponents[i].FindNextWindow();
268                buckets[i].resize(((size_t) 1) << (exponents[i].windowSize-1), Identity());
269        }
270
271        unsigned int expBitPosition = 0;
272        Element g = base;
273        bool notDone = true;
274
275        while (notDone)
276        {
277                notDone = false;
278                for (i=0; i<expCount; i++)
279                {
280                        if (!exponents[i].finished && expBitPosition == exponents[i].windowBegin)
281                        {
282                                Element &bucket = buckets[i][exponents[i].expWindow/2];
283                                if (exponents[i].negateNext)
284                                        Accumulate(bucket, Inverse(g));
285                                else
286                                        Accumulate(bucket, g);
287                                exponents[i].FindNextWindow();
288                        }
289                        notDone = notDone || !exponents[i].finished;
290                }
291
292                if (notDone)
293                {
294                        g = Double(g);
295                        expBitPosition++;
296                }
297        }
298
299        for (i=0; i<expCount; i++)
300        {
301                Element &r = *results++;
302                r = buckets[i][buckets[i].size()-1];
303                if (buckets[i].size() > 1)
304                {
305                        for (int j = (int)buckets[i].size()-2; j >= 1; j--)
306                        {
307                                Accumulate(buckets[i][j], buckets[i][j+1]);
308                                Accumulate(r, buckets[i][j]);
309                        }
310                        Accumulate(buckets[i][0], buckets[i][1]);
311                        r = Add(Double(r), buckets[i][0]);
312                }
313        }
314}
315
316template <class T> T AbstractRing<T>::Exponentiate(const Element &base, const Integer &exponent) const
317{
318        Element result;
319        SimultaneousExponentiate(&result, base, &exponent, 1);
320        return result;
321}
322
323template <class T> T AbstractRing<T>::CascadeExponentiate(const Element &x, const Integer &e1, const Element &y, const Integer &e2) const
324{
325        return MultiplicativeGroup().AbstractGroup<T>::CascadeScalarMultiply(x, e1, y, e2);
326}
327
328template <class Element, class Iterator> Element GeneralCascadeExponentiation(const AbstractRing<Element> &ring, Iterator begin, Iterator end)
329{
330        return GeneralCascadeMultiplication<Element>(ring.MultiplicativeGroup(), begin, end);
331}
332
333template <class T>
334void AbstractRing<T>::SimultaneousExponentiate(T *results, const T &base, const Integer *exponents, unsigned int expCount) const
335{
336        MultiplicativeGroup().AbstractGroup<T>::SimultaneousMultiply(results, base, exponents, expCount);
337}
338
339NAMESPACE_END
340
341#endif
Note: See TracBrowser for help on using the repository browser.