chiark / gitweb /
math/mpmont.c: Factor out the computational core of the algorithm.
[catacomb] / math / mpmont.h
1 /* -*-c-*-
2  *
3  * Montgomery reduction
4  *
5  * (c) 1999 Straylight/Edgeware
6  */
7
8 /*----- Licensing notice --------------------------------------------------*
9  *
10  * This file is part of Catacomb.
11  *
12  * Catacomb is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU Library General Public License as
14  * published by the Free Software Foundation; either version 2 of the
15  * License, or (at your option) any later version.
16  *
17  * Catacomb is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20  * GNU Library General Public License for more details.
21  *
22  * You should have received a copy of the GNU Library General Public
23  * License along with Catacomb; if not, write to the Free
24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25  * MA 02111-1307, USA.
26  */
27
28 #ifndef CATACOMB_MPMONT_H
29 #define CATACOMB_MPMONT_H
30
31 #ifdef __cplusplus
32   extern "C" {
33 #endif
34
35 /*----- Header files ------------------------------------------------------*/
36
37 #ifndef CATACOMB_MP_H
38 #  include "mp.h"
39 #endif
40
41 /*----- Notes on Montgomery reduction -------------------------------------*
42  *
43  * Given a little bit of precomputation, Montgomery reduction enables modular
44  * reductions of products to be calculated rather rapidly, without recourse
45  * to annoying things like division.
46  *
47  * Before starting, you need to do a little work.  In particular, the
48  * following things need to be worked out:
49  *
50  *   * %$m$%, which is the modulus you'll be working with.  This must be odd,
51  *     otherwise the whole thing doesn't work.  You're better off using
52  *     Barrett reduction if your modulus might be even.
53  *
54  *   * %$b$%, the radix of the number system you're in (here, it's
55  *     @MPW_MAX + 1@).
56  *
57  *   * %$m' = -m^{-1} \bmod b$%, a useful number for the reduction step.
58  *     (This means that the modulus mustn't be even.  This shouldn't be a
59  *     problem.)
60  *
61  *   * %$R = b^n > m > b^{n - 1}$%, or at least %$\log_2 R$%.
62  *
63  *   * %$R \bmod m$% and %$R^2 \bmod m$%, which are useful when doing
64  *     calculations such as exponentiation.
65  *
66  * Suppose that %$0 \le a_i \le (b^n + b^i - 1) m$% with %$a_i \equiv {}$%
67  * %$0 \pmod{b^i}$%.  Let %$w_i = m' a_i/b^i \bmod b$%, and set %$a_{i+1} =
68  * a_i + b^i w_i m$%.  Then obviously %$a_{i+1} \equiv {} $% %$a_i
69  * \pmod{m}$%, and less obviously %$a_{i+1}/b^i \equiv a_i/b^i + {}$% %$m m'
70  * a_i/b^i \equiv 0 \pmod{b}$% so %$a_{i+1} \equiv 0 \pmod{b^{i+1}}$%.
71  * Finally, we can bound %$a_{i+1} \le {}$% %$a_i + b^i (b - 1) m = {}$%
72  * %$a_i + (b^{i+1} - b^i) m \le (b^n + b^{i+1} - 1) m$%.  As a result, if
73  * we're given some %a_0%, we can calculate %$a_n \equiv 0 \pmod{R}$%, with
74  * $%a_n \equiv a_0 \pmod{n}$%, i.e., %$a_n/R \equiv a_0 R^{-1} \pmod{m}$%;
75  * furthermore, if %$0 \le a_0 < m + b^n%$ then %$0 \le a_n/R < 2 m$%, so a
76  * fully reduced result can be obtained with a single conditional
77  * subtraction.
78  *
79  * The result of reduing %$a$% is then %$a R^{-1}$% \bmod m$%.  This is
80  * actually rather useful for reducing products, if we run an extra factor of
81  * %$R$% through the calculation: the result of reducing the product of
82  * %$(x R)(y R) = x y R^2$% is then %$x y R \bmod m$%, preserving the running
83  * factor.  Thanks to distributivity, additions and subtractions can be
84  * performed on numbers in this form -- the extra factor of %$R$% just runs
85  * through all the calculations until it's finally stripped out by a final
86  * reduction operation.
87  */
88
89 /*----- Data structures ---------------------------------------------------*/
90
91 /* --- A Montgomery reduction context --- */
92
93 typedef struct mpmont {
94   mp *m;                                /* Modulus */
95   mp *mi;                               /* %$-m^{-1} \bmod R$% */
96   size_t n;                             /* %$\log_b R$% */
97   mp *r, *r2;                           /* %$R \bmod m$%, %$R^2 \bmod m$% */
98 } mpmont;
99
100 /*----- Functions provided ------------------------------------------------*/
101
102 /* --- @mpmont_create@ --- *
103  *
104  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
105  *              @mp *m@ = modulus to use
106  *
107  * Returns:     Zero on success, nonzero on error.
108  *
109  * Use:         Initializes a Montgomery reduction context ready for use.
110  *              The argument @m@ must be a positive odd integer.
111  */
112
113 extern int mpmont_create(mpmont */*mm*/, mp */*m*/);
114
115 /* --- @mpmont_destroy@ --- *
116  *
117  * Arguments:   @mpmont *mm@ = pointer to a Montgomery reduction context
118  *
119  * Returns:     ---
120  *
121  * Use:         Disposes of a context when it's no longer of any use to
122  *              anyone.
123  */
124
125 extern void mpmont_destroy(mpmont */*mm*/);
126
127 /* --- @mpmont_reduce@ --- *
128  *
129  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
130  *              @mp *d@ = destination
131  *              @mp *a@ = source, assumed positive
132  *
133  * Returns:     Result, %$a R^{-1} \bmod m$%.
134  */
135
136 extern mp *mpmont_reduce(mpmont */*mm*/, mp */*d*/, mp */*a*/);
137
138 /* --- @mpmont_mul@ --- *
139  *
140  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
141  *              @mp *d@ = destination
142  *              @mp *a, *b@ = sources, assumed positive
143  *
144  * Returns:     Result, %$a b R^{-1} \bmod m$%.
145  */
146
147 extern mp *mpmont_mul(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*b*/);
148
149 /* --- @mpmont_expr@ --- *
150  *
151  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
152  *              @mp *d@ = fake destination
153  *              @mp *a@ = base
154  *              @mp *e@ = exponent
155  *
156  * Returns:     Result, %$(a R^{-1})^e R \bmod m$%.  This is useful if
157  *              further modular arithmetic is to be performed on the result.
158  */
159
160 extern mp *mpmont_expr(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*e*/);
161
162 /* --- @mpmont_exp@ --- *
163  *
164  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
165  *              @mp *d@ = fake destination
166  *              @mp *a@ = base
167  *              @mp *e@ = exponent
168  *
169  * Returns:     Result, %$a^e \bmod m$%.
170  */
171
172 extern mp *mpmont_exp(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*e*/);
173
174 /* --- @mpmont_mexpr@ --- *
175  *
176  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
177  *              @mp *d@ = fake destination
178  *              @const mp_expfactor *f@ = pointer to array of factors
179  *              @size_t n@ = number of factors supplied
180  *
181  * Returns:     If the bases are %$g_0, g_1, \ldots, g_{n-1}$% and the
182  *              exponents are %$e_0, e_1, \ldots, e_{n-1}$% then the result
183  *              is:
184  *
185  *              %$g_0^{e_0} g_1^{e_1} \ldots g_{n-1}^{e_{n-1}} \bmod m$%
186  *
187  *
188  *              except that the %$g_i$% and result are in Montgomery form.
189  */
190
191 extern mp *mpmont_mexpr(mpmont */*mm*/, mp */*d*/,
192                         const mp_expfactor */*f*/, size_t /*n*/);
193
194 /* --- @mpmont_mexp@ --- *
195  *
196  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
197  *              @mp *d@ = fake destination
198  *              @const mp_expfactor *f@ = pointer to array of factors
199  *              @size_t n@ = number of factors supplied
200  *
201  * Returns:     Product of bases raised to exponents, all mod @m@.
202  *
203  * Use:         Convenient interface over @mpmont_mexpr@.
204  */
205
206 extern mp *mpmont_mexp(mpmont */*mm*/, mp */*d*/,
207                        const mp_expfactor */*f*/, size_t /*n*/);
208
209 /*----- That's all, folks -------------------------------------------------*/
210
211 #ifdef __cplusplus
212   }
213 #endif
214
215 #endif