bn_mp_montgomery_reduce.c 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118
  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, tomstdenis@gmail.com, http://math.libtomcrypt.com
  16. */
  17. /* computes xR**-1 == x (mod N) via Montgomery Reduction */
  18. int
  19. mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
  20. {
  21. int ix, res, digs;
  22. mp_digit mu;
  23. /* can the fast reduction [comba] method be used?
  24. *
  25. * Note that unlike in mul you're safely allowed *less*
  26. * than the available columns [255 per default] since carries
  27. * are fixed up in the inner loop.
  28. */
  29. digs = n->used * 2 + 1;
  30. if ((digs < MP_WARRAY) &&
  31. n->used <
  32. (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
  33. return fast_mp_montgomery_reduce (x, n, rho);
  34. }
  35. /* grow the input as required */
  36. if (x->alloc < digs) {
  37. if ((res = mp_grow (x, digs)) != MP_OKAY) {
  38. return res;
  39. }
  40. }
  41. x->used = digs;
  42. for (ix = 0; ix < n->used; ix++) {
  43. /* mu = ai * rho mod b
  44. *
  45. * The value of rho must be precalculated via
  46. * montgomery_setup() such that
  47. * it equals -1/n0 mod b this allows the
  48. * following inner loop to reduce the
  49. * input one digit at a time
  50. */
  51. mu = (mp_digit) (((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK);
  52. /* a = a + mu * m * b**i */
  53. {
  54. register int iy;
  55. register mp_digit *tmpn, *tmpx, u;
  56. register mp_word r;
  57. /* alias for digits of the modulus */
  58. tmpn = n->dp;
  59. /* alias for the digits of x [the input] */
  60. tmpx = x->dp + ix;
  61. /* set the carry to zero */
  62. u = 0;
  63. /* Multiply and add in place */
  64. for (iy = 0; iy < n->used; iy++) {
  65. /* compute product and sum */
  66. r = ((mp_word)mu) * ((mp_word)*tmpn++) +
  67. ((mp_word) u) + ((mp_word) * tmpx);
  68. /* get carry */
  69. u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
  70. /* fix digit */
  71. *tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK));
  72. }
  73. /* At this point the ix'th digit of x should be zero */
  74. /* propagate carries upwards as required*/
  75. while (u) {
  76. *tmpx += u;
  77. u = *tmpx >> DIGIT_BIT;
  78. *tmpx++ &= MP_MASK;
  79. }
  80. }
  81. }
  82. /* at this point the n.used'th least
  83. * significant digits of x are all zero
  84. * which means we can shift x to the
  85. * right by n.used digits and the
  86. * residue is unchanged.
  87. */
  88. /* x = x/b**n.used */
  89. mp_clamp(x);
  90. mp_rshd (x, n->used);
  91. /* if x >= n then x = x - n */
  92. if (mp_cmp_mag (x, n) != MP_LT) {
  93. return s_mp_sub (x, n, x);
  94. }
  95. return MP_OKAY;
  96. }
  97. #endif
  98. /* $Source: /cvs/libtom/libtommath/bn_mp_montgomery_reduce.c,v $ */
  99. /* $Revision: 1.3 $ */
  100. /* $Date: 2006/03/31 14:18:44 $ */