michael@0: Modular Reduction
michael@0: 
michael@0: Usually, modular reduction is accomplished by long division, using the
michael@0: mp_div() or mp_mod() functions.  However, when performing modular
michael@0: exponentiation, you spend a lot of time reducing by the same modulus
michael@0: again and again.  For this purpose, doing a full division for each
michael@0: multiplication is quite inefficient.
michael@0: 
michael@0: For this reason, the mp_exptmod() function does not perform modular
michael@0: reductions in the usual way, but instead takes advantage of an
michael@0: algorithm due to Barrett, as described by Menezes, Oorschot and
michael@0: VanStone in their book _Handbook of Applied Cryptography_, published
michael@0: by the CRC Press (see Chapter 14 for details).  This method reduces
michael@0: most of the computation of reduction to efficient shifting and masking
michael@0: operations, and avoids the multiple-precision division entirely.
michael@0: 
michael@0: Here is a brief synopsis of Barrett reduction, as it is implemented in
michael@0: this library.
michael@0: 
michael@0: Let b denote the radix of the computation (one more than the maximum
michael@0: value that can be denoted by an mp_digit).  Let m be the modulus, and
michael@0: let k be the number of significant digits of m.  Let x be the value to
michael@0: be reduced modulo m.  By the Division Theorem, there exist unique
michael@0: integers Q and R such that:
michael@0: 
michael@0: 	 x = Qm + R, 0 <= R < m
michael@0: 
michael@0: Barrett reduction takes advantage of the fact that you can easily
michael@0: approximate Q to within two, given a value M such that:
michael@0: 
michael@0: 	                  2k
michael@0: 	                 b
michael@0: 	    M = floor( ----- )
michael@0: 	                 m
michael@0: 
michael@0: Computation of M requires a full-precision division step, so if you
michael@0: are only doing a single reduction by m, you gain no advantage.
michael@0: However, when multiple reductions by the same m are required, this
michael@0: division need only be done once, beforehand.  Using this, we can use
michael@0: the following equation to compute Q', an approximation of Q:
michael@0: 
michael@0:                      x
michael@0:             floor( ------ ) M
michael@0:                       k-1
michael@0:                      b
michael@0: Q' = floor( ----------------- )
michael@0:                     k+1
michael@0:                    b
michael@0: 
michael@0: The divisions by b^(k-1) and b^(k+1) and the floor() functions can be
michael@0: efficiently implemented with shifts and masks, leaving only a single
michael@0: multiplication to be performed to get this approximation.  It can be
michael@0: shown that Q - 2 <= Q' <= Q, so in the worst case, we can get out with
michael@0: two additional subtractions to bring the value into line with the
michael@0: actual value of Q.
michael@0: 
michael@0: Once we've got Q', we basically multiply that by m and subtract from
michael@0: x, yielding:
michael@0: 
michael@0:    x - Q'm = Qm + R - Q'm
michael@0: 
michael@0: Since we know the constraint on Q', this is one of:
michael@0: 
michael@0:       R
michael@0:       m + R
michael@0:       2m + R
michael@0: 
michael@0: Since R < m by the Division Theorem, we can simply subtract off m
michael@0: until we get a value in the correct range, which will happen with no
michael@0: more than 2 subtractions:
michael@0: 
michael@0:      v = x - Q'm
michael@0: 
michael@0:      while(v >= m)
michael@0:        v = v - m
michael@0:      endwhile
michael@0: 
michael@0: 
michael@0: In random performance trials, modular exponentiation using this method
michael@0: of reduction gave around a 40% speedup over using the division for
michael@0: reduction.
michael@0: 
michael@0: ------------------------------------------------------------------
michael@0:  This Source Code Form is subject to the terms of the Mozilla Public
michael@0:  # License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0:  # file, You can obtain one at http://mozilla.org/MPL/2.0/.