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