X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a..HEAD:/math/mpmont.h diff --git a/math/mpmont.h b/math/mpmont.h index 9973ec3f..cd1d3999 100644 --- a/math/mpmont.h +++ b/math/mpmont.h @@ -54,22 +54,36 @@ * * %$b$%, the radix of the number system you're in (here, it's * @MPW_MAX + 1@). * - * * %$-m^{-1} \bmod b$%, a useful number for the reduction step. (This - * means that the modulus mustn't be even. This shouldn't be a problem.) + * * %$m' = -m^{-1} \bmod b$%, a useful number for the reduction step. + * (This means that the modulus mustn't be even. This shouldn't be a + * problem.) * * * %$R = b^n > m > b^{n - 1}$%, or at least %$\log_2 R$%. * * * %$R \bmod m$% and %$R^2 \bmod m$%, which are useful when doing * calculations such as exponentiation. * - * The result of a Montgomery reduction of %$x$% is %$x R^{-1} \bmod m$%, - * which doesn't look ever-so useful. The trick is to initially apply a - * factor of %$R$% to all of your numbers so that when you multiply and - * perform a Montgomery reduction you get %$(x R \cdot y R) R^{-1} \bmod m$%, - * which is just %$x y R \bmod m$%. Thanks to distributivity, even additions - * and subtractions can be performed on numbers in this form -- the extra - * factor of %$R$% just runs through all the calculations until it's finally - * stripped out by a final reduction operation. + * Suppose that %$0 \le a_i \le (b^n + b^i - 1) m$% with %$a_i \equiv {}$% + * %$0 \pmod{b^i}$%. Let %$w_i = m' a_i/b^i \bmod b$%, and set %$a_{i+1} = + * a_i + b^i w_i m$%. Then obviously %$a_{i+1} \equiv {} $% %$a_i + * \pmod{m}$%, and less obviously %$a_{i+1}/b^i \equiv a_i/b^i + {}$% %$m m' + * a_i/b^i \equiv 0 \pmod{b}$% so %$a_{i+1} \equiv 0 \pmod{b^{i+1}}$%. + * Finally, we can bound %$a_{i+1} \le {}$% %$a_i + b^i (b - 1) m = {}$% + * %$a_i + (b^{i+1} - b^i) m \le (b^n + b^{i+1} - 1) m$%. As a result, if + * we're given some %a_0%, we can calculate %$a_n \equiv 0 \pmod{R}$%, with + * $%a_n \equiv a_0 \pmod{n}$%, i.e., %$a_n/R \equiv a_0 R^{-1} \pmod{m}$%; + * furthermore, if %$0 \le a_0 < m + b^n%$ then %$0 \le a_n/R < 2 m$%, so a + * fully reduced result can be obtained with a single conditional + * subtraction. + * + * The result of reduing %$a$% is then %$a R^{-1}$% \bmod m$%. This is + * actually rather useful for reducing products, if we run an extra factor of + * %$R$% through the calculation: the result of reducing the product of + * %$(x R)(y R) = x y R^2$% is then %$x y R \bmod m$%, preserving the running + * factor. Thanks to distributivity, additions and subtractions can be + * performed on numbers in this form -- the extra factor of %$R$% just runs + * through all the calculations until it's finally stripped out by a final + * reduction operation. */ /*----- Data structures ---------------------------------------------------*/ @@ -112,29 +126,29 @@ extern void mpmont_destroy(mpmont */*mm*/); /* --- @mpmont_reduce@ --- * * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = destination * @mp *a@ = source, assumed positive * * Returns: Result, %$a R^{-1} \bmod m$%. */ -extern mp *mpmont_reduce(mpmont */*mm*/, mp */*d*/, mp */*a*/); +extern mp *mpmont_reduce(const mpmont */*mm*/, mp */*d*/, mp */*a*/); /* --- @mpmont_mul@ --- * * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = destination * @mp *a, *b@ = sources, assumed positive * * Returns: Result, %$a b R^{-1} \bmod m$%. */ -extern mp *mpmont_mul(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*b*/); +extern mp *mpmont_mul(const mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*b*/); /* --- @mpmont_expr@ --- * * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = fake destination * @mp *a@ = base * @mp *e@ = exponent @@ -143,11 +157,12 @@ extern mp *mpmont_mul(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*b*/); * further modular arithmetic is to be performed on the result. */ -extern mp *mpmont_expr(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*e*/); +extern mp *mpmont_expr(const mpmont */*mm*/, mp */*d*/, + mp */*a*/, mp */*e*/); /* --- @mpmont_exp@ --- * * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = fake destination * @mp *a@ = base * @mp *e@ = exponent @@ -155,11 +170,11 @@ extern mp *mpmont_expr(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*e*/); * Returns: Result, %$a^e \bmod m$%. */ -extern mp *mpmont_exp(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*e*/); +extern mp *mpmont_exp(const mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*e*/); /* --- @mpmont_mexpr@ --- * * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = fake destination * @const mp_expfactor *f@ = pointer to array of factors * @size_t n@ = number of factors supplied @@ -174,12 +189,12 @@ extern mp *mpmont_exp(mpmont */*mm*/, mp */*d*/, mp */*a*/, mp */*e*/); * except that the %$g_i$% and result are in Montgomery form. */ -extern mp *mpmont_mexpr(mpmont */*mm*/, mp */*d*/, +extern mp *mpmont_mexpr(const mpmont */*mm*/, mp */*d*/, const mp_expfactor */*f*/, size_t /*n*/); /* --- @mpmont_mexp@ --- * * - * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = fake destination * @const mp_expfactor *f@ = pointer to array of factors * @size_t n@ = number of factors supplied @@ -189,7 +204,7 @@ extern mp *mpmont_mexpr(mpmont */*mm*/, mp */*d*/, * Use: Convenient interface over @mpmont_mexpr@. */ -extern mp *mpmont_mexp(mpmont */*mm*/, mp */*d*/, +extern mp *mpmont_mexp(const mpmont */*mm*/, mp */*d*/, const mp_expfactor */*f*/, size_t /*n*/); /*----- That's all, folks -------------------------------------------------*/