monotone

monotone Mtn Source Tree

Root/cryptopp/algebra.cpp

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

Archive Download this file

Branches

Tags

Quick Links:     www.monotone.ca    -     Downloads    -     Documentation    -     Wiki    -     Code Forge    -     Build Status