My Rough/Slow Version Of Daniel J. Bernstein’s Elliptic Curve25519 Used For ECC-DHKE-ElGamal Pub/Pri Key Enc/Dec In Py & JS

Edit: I need to verify the math on the code below and optimize it as much as possible!

The message being encrypted by the public key and decrypted by the private key is [71, 73], as you can see below, in both Python and JavaScript.

Output:

point=[31338, 27137639241784266545116051737587625422240160049866585375429726394143891718199] ((y^2)%p)=508712386423338 == curve(x)=508712386423338
pub=[53977350860898987307302427493695333251161207860469554940047122752933735592647, 26210366092361168208973994509836336790598036553932356811338803611050445452207]
& [1445723211146282456267893060679523547963682350054864882998693087748930328718927, 2685102703073507930430754637571136937637097240643254035250750702470189871650545]
mesg=[71, 73]
time=254 ms

JS:

<script src="jsbn.js"></script>
<script src="jsbn2.js"></script>
<script>
var c = new BigInteger("486662", 10);
var p = new BigInteger("57896044618658097711785492504343953926634992332820282019728792003956564819949", 10);
var two = new BigInteger("2", 10);
var three = new BigInteger("3", 10);

function safe_sub(a, p)
{
	if (a.compareTo(BigInteger.ZERO) < 0)
	{
		var r = a.abs().divide(p);
		a = a.add(r.multiply(p));
		if (a.compareTo(BigInteger.ZERO) < 0)
		{
			a = a.add(p);
		}
	}
	return a;
}

function point_add(q, r, p)
{
	var ys = safe_sub(q[1].subtract(r[1]), p);
	var xs = safe_sub(q[0].subtract(r[0]), p);
	var s = ys.multiply(xs.modInverse(p)).mod(p);
	var rx = s.modPow(two, p);
	rx = safe_sub(rx.subtract(q[0]), p);
	rx = safe_sub(rx.subtract(r[0]), p);
	var rp = safe_sub(q[0].subtract(rx), p);
	var ry = safe_sub(s.multiply(rp).mod(p).subtract(q[1]), p);
	return [rx, ry];
}

function point_dub(a, q, p)
{
	var ttt = three.multiply(q[0].modPow(two, p)).mod(p);
	var tt = two.multiply(q[1]).mod(p);
	var s = ttt.add(a).multiply(tt.modInverse(p)).mod(p);
	tt = two.multiply(q[0]).mod(p);
	var rx = safe_sub(s.modPow(two, p).subtract(tt), p);
	var rp = safe_sub(q[0].subtract(rx), p);
	var ry = safe_sub(s.multiply(rp).mod(p).subtract(q[1]), p);
	return [rx, ry];
}

function point_mul(n, q, p)
{
	//if(this.isInfinity()) return this;
	//if(k.signum() == 0) return this.curve.getInfinity();
	
	var e = n;
	var h = e.multiply(three);
	
	var neg = [q[0], safe_sub(q[1].negate(), p)];
	var R = q;
	
	var i;
	for (i = (h.bitLength() - 2); i > 0; --i)
	{
		R = point_dub(BigInteger.ONE, R, p);
		
		var hBit = h.testBit(i);
		var eBit = e.testBit(i);
		
		if (hBit != eBit)
		{
			R = point_add(hBit ? q : neg, R, p);
		}
	}
	
	return R;
}

function ord_r(r, n)
{
	var k = three;
	while (1)
	{
		if (n.modPow(k, r).equals(BigInteger.ONE))
		{
			return k;
		}
		k = k.add(BigInteger.ONE);
	}
}

function ifrexp(x)
{
	var e = BigInteger.ZERO;
	while (x.mod(two).equals(BigInteger.ZERO))
	{
		x = x.divide(two);
		e = e.add(BigInteger.ONE);
	}
	return [x, e];
}

function tonelli(a, p)
{
	var pmo = p.subtract(BigInteger.ONE);
	if (a.modPow(pmo.divide(two), p).equals(pmo))
	{
		return -1;
	}
	var t = ifrexp(pmo);
	var s = t[0], e = t[1];
	var n = two;
	while (n.compareTo(p) < 0)
	{
		if (n.modPow(pmo.divide(two), p).equals(pmo))
		{
			break;
		}
		n = n.add(BigInteger.ONE);
	}
	var x = a.modPow(s.add(BigInteger.ONE).divide(two), p);
	var b = a.modPow(s, p);
	var g = n.modPow(s, p);
	var r = e;
	while (1)
	{
		var m = BigInteger.ZERO;
		while (m.compareTo(r) < 0)
		{
			if (ord_r(p, b).equals(two.pow(m)))
			{
				break;
			}
			if (m.add(BigInteger.ONE).equals(r))
			{
				break;
			}
			m = m.add(BigInteger.ONE);
		}
		if (m.equals(BigInteger.ZERO))
		{
			return x;
		}
		var pmi = two.pow(r.subtract(m).subtract(BigInteger.ONE));
		x = x.multiply(g.modPow(pmi, p)).mod(p);
		pmi = two.pow(r.subtract(m));
		g = g.modPow(pmi, p);
		b = b.multiply(g).mod(p);
		if (b.equals(BigInteger.ONE))
		{
			return x;
		}
		r = m;
	}
	return -1;
}

function curve_25519(x, p)
{
	var y2 = x.modPow(three, p).add(c.multiply(x.modPow(two, p))).add(x).mod(p);
	return tonelli(y2, p);
}

var i = new BigInteger("31337", 10);
var pnt = null
while (pnt == null)
{
	var q = [i, curve_25519(i, p)];
	if (q[1] == -1)
	{
		i = i.add(BigInteger.ONE);
		continue;
	}
	var x = q[0].modPow(three, p).add(c.multiply(q[0].modPow(two, p))).add(q[0]).mod(p);
	var y = q[1].modPow(two, p);
	if (x.equals(y))
	{
		pnt = q;
		document.write("point=["+pnt[0].toString(10)+", "+pnt[1].toString(10)+"] ((y^2)%p)="+y.toString(10)+" == "+"curve(x)="+x.toString(10)+"<br>");
	}
	i = i.add(BigInteger.ONE);
}

var m = new BigInteger("71", 10);
var n = new BigInteger("73", 10);

var nta = new Date().getTime();

var a = new BigInteger("23002347587565544268625339214141417360725798840824195817183950850507406831337", 10);
var aG = point_mul(a, pnt, p);

var b = new BigInteger("33951982198225404751798578318443600923434009233142787987410481710685782131337", 10);
var bG = point_mul(b, pnt, p);
var baG = point_mul(b, aG, p);

var mbaG = [m.multiply(baG[0]), n.multiply(baG[1])];

var abG = point_mul(a, bG, p);

var ntb = new Date().getTime();

document.write("pub=["+bG[0]+", "+bG[1]+"]<br> & ["+mbaG[0]+", "+mbaG[1]+"]<br>"+"mesg=["+mbaG[0].divide(abG[0]).toString(10)+", "+mbaG[1].divide(abG[1]).toString(10)+"]<br>"+"time="+(ntb-nta)+" ms<br>");
</script>

Output:

('point', [31338, 27137639241784266545116051737587625422240160049866585375429726394143891718199L], 'y^2%p=', 508712386423338L, '==', 'curve(x)', 508712386423338L)
pub=[53977350860898987307302427493695333251161207860469554940047122752933735592647L, 26210366092361168208973994509836336790598036553932356811338803611050445452207L]
 & [1445723211146282456267893060679523547963682350054864882998693087748930328718927L, 2685102703073507930430754637571136937637097240643254035250750702470189871650545L]
mesg=[71, 73]
time=145.143032074 ms

Py:

import time

p = (pow(2, 255) - 19)

def egcd(a, b):
	if (a == 0):
		return (b, 0, 1)
	else:
		(g, y, x) = egcd(b % a, a)
		return (g, x - (b / a) * y, y)

def inv_mod(a, m):
	(g, x, y) = egcd(a, m)
	if (g != 1):
		return -1
	else:
		return (x % m)

def pow_mod(b, e, m):
	r = 1
	b = (b % m)
	while (e > 0):
		if ((e % 2) == 1):
			r = ((r * b) % m)
		e = (e / 2)
		b = ((b * b) % m)
	return r

def safe_sub(a, p):
	if (a < 0):
		a += ((abs(a) / p) * p)
		if (a < 0):
			a += p
	return a

def point_add(q, r, p):
	ys = safe_sub(q[1] - r[1], p)
	xs = safe_sub(q[0] - r[0], p)
	s = ((ys * inv_mod(xs, p)) % p)
	rx = pow_mod(s, 2, p)
	rx = safe_sub(rx - q[0], p)
	rx = safe_sub(rx - r[0], p)
	rp = safe_sub(q[0] - rx, p)
	ry = safe_sub(((s * rp) % p) - q[1], p)
	return [rx, ry]

def point_dub(a, q, p):
	ttt = ((3 * pow_mod(q[0], 2, p)) % p)
	tt = ((2 * q[1]) % p)
	s = (((ttt + a) * inv_mod(tt, p)) % p)
	tt = ((2 * q[0]) % p)
	rx = safe_sub(pow_mod(s, 2, p) - tt, p)
	rp = safe_sub(q[0] - rx, p)
	ry = safe_sub(((s * rp) % p) - q[1], p)
	return [rx, ry]

def zpoint_mul(n, q, p):
	r = q; m = 1
	h = [[m, r]]; l = 1
	while (m < n):
		t = (m * 2)
		if (t <= n):
			r = point_dub(1, r, p); m = t
			h.append([m, r]); l += 1
		else:
			x = (l - 1)
			while (x > -1):
				while (m < n):
					t = (m + h[x][0])
					if (t <= n):
						r = point_add(h[x][1], r, p); m = t
					else:
						break
				x -= 1
	return r

def bitleng(a):
	b = 0
	while (a > 0):
		a = (a / 2)
		b += 1
	return b

def testbit(a, b):
	if ((a & pow(2, b)) > 0):
		return 1
	return 0

def point_mul(n, q, p):
	#if(this.isInfinity()) return this;
	#if(k.signum() == 0) return this.curve.getInfinity();
	
	e = n
	h = (e * 3)
	
	neg = [q[0], safe_sub(-1 * q[1], p)]
	R = q
	
	i = (bitleng(h) - 2)
	while (i > 0):
		R = point_dub(1, R, p)
		
		hBit = testbit(h, i)
		eBit = testbit(e, i)
		
		if (hBit != eBit):
			if (hBit):
				R = point_add(q, R, p)
			else:
				R = point_add(neg, R, p)
		
		i -= 1
	
	return R

def ord_r(r, n):
	k = 3
	while (1):
		if (pow_mod(n, k, r) == 1):
			return k
		k += 1

def ifrexp(x):
	e = 0
	while ((x % 2) == 0):
		x /= 2
		e += 1
	return (x, e)

def tonelli(a, p):
	if (pow_mod(a, (p - 1) / 2, p) == (p - 1)):
		raise ValueError("no sqrt possible")
	(s, e) = ifrexp(p - 1)
	n = 2
	while (n < p):
		if (pow_mod(n, (p - 1) / 2, p) == (p - 1)):
			break
		n += 1
	x = pow_mod(a, (s + 1) / 2, p)
	b = pow_mod(a, s, p)
	g = pow_mod(n, s, p)
	r = e
	while (1):
		m = 0
		while (m < r):
			if (ord_r(p, b) == pow(2, m)):
				break
			if ((m + 1) == r):
				break
			m += 1
		if (m == 0):
			return x
		x = ((x * pow_mod(g, pow(2, (r - m - 1)), p)) % p)
		g = pow_mod(g, pow(2, (r - m)), p)
		b = ((b * g) % p)
		if (b == 1):
			return x
		r = m
	return -1

def curve_25519(x, p):
	y2 = ((pow_mod(x, 3, p) + (486662 * pow_mod(x, 2, p)) + x) % p)
	return tonelli(y2, p)

i = 31337
pnt = None
while (pnt == None):
	try:
		q = [i, curve_25519(i, p)]
	except:
		i += 1
		continue
	x = ((pow_mod(q[0], 3, p) + (486662 * pow_mod(q[0], 2, p)) + q[0]) % p)
	y = pow_mod(q[1], 2, p)
	if (x == y):
		print("point",q,"y^2%p=",y,"==","curve(x)",x)
		pnt = q
	i += 1

m = 71; n = 73

x = time.time()

a = 23002347587565544268625339214141417360725798840824195817183950850507406831337
aG = point_mul(a, pnt, p)

b = 33951982198225404751798578318443600923434009233142787987410481710685782131337
bG = point_mul(b, pnt, p)
baG = point_mul(b, aG, p)

mbaG = [m * baG[0], n * baG[1]]

abG = point_mul(a, bG, p)

y = time.time()

print("pub=%s\n & %s\nmesg=[%d, %d]\ntime=%s ms" % (bG, mbaG, mbaG[0] / abG[0], mbaG[1] / abG[1], (y - x) * 1000))

Advertisements
My Rough/Slow Version Of Daniel J. Bernstein’s Elliptic Curve25519 Used For ECC-DHKE-ElGamal Pub/Pri Key Enc/Dec In Py & JS

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s