/* Script for pari/gp. Run as gp -q ecc-ref.gp */

out(apriv, A, bpriv, B, S) = print(	\
  "/* a_s */ \"", apriv, "\",\n",	\
  "/* a_x */ \"", component(A[1], 2), "\",\n",	\
  "/* a_y */ \"", component(A[2], 2), "\",\n",	\
  "/* b_s */ \"", bpriv, "\",\n",			\
  "/* b_x */ \"", component(B[1], 2), "\",\n",	\
  "/* b_y */ \"", component(B[2], 2), "\",\n",	\
  "/* s_x */ \"", component(S[1], 2), "\",\n",	\
  "/* s_y */ \"", component(S[2], 2), "\",");

p192 = 2^192 - 2^64 - 1;
b192 = 2455155546008943817740293915197451784769108058161191238065;
g = Mod([602046282375688656758213480587526111916698976636884684818, \
	    174050332293622031404857552280219410364023488927386650641], p192);
secp192 = ellinit(Mod([0,0,0,-3, b192], p192));
q = 6277101735386680763835789423176059013767194773182842284081;
if (ellpow(secp192, g, q) != [0], error("secp192 parameter error"));

a = 1+random(q-1);
b = 1+random(q-1);
A = ellpow(secp192, g, a);
B = ellpow(secp192, g, b);
S = ellpow(secp192, A, b);
if (S != ellpow(secp192, B, a), error("secp192 dh error"));
print("secp192");
out(a, A, b, B, S);

p224 = 2^224 - 2^96 + 1;
b224 = 18958286285566608000408668544493926415504680968679321075787234672564;
g = Mod([19277929113566293071110308034699488026831934219452440156649784352033,\
	 19926808758034470970197974370888749184205991990603949537637343198772], p224);

secp224 = ellinit(Mod([0,0,0,-3, b224], p224));
q = 26959946667150639794667015087019625940457807714424391721682722368061;
if (ellpow(secp224, g, q) != [0], error("secp224 parameter error"));

a = 1+random(q-1);
b = 1+random(q-1);
A = ellpow(secp224, g, a);
B = ellpow(secp224, g, b);
S = ellpow(secp224, A, b);
if (S != ellpow(secp224, B, a), error("secp224 dh error"));
print("secp224");
out(a, A, b, B, S);

p256 = 2^256 - 2^224 + 2^192 + 2^96 - 1;
b256 = 41058363725152142129326129780047268409114441015993725554835256314039467401291;
g = Mod([48439561293906451759052585252797914202762949526041747995844080717082404635286,\
	 36134250956749795798585127919587881956611106672985015071877198253568414405109], p256);

secp256 = ellinit(Mod([0,0,0,-3, b256], p256));
q = 115792089210356248762697446949407573529996955224135760342422259061068512044369;
if (ellpow(secp256, g, q) != [0], error("secp256 parameter error"));

a = 1+random(q-1);
b = 1+random(q-1);
A = ellpow(secp256, g, a);
B = ellpow(secp256, g, b);
S = ellpow(secp256, A, b);
if (S != ellpow(secp256, B, a), error("secp256 dh error"));
print("secp256");
out(a, A, b, B, S);

p384 = 2^384 - 2^128 - 2^96 + 2^32 - 1;
b384 = 27580193559959705877849011840389048093056905856361568521428707301988689241309860865136260764883745107765439761230575;
g = Mod([26247035095799689268623156744566981891852923491109213387815615900925518854738050089022388053975719786650872476732087,\
	 8325710961489029985546751289520108179287853048861315594709205902480503199884419224438643760392947333078086511627871], p384);

secp384 = ellinit(Mod([0,0,0,-3, b384], p384));
q = 39402006196394479212279040100143613805079739270465446667946905279627659399113263569398956308152294913554433653942643;
if (ellpow(secp384, g, q) != [0], error("secp384 parameter error"));

a = 1+random(q-1);
b = 1+random(q-1);
A = ellpow(secp384, g, a);
B = ellpow(secp384, g, b);
S = ellpow(secp384, A, b);
if (S != ellpow(secp384, B, a), error("secp384 dh error"));
print("secp384");
out(a, A, b, B, S);

p521 = 2^521 - 1;
b521 = 1093849038073734274511112390766805569936207598951683748994586394495953116150735016013708737573759623248592132296706313309438452531591012912142327488478985984;
g = Mod([2661740802050217063228768716723360960729859168756973147706671368418802944996427808491545080627771902352094241225065558662157113545570916814161637315895999846,\
	 3757180025770020463545507224491183603594455134769762486694567779615544477440556316691234405012945539562144444537289428522585666729196580810124344277578376784], p521);

secp521 = ellinit(Mod([0,0,0,-3, b521], p521));
q = 6864797660130609714981900799081393217269435300143305409394463459185543183397655394245057746333217197532963996371363321113864768612440380340372808892707005449;
if (ellpow(secp521, g, q) != [0], error("secp521 parameter error"));

a = 1+random(q-1);
b = 1+random(q-1);
A = ellpow(secp521, g, a);
B = ellpow(secp521, g, b);
S = ellpow(secp521, A, b);
if (S != ellpow(secp521, B, a), error("secp521 dh error"));
print("secp521");
out(a, A, b, B, S);

p25519 = 2^255 - 19;
b25519 = 486662;
x = Mod(9, p25519);
y = sqrt(x^3 + b25519*x^2 + x);
g = [x, y];

curve25519 = ellinit(Mod([0, b25519, 0, 1, 0], p25519));
q = 2^252 + 27742317777372353535851937790883648493;
if (ellpow(curve25519, g, q) != [0], error("curve25519 parameter error"));

a = 1+random(q-1);
b = 1+random(q-1);
A = ellpow(curve25519, g, a);
B = ellpow(curve25519, g, b);
S = ellpow(curve25519, A, b);
if (S != ellpow(curve25519, B, a), error("curve25519 dh error"));
print("curve25519");
out(a, A, b, B, S);

/* Convert point on curve25519 to a point on the twisted edwards curve */
beta = -sqrt(Mod(-486664, p25519));
ed25519(p) = [p[1] * beta / p[2], (p[1] - 1) / (p[1] + 1)];

Ae = ed25519(A);
Be = ed25519(B);
Se = ed25519(S);
print("ed25519");
out(a, Ae, b, Be, Se);

quit