iPXE
bigint.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012 Michael Brown <mbrown@fensystems.co.uk>.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License as
6  * published by the Free Software Foundation; either version 2 of the
7  * License, or any later version.
8  *
9  * This program is distributed in the hope that it will be useful, but
10  * WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12  * General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17  * 02110-1301, USA.
18  *
19  * You can also choose to distribute this program under the terms of
20  * the Unmodified Binary Distribution Licence (as given in the file
21  * COPYING.UBDL), provided that you have satisfied its requirements.
22  */
23 
24 FILE_LICENCE ( GPL2_OR_LATER_OR_UBDL );
25 
26 #include <stdint.h>
27 #include <string.h>
28 #include <assert.h>
29 #include <stdio.h>
30 #include <ipxe/bigint.h>
31 
32 /** @file
33  *
34  * Big integer support
35  */
36 
37 /** Minimum number of least significant bytes included in transcription */
38 #define BIGINT_NTOA_LSB_MIN 16
39 
40 /**
41  * Transcribe big integer (for debugging)
42  *
43  * @v value0 Element 0 of big integer to be transcribed
44  * @v size Number of elements
45  * @ret string Big integer in string form (may be abbreviated)
46  */
47 const char * bigint_ntoa_raw ( const bigint_element_t *value0,
48  unsigned int size ) {
49  const bigint_t ( size ) __attribute__ (( may_alias ))
50  *value = ( ( const void * ) value0 );
52  static char buf[256];
53  unsigned int count;
54  int threshold;
55  uint8_t byte;
56  char *tmp;
57  int i;
58 
59  /* Calculate abbreviation threshold, if any */
60  count = ( size * sizeof ( element ) );
61  threshold = count;
62  if ( count >= ( ( sizeof ( buf ) - 1 /* NUL */ ) / 2 /* "xx" */ ) ) {
63  threshold -= ( ( sizeof ( buf ) - 3 /* "..." */
64  - ( BIGINT_NTOA_LSB_MIN * 2 /* "xx" */ )
65  - 1 /* NUL */ ) / 2 /* "xx" */ );
66  }
67 
68  /* Transcribe bytes, abbreviating with "..." if necessary */
69  for ( tmp = buf, i = ( count - 1 ) ; i >= 0 ; i-- ) {
70  element = value->element[ i / sizeof ( element ) ];
71  byte = ( element >> ( 8 * ( i % sizeof ( element ) ) ) );
72  tmp += sprintf ( tmp, "%02x", byte );
73  if ( i == threshold ) {
74  tmp += sprintf ( tmp, "..." );
76  }
77  }
78  assert ( tmp < ( buf + sizeof ( buf ) ) );
79 
80  return buf;
81 }
82 
83 /**
84  * Conditionally swap big integers (in constant time)
85  *
86  * @v first0 Element 0 of big integer to be conditionally swapped
87  * @v second0 Element 0 of big integer to be conditionally swapped
88  * @v size Number of elements in big integers
89  * @v swap Swap first and second big integers
90  */
92  unsigned int size, int swap ) {
93  bigint_element_t mask;
95  unsigned int i;
96 
97  /* Construct mask */
98  mask = ( ( bigint_element_t ) ( ! swap ) - 1 );
99 
100  /* Conditionally swap elements */
101  for ( i = 0 ; i < size ; i++ ) {
102  xor = ( mask & ( first0[i] ^ second0[i] ) );
103  first0[i] ^= xor;
104  second0[i] ^= xor;
105  }
106 }
107 
108 /**
109  * Multiply big integers
110  *
111  * @v multiplicand0 Element 0 of big integer to be multiplied
112  * @v multiplicand_size Number of elements in multiplicand
113  * @v multiplier0 Element 0 of big integer to be multiplied
114  * @v multiplier_size Number of elements in multiplier
115  * @v result0 Element 0 of big integer to hold result
116  */
117 void bigint_multiply_raw ( const bigint_element_t *multiplicand0,
118  unsigned int multiplicand_size,
119  const bigint_element_t *multiplier0,
120  unsigned int multiplier_size,
121  bigint_element_t *result0 ) {
122  unsigned int result_size = ( multiplicand_size + multiplier_size );
123  const bigint_t ( multiplicand_size ) __attribute__ (( may_alias ))
124  *multiplicand = ( ( const void * ) multiplicand0 );
125  const bigint_t ( multiplier_size ) __attribute__ (( may_alias ))
126  *multiplier = ( ( const void * ) multiplier0 );
127  bigint_t ( result_size ) __attribute__ (( may_alias ))
128  *result = ( ( void * ) result0 );
129  bigint_element_t multiplicand_element;
130  const bigint_element_t *multiplier_element;
131  bigint_element_t *result_element;
132  bigint_element_t carry_element;
133  unsigned int i;
134  unsigned int j;
135 
136  /* Zero required portion of result
137  *
138  * All elements beyond the length of the multiplier will be
139  * written before they are read, and so do not need to be
140  * zeroed in advance.
141  */
142  memset ( result, 0, sizeof ( *multiplier ) );
143 
144  /* Multiply integers one element at a time, adding the low
145  * half of the double-element product directly into the
146  * result, and maintaining a running single-element carry.
147  *
148  * The running carry can never overflow beyond a single
149  * element. At each step, the calculation we perform is:
150  *
151  * carry:result[i+j] := ( ( multiplicand[i] * multiplier[j] )
152  * + result[i+j] + carry )
153  *
154  * The maximum value (for n-bit elements) is therefore:
155  *
156  * (2^n - 1)*(2^n - 1) + (2^n - 1) + (2^n - 1) = 2^(2n) - 1
157  *
158  * This is precisely the maximum value for a 2n-bit integer,
159  * and so the carry out remains within the range of an n-bit
160  * integer, i.e. a single element.
161  */
162  for ( i = 0 ; i < multiplicand_size ; i++ ) {
163  multiplicand_element = multiplicand->element[i];
164  multiplier_element = &multiplier->element[0];
165  result_element = &result->element[i];
166  carry_element = 0;
167  for ( j = 0 ; j < multiplier_size ; j++ ) {
168  bigint_multiply_one ( multiplicand_element,
169  *(multiplier_element++),
170  result_element++,
171  &carry_element );
172  }
173  *result_element = carry_element;
174  }
175 }
176 
177 /**
178  * Reduce big integer R^2 modulo N
179  *
180  * @v modulus0 Element 0 of big integer modulus
181  * @v result0 Element 0 of big integer to hold result
182  * @v size Number of elements in modulus and result
183  *
184  * Reduce the value R^2 modulo N, where R=2^n and n is the number of
185  * bits in the representation of the modulus N, including any leading
186  * zero bits.
187  */
188 void bigint_reduce_raw ( const bigint_element_t *modulus0,
189  bigint_element_t *result0, unsigned int size ) {
190  const bigint_t ( size ) __attribute__ (( may_alias ))
191  *modulus = ( ( const void * ) modulus0 );
192  bigint_t ( size ) __attribute__ (( may_alias ))
193  *result = ( ( void * ) result0 );
194  const unsigned int width = ( 8 * sizeof ( bigint_element_t ) );
195  unsigned int shift;
196  int max;
197  int sign;
198  int msb;
199  int carry;
200 
201  /* We have the constants:
202  *
203  * N = modulus
204  *
205  * n = number of bits in the modulus (including any leading zeros)
206  *
207  * R = 2^n
208  *
209  * Let r be the extension of the n-bit result register by a
210  * separate two's complement sign bit, such that -R <= r < R,
211  * and define:
212  *
213  * x = r * 2^k
214  *
215  * as the value being reduced modulo N, where k is a
216  * non-negative integer bit shift.
217  *
218  * We want to reduce the initial value R^2=2^(2n), which we
219  * may trivially represent using r=1 and k=2n.
220  *
221  * We then iterate over decrementing k, maintaining the loop
222  * invariant:
223  *
224  * -N <= r < N
225  *
226  * On each iteration we must first double r, to compensate for
227  * having decremented k:
228  *
229  * k' = k - 1
230  *
231  * r' = 2r
232  *
233  * x = r * 2^k = 2r * 2^(k-1) = r' * 2^k'
234  *
235  * Note that doubling the n-bit result register will create a
236  * value of n+1 bits: this extra bit needs to be handled
237  * separately during the calculation.
238  *
239  * We then subtract N (if r is currently non-negative) or add
240  * N (if r is currently negative) to restore the loop
241  * invariant:
242  *
243  * 0 <= r < N => r" = 2r - N => -N <= r" < N
244  * -N <= r < 0 => r" = 2r + N => -N <= r" < N
245  *
246  * Note that since N may use all n bits, the most significant
247  * bit of the n-bit result register is not a valid two's
248  * complement sign bit for r: the extra sign bit therefore
249  * also needs to be handled separately.
250  *
251  * Once we reach k=0, we have x=r and therefore:
252  *
253  * -N <= x < N
254  *
255  * After this last loop iteration (with k=0), we may need to
256  * add a single multiple of N to ensure that x is positive,
257  * i.e. lies within the range 0 <= x < N.
258  *
259  * Since neither the modulus nor the value R^2 are secret, we
260  * may elide approximately half of the total number of
261  * iterations by constructing the initial representation of
262  * R^2 as r=2^m and k=2n-m (for some m such that 2^m < N).
263  */
264 
265  /* Initialise x=R^2 */
266  memset ( result, 0, sizeof ( *result ) );
267  max = ( bigint_max_set_bit ( modulus ) - 2 );
268  if ( max < 0 ) {
269  /* Degenerate case of N=0 or N=1: return a zero result */
270  return;
271  }
273  shift = ( ( 2 * size * width ) - max );
274  sign = 0;
275 
276  /* Iterate as described above */
277  while ( shift-- ) {
278 
279  /* Calculate 2r, storing extra bit separately */
280  msb = bigint_shl ( result );
281 
282  /* Add or subtract N according to current sign */
283  if ( sign ) {
284  carry = bigint_add ( modulus, result );
285  } else {
286  carry = bigint_subtract ( modulus, result );
287  }
288 
289  /* Calculate new sign of result
290  *
291  * We know the result lies in the range -N <= r < N
292  * and so the tuple (old sign, msb, carry) cannot ever
293  * take the values (0, 1, 0) or (1, 0, 0). We can
294  * therefore treat these as don't-care inputs, which
295  * allows us to simplify the boolean expression by
296  * ignoring the old sign completely.
297  */
298  assert ( ( sign == msb ) || carry );
299  sign = ( msb ^ carry );
300  }
301 
302  /* Add N to make result positive if necessary */
303  if ( sign )
304  bigint_add ( modulus, result );
305 
306  /* Sanity check */
307  assert ( ! bigint_is_geq ( result, modulus ) );
308 }
309 
310 /**
311  * Compute inverse of odd big integer modulo any power of two
312  *
313  * @v invertend0 Element 0 of odd big integer to be inverted
314  * @v inverse0 Element 0 of big integer to hold result
315  * @v size Number of elements in invertend and result
316  */
317 void bigint_mod_invert_raw ( const bigint_element_t *invertend0,
318  bigint_element_t *inverse0, unsigned int size ) {
319  const bigint_t ( size ) __attribute__ (( may_alias ))
320  *invertend = ( ( const void * ) invertend0 );
321  bigint_t ( size ) __attribute__ (( may_alias ))
322  *inverse = ( ( void * ) inverse0 );
323  bigint_element_t accum;
325  unsigned int i;
326 
327  /* Sanity check */
328  assert ( bigint_bit_is_set ( invertend, 0 ) );
329 
330  /* Initialise output */
331  memset ( inverse, 0xff, sizeof ( *inverse ) );
332 
333  /* Compute inverse modulo 2^(width)
334  *
335  * This method is a lightly modified version of the pseudocode
336  * presented in "A New Algorithm for Inversion mod p^k (Koç,
337  * 2017)".
338  *
339  * Each inner loop iteration calculates one bit of the
340  * inverse. The residue value is the two's complement
341  * negation of the value "b" as used by Koç, to allow for
342  * division by two using a logical right shift (since we have
343  * no arithmetic right shift operation for big integers).
344  *
345  * The residue is stored in the as-yet uncalculated portion of
346  * the inverse. The size of the residue therefore decreases
347  * by one element for each outer loop iteration. Trivial
348  * inspection of the algorithm shows that any higher bits
349  * could not contribute to the eventual output value, and so
350  * we may safely reuse storage this way.
351  *
352  * Due to the suffix property of inverses mod 2^k, the result
353  * represents the least significant bits of the inverse modulo
354  * an arbitrarily large 2^k.
355  */
356  for ( i = size ; i > 0 ; i-- ) {
357  const bigint_t ( i ) __attribute__ (( may_alias ))
358  *addend = ( ( const void * ) invertend );
359  bigint_t ( i ) __attribute__ (( may_alias ))
360  *residue = ( ( void * ) inverse );
361 
362  /* Calculate one element's worth of inverse bits */
363  for ( accum = 0, bit = 1 ; bit ; bit <<= 1 ) {
364  if ( bigint_bit_is_set ( residue, 0 ) ) {
365  accum |= bit;
366  bigint_add ( addend, residue );
367  }
368  bigint_shr ( residue );
369  }
370 
371  /* Store in the element no longer required to hold residue */
372  inverse->element[ i - 1 ] = accum;
373  }
374 
375  /* Correct order of inverse elements */
376  for ( i = 0 ; i < ( size / 2 ) ; i++ ) {
377  accum = inverse->element[i];
378  inverse->element[i] = inverse->element[ size - 1 - i ];
379  inverse->element[ size - 1 - i ] = accum;
380  }
381 }
382 
383 /**
384  * Perform relaxed Montgomery reduction (REDC) of a big integer
385  *
386  * @v modulus0 Element 0 of big integer odd modulus
387  * @v value0 Element 0 of big integer to be reduced
388  * @v result0 Element 0 of big integer to hold result
389  * @v size Number of elements in modulus and result
390  * @ret carry Carry out
391  *
392  * The value to be reduced will be made divisible by the size of the
393  * modulus while retaining its residue class (i.e. multiples of the
394  * modulus will be added until the low half of the value is zero).
395  *
396  * The result may be expressed as
397  *
398  * tR = x + mN
399  *
400  * where x is the input value, N is the modulus, R=2^n (where n is the
401  * number of bits in the representation of the modulus, including any
402  * leading zero bits), and m is the number of multiples of the modulus
403  * added to make the result tR divisible by R.
404  *
405  * The maximum addend is mN <= (R-1)*N (and such an m can be proven to
406  * exist since N is limited to being odd and therefore coprime to R).
407  *
408  * Since the result of this addition is one bit larger than the input
409  * value, a carry out bit is also returned. The caller may be able to
410  * prove that the carry out is always zero, in which case it may be
411  * safely ignored.
412  *
413  * The upper half of the output value (i.e. t) will also be copied to
414  * the result pointer. It is permissible for the result pointer to
415  * overlap the lower half of the input value.
416  *
417  * External knowledge of constraints on the modulus and the input
418  * value may be used to prove constraints on the result. The
419  * constraint on the modulus may be generally expressed as
420  *
421  * R > kN
422  *
423  * for some positive integer k. The value k=1 is allowed, and simply
424  * expresses that the modulus fits within the number of bits in its
425  * own representation.
426  *
427  * For classic Montgomery reduction, we have k=1, i.e. R > N and a
428  * separate constraint that the input value is in the range x < RN.
429  * This gives the result constraint
430  *
431  * tR < RN + (R-1)N
432  * < 2RN - N
433  * < 2RN
434  * t < 2N
435  *
436  * A single subtraction of the modulus may therefore be required to
437  * bring it into the range t < N.
438  *
439  * When the input value is known to be a product of two integers A and
440  * B, with A < aN and B < bN, we get the result constraint
441  *
442  * tR < abN^2 + (R-1)N
443  * < (ab/k)RN + RN - N
444  * < (1 + ab/k)RN
445  * t < (1 + ab/k)N
446  *
447  * If we have k=a=b=1, i.e. R > N with A < N and B < N, then the
448  * result is in the range t < 2N and may require a single subtraction
449  * of the modulus to bring it into the range t < N so that it may be
450  * used as an input on a subsequent iteration.
451  *
452  * If we have k=4 and a=b=2, i.e. R > 4N with A < 2N and B < 2N, then
453  * the result is in the range t < 2N and may immediately be used as an
454  * input on a subsequent iteration, without requiring a subtraction.
455  *
456  * Larger values of k may be used to allow for larger values of a and
457  * b, which can be useful to elide intermediate reductions in a
458  * calculation chain that involves additions and subtractions between
459  * multiplications (as used in elliptic curve point addition, for
460  * example). As a general rule: each intermediate addition or
461  * subtraction will require k to be doubled.
462  *
463  * When the input value is known to be a single integer A, with A < aN
464  * (as used when converting out of Montgomery form), we get the result
465  * constraint
466  *
467  * tR < aN + (R-1)N
468  * < RN + (a-1)N
469  *
470  * If we have a=1, i.e. A < N, then the constraint becomes
471  *
472  * tR < RN
473  * t < N
474  *
475  * and so the result is immediately in the range t < N with no
476  * subtraction of the modulus required.
477  *
478  * For any larger value of a, the result value t=N becomes possible.
479  * Additional external knowledge may potentially be used to prove that
480  * t=N cannot occur. For example: if the caller is performing modular
481  * exponentiation with a prime modulus (or, more generally, a modulus
482  * that is coprime to the base), then there is no way for a non-zero
483  * base value to end up producing an exact multiple of the modulus.
484  * If t=N cannot be disproved, then conversion out of Montgomery form
485  * may require an additional subtraction of the modulus.
486  */
489  bigint_element_t *result0,
490  unsigned int size ) {
491  const bigint_t ( size ) __attribute__ (( may_alias ))
492  *modulus = ( ( const void * ) modulus0 );
493  union {
494  bigint_t ( size * 2 ) full;
495  struct {
496  bigint_t ( size ) low;
497  bigint_t ( size ) high;
498  } __attribute__ (( packed ));
499  } __attribute__ (( may_alias )) *value = ( ( void * ) value0 );
500  bigint_t ( size ) __attribute__ (( may_alias ))
501  *result = ( ( void * ) result0 );
502  static bigint_t ( 1 ) cached;
503  static bigint_t ( 1 ) negmodinv;
504  bigint_element_t multiple;
506  unsigned int i;
507  unsigned int j;
508  int overflow;
509 
510  /* Sanity checks */
511  assert ( bigint_bit_is_set ( modulus, 0 ) );
512 
513  /* Calculate inverse (or use cached version) */
514  if ( cached.element[0] != modulus->element[0] ) {
515  bigint_mod_invert ( modulus, &negmodinv );
516  negmodinv.element[0] = -negmodinv.element[0];
517  cached.element[0] = modulus->element[0];
518  }
519 
520  /* Perform multiprecision Montgomery reduction */
521  for ( i = 0 ; i < size ; i++ ) {
522 
523  /* Determine scalar multiple for this round */
524  multiple = ( value->low.element[i] * negmodinv.element[0] );
525 
526  /* Multiply value to make it divisible by 2^(width*i) */
527  carry = 0;
528  for ( j = 0 ; j < size ; j++ ) {
529  bigint_multiply_one ( multiple, modulus->element[j],
530  &value->full.element[ i + j ],
531  &carry );
532  }
533 
534  /* Since value is now divisible by 2^(width*i), we
535  * know that the current low element must have been
536  * zeroed.
537  */
538  assert ( value->low.element[i] == 0 );
539 
540  /* Store the multiplication carry out in the result,
541  * avoiding the need to immediately propagate the
542  * carry through the remaining elements.
543  */
544  result->element[i] = carry;
545  }
546 
547  /* Add the accumulated carries */
548  overflow = bigint_add ( result, &value->high );
549 
550  /* Copy to result buffer */
551  bigint_copy ( &value->high, result );
552 
553  return overflow;
554 }
555 
556 /**
557  * Perform classic Montgomery reduction (REDC) of a big integer
558  *
559  * @v modulus0 Element 0 of big integer odd modulus
560  * @v value0 Element 0 of big integer to be reduced
561  * @v result0 Element 0 of big integer to hold result
562  * @v size Number of elements in modulus and result
563  */
566  bigint_element_t *result0,
567  unsigned int size ) {
568  const bigint_t ( size ) __attribute__ (( may_alias ))
569  *modulus = ( ( const void * ) modulus0 );
570  union {
571  bigint_t ( size * 2 ) full;
572  struct {
573  bigint_t ( size ) low;
574  bigint_t ( size ) high;
575  } __attribute__ (( packed ));
576  } __attribute__ (( may_alias )) *value = ( ( void * ) value0 );
577  bigint_t ( size ) __attribute__ (( may_alias ))
578  *result = ( ( void * ) result0 );
579  int overflow;
580  int underflow;
581 
582  /* Sanity check */
583  assert ( ! bigint_is_geq ( &value->high, modulus ) );
584 
585  /* Perform relaxed Montgomery reduction */
586  overflow = bigint_montgomery_relaxed ( modulus, &value->full, result );
587 
588  /* Conditionally subtract the modulus once */
589  underflow = bigint_subtract ( modulus, result );
590  bigint_swap ( result, &value->high, ( underflow & ~overflow ) );
591 
592  /* Sanity check */
593  assert ( ! bigint_is_geq ( result, modulus ) );
594 }
595 
596 /**
597  * Perform generalised exponentiation via a Montgomery ladder
598  *
599  * @v result0 Element 0 of result (initialised to identity element)
600  * @v multiple0 Element 0 of multiple (initialised to generator)
601  * @v size Number of elements in result and multiple
602  * @v exponent0 Element 0 of exponent
603  * @v exponent_size Number of elements in exponent
604  * @v op Montgomery ladder commutative operation
605  * @v ctx Operation context (if needed)
606  * @v tmp Temporary working space (if needed)
607  *
608  * The Montgomery ladder may be used to perform any operation that is
609  * isomorphic to exponentiation, i.e. to compute the result
610  *
611  * r = g^e = g * g * g * g * .... * g
612  *
613  * for an arbitrary commutative operation "*", generator "g" and
614  * exponent "e".
615  *
616  * The result "r" is computed in constant time (assuming that the
617  * underlying operation is constant time) in k steps, where k is the
618  * number of bits in the big integer representation of the exponent.
619  *
620  * The result "r" must be initialised to the operation's identity
621  * element, and the multiple must be initialised to the generator "g".
622  * On exit, the multiple will contain
623  *
624  * m = r * g = g^(e+1)
625  *
626  * Note that the terminology used here refers to exponentiation
627  * defined as repeated multiplication, but that the ladder may equally
628  * well be used to perform any isomorphic operation (such as
629  * multiplication defined as repeated addition).
630  */
632  bigint_element_t *multiple0, unsigned int size,
633  const bigint_element_t *exponent0,
634  unsigned int exponent_size, bigint_ladder_op_t *op,
635  const void *ctx, void *tmp ) {
636  bigint_t ( size ) __attribute__ (( may_alias ))
637  *result = ( ( void * ) result0 );
638  bigint_t ( size ) __attribute__ (( may_alias ))
639  *multiple = ( ( void * ) multiple0 );
640  const bigint_t ( exponent_size ) __attribute__ (( may_alias ))
641  *exponent = ( ( const void * ) exponent0 );
642  const unsigned int width = ( 8 * sizeof ( bigint_element_t ) );
643  unsigned int bit = ( exponent_size * width );
644  int toggle = 0;
645  int previous;
646 
647  /* We have two physical storage locations: Rres (the "result"
648  * register) and Rmul (the "multiple" register).
649  *
650  * On entry we have:
651  *
652  * Rres = g^0 = 1 (the identity element)
653  * Rmul = g (the generator)
654  *
655  * For calculation purposes, define two virtual registers R[0]
656  * and R[1], mapped to the underlying physical storage
657  * locations via a boolean toggle state "t":
658  *
659  * R[t] -> Rres
660  * R[~t] -> Rmul
661  *
662  * Define the initial toggle state to be t=0, so that we have:
663  *
664  * R[0] = g^0 = 1 (the identity element)
665  * R[1] = g (the generator)
666  *
667  * The Montgomery ladder then iterates over each bit "b" of
668  * the exponent "e" (from MSB down to LSB), computing:
669  *
670  * R[1] := R[0] * R[1], R[0] := R[0] * R[0] (if b=0)
671  * R[0] := R[1] * R[0], R[1] := R[1] * R[1] (if b=1)
672  *
673  * i.e.
674  *
675  * R[~b] := R[b] * R[~b], R[b] := R[b] * R[b]
676  *
677  * To avoid data-dependent memory access patterns, we
678  * implement this as:
679  *
680  * Rmul := Rres * Rmul
681  * Rres := Rres * Rres
682  *
683  * and update the toggle state "t" (via constant-time
684  * conditional swaps) as needed to ensure that R[b] always
685  * maps to Rres and R[~b] always maps to Rmul.
686  */
687  while ( bit-- ) {
688 
689  /* Update toggle state t := b
690  *
691  * If the new toggle state is not equal to the old
692  * toggle state, then we must swap R[0] and R[1] (in
693  * constant time).
694  */
695  previous = toggle;
696  toggle = bigint_bit_is_set ( exponent, bit );
697  bigint_swap ( result, multiple, ( toggle ^ previous ) );
698 
699  /* Calculate Rmul = R[~b] := R[b] * R[~b] = Rres * Rmul */
700  op ( result->element, multiple->element, size, ctx, tmp );
701 
702  /* Calculate Rres = R[b] := R[b] * R[b] = Rres * Rres */
703  op ( result->element, result->element, size, ctx, tmp );
704  }
705 
706  /* Reset toggle state to zero */
707  bigint_swap ( result, multiple, toggle );
708 }
709 
710 /**
711  * Perform modular multiplication as part of a Montgomery ladder
712  *
713  * @v operand Element 0 of first input operand (may overlap result)
714  * @v result Element 0 of second input operand and result
715  * @v size Number of elements in operands and result
716  * @v ctx Operation context (odd modulus, or NULL)
717  * @v tmp Temporary working space
718  */
719 void bigint_mod_exp_ladder ( const bigint_element_t *multiplier0,
720  bigint_element_t *result0, unsigned int size,
721  const void *ctx, void *tmp ) {
722  const bigint_t ( size ) __attribute__ (( may_alias ))
723  *multiplier = ( ( const void * ) multiplier0 );
724  bigint_t ( size ) __attribute__ (( may_alias ))
725  *result = ( ( void * ) result0 );
726  const bigint_t ( size ) __attribute__ (( may_alias ))
727  *modulus = ( ( const void * ) ctx );
728  bigint_t ( size * 2 ) __attribute__ (( may_alias ))
729  *product = ( ( void * ) tmp );
730 
731  /* Multiply and reduce */
733  if ( modulus ) {
734  bigint_montgomery ( modulus, product, result );
735  } else {
737  }
738 }
739 
740 /**
741  * Perform modular exponentiation of big integers
742  *
743  * @v base0 Element 0 of big integer base
744  * @v modulus0 Element 0 of big integer modulus
745  * @v exponent0 Element 0 of big integer exponent
746  * @v result0 Element 0 of big integer to hold result
747  * @v size Number of elements in base, modulus, and result
748  * @v exponent_size Number of elements in exponent
749  * @v tmp Temporary working space
750  */
752  const bigint_element_t *modulus0,
753  const bigint_element_t *exponent0,
754  bigint_element_t *result0,
755  unsigned int size, unsigned int exponent_size,
756  void *tmp ) {
757  const bigint_t ( size ) __attribute__ (( may_alias )) *base =
758  ( ( const void * ) base0 );
759  const bigint_t ( size ) __attribute__ (( may_alias )) *modulus =
760  ( ( const void * ) modulus0 );
761  const bigint_t ( exponent_size ) __attribute__ (( may_alias ))
762  *exponent = ( ( const void * ) exponent0 );
763  bigint_t ( size ) __attribute__ (( may_alias )) *result =
764  ( ( void * ) result0 );
765  const unsigned int width = ( 8 * sizeof ( bigint_element_t ) );
766  struct {
767  struct {
768  bigint_t ( size ) modulus;
769  bigint_t ( size ) stash;
770  };
771  union {
772  bigint_t ( 2 * size ) full;
773  bigint_t ( size ) low;
774  } product;
775  } *temp = tmp;
776  const uint8_t one[1] = { 1 };
777  bigint_element_t submask;
778  unsigned int subsize;
779  unsigned int scale;
780  unsigned int i;
781 
782  /* Sanity check */
783  assert ( sizeof ( *temp ) == bigint_mod_exp_tmp_len ( modulus ) );
784 
785  /* Handle degenerate case of zero modulus */
786  if ( ! bigint_max_set_bit ( modulus ) ) {
787  memset ( result, 0, sizeof ( *result ) );
788  return;
789  }
790 
791  /* Factor modulus as (N * 2^scale) where N is odd */
792  bigint_copy ( modulus, &temp->modulus );
793  for ( scale = 0 ; ( ! bigint_bit_is_set ( &temp->modulus, 0 ) ) ;
794  scale++ ) {
795  bigint_shr ( &temp->modulus );
796  }
797  subsize = ( ( scale + width - 1 ) / width );
798  submask = ( ( 1UL << ( scale % width ) ) - 1 );
799  if ( ! submask )
800  submask = ~submask;
801 
802  /* Calculate (R^2 mod N) */
803  bigint_reduce ( &temp->modulus, &temp->stash );
804 
805  /* Initialise result = Montgomery(1, R^2 mod N) */
806  bigint_grow ( &temp->stash, &temp->product.full );
807  bigint_montgomery ( &temp->modulus, &temp->product.full, result );
808 
809  /* Convert base into Montgomery form */
810  bigint_multiply ( base, &temp->stash, &temp->product.full );
811  bigint_montgomery ( &temp->modulus, &temp->product.full,
812  &temp->stash );
813 
814  /* Calculate x1 = base^exponent modulo N */
815  bigint_ladder ( result, &temp->stash, exponent, bigint_mod_exp_ladder,
816  &temp->modulus, &temp->product );
817 
818  /* Convert back out of Montgomery form */
819  bigint_grow ( result, &temp->product.full );
820  bigint_montgomery_relaxed ( &temp->modulus, &temp->product.full,
821  result );
822 
823  /* Handle even moduli via Garner's algorithm */
824  if ( subsize ) {
825  const bigint_t ( subsize ) __attribute__ (( may_alias ))
826  *subbase = ( ( const void * ) base );
827  bigint_t ( subsize ) __attribute__ (( may_alias ))
828  *submodulus = ( ( void * ) &temp->modulus );
829  bigint_t ( subsize ) __attribute__ (( may_alias ))
830  *substash = ( ( void * ) &temp->stash );
831  bigint_t ( subsize ) __attribute__ (( may_alias ))
832  *subresult = ( ( void * ) result );
833  union {
834  bigint_t ( 2 * subsize ) full;
835  bigint_t ( subsize ) low;
836  } __attribute__ (( may_alias ))
837  *subproduct = ( ( void * ) &temp->product.full );
838 
839  /* Calculate x2 = base^exponent modulo 2^k */
840  bigint_init ( substash, one, sizeof ( one ) );
841  bigint_copy ( subbase, submodulus );
842  bigint_ladder ( substash, submodulus, exponent,
843  bigint_mod_exp_ladder, NULL, subproduct );
844 
845  /* Reconstruct N */
846  bigint_copy ( modulus, &temp->modulus );
847  for ( i = 0 ; i < scale ; i++ )
848  bigint_shr ( &temp->modulus );
849 
850  /* Calculate N^-1 modulo 2^k */
851  bigint_mod_invert ( submodulus, &subproduct->low );
852  bigint_copy ( &subproduct->low, submodulus );
853 
854  /* Calculate y = (x2 - x1) * N^-1 modulo 2^k */
855  bigint_subtract ( subresult, substash );
856  bigint_multiply ( substash, submodulus, &subproduct->full );
857  subproduct->low.element[ subsize - 1 ] &= submask;
858  bigint_grow ( &subproduct->low, &temp->stash );
859 
860  /* Reconstruct N */
861  bigint_mod_invert ( submodulus, &subproduct->low );
862  bigint_copy ( &subproduct->low, submodulus );
863 
864  /* Calculate x = x1 + N * y */
865  bigint_multiply ( &temp->modulus, &temp->stash,
866  &temp->product.full );
867  bigint_add ( &temp->product.low, result );
868  }
869 }
#define __attribute__(x)
Definition: compiler.h:10
uint32_t base
Base.
Definition: librm.h:252
void bigint_multiply_raw(const bigint_element_t *multiplicand0, unsigned int multiplicand_size, const bigint_element_t *multiplier0, unsigned int multiplier_size, bigint_element_t *result0)
Multiply big integers.
Definition: bigint.c:117
uint32_t low
Low 16 bits of address.
Definition: myson.h:19
void bigint_mod_exp_raw(const bigint_element_t *base0, const bigint_element_t *modulus0, const bigint_element_t *exponent0, bigint_element_t *result0, unsigned int size, unsigned int exponent_size, void *tmp)
Perform modular exponentiation of big integers.
Definition: bigint.c:751
#define max(x, y)
Definition: ath.h:39
#define bigint_max_set_bit(value)
Find highest bit set in big integer.
Definition: bigint.h:198
int carry
Definition: bigint.h:65
void bigint_ladder_raw(bigint_element_t *result0, bigint_element_t *multiple0, unsigned int size, const bigint_element_t *exponent0, unsigned int exponent_size, bigint_ladder_op_t *op, const void *ctx, void *tmp)
Perform generalised exponentiation via a Montgomery ladder.
Definition: bigint.c:631
#define sprintf(buf, fmt,...)
Write a formatted string to a buffer.
Definition: stdio.h:36
uint8_t size
Entry size (in 32-bit words)
Definition: ena.h:16
#define bigint_grow(source, dest)
Grow big integer.
Definition: bigint.h:208
#define bigint_init(value, data, len)
Initialise big integer.
Definition: bigint.h:61
static unsigned int unsigned int bit
Definition: bigint.h:391
void bigint_multiply_one(const bigint_element_t multiplicand, const bigint_element_t multiplier, bigint_element_t *result, bigint_element_t *carry)
struct golan_eq_context ctx
Definition: CIB_PRM.h:28
#define bigint_montgomery_relaxed(modulus, value, result)
Perform relaxed Montgomery reduction (REDC) of a big integer.
Definition: bigint.h:299
value element[index]
Definition: bigint.h:397
static uint32_t * value0
Definition: bigint.h:58
Big integer support.
static u32 xor(u32 a, u32 b)
Definition: tlan.h:457
int bigint_montgomery_relaxed_raw(const bigint_element_t *modulus0, bigint_element_t *value0, bigint_element_t *result0, unsigned int size)
Perform relaxed Montgomery reduction (REDC) of a big integer.
Definition: bigint.c:487
unsigned long tmp
Definition: linux_pci.h:63
#define bigint_is_geq(value, reference)
Compare big integers.
Definition: bigint.h:144
#define bigint_mod_invert(invertend, inverse)
Compute inverse of odd big integer modulo any power of two.
Definition: bigint.h:285
Assertions.
assert((readw(&hdr->flags) &(GTF_reading|GTF_writing))==0)
#define bigint_shr(value)
Shift big integer right.
Definition: bigint.h:121
#define bigint_shrink(source, dest)
Shrink big integer.
Definition: bigint.h:221
pseudo_bit_t value[0x00020]
Definition: arbel.h:13
#define bigint_copy(source, dest)
Copy big integer.
Definition: bigint.h:234
uint32_t bigint_element_t
Element of a big integer.
Definition: bigint.h:15
uint16_t count
Number of entries.
Definition: ena.h:22
uint32_t high
High 32 bits of address.
Definition: myson.h:20
unsigned char uint8_t
Definition: stdint.h:10
#define bigint_swap(first, second, swap)
Conditionally swap big integers (in constant time)
Definition: bigint.h:246
unsigned char byte
Definition: smc9000.h:38
void bigint_montgomery_raw(const bigint_element_t *modulus0, bigint_element_t *value0, bigint_element_t *result0, unsigned int size)
Perform classic Montgomery reduction (REDC) of a big integer.
Definition: bigint.c:564
uint16_t result
Definition: hyperv.h:33
#define bigint_shl(value)
Shift big integer left.
Definition: bigint.h:110
FILE_LICENCE(GPL2_OR_LATER_OR_UBDL)
static uint16_t struct vmbus_xfer_pages_operations * op
Definition: netvsc.h:327
#define bigint_montgomery(modulus, value, result)
Perform classic Montgomery reduction (REDC) of a big integer.
Definition: bigint.h:313
static const uint32_t multiplier
Port multiplier number.
Definition: bigint.h:330
void bigint_swap_raw(bigint_element_t *first0, bigint_element_t *second0, unsigned int size, int swap)
Conditionally swap big integers (in constant time)
Definition: bigint.c:91
#define bigint_multiply(multiplicand, multiplier, result)
Multiply big integers.
Definition: bigint.h:259
#define bigint_reduce(modulus, result)
Reduce big integer R^2 modulo N.
Definition: bigint.h:273
uint8_t product
Product string.
Definition: smbios.h:16
#define bigint_ladder(result, multiple, exponent, op, ctx, tmp)
Perform generalised exponentiation via a Montgomery ladder.
Definition: bigint.h:329
void() bigint_ladder_op_t(const bigint_element_t *operand0, bigint_element_t *result0, unsigned int size, const void *ctx, void *tmp)
A big integer Montgomery ladder commutative operation.
Definition: bigint.h:377
const char * bigint_ntoa_raw(const bigint_element_t *value0, unsigned int size)
Transcribe big integer (for debugging)
Definition: bigint.c:47
#define bigint_bit_is_set(value, bit)
Test if bit is set in big integer.
Definition: bigint.h:178
#define bigint_subtract(subtrahend, value)
Subtract big integers.
Definition: bigint.h:98
void bigint_mod_invert_raw(const bigint_element_t *invertend0, bigint_element_t *inverse0, unsigned int size)
Compute inverse of odd big integer modulo any power of two.
Definition: bigint.c:317
#define bigint_mod_exp_tmp_len(modulus)
Calculate temporary working space required for moduluar exponentiation.
Definition: bigint.h:360
#define bigint_add(addend, value)
Add big integers.
Definition: bigint.h:86
#define NULL
NULL pointer (VOID *)
Definition: Base.h:321
String functions.
#define bigint_set_bit(value, bit)
Set bit in big integer.
Definition: bigint.h:155
void bigint_reduce_raw(const bigint_element_t *modulus0, bigint_element_t *result0, unsigned int size)
Reduce big integer R^2 modulo N.
Definition: bigint.c:188
typedef bigint_t(X25519_SIZE) x25519_t
An X25519 unsigned big integer used in internal calculations.
#define BIGINT_NTOA_LSB_MIN
Minimum number of least significant bytes included in transcription.
Definition: bigint.c:38
void * memset(void *dest, int character, size_t len) __nonnull
void bigint_mod_exp_ladder(const bigint_element_t *multiplier0, bigint_element_t *result0, unsigned int size, const void *ctx, void *tmp)
Perform modular multiplication as part of a Montgomery ladder.
Definition: bigint.c:719