Crypto-related Codes in Sage

Code Crypto Sage

RSA Code

Basis Algorithm

#https://www.math.ucdavis.edu/~anne/FQ2010/Number_Theory_RSA.html
def rsa_kg(bits):
    # only prove correctness up to 1024 bits
    proof = (bits <= 1024)
    p = next_prime(ZZ.random_element(2**(bits//2+1)), proof=proof)
    q = next_prime(ZZ.random_element(2**(bits//2+1)), proof=proof)
    n = p*q
    phi_n = (p-1)*(q-1)
    while True:
        e = ZZ.random_element(1,phi_n)
        if gcd(e,phi_n) == 1: break
    d = lift(Mod(e,phi_n)^(-1))
    return e, d, n 
def encrypt(m, e, n):
    return lift(Mod(m,n)^e)
def decrypt(c, d, n):
    return lift(Mod(c,n)^d) 
def encode(s):
    s = str(s)
    return sum(ord(s[i])*256^i for i in range(len(s)))
def decode(n):
    n = Integer(n)
    v = []
    while n != 0:
        v.append(chr(n % 256))
        n //= 256     # this replaces n by floor(n/256)
    return ''.join(v) 
e,d,n = rsa_kg(1024) 
m = encode('Meet me at 4.'); m 
c = encrypt(m,e,n); c
M = decrypt(c,d,n); M 
decode(M) 

Common Modulus Attack

Given e, d, N and e1, compute d1 such that d1 is a valid decrypting exponent for the public key (e1 , N )

k1 = gcd(e1 , e*d − 1)
k2 = (e*d − 1) / k1
d1 = xgcd(e1, k2) [1]

Small d attack

攻击成功的条件：q < p < 2q and that d < n^{1/4}/3

0、生成弱秘钥

def rsa_weak_kg(bits):
    # only prove correctness up to 1024 bits
    proof = (bits <= 1024)
    p = next_prime(ZZ.random_element(2**(bits//2+1)), proof=proof)
    q = next_prime(ZZ.random_element(2**(bits//2+1)), proof=proof)
    n = p*q
    phi_n = (p-1)*(q-1)
    while True:
        d = ZZ.random_element(2**(bits//4))
        if gcd(d,phi_n) == 1 and 36*pow(d, 4) < n: break
    e = lift(Mod(d,phi_n)^(-1))
    return e, d, n 

1、算e/n的连分数

fra = (e/n).continued_fraction()

2、依次选择fra对e/N的逼近，记为k/d，利用这个k和d进行N的分解

# Define f(x)= (x - p)(x - q) = 0
# x = var('x')
# a , b, c = 1, 1, 1
# qe = (a*x^2 + b*x + c == 0)
# Trials
def small_d_attack(fra):
  for i in range(1, len(fra)):
    k = fra.numerator(i)
    d = fra.denominator(i)
    if k != 0:
        phi = (d*e - 1)/k
        a = 1
        b = -(n - phi + 1)
        c = n
        # Define f(x)= (x - p)(x - q) = 0
        x = var('x')
        qe = (a*x^2 + b*x + c == 0)
        sol = solve(qe, x)
        p = sol[0].right()
        q = sol[1].right()
        if p.is_integer() and q.is_integer() and q*p == n:
            print sol[0].right(), sol[1].right()
            print "We found the decrypting exponent.", d
            break

Discrete Logarithm

Sage下Dlog的求解命令

F.<a>=GF(2^16+1)
g=F.multiplicative_generator()
u=g**12345
discrete_log(u,g) 
# discrete_log_rho(u,g) can't work for composite order group.

DLOG实例生成

# p is a prime and q is (p-1)/2 and also a prime. Thus Zq* is a subgroup of Zp* with prime order.
def dlog_gen(n):
    p = next_prime(n)
    while not is_prime( floor((p-1)/2) ):
        p = next_prime(p)
    x = randint(1,p-1)
    y = randint(1,p-1)
    g = x*x % p
    h = y*y % p
    return [p,floor( (p-1)/2 ),g,h]

Baby-Step Giant-Step DLOG

# Generate a random DLOG instance.
def dlog_gen(bits):
    p = next_prime(ZZ.random_element(2**bits))
    q = floor((p-1)/2)
    while not is_prime(q):
        p = next_prime(p)
        q = floor((p-1)/2)
    g = 1
    while g == 1:
        x = randint(1, p-1)
        g = x*x % p
    return [p,q,g]
# Baby Step Giant Step DLP problem a = g**x mod n
# Input: g,a,n where g is the generator, a = g^x mod n,
def dlog_bgs(a, g, n):
    s = floor(sqrt(n))
    A = []
    B = []
    for r in range(0,s):
        value = a*(g^r) % n
        A.append(value)
    for t in range(1,s+1):
        value = g^(t*s) % n
        B.append(value)
    #print A
    #print B
    x1,x2 =0,0
    for u in A:
        for v in B:
            if u == v:
                x1 = A.index(u)
                x2 = B.index(v)
                #print x1,x2
                break
    #print 'the value of x is ', ((x2+1)*s - x1) % n # Answer
    return ((x2+1)*s - x1) % n

Pollard rho discrete logarithm

# A simple Pollard rho discrete logarithm
# implementation and has some limitations:
# 1. p must be a prime that equals 2q + 1
# 2. q must be a prime, too
# 3. g generates a sub group with order q
# 4. a belongs to <g>, the subgroup generated by g
# these four limitations made this program simpler
# x = log_g(a) in Zp
# assert: classify(1) != 1
def classify(x):
    # 3n + 2 -> S0
    # 3n     -> S1
    # 3n + 1 -> S2
    return (x + 1) % 3
def succssor(x, s, t, p, q, g, a):
    c = classify(x)
    if c == 0:
        return (g * x) % p, (s + 1) % q, t
    elif c == 1:
        return (x * x) % p, (2 * s) % q, (2 * t) % q
    else: # c == 2
        return (a * x) % p, s, (t + 1) % q
def dlog_rho(a, g, p, q):
    xa, sa, ta = 1, 0, 0
    xb, sb, tb = succssor(1, 0, 0, p, q, g, a)
    while xa != xb:
          xa, sa, ta = succssor(xa, sa, ta, p, q, g, a)
          xb, sb, tb = succssor(xb, sb, tb, p, q, g, a)
          xb, sb, tb = succssor(xb, sb, tb, p, q, g, a)
    s, t = sa - sb, tb - ta
    if s == 0:
        return 'fail'
    if s < 0:
        s = s + q
    if t < 0:
        t = t + q
    res = xgcd(t, q)[1]
    if res < 0 :
        res = res + q
    res = (res * s) % q
    return res
# How to use in Sage and compare with discrete_log and discrete_log_rho
bits = 30
p,q,g = dlog_gen(30)
x = 120
a = Mod(g,p)^x
time discrete_log(a, g, p)
time discrete_log_rho(a, g, q)
#dlog_bgs(g, a, p)
#dlog_rho(g, a, p, q)

WStein 2017 Elementary Number Theory

Chapter 6 ECC

#Code from WStein 2017 Elementary Number Theory
from sage.libs.libecm import ecmfactor
ecmfactor(2^128+1,1000,sigma=227140902)
# Pollard p-1 method 
def pollard(N, B=10^5, stop=10):
    m = prod([p^int(math.log(B)/math.log(p)) for p in prime_range(B+1)])
    for a in [2..stop]:
         x = (Mod(a,N)^m - 1).lift()
         if x == 0: continue
         g = gcd(x, N)
         if g!=1 or g != N: return g
    return 1
# Lenstra's EC method
def ecm(N, B=10^3, trials=10):
    m = prod([p^int(math.log(B)/math.log(p)) for p in prime_range(B+1)])
    R = Integers(N)
    # Make Sage think that R is a field:
    R.is_field = lambda : True
    for _ in range(trials):
        while True:
            a = R.random_element()
            if gcd(4*a.lift()^3 + 27, N) == 1: break
        try:
            m * EllipticCurve([a, 1])([0,1])
        except ZeroDivisionError, msg:
            # msg: "Inverse of <int> does not exist"
            return gcd(Integer(str(msg).split()[2]), N)
    return 1
#set_random_seed(2)
#ecm(5959, B=20)
#ecm(next_prime(10^20)*next_prime(10^7), B=10^3)