Crypto++  5.6.3
Free C++ class library of cryptographic schemes
algebra.cpp
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 
13 NAMESPACE_BEGIN(CryptoPP)
14 
15 template <class T> const T& AbstractGroup<T>::Double(const Element &a) const
16 {
17  return this->Add(a, a);
18 }
19 
20 template <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 
27 template <class T> T& AbstractGroup<T>::Accumulate(Element &a, const Element &b) const
28 {
29  return a = this->Add(a, b);
30 }
31 
32 template <class T> T& AbstractGroup<T>::Reduce(Element &a, const Element &b) const
33 {
34  return a = this->Subtract(a, b);
35 }
36 
37 template <class T> const T& AbstractRing<T>::Square(const Element &a) const
38 {
39  return this->Multiply(a, a);
40 }
41 
42 template <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 
49 template <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 
56 template <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 
70 template <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 
90 template <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 
97 template <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 
172 template <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 
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 
255 template <class T>
256 void 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  assert(expBegin->NotNegative());
266  exponents.push_back(WindowSlider(*expBegin++, InversionIsFast(), 0));
267  exponents[i].FindNextWindow();
268  buckets[i].resize(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 
316 template <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 
323 template <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 
328 template <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 
333 template <class T>
334 void 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 
339 NAMESPACE_END
340 
341 #endif
bool GetBit(size_t i) const
Provides the i-th bit of the Integer.
Definition: integer.cpp:2971
unsigned int WordCount() const
Determines the number of words required to represent the Integer.
Definition: integer.cpp:3191
Abstract Euclidean Domain.
Definition: algebra.h:148
Classes for performing mathematics over different fields.
Quotient Ring.
Definition: algebra.h:228
static const Integer & One()
Integer representing 1.
Definition: integer.cpp:2944
Abstract Ring.
Definition: algebra.h:52
unsigned int BitCount() const
Determines the number of bits required to represent the Integer.
Definition: integer.cpp:3205
Multiple precision integer with arithmetic operations.
Definition: integer.h:31
static void Divide(Integer &r, Integer &q, const Integer &a, const Integer &d)
calculate r and q such that (a == d*q + r) && (0 <= r < abs(d))
Definition: integer.cpp:3899
Abstract Group.
Definition: algebra.h:27
const T & STDMAX(const T &a, const T &b)
Replacement function for std::max.
Definition: misc.h:407
Crypto++ library namespace.
bool NotNegative() const
Determines if the Integer is non-negative.
Definition: integer.h:321