viff

changeset 1302:7fd61e465bc6

Added implementation of Montgomery Exponentiation.
author Marc X. Makkes
date Mon, 19 Oct 2009 11:03:18 +0200
parents c44e2e1a9279
children deb7659465d2
files viff/montgomery_exponentiation.py
diffstat 1 files changed, 166 insertions(+), 0 deletions(-) [+]
line diff
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/viff/montgomery_exponentiation.py	Mon Oct 19 11:03:18 2009 +0200
     1.3 @@ -0,0 +1,166 @@
     1.4 +# Copyright 2009 VIFF Development Team.
     1.5 +#
     1.6 +# This file is part of VIFF, the Virtual Ideal Functionality Framework.
     1.7 +#
     1.8 +# VIFF is free software: you can redistribute it and/or modify it
     1.9 +# under the terms of the GNU Lesser General Public License (LGPL) as
    1.10 +# published by the Free Software Foundation, either version 3 of the
    1.11 +# License, or (at your option) any later version.
    1.12 +#
    1.13 +# VIFF is distributed in the hope that it will be useful, but WITHOUT
    1.14 +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.15 +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
    1.16 +# Public License for more details.
    1.17 +#
    1.18 +# You should have received a copy of the GNU Lesser General Public
    1.19 +# License along with VIFF. If not, see <http://www.gnu.org/licenses/>.
    1.20 +
    1.21 +
    1.22 +def sizeinbits(a):
    1.23 +    acc = 0
    1.24 +    while a != 0:
    1.25 +        acc += 1
    1.26 +        a = a >> 1
    1.27 +    return acc
    1.28 +
    1.29 +
    1.30 +def calc_r(n):
    1.31 +    """Compute r."""
    1.32 +    n_size = sizeinbits(n)
    1.33 +    r = 2**(n_size+1)
    1.34 +
    1.35 +    while gcd(r, n) != 1:
    1.36 +        r = r << 1
    1.37 +		
    1.38 +    return r
    1.39 +
    1.40 +
    1.41 +def calc_r_r_inv(n):
    1.42 +    """Compute r and r inverse."""
    1.43 +    assert (n % 2) == 1, "n must be odd."
    1.44 +    n_size = sizeinbits(n)
    1.45 +    r = 2**(n_size+1)
    1.46 +
    1.47 +    while gcd(r, n) != 1:
    1.48 +        r = r << 1
    1.49 +		
    1.50 +    return r, inversemod(r, n)
    1.51 +
    1.52 +
    1.53 +def calc_np(n, r):
    1.54 +    """Compute n'."""
    1.55 +    assert (n % 2) == 1, "n must be odd."
    1.56 +    n_prime = inversemod(n, r)
    1.57 +    return r - n_prime 
    1.58 +
    1.59 +
    1.60 +def inversemod(a, n):
    1.61 +    g, x, y = xgcd(a, n)
    1.62 +    if g != 1:
    1.63 +        raise ZeroDivisionError, (a,n)
    1.64 +    assert g == 1, "a must be coprime to n."
    1.65 +    return x%n
    1.66 +
    1.67 +
    1.68 +def gcd(a, b):                                       
    1.69 +    if a < 0:  a = -a
    1.70 +    if b < 0:  b = -b
    1.71 +    if a == 0: return b
    1.72 +    if b == 0: return a
    1.73 +    while b != 0: 
    1.74 +        (a, b) = (b, a%b)
    1.75 +    return a
    1.76 +
    1.77 +
    1.78 +def xgcd(a, b):
    1.79 +    if a == 0 and b == 0: return (0, 0, 1)
    1.80 +    if a == 0: return (abs(b), 0, b/abs(b))
    1.81 +    if b == 0: return (abs(a), a/abs(a), 0)
    1.82 +    x_sign = 1; y_sign = 1
    1.83 +    if a < 0: a = -a; x_sign = -1
    1.84 +    if b < 0: b = -b; y_sign = -1
    1.85 +    x = 1; y = 0; r = 0; s = 1
    1.86 +    while b != 0:
    1.87 +        (c, q) = (a%b, a/b)
    1.88 +        (a, b, r, s, x, y) = (b, c, x-q*r, y-q*s, r, s)
    1.89 +    return (a, x*x_sign, y*y_sign)
    1.90 +
    1.91 +
    1.92 +def montgomery_exponentiation_reduction(a, r , n ):
    1.93 +    return (a * r) % n
    1.94 +
    1.95 +
    1.96 +def montgomery_product(a, b, n_prime, size_of_r, r, n):
    1.97 +    t = a * b 
    1.98 +    m = (t * n_prime) & r -1
    1.99 +    u = (t + m * n ) >> size_of_r - 1
   1.100 +    if u >= n:
   1.101 +        return u -n 
   1.102 +    return u
   1.103 +
   1.104 +
   1.105 +def montgomery_exponentiation(a, x, n, n_prime, r):
   1.106 +    ah = (a * r) % n 
   1.107 +    xh = r % n 
   1.108 +    x_s = sizeinbits(x) - 1
   1.109 +    px = 2**x_s
   1.110 +    size_of_r = sizeinbits(r)
   1.111 +    while px != 0:
   1.112 +        t = xh * xh 
   1.113 +        m = (t * n_prime) & r -1
   1.114 +        u = (t + m * n ) >> size_of_r - 1
   1.115 +        if u >= n:
   1.116 +            xh = u - n 
   1.117 +        else:
   1.118 +            xh = u
   1.119 +        if (px & x) > 0:
   1.120 +            t = ah * xh 
   1.121 +            m = (t * n_prime) & r - 1
   1.122 +            u = (t + m * n ) >> size_of_r - 1
   1.123 +            if u >= n:
   1.124 +                xh =  u -n 
   1.125 +            else:
   1.126 +                xh = u
   1.127 +	px = px >> 1
   1.128 +
   1.129 +    m = (xh * n_prime) & r - 1
   1.130 +    u = (xh + m * n ) >> size_of_r - 1
   1.131 +    if u >= n:
   1.132 +        return u - n 
   1.133 +    return u
   1.134 +
   1.135 +def montgomery_exponentiation2(a, x, n, n_prime,  r):
   1.136 +    ah = (a * r) % n 
   1.137 +    xh = r % n 
   1.138 +    x_s = sizeinbits(x) - 1
   1.139 +    px = 2**x_s
   1.140 +    size_of_r = sizeinbits(r)
   1.141 +    while px != 0:
   1.142 +        xh = montgomery_product(xh, xh, n_prime, size_of_r, r, n)
   1.143 +        if (px & x) > 0:
   1.144 +            xh = montgomery_product(ah, xh, n_prime, size_of_r, r, n)
   1.145 +	px = px >> 1
   1.146 +
   1.147 +    x  = montgomery_product(xh, 1, n_prime, size_of_r, r, n)
   1.148 +    return x
   1.149 +
   1.150 +
   1.151 +
   1.152 +
   1.153 +# n = 2328734783
   1.154 +# r, r_inv = calc_r(n)
   1.155 +# n_prime = calc_np(n, r) 
   1.156 +
   1.157 +if __name__ ==  '__main__':
   1.158 +    n = 2328734783
   1.159 +    r, r_inv = calc_r_r_inv(n)
   1.160 +    n_prime = calc_np(n, r) 
   1.161 +    a = 2987
   1.162 +    x = 1093874
   1.163 +    print montgomery_exponentiation(a, x, n, n_prime, r)
   1.164 +    print pow(a, x, n)
   1.165 +    print a**x % n
   1.166 +
   1.167 +
   1.168 +
   1.169 +