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))

Leave a comment