(* curverecords

   Recorded data for Apple ECC curves.
   
   R. Crandall
   3 Apr 2001
   

 *)

pointQ[x_] := (JacobiSymbol[x^3 + c x^2 + a x + b, p] > -1);

(* Next, binary expansion for very old M'ca versions,
   otherwise use IntegerDigits[.,2]. *)
bitList[k_] := Block[{li = {}, j = k},
	While[j > 0,
	    li = Append[li, Mod[j,2]];
	    j = Floor[j/2];
	];
	Return[Reverse[li]];
	];
	
ellinv[n_] := PowerMod[n,-1,p];
(* Next, obtain actual x,y coords via normalization:
   {x,y,z} := {X/Z^2, Y/Z^3, 1}. *)
normalize[pt_] := Block[{z,z2,z3},
		If[pt[[3]] == 0, Return[pt]];
		z = ellinv[pt[[3]]];
		z2 = Mod[z^2,p];
		z3 = Mod[z z2,p];
		Return[{Mod[pt[[1]] z2, p], Mod[pt[[2]] z3, p], 1}];
		];

ellneg[pt_] := Mod[pt * {1,-1,1}, p];
ellsub[pt1_, pt2_] := elladd[pt1, ellneg[pt2]];
elldouble[pt_] := Block[{x,y,z,m,y2,s},
	x = pt[[1]]; y = pt[[2]]; z = pt[[3]];
	If[(y==0) || (z==0), Return[{1,1,0}]];
	m = Mod[3 x^2 + a Mod[Mod[z^2,p]^2,p],p];
	z = Mod[2 y z, p];
	y2 = Mod[y^2,p];
	s = Mod[4 x y2,p]; 
	x = Mod[m^2 - 2s,p];
	y = Mod[m(s - x) - 8 y2^2,p];
	Return[{x,y,z}];
];

elladd[pt0_, pt1_] := Block[
	{x0,y0,z0,x1,y1,z1,
	t1,t2,t3,t4,t5,t6,t7},
	x0 = pt0[[1]]; y0 = pt0[[2]]; z0 = pt0[[3]];
	x1 = pt1[[1]]; y1 = pt1[[2]]; z1 = pt1[[3]];
	If[z0 == 0, Return[pt1]];
	If[z1 == 0, Return[pt0]];

	t1 = x0;
	t2 = y0;
	t3 = z0;
	t4 = x1;
	t5 = y1;
	If[(z1 != 1),
		t6 = z1;
		t7 = Mod[t6^2, p];
		t1 = Mod[t1 t7, p];
		t7 = Mod[t6 t7, p];
		t2 = Mod[t2 t7, p];
	];
	t7 = Mod[t3^2, p];
	t4 = Mod[t4 t7, p];
	t7 = Mod[t3 t7, p];
	t5 = Mod[t5 t7, p];
	t4 = Mod[t1-t4, p];
	t5 = Mod[t2 - t5, p];
	If[t4 == 0, If[t5 == 0,
				    Return[elldouble[pt0]],
	   				Return[{1,1,0}]
	   			]
	];
	t1 = Mod[2t1 - t4,p];
	t2 = Mod[2t2 - t5, p];
	If[z1 != 1, t3 = Mod[t3 t6, p]];
	t3 = Mod[t3 t4, p];
	t7 = Mod[t4^2, p];
	t4 = Mod[t4 t7, p];
	t7 = Mod[t1 t7, p];
	t1 = Mod[t5^2, p];
	t1 = Mod[t1-t7, p];
	t7 = Mod[t7 - 2t1, p];
	t5 = Mod[t5 t7, p];
	t4 = Mod[t2 t4, p];
	t2 = Mod[t5-t4, p];
	If[EvenQ[t2], t2 = t2/2, t2 = (p+t2)/2];
	Return[{t1, t2, t3}];
];
		
(* Next, elliptic-multiply a normalized pt by k. *)
elliptic[pt_, k_] := Block[{pt2, hh, kk, hb, kb, lenh, lenk},
	If[k==0, Return[{1,1,0}]];
	hh = Reverse[bitList[3k]];
	kk = Reverse[bitList[k]];
	pt2 = pt;
	lenh = Length[hh];
	lenk = Length[kk];
	Do[
		pt2 = elldouble[pt2];
		hb = hh[[b]];
		If[b <= lenk, kb = kk[[b]], kb = 0];
		If[{hb,kb} == {1,0},
			pt2 = elladd[pt2, pt],
			If[{hb, kb} == {0,1},
			pt2 = ellsub[pt2, pt]]
		]
	   ,{b, lenh-1, 2,-1}
	 ];
	Return[pt2];
];

(* Next, provide point-finding functions. *)

(* Next, perform (a + b w)^n (mod p), where pair = {a,b}, w2 = w^2. *)
pow[pair_, w2_, n_, p_] := Block[{bitlist, z},
    bitlist = bitList[n];
    z = pair;
	Do[	
	   zi = Mod[z[[2]]^2,p];
	   z = {Mod[z[[1]]^2 + w2 zi, p], Mod[2 z[[1]] z[[2]], p]};   
	   If[bitlist[[q]] == 1,
	       zi = Mod[pair[[2]] z[[2]], p];
	   	   z = {Mod[pair[[1]] z[[1]] + w2 zi, p],
	   	         Mod[pair[[1]] z[[2]] + pair[[2]] z[[1]], p]};
	   ],
	   {q,2,Length[bitlist]}
    ];
    Return[z]
];


sqrt[x_, p_] := Module[{t, b, w2},
    If[Mod[x,p] == 0, Return[0]];
	If[Mod[p,4] == 3, Return[PowerMod[x, (p+1)/4, p]]];
	If[Mod[p,8] == 5,
		b = PowerMod[x, (p-1)/4, p];
		If[b==1, Return[PowerMod[x, (p+3)/8, p]],
			Return[Mod[2x PowerMod[4x, (p-5)/8,p],p]]
		]
	];
	t = 2;
    While[True,
      w2 = Mod[t^2 - x, p];
      If[JacobiSymbol[w2,p] == -1, Break[]];
      ++t
    ];
    (* Next, raise (t + Sqrt[w2])^((p+1)/2). *)
    t = pow[{t,1},w2, (p+1)/2, p];
    Return[t[[1]]];
    ];

findpoint[start_] := Block[{x = start, y, s},
	While[True,
	   s = Mod[x(Mod[x^2+a,p]) + b, p];
	   y = sqrt[s, p];
	   If[Mod[y^2, p] == s, Break[]];
	   ++x;
	];
	Return[{x, y, 1}]
];

report[a_] := Module[{ia = IntegerDigits[a,65536]},    
      Prepend[Reverse[ia], Length[ia]]
      ];

report[a_] := Module[{ia = IntegerDigits[a,65536]},    
      Prepend[Reverse[ia], Length[ia]]
      ];
    
    
(* Case of Weierstrass/feemod curve. *)
p = 2^127 + 57675
report[p]
r = 512000; s = 512001;

a = Mod[-3 r s^3, p]
report[a]
b = Mod[-2 r s^5, p]
report[b]
pt = findpoint[3];
pt
plusOrd = 170141183460469231756943134065055014407
report[plusOrd]
PrineQ[plusOrd]
minusOrd = 170141183460469231706431473366713312401
report[minusOrd]
PrimeQ[minusOrd]
elliptic[pt, plusOrd]
elliptic[pt, minusOrd]

\( (*\ Case\ of\ Weierstrass/gen . \ mod\ \(curve . \)\ *) 

    p\  = \ 1654338658923174831024422729553880293604080853451; 

    Mod[p, 4]\)Length[IntegerDigits[p, 2]]report[p]\(PrimeQ[p]

\)\(a\  = \ \(-152\); 
report[a]\)b  =  Mod[722,  p]report[b]ptplus  =  
      findpoint[1245904487553815885170631576005220733978383542270]ptminus  =  
      findpoint[1173563507729187954550227059395955904200719019884]\(plusOrd\  = \ \ 1654338658923174831024425147405519522862430265804; 

    report[plusOrd]\)PrimeQ[plusOrd]minusOrd  =  2  p + 2  -  plusOrdreport[minusOrd]\(PrimeQ[minusOrd]
\)pt2  =  
      elliptic[ptplus,  plusOrd/(2^2 * 23 * 359 * 479 * 102107)]
        pt3  =  elliptic[ptminus,  minusOrd/(2^2 * 5^2 * 17^2)]

(* Case of Weierstrass/feemod curve. *)
p = 2^160 + 5875
report[p]
PrimeQ[p]
r = 512; s = 513;

a = Mod[-3 r s^3, p]
report[a]
b = Mod[2 r s^5, p]
report[b]
pt = findpoint[3];
pt
plusOrd =  1461501637330902918203687223801810245920805144027
report[plusOrd]
PrimeQ[plusOrd]
minusOrd = 1461501637330902918203682441630755793391059953677
report[minusOrd]
PrimeQ[minusOrd]

elliptic[pt, plusOrd]
elliptic[pt, minusOrd]

\( (*\ Case\ of\ NIST\ P - 192. \ *) 

    p\  = \ 6277101735386680763835789423207666416083908700390324961279; 

    Mod[p, 4]\)Length[IntegerDigits[p, 2]]report[p]\(PrimeQ[p]

\)\(a\  = \ \(-3\); 
report[a]\)b  =  
      Mod[-2455155546008943817740293915197451784769108058161191238065,  
        p]report[b]\(pt\  = \ findpoint[3]; 
pt\)plusOrd  =   
      6277101735386680763835789423176059013767194773182842284081report[plusOrd]PrimeQ[plusOrd]minusOrd  =  2  p + 2  -  plusOrdreport[minusOrd]PrimeQ[minusOrd]elliptic[pt,  plusOrd]pt2  =  elliptic[pt,  23]pt  =  elliptic[pt2,  minusOrd/23]report[minusOrd/23]

normalize[pt2]

report[%[[1]]]

minusOrdplusOrd

IntegerDigits[13, 2]

minusOrd/23\(	report[minusOrd/23]\)

\(ord\  = \ 1024120625531724089187207582052247831; 

    Floor[4^Length[IntegerDigits[plusOrd, \ 2]]/plusOrd]\)

report[%]

IntegerDigits[
      6277101735386680763835789423207666416083908700390324961279, 2]

abJacobiSymbol[3^3  +  a * 3  -  b,  p]

minusOrd

\(o\  = 272917466755942641905903887966924948626114027286861201673; 

    Floor[4^Length[IntegerDigits[o, 2]]/o]\)

report[%]

IntegerDigits[1654338658923174831024422729553880293604080853451, 2]