Benutzer:Georg-Johann/Skripte

Extended Euclidean Algorithm

#!/usr/bin/env python

def xgcd (a, b):
	"""Return a triple (g,x,y) such that
	   * g = gcd (a,b) >= 0
	   * g = a*x + b*y
	   * gcd (0,0) := 0
	   Usage:
		(g, x, y) = xgcd (105, 170)
	"""
	(r0, r1) = (a, b)
	(s0, s1) = (1, 0)
	(t0, t1) = (0, 1)

	while r1:
		q = r0 // r1 # Floor division.
		(r0, r1) = (r1, r0 - q * r1)
		(s0, s1) = (s1, s0 - q * s1)
		(t0, t1) = (t1, t0 - q * t1)
	return (r0, s0, t0) if r0 > 0 else (-r0, -s0, -t0)

Arithmetic on Elliptic Curve mod n

Arithmetic in short Weierstraß form.

def EC_has_point (n, a, b, P):
	""" Test whether  P in E(a,b) mod n.  E(a,b): y^2 = x^3 + ax + b.
	    0 is represented as empty list ().
	"""
	if len(P) == 0:
		return True
	assert len(P) == 2
	(x,y) = P
	x2 = (x * x) % n
	y2 = (y * y) % n
	x3_ax_b = ((x2 + a) * x + b) % n
	return y2 == x3_ax_b

def ECadd (n, a, b, P, Q, check=False):
	""" Compute P + Q on E(a,b) mod n, E(a,b): y^2 = x^3 + ax + b.
	    P, Q in E(a,b).
	    0 is represented as empty list ().
	    Return a pair (R,d).
	    * If d = 1 then R is the sum of P and Q.
	    * If d > 1 then d is a divisor of n and R is ().
	"""
	if check:
		assert EC_has_point (n, a, b, P)
		assert EC_has_point (n, a, b, Q)
		
	if len(Q) == 0:		return (P,1)	# Q = 0: return P
	if len(P) == 0:		return (Q,1)	# P = 0: return Q

	(xQ,yQ) = Q
	(xP,yP) = P

	if (xQ - xP) % n == 0:
		if (yQ + yP) % n == 0:
			# Have P = -Q.
			return ((),1)
		# Have P = Q.
		(g,_,inv_2yP) = xgcd (n, 2 * yP)
		if g != 1:
			# 2*yP is not invertible mod n, found a divisor of n.
			return ((),g)
		# s = (3 xP^2 + a) / (2 yP)
		s = (3 * xP ** 2 + a) % n
		s = (s * inv_2yP) % n
	else:
		# Have P != Q != -P.
		(g,_,inv_xP_xQ) = xgcd (n, xP - xQ)
		if g != 1:
			# xP-xQ is not invertible mod n, found a divisor of n.
			return ((),g)
		s = ((yP - yQ) * inv_xP_xQ) % n

	xR = (s*s - xP - xQ) % n
	yR = (-yP - s * (xR - xP)) % n

	return ((xR,yR),1)

def ECmul (n, a, b, k, P, check=False):
	""" Compute k * P on E(a,b) mod n, E(a,b): y^2 = x^3 + ax + b.
	    k >= 0 in Z, P in E(a,b).
	    0 is represented as empty list ().
	    Return a pair (R,d).
	    * If d = 1 then R = k * P.
	    * If d > 1 then d is a divisor of n and R is ().
	"""
	# In the m-th round of the loop below, PP = 2^{m-1} * P.
	PP = P
	(kP,d) = ((),1)
	while True:
		if k & 1:	(kP,d) = ECadd (n, a, b, kP, PP, check)
		if d > 1:	return ((),d)
		if k < 2:	break
		k >>= 1
		(PP,d) = ECadd (n, a, b, PP, PP, check)
		if d > 1:	return ((),d)
	return (kP,1)