Clearmatics Libff  0.1
C++ library for Finite Fields and Elliptic Curves
fp.tcc
Go to the documentation of this file.
1 /** @file
2  *****************************************************************************
3  Implementation of arithmetic in the finite field F[p], for prime p of fixed
4  length.
5  *****************************************************************************
6  * @author This file is part of libff, developed by SCIPR Lab
7  * and contributors (see AUTHORS).
8  * @copyright MIT license (see LICENSE file)
9  *****************************************************************************/
10 
11 #ifndef FP_TCC_
12 #define FP_TCC_
13 #include <cassert>
14 #include <cmath>
15 #include <cstdlib>
16 #include <libff/algebra/fields/field_serialization.hpp>
17 #include <libff/algebra/fields/field_utils.hpp>
18 #include <libff/algebra/fields/fp_aux.tcc>
19 #include <limits>
20 
21 namespace libff
22 {
23 
24 template<mp_size_t n, const bigint<n> &modulus>
25 bool Fp_model<n, modulus>::s_initialized = false;
26 
27 template<mp_size_t n, const bigint<n> &modulus>
28 Fp_model<n, modulus> Fp_model<n, modulus>::s_zero;
29 
30 template<mp_size_t n, const bigint<n> &modulus>
31 Fp_model<n, modulus> Fp_model<n, modulus>::s_one;
32 
33 template<mp_size_t n, const bigint<n> &modulus>
34 void Fp_model<n, modulus>::static_init()
35 {
36  // This function may be entered several times, in particular where an
37  // individual field is shread between curves.
38  if (s_initialized) {
39  return;
40  }
41  // Initialize s_zero and s_one
42  s_zero.mont_repr.clear();
43 
44  s_one.mont_repr.data[0] = 1;
45  s_one.mul_reduce(Rsquared);
46 
47  s_initialized = true;
48 }
49 
50 template<mp_size_t n, const bigint<n> &modulus>
51 void Fp_model<n, modulus>::mul_reduce(const bigint<n> &other)
52 {
53  // stupid pre-processor tricks; beware
54 #if defined(__x86_64__) && defined(USE_ASM)
55  if (n == 3) { // Use asm-optimized Comba multiplication and reduction
56  mp_limb_t res[2 * n];
57  mp_limb_t c0, c1, c2;
58  COMBA_3_BY_3_MUL(c0, c1, c2, res, this->mont_repr.data, other.data);
59 
60  mp_limb_t k;
61  mp_limb_t tmp1, tmp2, tmp3;
62  REDUCE_6_LIMB_PRODUCT(k, tmp1, tmp2, tmp3, inv, res, modulus.data);
63 
64  // subtract t > mod
65  __asm__( // Preserve alignment
66  "/* check for overflow */ \n\t" //
67  MONT_CMP(16) //
68  MONT_CMP(8) //
69  MONT_CMP(0) //
70  "/* subtract mod if overflow */ \n\t" //
71  "subtract%=: \n\t" //
72  MONT_FIRSTSUB() //
73  MONT_NEXTSUB(8) //
74  MONT_NEXTSUB(16) //
75  "done%=: \n\t" //
76  :
77  : [tmp] "r"(res + n), [M] "r"(modulus.data)
78  : "cc", "memory", "%rax");
79  mpn_copyi(this->mont_repr.data, res + n, n);
80  } else if (n == 4) { // use asm-optimized "CIOS method"
81 
82  mp_limb_t tmp[n + 1];
83  mp_limb_t T0 = 0, T1 = 1, cy = 2, u = 3; // TODO: fix this
84 
85  __asm__( // Preserve alignment
86  MONT_PRECOMPUTE() //
87  MONT_FIRSTITER(1) //
88  MONT_FIRSTITER(2) //
89  MONT_FIRSTITER(3) //
90  MONT_FINALIZE(3) //
91  MONT_ITERFIRST(1) //
92  MONT_ITERITER(1, 1) //
93  MONT_ITERITER(1, 2) //
94  MONT_ITERITER(1, 3) //
95  MONT_FINALIZE(3) //
96  MONT_ITERFIRST(2) //
97  MONT_ITERITER(2, 1) //
98  MONT_ITERITER(2, 2) //
99  MONT_ITERITER(2, 3) //
100  MONT_FINALIZE(3) //
101  MONT_ITERFIRST(3) //
102  MONT_ITERITER(3, 1) //
103  MONT_ITERITER(3, 2) //
104  MONT_ITERITER(3, 3) //
105  MONT_FINALIZE(3) //
106  "/* check for overflow */ " //
107  "\n\t" //
108  MONT_CMP(24) //
109  MONT_CMP(16) //
110  MONT_CMP(8) //
111  MONT_CMP(0) //
112  "/* subtract mod if overflow */ " //
113  "\n\t" //
114  "subtract%=: " //
115  "\n\t" //
116  MONT_FIRSTSUB() //
117  MONT_NEXTSUB(8) //
118  MONT_NEXTSUB(16) //
119  MONT_NEXTSUB(24) //
120  "done%=: " //
121  " \n\t" //
122  :
123  : [tmp] "r"(tmp),
124  [A] "r"(this->mont_repr.data),
125  [B] "r"(other.data),
126  [inv] "r"(inv),
127  [M] "r"(modulus.data),
128  [T0] "r"(T0),
129  [T1] "r"(T1),
130  [cy] "r"(cy),
131  [u] "r"(u)
132  : "cc", "memory", "%rax", "%rdx");
133  mpn_copyi(this->mont_repr.data, tmp, n);
134  } else if (n == 5) {
135  // use asm-optimized "CIOS method"
136  mp_limb_t tmp[n + 1];
137  mp_limb_t T0 = 0, T1 = 1, cy = 2, u = 3; // TODO: fix this
138 
139  __asm__( // Preserve alignment
140  MONT_PRECOMPUTE() //
141  MONT_FIRSTITER(1) //
142  MONT_FIRSTITER(2) //
143  MONT_FIRSTITER(3) //
144  MONT_FIRSTITER(4) //
145  MONT_FINALIZE(4) //
146  MONT_ITERFIRST(1) //
147  MONT_ITERITER(1, 1) //
148  MONT_ITERITER(1, 2) //
149  MONT_ITERITER(1, 3) //
150  MONT_ITERITER(1, 4) //
151  MONT_FINALIZE(4) //
152  MONT_ITERFIRST(2) //
153  MONT_ITERITER(2, 1) //
154  MONT_ITERITER(2, 2) //
155  MONT_ITERITER(2, 3) //
156  MONT_ITERITER(2, 4) //
157  MONT_FINALIZE(4) //
158  MONT_ITERFIRST(3) //
159  MONT_ITERITER(3, 1) //
160  MONT_ITERITER(3, 2) //
161  MONT_ITERITER(3, 3) //
162  MONT_ITERITER(3, 4) //
163  MONT_FINALIZE(4) //
164  MONT_ITERFIRST(4) //
165  MONT_ITERITER(4, 1) //
166  MONT_ITERITER(4, 2) //
167  MONT_ITERITER(4, 3) //
168  MONT_ITERITER(4, 4) //
169  MONT_FINALIZE(4) //
170  "/* check for overflow */ " //
171  "\n\t" //
172  MONT_CMP(32) //
173  MONT_CMP(24) //
174  MONT_CMP(16) //
175  MONT_CMP(8) //
176  MONT_CMP(0) //
177  "/* subtract mod if overflow */ " //
178  " \n\t" //
179  "subtract%=: " //
180  " \n\t" //
181  MONT_FIRSTSUB() //
182  MONT_NEXTSUB(8) //
183  MONT_NEXTSUB(16) //
184  MONT_NEXTSUB(24) //
185  MONT_NEXTSUB( //
186  32) //
187  "done%=: " //
188  " " //
189  " \n\t" //
190  :
191  : [tmp] "r"(tmp),
192  [A] "r"(this->mont_repr.data),
193  [B] "r"(other.data),
194  [inv] "r"(inv),
195  [M] "r"(modulus.data),
196  [T0] "r"(T0),
197  [T1] "r"(T1),
198  [cy] "r"(cy),
199  [u] "r"(u)
200  : "cc", "memory", "%rax", "%rdx");
201  mpn_copyi(this->mont_repr.data, tmp, n);
202  } else
203 #endif
204  {
205  mp_limb_t res[2 * n];
206  mpn_mul_n(res, this->mont_repr.data, other.data, n);
207 
208  // The Montgomery reduction here is based on Algorithm 14.32 in
209  // Handbook of Applied Cryptography
210  // <http://cacr.uwaterloo.ca/hac/about/chap14.pdf>.
211  for (size_t i = 0; i < n; ++i) {
212  mp_limb_t k = inv * res[i];
213  // calculate res = res + k * mod * b^i
214  mp_limb_t carryout = mpn_addmul_1(res + i, modulus.data, n, k);
215  carryout = mpn_add_1(res + n + i, res + n + i, n - i, carryout);
216  assert(carryout == 0);
217  }
218 
219  if (mpn_cmp(res + n, modulus.data, n) >= 0) {
220  const mp_limb_t borrow =
221  mpn_sub(res + n, res + n, n, modulus.data, n);
222  assert(borrow == 0);
223  UNUSED(borrow);
224  }
225 
226  mpn_copyi(this->mont_repr.data, res + n, n);
227  }
228 }
229 
230 template<mp_size_t n, const bigint<n> &modulus>
231 Fp_model<n, modulus>::Fp_model(const bigint<n> &b)
232 {
233  mpn_copyi(this->mont_repr.data, Rsquared.data, n);
234  mul_reduce(b);
235 }
236 
237 template<mp_size_t n, const bigint<n> &modulus>
238 Fp_model<n, modulus>::Fp_model(const long x, const bool is_unsigned)
239 {
240  static_assert(
241  std::numeric_limits<mp_limb_t>::max() >=
242  static_cast<unsigned long>(std::numeric_limits<long>::max()),
243  "long won't fit in mp_limb_t");
244  if (is_unsigned || x >= 0) {
245  this->mont_repr.data[0] = (mp_limb_t)x;
246  } else {
247  const mp_limb_t borrow =
248  mpn_sub_1(this->mont_repr.data, modulus.data, n, (mp_limb_t)-x);
249  assert(borrow == 0);
250  UNUSED(borrow);
251  }
252 
253  mul_reduce(Rsquared);
254 }
255 
256 template<mp_size_t n, const bigint<n> &modulus>
257 void Fp_model<n, modulus>::set_ulong(const unsigned long x)
258 {
259  this->mont_repr.clear();
260  this->mont_repr.data[0] = x;
261  mul_reduce(Rsquared);
262 }
263 
264 template<mp_size_t n, const bigint<n> &modulus>
265 void Fp_model<n, modulus>::clear()
266 {
267  this->mont_repr.clear();
268 }
269 
270 template<mp_size_t n, const bigint<n> &modulus>
271 bigint<n> Fp_model<n, modulus>::as_bigint() const
272 {
273  bigint<n> one;
274  one.clear();
275  one.data[0] = 1;
276 
277  Fp_model<n, modulus> res(*this);
278  res.mul_reduce(one);
279 
280  return (res.mont_repr);
281 }
282 
283 template<mp_size_t n, const bigint<n> &modulus>
284 unsigned long Fp_model<n, modulus>::as_ulong() const
285 {
286  return this->as_bigint().as_ulong();
287 }
288 
289 template<mp_size_t n, const bigint<n> &modulus>
290 bool Fp_model<n, modulus>::operator==(const Fp_model &other) const
291 {
292  return (this->mont_repr == other.mont_repr);
293 }
294 
295 template<mp_size_t n, const bigint<n> &modulus>
296 bool Fp_model<n, modulus>::operator!=(const Fp_model &other) const
297 {
298  return (this->mont_repr != other.mont_repr);
299 }
300 
301 template<mp_size_t n, const bigint<n> &modulus>
302 bool Fp_model<n, modulus>::is_zero() const
303 {
304  // zero maps to zero
305  return (this->mont_repr.is_zero());
306 }
307 
308 template<mp_size_t n, const bigint<n> &modulus>
309 void Fp_model<n, modulus>::print() const
310 {
311  Fp_model<n, modulus> tmp;
312  tmp.mont_repr.data[0] = 1;
313  tmp.mul_reduce(this->mont_repr);
314 
315  tmp.mont_repr.print();
316 }
317 
318 template<mp_size_t n, const bigint<n> &modulus>
319 const Fp_model<n, modulus> &Fp_model<n, modulus>::zero()
320 {
321  assert(s_initialized);
322  return s_zero;
323 }
324 
325 template<mp_size_t n, const bigint<n> &modulus>
326 const Fp_model<n, modulus> &Fp_model<n, modulus>::one()
327 {
328  assert(s_initialized);
329  return s_one;
330 }
331 
332 template<mp_size_t n, const bigint<n> &modulus>
333 Fp_model<n, modulus> Fp_model<n, modulus>::geometric_generator()
334 {
335  Fp_model<n, modulus> res;
336  res.mont_repr.data[0] = 2;
337  res.mul_reduce(Rsquared);
338  return res;
339 }
340 
341 template<mp_size_t n, const bigint<n> &modulus>
342 Fp_model<n, modulus> Fp_model<n, modulus>::arithmetic_generator()
343 {
344  Fp_model<n, modulus> res;
345  res.mont_repr.data[0] = 1;
346  res.mul_reduce(Rsquared);
347  return res;
348 }
349 
350 template<mp_size_t n, const bigint<n> &modulus>
351 Fp_model<n, modulus> &Fp_model<n, modulus>::operator+=(
352  const Fp_model<n, modulus> &other)
353 {
354 #ifdef PROFILE_OP_COUNTS
355  this->add_cnt++;
356 #endif
357 #if defined(__x86_64__) && defined(USE_ASM)
358  if (n == 3) {
359  __asm__( // Preserve alignment
360  "/* perform bignum addition */ \n\t" //
361  ADD_FIRSTADD() //
362  ADD_NEXTADD(8) //
363  ADD_NEXTADD(16) //
364  "/* if overflow: subtract */ \n\t" //
365  "/* (tricky point: if A and B are in the range we do " //
366  "not need to do anything special for the possible " //
367  "carry flag) */ \n\t" //
368  "jc subtract%= \n\t" //
369  "/* check for overflow */ \n\t" //
370  ADD_CMP(16) //
371  ADD_CMP(8) //
372  ADD_CMP(0) //
373  "/* subtract mod if overflow */ \n\t" //
374  "subtract%=: \n\t" //
375  ADD_FIRSTSUB() //
376  ADD_NEXTSUB(8) //
377  ADD_NEXTSUB(16) //
378  "done%=: \n\t" //
379  :
380  : [A] "r"(this->mont_repr.data),
381  [B] "r"(other.mont_repr.data),
382  [mod] "r"(modulus.data)
383  : "cc", "memory", "%rax");
384  } else if (n == 4) {
385  __asm__( // Preserve alignment
386  "/* perform bignum addition */ \n\t" //
387  ADD_FIRSTADD() //
388  ADD_NEXTADD(8) //
389  ADD_NEXTADD(16) //
390  ADD_NEXTADD(24) //
391  "/* if overflow: subtract */ \n\t" //
392  "/* (tricky point: if A and B are in the range we do " //
393  "not need to do anything special for the possible " //
394  "carry flag) */ \n\t" //
395  "jc subtract%= \n\t" //
396  "/* check for overflow */ \n\t" //
397  ADD_CMP(24) //
398  ADD_CMP(16) //
399  ADD_CMP(8) //
400  ADD_CMP(0) //
401  "/* subtract mod if overflow */ \n\t" //
402  "subtract%=: \n\t" //
403  ADD_FIRSTSUB() //
404  ADD_NEXTSUB(8) //
405  ADD_NEXTSUB(16) //
406  ADD_NEXTSUB(24) //
407  "done%=: \n\t" //
408  :
409  : [A] "r"(this->mont_repr.data),
410  [B] "r"(other.mont_repr.data),
411  [mod] "r"(modulus.data)
412  : "cc", "memory", "%rax");
413  } else if (n == 5) {
414  __asm__( // Preserve alignment
415  "/* perform bignum addition */ \n\t" //
416  ADD_FIRSTADD() //
417  ADD_NEXTADD(8) //
418  ADD_NEXTADD(16) //
419  ADD_NEXTADD(24) //
420  ADD_NEXTADD(32) //
421  "/* if overflow: subtract */ \n\t" //
422  "/* (tricky point: if A and B are in the range we do " //
423  "not need to do anything special for the possible " //
424  "carry flag) */ \n\t" //
425  "jc subtract%= \n\t" //
426  "/* check for overflow */ \n\t" //
427  ADD_CMP(32) //
428  ADD_CMP(24) //
429  ADD_CMP(16) //
430  ADD_CMP(8) //
431  ADD_CMP(0) //
432  "/* subtract mod if overflow */ \n\t" //
433  "subtract%=: \n\t" //
434  ADD_FIRSTSUB() //
435  ADD_NEXTSUB(8) //
436  ADD_NEXTSUB(16) //
437  ADD_NEXTSUB(24) //
438  ADD_NEXTSUB(32) //
439  "done%=: " //
440  " \n\t" //
441  :
442  : [A] "r"(this->mont_repr.data),
443  [B] "r"(other.mont_repr.data),
444  [mod] "r"(modulus.data)
445  : "cc", "memory", "%rax");
446  } else
447 #endif
448  {
449  mp_limb_t scratch[n + 1];
450  const mp_limb_t carry =
451  mpn_add_n(scratch, this->mont_repr.data, other.mont_repr.data, n);
452  scratch[n] = carry;
453 
454  if (carry || mpn_cmp(scratch, modulus.data, n) >= 0) {
455  const mp_limb_t borrow =
456  mpn_sub(scratch, scratch, n + 1, modulus.data, n);
457  assert(borrow == 0);
458  UNUSED(borrow);
459  }
460 
461  mpn_copyi(this->mont_repr.data, scratch, n);
462  }
463 
464  return *this;
465 }
466 
467 template<mp_size_t n, const bigint<n> &modulus>
468 Fp_model<n, modulus> &Fp_model<n, modulus>::operator-=(
469  const Fp_model<n, modulus> &other)
470 {
471 #ifdef PROFILE_OP_COUNTS
472  this->sub_cnt++;
473 #endif
474 #if defined(__x86_64__) && defined(USE_ASM)
475  if (n == 3) {
476  __asm__( // Preserve alignment
477  SUB_FIRSTSUB() //
478  SUB_NEXTSUB(8) //
479  SUB_NEXTSUB(16) //
480  "jnc done%=\n\t" //
481  SUB_FIRSTADD() //
482  SUB_NEXTADD(8) //
483  SUB_NEXTADD(16) //
484  "done%=:\n\t" //
485  :
486  : [A] "r"(this->mont_repr.data),
487  [B] "r"(other.mont_repr.data),
488  [mod] "r"(modulus.data)
489  : "cc", "memory", "%rax");
490  } else if (n == 4) {
491  __asm__( // Preserve alignment
492  SUB_FIRSTSUB() //
493  SUB_NEXTSUB(8) //
494  SUB_NEXTSUB(16) //
495  SUB_NEXTSUB(24) //
496  "jnc done%=\n\t" //
497  SUB_FIRSTADD() //
498  SUB_NEXTADD(8) //
499  SUB_NEXTADD(16) //
500  SUB_NEXTADD(24) //
501  "done%=:\n\t" //
502  :
503  : [A] "r"(this->mont_repr.data),
504  [B] "r"(other.mont_repr.data),
505  [mod] "r"(modulus.data)
506  : "cc", "memory", "%rax");
507  } else if (n == 5) {
508  __asm__( // Preserve alignment
509  SUB_FIRSTSUB() //
510  SUB_NEXTSUB(8) //
511  SUB_NEXTSUB(16) //
512  SUB_NEXTSUB(24) //
513  SUB_NEXTSUB(32) //
514  "jnc done%=\n\t" //
515  SUB_FIRSTADD() //
516  SUB_NEXTADD(8) //
517  SUB_NEXTADD(16) //
518  SUB_NEXTADD(24) //
519  SUB_NEXTADD(32) //
520  "done%=:\n\t" //
521  :
522  : [A] "r"(this->mont_repr.data),
523  [B] "r"(other.mont_repr.data),
524  [mod] "r"(modulus.data)
525  : "cc", "memory", "%rax");
526  } else
527 #endif
528  {
529  mp_limb_t scratch[n + 1];
530  if (mpn_cmp(this->mont_repr.data, other.mont_repr.data, n) < 0) {
531  const mp_limb_t carry =
532  mpn_add_n(scratch, this->mont_repr.data, modulus.data, n);
533  scratch[n] = carry;
534  } else {
535  mpn_copyi(scratch, this->mont_repr.data, n);
536  scratch[n] = 0;
537  }
538 
539  const mp_limb_t borrow =
540  mpn_sub(scratch, scratch, n + 1, other.mont_repr.data, n);
541  assert(borrow == 0);
542  UNUSED(borrow);
543 
544  mpn_copyi(this->mont_repr.data, scratch, n);
545  }
546  return *this;
547 }
548 
549 template<mp_size_t n, const bigint<n> &modulus>
550 Fp_model<n, modulus> &Fp_model<n, modulus>::operator*=(
551  const Fp_model<n, modulus> &other)
552 {
553 #ifdef PROFILE_OP_COUNTS
554  this->mul_cnt++;
555 #endif
556 
557  mul_reduce(other.mont_repr);
558  return *this;
559 }
560 
561 template<mp_size_t n, const bigint<n> &modulus>
562 Fp_model<n, modulus> &Fp_model<n, modulus>::operator^=(const unsigned long pow)
563 {
564  (*this) = power<Fp_model<n, modulus>>(*this, pow);
565  return (*this);
566 }
567 
568 template<mp_size_t n, const bigint<n> &modulus>
569 template<mp_size_t m>
570 Fp_model<n, modulus> &Fp_model<n, modulus>::operator^=(const bigint<m> &pow)
571 {
572  (*this) = power<Fp_model<n, modulus>, m>(*this, pow);
573  return (*this);
574 }
575 
576 template<mp_size_t n, const bigint<n> &modulus>
577 Fp_model<n, modulus> Fp_model<n, modulus>::operator+(
578  const Fp_model<n, modulus> &other) const
579 {
580  Fp_model<n, modulus> r(*this);
581  return (r += other);
582 }
583 
584 template<mp_size_t n, const bigint<n> &modulus>
585 Fp_model<n, modulus> Fp_model<n, modulus>::operator-(
586  const Fp_model<n, modulus> &other) const
587 {
588  Fp_model<n, modulus> r(*this);
589  return (r -= other);
590 }
591 
592 template<mp_size_t n, const bigint<n> &modulus>
593 Fp_model<n, modulus> Fp_model<n, modulus>::operator*(
594  const Fp_model<n, modulus> &other) const
595 {
596  Fp_model<n, modulus> r(*this);
597  return (r *= other);
598 }
599 
600 template<mp_size_t n, const bigint<n> &modulus>
601 Fp_model<n, modulus> Fp_model<n, modulus>::operator^(
602  const unsigned long pow) const
603 {
604  Fp_model<n, modulus> r(*this);
605  return (r ^= pow);
606 }
607 
608 template<mp_size_t n, const bigint<n> &modulus>
609 template<mp_size_t m>
610 Fp_model<n, modulus> Fp_model<n, modulus>::operator^(const bigint<m> &pow) const
611 {
612  Fp_model<n, modulus> r(*this);
613  return (r ^= pow);
614 }
615 
616 template<mp_size_t n, const bigint<n> &modulus>
617 Fp_model<n, modulus> Fp_model<n, modulus>::operator-() const
618 {
619 #ifdef PROFILE_OP_COUNTS
620  this->sub_cnt++;
621 #endif
622 
623  if (this->is_zero()) {
624  return (*this);
625  } else {
626  Fp_model<n, modulus> r;
627  mpn_sub_n(r.mont_repr.data, modulus.data, this->mont_repr.data, n);
628  return r;
629  }
630 }
631 
632 template<mp_size_t n, const bigint<n> &modulus>
633 Fp_model<n, modulus> Fp_model<n, modulus>::squared() const
634 {
635 #ifdef PROFILE_OP_COUNTS
636  this->sqr_cnt++;
637  // zero out the upcoming mul
638  this->mul_cnt--;
639 #endif
640  // stupid pre-processor tricks; beware
641 #if defined(__x86_64__) && defined(USE_ASM)
642  // use asm-optimized Comba squaring
643  if (n == 3) {
644  mp_limb_t res[2 * n];
645  mp_limb_t c0, c1, c2;
646  COMBA_3_BY_3_SQR(c0, c1, c2, res, this->mont_repr.data);
647 
648  mp_limb_t k;
649  mp_limb_t tmp1, tmp2, tmp3;
650  REDUCE_6_LIMB_PRODUCT(k, tmp1, tmp2, tmp3, inv, res, modulus.data);
651 
652  // subtract t > mod
653  __asm__ volatile( // Preserve alignment
654  "/* check for overflow */ \n\t" //
655  MONT_CMP(16) //
656  MONT_CMP(8) //
657  MONT_CMP(0) //
658  "/* subtract mod if overflow */ \n\t" //
659  "subtract%=: \n\t" //
660  MONT_FIRSTSUB() //
661  MONT_NEXTSUB(8) //
662  MONT_NEXTSUB(16) //
663  "done%=: \n\t" //
664  :
665  : [tmp] "r"(res + n), [M] "r"(modulus.data)
666  : "cc", "memory", "%rax");
667 
668  Fp_model<n, modulus> r;
669  mpn_copyi(r.mont_repr.data, res + n, n);
670  return r;
671  } else
672 #endif
673  {
674  Fp_model<n, modulus> r(*this);
675  return (r *= r);
676  }
677 }
678 
679 template<mp_size_t n, const bigint<n> &modulus>
680 Fp_model<n, modulus> &Fp_model<n, modulus>::invert()
681 {
682 #ifdef PROFILE_OP_COUNTS
683  this->inv_cnt++;
684 #endif
685 
686  assert(!this->is_zero());
687 
688  // gp should have room for vn = n limbs
689  bigint<n> g;
690 
691  mp_limb_t s[n + 1]; /* sp should have room for vn+1 limbs */
692  mp_size_t sn;
693 
694  // both source operands are destroyed by mpn_gcdext
695  bigint<n> v = modulus;
696 
697  // computes gcd(u, v) = g = u*s + v*t, so s*u will be 1 (mod v)
698  const mp_size_t gn =
699  mpn_gcdext(g.data, s, &sn, this->mont_repr.data, n, v.data, n);
700  assert(gn == 1 && g.data[0] == 1); // inverse exists
701  UNUSED(gn);
702 
703  // division result fits into q, as sn <= n+1
704  // sn < 0 indicates negative sn; will fix up later
705  mp_limb_t q;
706 
707  if (std::abs(sn) >= n) {
708  // if sn could require modulus reduction, do it here
709  mpn_tdiv_qr(
710  &q, this->mont_repr.data, 0, s, std::abs(sn), modulus.data, n);
711  } else {
712  // otherwise just copy it over
713  mpn_zero(this->mont_repr.data, n);
714  mpn_copyi(this->mont_repr.data, s, std::abs(sn));
715  }
716 
717  // fix up the negative sn
718  if (sn < 0) {
719  const mp_limb_t borrow = mpn_sub_n(
720  this->mont_repr.data, modulus.data, this->mont_repr.data, n);
721  assert(borrow == 0);
722  UNUSED(borrow);
723  }
724 
725  mul_reduce(Rcubed);
726  return *this;
727 }
728 
729 template<mp_size_t n, const bigint<n> &modulus>
730 Fp_model<n, modulus> Fp_model<n, modulus>::inverse() const
731 {
732  Fp_model<n, modulus> r(*this);
733  return (r.invert());
734 }
735 
736 template<mp_size_t n, const bigint<n> &modulus>
737 Fp_model<n, modulus> Fp_model<n, modulus>::random_element()
738 {
739  // note that as Montgomery representation is a bijection then
740  // selecting a random element of {xR} is the same as selecting a
741  // random element of {x}
742  Fp_model<n, modulus> r;
743  do {
744  r.mont_repr.randomize();
745 
746  // clear all bits higher than MSB of modulus
747  size_t bitno = GMP_NUMB_BITS * n - 1;
748  while (modulus.test_bit(bitno) == false) {
749  const std::size_t part = bitno / GMP_NUMB_BITS;
750  const std::size_t bit = bitno - (GMP_NUMB_BITS * part);
751 
752  r.mont_repr.data[part] &= ~(1ul << bit);
753 
754  bitno--;
755  }
756  }
757  // if r.data is still >= modulus -- repeat (rejection sampling)
758  while (mpn_cmp(r.mont_repr.data, modulus.data, n) >= 0);
759 
760  return r;
761 }
762 
763 template<mp_size_t n, const bigint<n> &modulus>
764 Fp_model<n, modulus> Fp_model<n, modulus>::sqrt() const
765 {
766  Fp_model<n, modulus> one = Fp_model<n, modulus>::one();
767 
768  size_t v = Fp_model<n, modulus>::s;
769  Fp_model<n, modulus> z = Fp_model<n, modulus>::nqr_to_t;
770  Fp_model<n, modulus> w = (*this) ^ Fp_model<n, modulus>::t_minus_1_over_2;
771  Fp_model<n, modulus> x = (*this) * w;
772  Fp_model<n, modulus> b = x * w; // b = (*this)^t
773 
774 #if DEBUG
775  // check if square with euler's criterion
776  Fp_model<n, modulus> check = b;
777  for (size_t i = 0; i < v - 1; ++i) {
778  check = check.squared();
779  }
780  if (check != one) {
781  assert(0);
782  }
783 #endif
784 
785  // compute square root with Tonelli--Shanks
786  // (does not terminate if not a square!)
787 
788  while (b != one) {
789  size_t m = 0;
790  Fp_model<n, modulus> b2m = b;
791  while (b2m != one) {
792  // invariant: b2m = b^(2^m) after entering this loop
793  b2m = b2m.squared();
794  m += 1;
795  }
796 
797  // w = z^2^(v-m-1)
798  int j = v - m - 1;
799  w = z;
800  while (j > 0) {
801  w = w.squared();
802  --j;
803  }
804 
805  z = w.squared();
806  b = b * z;
807  x = x * w;
808  v = m;
809  }
810 
811  return x;
812 }
813 
814 template<mp_size_t n, const bigint<n> &modulus>
815 std::ostream &operator<<(std::ostream &out, const Fp_model<n, modulus> &p)
816 {
817  field_write<DEFAULT_ENCODING, DEFAULT_FORM>(p, out);
818  return out;
819 }
820 
821 template<mp_size_t n, const bigint<n> &modulus>
822 std::istream &operator>>(std::istream &in, Fp_model<n, modulus> &p)
823 {
824  field_read<DEFAULT_ENCODING, DEFAULT_FORM>(p, in);
825  return in;
826 }
827 
828 } // namespace libff
829 
830 #endif // FP_TCC_
Generated on Thu Aug 18 2022 12:42:19 for Clearmatics Libff by   doxygen 1.8.17