282
|
1 #include <tommath.h> |
|
2 #ifdef BN_MP_MONTGOMERY_REDUCE_C |
|
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis |
|
4 * |
|
5 * LibTomMath is a library that provides multiple-precision |
|
6 * integer arithmetic as well as number theoretic functionality. |
|
7 * |
|
8 * The library was designed directly after the MPI library by |
|
9 * Michael Fromberger but has been written from scratch with |
|
10 * additional optimizations in place. |
|
11 * |
|
12 * The library is free for all purposes without any express |
|
13 * guarantee it works. |
|
14 * |
|
15 * Tom St Denis, [email protected], http://math.libtomcrypt.org |
|
16 */ |
|
17 |
|
18 /* computes xR**-1 == x (mod N) via Montgomery Reduction */ |
|
19 int |
|
20 mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho) |
|
21 { |
|
22 int ix, res, digs; |
|
23 mp_digit mu; |
|
24 |
|
25 /* can the fast reduction [comba] method be used? |
|
26 * |
|
27 * Note that unlike in mul you're safely allowed *less* |
|
28 * than the available columns [255 per default] since carries |
|
29 * are fixed up in the inner loop. |
|
30 */ |
|
31 digs = n->used * 2 + 1; |
|
32 if ((digs < MP_WARRAY) && |
|
33 n->used < |
|
34 (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { |
|
35 return fast_mp_montgomery_reduce (x, n, rho); |
|
36 } |
|
37 |
|
38 /* grow the input as required */ |
|
39 if (x->alloc < digs) { |
|
40 if ((res = mp_grow (x, digs)) != MP_OKAY) { |
|
41 return res; |
|
42 } |
|
43 } |
|
44 x->used = digs; |
|
45 |
|
46 for (ix = 0; ix < n->used; ix++) { |
|
47 /* mu = ai * rho mod b |
|
48 * |
|
49 * The value of rho must be precalculated via |
|
50 * montgomery_setup() such that |
|
51 * it equals -1/n0 mod b this allows the |
|
52 * following inner loop to reduce the |
|
53 * input one digit at a time |
|
54 */ |
|
55 mu = (mp_digit) (((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK); |
|
56 |
|
57 /* a = a + mu * m * b**i */ |
|
58 { |
|
59 register int iy; |
|
60 register mp_digit *tmpn, *tmpx, u; |
|
61 register mp_word r; |
|
62 |
|
63 /* alias for digits of the modulus */ |
|
64 tmpn = n->dp; |
|
65 |
|
66 /* alias for the digits of x [the input] */ |
|
67 tmpx = x->dp + ix; |
|
68 |
|
69 /* set the carry to zero */ |
|
70 u = 0; |
|
71 |
|
72 /* Multiply and add in place */ |
|
73 for (iy = 0; iy < n->used; iy++) { |
|
74 /* compute product and sum */ |
|
75 r = ((mp_word)mu) * ((mp_word)*tmpn++) + |
|
76 ((mp_word) u) + ((mp_word) * tmpx); |
|
77 |
|
78 /* get carry */ |
|
79 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT)); |
|
80 |
|
81 /* fix digit */ |
|
82 *tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK)); |
|
83 } |
|
84 /* At this point the ix'th digit of x should be zero */ |
|
85 |
|
86 |
|
87 /* propagate carries upwards as required*/ |
|
88 while (u) { |
|
89 *tmpx += u; |
|
90 u = *tmpx >> DIGIT_BIT; |
|
91 *tmpx++ &= MP_MASK; |
|
92 } |
|
93 } |
|
94 } |
|
95 |
|
96 /* at this point the n.used'th least |
|
97 * significant digits of x are all zero |
|
98 * which means we can shift x to the |
|
99 * right by n.used digits and the |
|
100 * residue is unchanged. |
|
101 */ |
|
102 |
|
103 /* x = x/b**n.used */ |
|
104 mp_clamp(x); |
|
105 mp_rshd (x, n->used); |
|
106 |
|
107 /* if x >= n then x = x - n */ |
|
108 if (mp_cmp_mag (x, n) != MP_LT) { |
|
109 return s_mp_sub (x, n, x); |
|
110 } |
|
111 |
|
112 return MP_OKAY; |
|
113 } |
|
114 #endif |