('(secret)', 'e=17211717500188099600568359980754843411619158773615107210604175251693280840607')

('(pub)', 'P', [15395204427582176828288140787842485523287609441875116301821962120815246523165L, 50972925492742860438904710940702337165443708937045978727200475339385656596802L])
('(pub)', 'Q', [54379753906375494870435841730857804485220282788952166778590221998932647140003L, 10542862410506568609831388869015874299228416047841839543681429954753256816461L])

Running first 3 prng outputs:
(1, '(pub)', 'rQ', [35085355371653225423586630859880842591953206407407058841980055259804158729190L, 17642611272178418467727837731595189739886826601909812633308152323681611088547L])
(2, '(pub)', 'rQ', [49239914602839976828302644790025511534554523065665058491590572885136905915179L, 8379389248907743775496938151971165397933498701999710225313910125914506554304L])
(3, '(pub)', 'rQ', [46183005535555809253746127552459776700939930656325949263490016267446108722393L, 39473869171813125431282259087083131250343325900631309529780555199432408899154L])

Predicting next prng output:
(4, '(pub)', 'sQ', [41709487842035263418257031106229715942424730101066821754384094668791902320545L, 47532264341329149594994883777108416347984877171111029279692175022666589952933L])
Running next prng output:
(4, '(pub)', 'rQ', [41709487842035263418257031106229715942424730101066821754384094668791902320545L, 47532264341329149594994883777108416347984877171111029279692175022666589952933L])

import sys
import random

def egcd(a, b):
	s = 0; news = 1
	r = b; newr = a
	while r != 0:
		quot = (newr / r)
		prev = s
		s = (news - (quot * prev))
		news = prev
		prev = r
		r = (newr - (quot * prev))
		newr = prev
	if (news < 0):
		news += b
	return news

def dub(p, a, b, n):
	# l = 3*x^2 + 2*a*x + 1 / 2*b*y
	xx = pow(p[0], 2, n)
	ya = (((3 * xx) + (2 * a * p[0]) + 1) % n)
	yd = ((2 * b * p[1]) % n)
	l = ((ya * egcd(yd, n)) % n)
	# xr = b*l^2 - a - 2*x
	ll = pow(l, 2, n)
	newx = (((b * ll) - a - (2 * p[0])) % n)
	# yr = ((3*x + a) * l) - b*l^3 - y
	yb = (((3 * p[0]) + a) % n)
	newy = (((yb * l) - (b * ll * l) - p[1]) % n)
	return [newx, newy]

def add(p, q, a, b, n):
	# l = (Qy - Py) / (Qx - Px)
	ya = ((q[1] - p[1]) % n)
	yd = ((q[0] - p[0]) % n)
	l = ((ya * egcd(yd, n)) % n)
	# xr = b*l^2 - a - Px - Qx
	ll = pow(l, 2, n)
	newx = (((b * ll) - a - p[0] - q[0]) % n)
	# yr = ((2*Px + Qx + a) * l) - b*l^3 - Py
	yb = (((2 * p[0]) + q[0] + a) % n)
	yc = ((b * ll * l) % n)
	newy = ((((yb * l) % n) - yc - p[1]) % n)
	return [newx, newy]

def mul(m, p, a, b, n):
	r = None
	while (m > 0):
		if ((m % 2) == 1):
			if (r == None):
				r = p
			else:
				r = add(r, p, a, b, n)
		p = dub(p, a, b, n)
		m = (m >> 1)
	return r

a = 486662
b = 1
n = 57896044618658097711785492504343953926634992332820282019728792003956564819949
try:
	p = [int(sys.argv[2]), int(sys.argv[3])]
except:
	p = [9, 43114425171068552920764898935933967039370386198203806730763910166200978582548]

# e is our secret key here
d = random.randint(3, n - 3)
e = random.randint(3, n - 3)
print("(secret)","e=%d"%e)
print("")

# let G = any base point
# calculate public point constants: P = edG & Q = dG
g = mul(d, p, a, b, n)
g = mul(e, g, a, b, n)
print("(pub)","P",g)
h = mul(d, p, a, b, n)
print("(pub)","Q",h)
print("")

# multiply random seed with the points: P' = rP & Q' = rQ
# these then equate to: P' = redG & Q' = rdG
# the 'key' to this is basing the next r' off of the output of P' and publishing the point info of Q' and repeat
r = random.randint(3, n - 3)
print("Running first 3 prng outputs:")
for x in range(1, 4):
	t = mul(r, g, a, b, n)
	u = mul(r, h, a, b, n)
	print(x,"(pub)","rQ",u)
	r = ((t[0] + t[1]) % n)
print("")

# to predict the next output then multiply: eQ' == erdQ == erdG == redG == P'
# you can now get r' and calculate the next P'' & Q''
print("Calculating next prng output (using e & latest Q):")
z = mul(e, [u[0], u[1]], a, b, n)
s = ((z[0] + z[1]) % n)
v = mul(s, z, a, b, n)
w = mul(s, h, a, b, n)
print(x+1,"(pub)","sQ",w)

# just to be sure
print("Running next prng output:")
for y in range(0, 1):
	t = mul(r, g, a, b, n)
	u = mul(r, h, a, b, n)
	print(x+1,"(pub)","rQ",u)
	r = ((t[0] + t[1]) % n)
print("")