chiark - git - mdw - catacomb/blob - 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