/rev/7fd61e465bc6 changeset: 1302:7fd61e465bc6 user: Marc X. Makkes date: Mon Oct 19 11:03:18 2009 +0200 summary: Added implementation of Montgomery Exponentiation.
diffstat: viff/montgomery_exponentiation.py | 166 +++++++++++++++++++++++++++++++++++++++++ 1 files changed, 166 insertions(+), 0 deletions(-) diffs (170 lines): diff -r c44e2e1a9279 -r 7fd61e465bc6 viff/montgomery_exponentiation.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/viff/montgomery_exponentiation.py Mon Oct 19 11:03:18 2009 +0200 @@ -0,0 +1,166 @@ +# Copyright 2009 VIFF Development Team. +# +# This file is part of VIFF, the Virtual Ideal Functionality Framework. +# +# VIFF is free software: you can redistribute it and/or modify it +# under the terms of the GNU Lesser General Public License (LGPL) as +# published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# VIFF is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General +# Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with VIFF. If not, see <http://www.gnu.org/licenses/>. + + +def sizeinbits(a): + acc = 0 + while a != 0: + acc += 1 + a = a >> 1 + return acc + + +def calc_r(n): + """Compute r.""" + n_size = sizeinbits(n) + r = 2**(n_size+1) + + while gcd(r, n) != 1: + r = r << 1 + + return r + + +def calc_r_r_inv(n): + """Compute r and r inverse.""" + assert (n % 2) == 1, "n must be odd." + n_size = sizeinbits(n) + r = 2**(n_size+1) + + while gcd(r, n) != 1: + r = r << 1 + + return r, inversemod(r, n) + + +def calc_np(n, r): + """Compute n'.""" + assert (n % 2) == 1, "n must be odd." + n_prime = inversemod(n, r) + return r - n_prime + + +def inversemod(a, n): + g, x, y = xgcd(a, n) + if g != 1: + raise ZeroDivisionError, (a,n) + assert g == 1, "a must be coprime to n." + return x%n + + +def gcd(a, b): + if a < 0: a = -a + if b < 0: b = -b + if a == 0: return b + if b == 0: return a + while b != 0: + (a, b) = (b, a%b) + return a + + +def xgcd(a, b): + if a == 0 and b == 0: return (0, 0, 1) + if a == 0: return (abs(b), 0, b/abs(b)) + if b == 0: return (abs(a), a/abs(a), 0) + x_sign = 1; y_sign = 1 + if a < 0: a = -a; x_sign = -1 + if b < 0: b = -b; y_sign = -1 + x = 1; y = 0; r = 0; s = 1 + while b != 0: + (c, q) = (a%b, a/b) + (a, b, r, s, x, y) = (b, c, x-q*r, y-q*s, r, s) + return (a, x*x_sign, y*y_sign) + + +def montgomery_exponentiation_reduction(a, r , n ): + return (a * r) % n + + +def montgomery_product(a, b, n_prime, size_of_r, r, n): + t = a * b + m = (t * n_prime) & r -1 + u = (t + m * n ) >> size_of_r - 1 + if u >= n: + return u -n + return u + + +def montgomery_exponentiation(a, x, n, n_prime, r): + ah = (a * r) % n + xh = r % n + x_s = sizeinbits(x) - 1 + px = 2**x_s + size_of_r = sizeinbits(r) + while px != 0: + t = xh * xh + m = (t * n_prime) & r -1 + u = (t + m * n ) >> size_of_r - 1 + if u >= n: + xh = u - n + else: + xh = u + if (px & x) > 0: + t = ah * xh + m = (t * n_prime) & r - 1 + u = (t + m * n ) >> size_of_r - 1 + if u >= n: + xh = u -n + else: + xh = u + px = px >> 1 + + m = (xh * n_prime) & r - 1 + u = (xh + m * n ) >> size_of_r - 1 + if u >= n: + return u - n + return u + +def montgomery_exponentiation2(a, x, n, n_prime, r): + ah = (a * r) % n + xh = r % n + x_s = sizeinbits(x) - 1 + px = 2**x_s + size_of_r = sizeinbits(r) + while px != 0: + xh = montgomery_product(xh, xh, n_prime, size_of_r, r, n) + if (px & x) > 0: + xh = montgomery_product(ah, xh, n_prime, size_of_r, r, n) + px = px >> 1 + + x = montgomery_product(xh, 1, n_prime, size_of_r, r, n) + return x + + + + +# n = 2328734783 +# r, r_inv = calc_r(n) +# n_prime = calc_np(n, r) + +if __name__ == '__main__': + n = 2328734783 + r, r_inv = calc_r_r_inv(n) + n_prime = calc_np(n, r) + a = 2987 + x = 1093874 + print montgomery_exponentiation(a, x, n, n_prime, r) + print pow(a, x, n) + print a**x % n + + + + _______________________________________________ viff-commits mailing list [email protected] http://lists.viff.dk/listinfo.cgi/viff-commits-viff.dk
