"Fossies" - the Fresh Open Source Software Archive 
As a special service "Fossies" has tried to format the requested text file into HTML format (style:
standard) with prefixed line numbers.
Alternatively you can here
view or
download the uninterpreted source code file.
1 _grp_ell_ops := rec();
2 _grp_ell_ops.Print :=function(r)
3 if Length(r.L)= 2 then
4 Print("EllipticCurve over ",r.F," defined by \n y^2 = x^3 ");
5 if r.L[1]<> 0 then Print("+ ",r.L[1],"*x "); fi;
6 if r.L[2]<> 0 then Print("+ ",r.L[2]); fi;
7 Print("\n");
8 else Print("EllipticCurve over ",r.F," defined by \n y^2 ");
9 if r.L[1]<> 0 then Print("+ ",r.L[1],"*x*y "); fi;
10 if r.L[2]<> 0 then Print("+ ",r.L[2],"*y "); fi;
11 Print(" = x^3 ");
12 if r.L[3]<> 0 then Print("+ ",r.L[3],"*x^2 ");fi;
13 if r.L[4]<> 0 then Print("+ ",r.L[4],"*x "); fi;
14 if r.L[5]<> 0 then Print("+ ",r.L[5]); fi;
15 Print("\n");
16 fi;
17 end;
18
19 _grp_ell_ops.\=:= function(EC1,EC2)
20 local L1,L2;
21 if EC1.F<>EC2.F then return FALSE;
22 fi;
23 L1:=EC1.L;
24 L2:=EC2.L;
25 if Length(L1) < Length(L2) then
26 if L2[1] <> 0 then return FALSE;
27 elif L2[2] <> 0 then return FALSE;
28 elif L2[3] <>0 then return FALSE;
29 fi;
30 fi;
31 return L1 = L2;
32 end;
33
34 _elt_grp_ell_ops:= rec();
35
36
37
38 _NewE:=function()
39 ReadLib("elliptic");
40 Print("\n");
41 end;
42
43 _Bin := function (n)
44 local i,L;
45 i:= 0;
46 L:= [];
47 while n>0 do
48 L[i+1] := n mod 2;
49 n := Div(n,2);
50 i:= i+1;
51 od;
52 return L;
53 end;
54
55 ## bar from overwriting global X which belongs to constant ZX
56 PX:= function(P)
57 if P = Zero(Parent(P)) then return INFTY;
58 else
59 return P.x;
60 fi;
61 end;
62
63 PY:= function(P)
64 if P = Zero(Parent(P)) then return INFTY;
65 else
66 return P.y;
67 fi;
68 end;
69
70 _WNF:= function(F,L)
71 local a1,a3,a2,a4,a6,b2,b4,b6;
72 a1:= L[1];
73 a3:= L[2];
74 a2:= L[3];
75 a4:= L[4];
76 a6:= L[5];
77 b2:=a1^2+4*a2;
78 b4:= 2*a4+a1*a3;
79 b6:= a3^2+4*a6;
80 return [b2,b4,b6];
81 end;
82
83 _WeierstrassNormalForm := function(F,L)
84 local a1,a2,a3,a4,a6,b2,b4,b6,A,B,c4,c6;
85 if Length(L) = 2 then return [F,L]; fi;
86 if Characteristic(F) <5 then
87 return Error("Characteristic of base field is 2 or 3, cannot create WeierstrassNormalForm\n");
88 fi;
89 a1:= L[1];
90 a3:= L[2];
91 a2:= L[3];
92 a4:= L[4];
93 a6:= L[5];
94 b2:=a1^2+4*a2;
95 b4:= 2*a4+a1*a3;
96 b6:= a3^2+4*a6;
97 c4:= b2^2-24*b4;
98 c6:= -b2^3+36*b2*b4-216*b6;
99 A:= -27*c4;
100 B:= -54*c6;
101 return [F,[Coerce(F,A),Coerce(F,B)]];
102 end;
103
104 _Discriminant:= function(F,L)
105 local b2,b4,b6,b8,a1,a3,a2,a4,a6;
106 if Length(L) = 5 then
107 a1:= L[1];
108 a3:= L[2];
109 a2:= L[3];
110 a4:= L[4];
111 a6:= L[5];
112 b2:= a1^2+4*a2;
113 b4:= 2*a4 +a1*a3;
114 b6:= a3^2 +4*a6;
115 b8:= (a1^2)*a6 +4*a2*a6-a1*a3*a4+a2*(a3^2)-a4^2;
116 return (-(b2^2)*b8-8*(b4^3)-27*(b6^2)+9*b2*b4*b6);
117 else
118 return (-16*(4*(L[1]^3)+27*(L[2]^2)));
119 fi;
120 end;
121
122
123
124 _jInv:= function(F,L)
125 local EWei,LWei,c4;
126 if Characteristic(F)>3 then
127 if Length(L) = 2 then
128 LWei:=L;
129 else
130 EWei:= _WeierstrassNormalForm(F,L);
131 LWei:= EWei[2];
132 fi;
133 return 1728*((4*(LWei[1]^3))/(4*(LWei[1]^3)+27*(LWei[2]^2)));
134 else
135 LWei:= _WNF(F,L);
136 c4:= LWei[1]^2-24*LWei[2];
137 return c4^3/_Discriminant(F,L);
138 fi;
139 end;
140
141 _Discriminant2:= function(EC)
142 return _Discriminant(EC.F,EC.L);
143 end;
144
145
146
147
148 IsSupersingular:= function(E)
149 local f,eltseq,p,BR,PA,E2,m,t;
150 BR:= BaseRing(E);
151 p:= Characteristic(BR);
152 if p = 2 then
153 #if E.L[3]=0 and E.L[4]=0 and E.L[5] = 0 then
154 return jInvariant(E) = 0;
155 else
156 t:=Size(BR)+1-Size(E);
157 return t mod p = 0;
158 # PA:= PolynomialAlgebra(BR);
159 # if Length(E.L) = 2 then
160 # E2:= E;
161 # else E2:=ShortWeierstrassNormalForm(E);
162 # fi;
163 # m:= Coerce(Integers(),((p-1)/2));
164 # f:= Coerce(PA,[E2.L[2],E2.L[1],0,1])^m;
165 # eltseq:= Eltseq(f);
166 # if eltseq[p] = 0 then return true;
167 # else return false;
168 # fi;
169 fi;
170 end;
171
172
173 jInvariant:= function(EC)
174 local L,F,a1,a3,a2,a4,a6,b2,b4,b6,b8,d;
175 L:=Coefficients(EC);
176 F:=BaseRing(EC);
177 if Length(L)=2 then
178 a1:=0;
179 a3:=0;
180 a2:=0;
181 a4:=L[1];
182 a6:=L[2];
183 else a1:=L[1];
184 a3:=L[2];
185 a2:=L[3];
186 a4:=L[4];
187 a6:=L[5];
188 fi;
189 b2:= a1^2+4*a2;
190 b4:= 2*a4 +a1*a3;
191 b6:= a3^2 +4*a6;
192 b8:= (a1^2)*a6 +4*a2*a6-a1*a3*a4+a2*(a3^2)-a4^2;
193 d:= (-(b2^2)*b8-8*(b4^3)-27*(b6^2)+9*b2*b4*b6);
194 return (b2^2-24*b4)^3/d;
195 end;
196
197
198
199 _PointMult := function(n,P)
200 local Li,M,Q,i,j,y1,L;
201 if P.x = "INFTY POINT" then return P; fi;
202 if n = 0 then
203 return rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops);
204 fi;
205 if n <0 then
206 n := -n;
207 Li:= P.parent.L;
208 if Length(Li)=2 then
209 P := rec (x := P.x, y := -P.y, parent := P.parent, type := elt-grp^ell, operations:=_elt_grp_ell_ops);
210 else
211 y1:=(-Li[1]*P.x-Li[2]-P.y);
212 P := rec (x := P.x, y := y1, parent := P.parent, type := elt-grp^ell,operations:=_elt_grp_ell_ops);
213 fi;
214 fi;
215 L:= _Bin(n);
216 M:= [P];
217 j := 1;
218 for i in [2..Length(L)] do
219 M[i]:= M[i-1]+M[i-1];
220 od;
221 while L[j] = 0 do
222 j := j+1;
223 od;
224 Q := M[j];
225
226 for i in [(j+1)..Length(L)] do
227 if L[i] = 1 then
228 Q := M[i]+Q;
229 fi;
230 od;
231 return Q;
232 end;
233
234 _elt_grp_ell_ops.\*:= function(n,P) return _PointMult(n,P); end;
235
236 _Parent:= function(P)
237 return P.parent;
238 end;
239
240 _GetPoly:= function(F,L)
241 local kx,kxy,f,i,x,M;
242 kx := PolynomialAlgebra (F);
243 kxy:= PolynomialAlgebra (kx);
244 x:= Coerce(kxy,kx.1);
245 M:=[];
246 for i in [1..Length(L)] do
247 M[i]:=Coerce(kxy,L[i]);
248 od;
249 if Length(L) = 2 then
250 f:= kxy.1^2 - x^3 -(M[1])* x -M[2];
251 else
252 f := (kxy.1^2) +M[1]*(kx.1)*(kxy.1) + M[2]*kxy.1- (kx.1^3) - M[3]*(kx.1^2) - M[4]*kx.1 -M[5];
253 fi;
254 return f;
255 end;
256
257 _GetFunf:= function(F,L)
258 return FunctionField(_GetPoly(F,L));
259 end;
260
261 EllipticFunctionField := function(EC)
262 local L,kx2,kxy2,f;
263 L:= EC.L;
264 kx2:= RationalFunctionField(BaseRing(EC));
265 kxy2:= PolynomialAlgebra(kx2);
266 f:= _GetPoly(BaseRing(EC),L);
267 f:= Coerce(kxy2,f);
268 return FunctionField(f);
269 end;
270
271 _OldFunf:= function(EC)
272 local f;
273 f:= _GetPoly(BaseRing(EC),EC.L);
274 return FunctionField(f);
275 end;
276
277
278 _Zero:= function(EC)
279 return rec(x := "INFTY POINT",parent:= EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
280 end;
281
282 _AddShort := function(P,Q)
283 local x3, y3,m;
284 if P.x=Q.x then
285 if P.y <> Q.y then
286 return rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops);
287 elif (P.y = Q.y and P.y = 0) then
288 return rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops);
289 fi;
290 fi;
291 if P=Q then
292 m:= Coerce(P.parent.F,(3*(P.x)^2+P.parent.L[1])/(2*P.y));
293 x3 := m^2 -2*P.x;
294 y3 := m*(P.x-x3)-P.y;
295 else
296 m := Coerce(P.parent.F,(Q.y-P.y)/ (Q.x - P.x));
297 x3 := m^2 - P.x -Q.x;
298 y3 := m*(P.x -x3) - P.y;
299 fi;
300
301 return rec(x:= Coerce(P.parent.F,x3),y := Coerce(P.parent.F,y3), parent := P.parent, type := elt-grp^ell, operations :=_elt_grp_ell_ops);
302 end;
303
304 _AddLong := function(P,Q)
305 local l, v,L,x1,x2,y1,y2,x3,y3;
306 x1 := P.x;
307 x2 := Q.x;
308 y1 := P.y;
309 y2 := Q.y;
310 L:= P.parent.L;
311
312 if x1 = x2 then
313 if (y1 + y2 + L[1]*x2+L[2] = 0) then
314 return rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops);
315 fi;
316 l := (3*(x1^2) +2*L[3]*x1+L[4]-L[1]*y1)/(2*y1+L[1]*x1+L[2]);
317 v := (-(x1^3) +L[4]*x1+2*L[5]-L[2]*y1)/(2*y1+L[1]*x1+L[2]);
318 else
319 l := (y2-y1)/(x2-x1);
320 v := (y1*x2-y2*x1)/(x2-x1);
321 fi;
322 x3 := l^2+L[1]*l-L[3]-x1-x2;
323 y3 := -(l+L[1])*x3 -v -L[2];
324 return rec(x:= Coerce(P.parent.F,x3),y := Coerce(P.parent.F,y3), parent := P.parent, type := elt-grp^ell, operations :=_elt_grp_ell_ops);
325 end;
326
327
328
329 _elt_grp_ell_ops.\+:= function(P,Q)
330 local x3, y3,m;
331 if ((P.parent) <> (Q.parent)) then
332 Print("Cannot add Points on different curves\n");
333 return;
334 elif P.x = "INFTY POINT" then return Q;
335 elif Q.x = "INFTY POINT" then return P;
336 fi;
337 if Length(Q.parent.L) =2 then return _AddShort(P,Q);
338 else return _AddLong(P,Q);
339 fi;
340 end;
341
342 Subcurve:=function(E1,E2)
343 if E1.L<> E2.L then
344 return false;
345 fi;
346 return IsSubset(E1.F,E2.F);
347 end;
348
349
350 _elt_grp_ell_ops.\=:= function(P,Q)
351 local E1,E2;
352 E1:=Q.parent;
353 E2:=P.parent;
354 if E1 <> E2 then
355 if not Subcurve(E1,E2) and not Subcurve(E2,E1) then
356 return FALSE;
357 fi;
358 fi;
359 if Q.x = "INFTY POINT" then
360 if P.x <> "INFTY POINT" then
361 return FALSE;
362 else return TRUE;
363 fi;
364 fi;
365 if P.x = "INFTY POINT" then return FALSE;
366 else return (Q.x = P.x and Q.y = P.y);
367 fi;
368 end;
369
370
371
372
373 _elt_grp_ell_ops.Print :=function(r)
374 if r.x = "INFTY POINT" then
375 Print("INFTY POINT");
376 else Print("[",r.x,",",r.y,"]");
377 fi;
378 end;
379
380
381
382 _SwapVars := function (f)
383 local poly, listpol1,listpol2,kx,kxy,l,g;
384 poly := 0;
385 kxy := Parent(f);
386 kx := CoefficientRing(kxy);
387 listpol1 := List(f);
388 listpol2 := [];
389 for l in [1..Length(listpol1)] do
390 if listpol1[l] = 0 then
391 listpol2[l] := 0;
392 else listpol2[l] := List(listpol1[l]);
393 fi;
394 od;
395 for l in [1..Length(listpol1)] do
396 listpol2[l]:=Coerce(kxy,listpol2[l]);
397 od;
398 for l in [1..Length(listpol1)] do
399 poly := poly + (listpol2[l])*kx.1^(l-1);
400 od;
401 return poly;
402 end;
403
404
405
406
407
408 _Discriminant:= function(F,L)
409 local LWei,b8,a1,a3,a2,a4,a6,b2,b4,b6;
410 if Length(L)=2 then
411 return Coerce(F,-16*(4*L[1]^3+27*L[2]^2));
412 fi;
413 LWei:= _WNF(F,L);
414 a1:= L[1];
415 a3:= L[2];
416 a2:= L[3];
417 a4:= L[4];
418 a6:= L[5];
419 b2:=LWei[1];
420 b4:= LWei[2];
421 b6:= LWei[3];
422 b8:= a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2;
423 return Coerce(F,-b2^2*b8-8*b4^3-27*b6^2+9*b2*b4*b6);
424 end;
425
426 ExistsEllipticCurve:= function(p,L)
427 local i,F,Funf,g,t;
428 if (Length(L)<>2 and Length(L)<>5) then
429 return false;
430 else
431 t := Type(p);
432 if Is(t,fld) then
433 if IsRec(p) then F:= p.base;
434 else F:=p;
435 fi;
436 elif IsPrimePower(p) = FALSE then
437 return false;
438 else
439 F := FiniteField(p);
440 fi;
441 for i in [1..Length(L)] do
442 L[i] := Coerce(F,L[i]);
443 od;
444 if _Discriminant(F,L) = 0 then return false; fi;
445 Funf := _GetFunf(F,L);
446 g := Genus(Funf);
447 if g<>1 then return false;fi;
448 return rec(base:=true,ext1:=rec(F:= F,L:= L,size:=-1,type := grp^ell,operations :=_grp_ell_ops ));
449 fi;
450 end;
451
452
453 EllipticCurve := function (p,L)
454 local i,F,Funf,g,t;
455 if (Length(L)<>2 and Length(L)<>5) then
456 return Error("Length of defining list must be 2 or 5\n");
457 else
458 t := Type(p);
459 if Is(t,fld) then
460 if IsRec(p) then F:= p.base;
461 else F:=p;
462 fi;
463 elif IsPrimePower(p) = FALSE then
464 return Error("first argument must be a prime power\n");
465 else
466 F := FiniteField(p);
467 fi;
468 for i in [1..Length(L)] do
469 L[i] := Coerce(F,L[i]);
470 od;
471 if _Discriminant(F,L) = 0 then return Error("Discriminant of defining polynomial is 0\n"); fi;
472 Funf := _GetFunf(F,L);
473 g := Genus(Funf);
474 if g<>1 then return Error("Genus of FunctionField is ",g,".\n");fi;
475 return rec(F:= F,L:= L,size:=-1,type := grp^ell,operations :=_grp_ell_ops );
476 fi;
477 end;
478
479
480 ChangeBaseRing:=function(E,F)
481 local E2;
482 if Length(E.L)=2 then
483 E2:= EllipticCurve (F,[Coerce(F,E.L[1]),Coerce(F,E.L[2])]);
484 E2.size:=-1;
485 else E2:= EllipticCurve(F,[Coerce(F,E.L[1]),Coerce(F,E.L[2]),Coerce(F,E.L[3]),Coerce(F,E.L[4]),Coerce(F,E.L[5])]);
486 E2.size:=-1;
487 fi;
488 return E2;
489 end;
490
491 ShortWeierstrassNormalForm:= function(EC)
492 local L;
493 if Characteristic(EC.F) >3 then
494 L:=_WeierstrassNormalForm(EC.F,EC.L);
495 return EllipticCurve(L[1],L[2]);
496 else return Error("Characteristic of base field is 2 or 3, Short WNF cannot be computed \n");
497 fi;
498 end;
499
500
501 #WeierstrassNormalForm := function(EC)
502 #if Characteristic(EC.F) = 2 then
503 # Print("cannot create short WeierstrassNormalForm\n");
504 # return;
505 #elif Characteristic(EC.F) = 3 then
506 # Print("char=3, cannot create complete short WeierstrassNormalForm");
507 #fi;
508 #return EllipticCurve(_WeierstrassNormalForm(EC.F,EC.L)[1],_WeierstrassNormalForm(EC.F,EC.L)[2]);
509 #end;
510
511
512 _EvalPlace := function(r,g)
513 local kx,kxy,l,lis,i;
514 kxy := Parent(DefiningPolynomial(Parent(g)));
515 kx := CoefficientRing(kxy);
516 l := List(g);
517 lis := [];
518 for i in [1.. Length(l)] do
519 lis[i] := Evaluate(l[i],r);
520 od;
521 return Coerce(kx, lis);
522 end;
523
524
525
526 _elt_grp_ell_ops.\-:= function(P,Q)
527 local S;
528 if Q.x = "INFTY POINT" then
529 S:= Copy(Q);
530 else S := _PointMult(-1,Q);
531 fi;
532 return P+S;
533 end;
534
535 Point:= function(pl,EC)
536 local L,Funf,kx,kxy,gen,g,r1,r2,i,lis,k,poly,eval1,eval2,f,x,y;
537 if Degree(pl)<>1 then return Error("Degree of given place is not 1."); fi;
538 kx := PolynomialAlgebra(EC.F);
539 kxy := PolynomialAlgebra(kx);
540 if IsFinite(pl) = FALSE then
541 return rec(x := "INFTY POINT",parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
542 else
543 g:=Generators(Ideal(pl));
544 r1:=Roots(Coerce(kx,g[1]))[1][1];
545 lis :=List(g[2]);
546 eval1 := FALSE;
547 eval2 := FALSE;
548 f:=DefiningPolynomial(_OldFunf(EC));
549 for k in [1..Length(lis)] do
550 if (Degree(Denominator(lis[k])) > 0 and eval1 = FALSE) then
551 eval1 := TRUE;
552 poly := _SwapVars(f);
553
554 elif (Degree(lis[k])>0 and eval2 = FALSE) then
555 eval2 := TRUE;
556 poly := _EvalPlace(r1,g[2]);
557 fi;
558 od;
559 if eval1 = TRUE then
560 r2 :=Roots(Evaluate (poly,r1))[1][1];
561 elif eval2 = TRUE then
562 r2 := Roots(poly)[1][1];
563 else
564 r2 :=Roots(Coerce(kxy,lis))[1][1];
565 fi;
566 x:=r1;
567 y:=Coerce(CoefficientRing(kx),r2);
568 if not Evaluate(Evaluate(f,y),x) = Element(EC.F,0) then
569 return Error("Given place does not belong to given curve");
570 fi;
571 return rec(x:=x,y:=y,parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
572 fi;
573 end;
574
575 Points := function (EC)
576 local L,Funf,kx,kxy,gen,p,g,r1,r2,i,retL,lis,k,poly,eval1,eval2;
577 if (Type(EC) <> grp^ell) then
578 Print(Type(EC),"\n");
579 Print("usage : Points (<grp^ell>);\n");
580 return;
581 fi;
582 Funf := _OldFunf(EC);
583 p:=Places(Funf,1);
584 kx := PolynomialAlgebra(EC.F);
585 kxy := PolynomialAlgebra(kx);
586 retL := [];
587 for i in [1..Length(p)] do
588 if IsFinite(p[i]) = FALSE then
589 retL[i] :=rec(x := "INFTY POINT",parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
590 else
591 g:=Generators(Ideal(p[i]));
592 r1:=Roots(Coerce(kx,g[1]))[1][1];
593 lis :=List(g[2]);
594 eval1 := FALSE;
595 eval2 := FALSE;
596 for k in [1..Length(lis)] do
597 if (Degree(Denominator(lis[k])) > 0 and eval1 = FALSE) then
598 eval1 := TRUE;
599 poly := _SwapVars(DefiningPolynomial(Funf));
600
601 elif (Degree(lis[k])>0 and eval2 = FALSE) then
602 eval2 := TRUE;
603 poly := _EvalPlace(r1,g[2]);
604 fi;
605 od;
606 if eval1 = TRUE then
607 r2 :=Roots(Evaluate (poly,r1))[1][1];
608 elif eval2 = TRUE then
609 r2 := Roots(poly)[1][1];
610 else
611 r2 :=Roots(Coerce(kxy,lis))[1][1];
612 fi;
613 retL[i]:= rec(x:=r1,y:=Coerce(CoefficientRing(kx),r2),parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
614 fi;
615 od;
616 return retL;
617 end;
618
619 _Place:=function(P)
620 local f1,f2,Funf,par,x,y,G,o,Id,kx,f,kxy,Pl,pl,i,kx2,kxy2;
621 par := Parent(P);
622 Funf:=_OldFunf(par);
623 x:= P.x;
624 if x= "INFTY POINT" then
625 o:=EquationOrderInfinite(Funf);
626 G:= Decomposition(o);
627 return Place(G[1]);
628 fi;
629 y:=P.y;
630 kx:= PolynomialAlgebra(par.F);
631 kxy:= PolynomialAlgebra(kx);
632 kx2:= BaseField(Funf);
633 kxy2:= PolynomialAlgebra(kx2);
634 f1:= MinimalPolynomial(Coerce(par.F,x));
635 G:= Generator(Funf,1);
636 f2:= G-y;
637 f2:= Coerce(Funf,f2);
638 f1:= Coerce(CoefficientRing(Funf),f1);
639 o:= EquationOrderFinite(Funf);
640 IsMaximal(o);
641 Id:= Ideal(o,[f1,f2]);
642 if (IsPrime(Id)) then
643 pl:= Place(Id);
644 else f:= Factorization(Id);
645 Pl:= [];
646 for i in [1..Length(f)] do
647 Pl[i] := Place(f[i][1]);
648 if Point(Pl[i],par) =P then
649 pl:= Pl[i]; i:= Length(f)+1;
650 fi;
651 od;
652 fi;
653 return pl;
654 end;
655
656
657
658
659 RandomPoint:= function(EC)
660 local P,F,kx,r1,r2,g,lis,eval1,eval2,k,poly,Funf,kxy,finished;
661 Funf := _OldFunf(EC);
662 finished := FALSE;
663 while finished = FALSE do
664 P := RandomPlace(Funf,1);
665 P:= P.ext1;
666 kx := PolynomialAlgebra(EC.F);
667 kxy := PolynomialAlgebra(kx);
668 if IsFinite(P) then finished := TRUE;
669 g:=Generators(Ideal(P));
670 r1:=Roots(Coerce(kx,g[1]))[1][1];
671 lis :=List(g[2]);
672 eval1 := FALSE;
673 eval2 := FALSE;
674 for k in [1..Length(lis)] do
675 if (Degree(Denominator(lis[k])) > 0 and eval1 = FALSE) then
676 eval1 := TRUE;
677 poly := _SwapVars(DefiningPolynomial(Funf));
678
679 elif (Degree(lis[k])>0 and eval2 = FALSE) then
680 eval2 := TRUE;
681 poly := _EvalPlace(r1,g[2]);
682 fi;
683 od;
684 if eval1 = TRUE then
685 r2 :=Roots(Evaluate (poly,r1))[1][1];
686 elif eval2 = TRUE then
687 r2 := Roots(poly)[1][1];
688 else
689 r2 :=Roots(Coerce(kxy,lis))[1][1];
690 fi;
691 return rec(x:=r1,y:=Coerce(CoefficientRing(kx),r2),parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
692 fi;
693 od;
694 end;
695
696 RandomPointWithInf:= function(EC)
697 local P,F,kx,r1,r2,g,lis,eval1,eval2,k,poly,Funf,kxy;
698 Funf := _OldFunf(EC);
699 P := RandomPlace(Funf,1);
700 P:= P.ext1;
701 kx := PolynomialAlgebra(EC.F);
702 kxy := PolynomialAlgebra(kx);
703 if (IsFinite(P) = FALSE) then
704 return rec(x := "INFTY POINT",parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
705 else
706 g:=Generators(Ideal(P));
707 r1:=Roots(Coerce(kx,g[1]))[1][1];
708 lis :=List(g[2]);
709 eval1 := FALSE;
710 eval2 := FALSE;
711 for k in [1..Length(lis)] do
712 if (Degree(Denominator(lis[k])) > 0 and eval1 = FALSE) then
713 eval1 := TRUE;
714 poly := _SwapVars(DefiningPolynomial(Funf));
715
716 elif (Degree(lis[k])>0 and eval2 = FALSE) then
717 eval2 := TRUE;
718 poly := _EvalPlace(r1,g[2]);
719 fi;
720 od;
721 if eval1 = TRUE then
722 r2 :=Roots(Evaluate (poly,r1))[1][1];
723 elif eval2 = TRUE then
724 r2 := Roots(poly)[1][1];
725 else
726 r2 :=Roots(Coerce(kxy,lis))[1][1];
727 fi;
728 return rec(x:=r1,y:=Coerce(CoefficientRing(kx),r2),parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
729 fi;
730 end;
731
732
733
734
735
736
737
738 Subgroup:= function(P)
739 local i,L,stop,R;
740 L:= [];
741 stop := FALSE;
742 i := 1;
743 while stop = FALSE do
744 R := _PointMult(i,P);
745 L[i] := R;
746 if (R.x = "INFTY POINT") then
747 stop := TRUE;
748 fi;
749 i := i+1;
750 od;
751 #if GetVerbose("Elliptic") <> 1 then
752 #Print("Length: ",Length(L),"\n");
753 #fi;
754 return L;
755 end;
756
757 _IndPoints:= function(P1,P2)
758 local L1,L2,i,j;
759 if P1 = P2 then return FALSE; fi;
760 L1 := Subgroup(P1);
761 L2 := Subgroup(P2);
762 for i in [1..Length(L2)] do
763 if L2[i] in L1 then return FALSE;fi;
764 od;
765 return TRUE;
766 end;
767
768
769
770 _DisSub := function(P1,P2)
771 local sub1,sub2,i,j;
772 sub1:= Subgroup(P1);
773 sub2:= Subgroup(P2);
774 for i in [1..Length(sub1)-1] do
775 for j in [1..Length(sub2)-1] do
776 if sub1[i]= sub2[j] then return FALSE;
777 fi;
778 od;
779 od;
780 return TRUE;
781 end;
782
783
784
785
786
787 _Dep:= function(P1,P2)
788 local sub,i;
789 sub:= Subgroup(P2);
790 for i in [1..Length(sub)] do
791 if sub[i] = P1 then return TRUE;
792 fi;
793 od;
794 return FALSE;
795 end;
796
797
798
799
800 __Order2 := function(P)
801 local q,Q,L,par,m,j,k,boo1,boo2,R,l,neg,M,fact,newF,i,found,tried,testP,n,ok;
802 if P.x = "INFTY POINT" then
803 return 1;
804 fi;
805 par:= P.parent;
806 q := Size(par.F);
807 Q := _PointMult(q+1,P);
808 m:= (Ceiling(q^(1/4)))+1;
809 if m > 2^27 then return Error("Curve too large, cannot compute order of point."); fi;
810 L := [];
811 L[1]:= rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops);
812 for j in [1..2*m-1] do
813 L[j+1] := L[j] + P;
814 od;
815 found := FALSE;
816 tried := FALSE;
817 while (found = FALSE) do
818 if tried = FALSE then
819 k := -1;
820 fi;
821 boo1 := TRUE;
822 boo2 := TRUE;
823 neg := FALSE;
824 while (boo1) do
825 k := k+1;
826 R := Q + _PointMult(k*2*m,P);
827 l := 0;
828 while (boo2 and l<Length(L)) do
829 l := l+1;
830 if L[l] = R then
831 boo2 := FALSE;
832 boo1 := FALSE;
833 fi;
834 od;
835 od;
836 if (boo1= TRUE and boo2 = TRUE) then Print("noMatch found\n");fi;
837 j:= l-1;
838 M := q+1+2*m*k-j;
839 newF := TRUE;
840 while (newF = TRUE) do
841 newF := FALSE;
842 if M = 0 then
843 found := FALSE;
844 tried := TRUE;
845 else fact := Factorization(M);
846 testP := (_PointMult(M,P));
847 if testP.x <>"INFTY POINT" then Print("Error!,_Order\n");fi;
848 found := TRUE;
849 i := 0;
850 ok := TRUE;
851 while ok= TRUE and i <Length(fact) do
852 i := i+1;
853 if _PointMult(Floor(M/((fact[i][1]))),P) = rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops) then
854 M := Abs(Floor(M/((fact[i][1]))));
855 newF := TRUE;
856 ok := FALSE;
857 fi;
858 od;
859 fi;
860 od;
861 od;
862 return Abs(M);
863 end;
864
865 __Ordernor:= function(P)
866 local P2,i,j,L1,L2,p,s,Q,R,val,zer,E,found,m;
867 E:= Parent(P);
868 p:= Size(BaseRing(E));
869 s:= Ceiling(p^(1/4));
870 L1:=[P];
871 zer:= Zero(E);
872 if P = zer then
873 return 1;
874 fi;
875 for i in [2..s] do
876 L1[i]:= L1[i-1] +P;
877 if L1[i] = zer then return i;
878 fi;
879 od;
880 Q:= (2*s+1)*P;
881 R:= (p+1)*P;
882 found:=false;
883 j:=0;
884 while not found do
885 if j > 0 then
886 R:= R+Q;
887 fi;
888 for i in [1..Length(L1)] do
889 if R = L1[i] then
890 found:=true;
891 val :=i;
892 i:= Length(L1);
893 elif R = -L1[i] then
894 found:=true;
895 val := -i;
896 i:= Length(L1);
897 fi;
898 od;
899 j:=j+1;
900 od;
901 m:= p+1+(2*s+1)*(j-1)-val;
902 for i in [s..m] do
903 if m mod i = 0 then
904 if i*P = zer then
905 #Print(P,i,"\n");
906 return i;
907 fi;
908 fi;
909 od;
910 #Print(P,m,"\n");
911 return m;
912 end;
913
914
915
916 __Ordernor2 := function(P)
917 local q,Q,L,par,m,j,k,boo1,boo2,R,l,neg,M,fact,newF,i,found,tried,testP,n,ok;
918 if P.x = "INFTY POINT" then
919 return 1;
920 fi;
921 par:= P.parent;
922 q := Size(par.F);
923 Q := _PointMult(q+1,P);
924 m:= (Ceiling(q^(1/4)))+1;
925 if m >=2^28 then return Error("Curve too large, cannot compute order of point.\n"); fi;
926 L := [];
927 L[1]:= rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops);
928 for j in [1..m-1] do
929 L[j+1] := L[j] + P;
930 od;
931 if Characteristic(par.F) =2 then
932 for j in [m..2*m-1] do
933 L[j+1] := L[j] + P;
934 od;
935 fi;
936 found := FALSE;
937 tried := FALSE;
938 while (found = FALSE) do
939 if tried = FALSE then
940 k := -m-1;
941 fi;
942 boo1 := TRUE;
943 boo2 := TRUE;
944 neg := FALSE;
945 while (boo1 ) do
946 k := k+1;
947 R := Q + _PointMult(k*2*m,P);
948 l := 0;
949 while (boo2 and l<Length(L)) do
950 l := l+1;
951 if R = _PointMult(-1,L[l]) then
952 boo2 := FALSE;
953 boo1 := FALSE;
954 neg := TRUE;
955 elif L[l] = R then
956 boo2 := FALSE;
957 boo1 := FALSE;
958 fi;
959
960 od;
961 od;
962 if (boo1= TRUE and boo2 = TRUE) then Print("noMatch found\n");fi;
963 j:= l-1;
964 if (Characteristic(par.F) =2) then
965 M := q+1+Abs(2*m*k);
966 else M:= q+1+2*m*k;
967 fi;
968 if (neg and (Characteristic(par.F) <>2))then M := M+j;
969 else M := M-j;
970 fi;
971
972 newF := TRUE;
973 while (newF = TRUE) do
974 newF := FALSE;
975 if M = 0 then
976 found := FALSE;
977 tried := TRUE;
978 else fact := Factorization(M);
979 testP := (_PointMult(M,P));
980 if testP.x <>"INFTY POINT" then Print("Error!,_Order\n");fi;
981 found := TRUE;
982 i := 0;
983 ok := TRUE;
984 while ok= TRUE and i <Length(fact) do
985 i := i+1;
986 if _PointMult(Floor(M/((fact[i][1]))),P) = rec(x := "INFTY POINT",parent:= P.parent,type := elt-grp^ell,operations := _elt_grp_ell_ops) then
987 M := Abs(Floor(M/((fact[i][1]))));
988 newF := TRUE;
989 ok := FALSE;
990 fi;
991 od;
992 fi;
993 od;
994 od;
995 return Abs(M);
996 end;
997
998 _Order:= function(P)
999 local par;
1000 par := P.parent;
1001 if Characteristic(par.F) =2 then return __Order2(P);
1002 else return __Ordernor(P);
1003 fi;
1004 end;
1005
1006 _Size2:=function(EC)
1007 local found,P,N,q,o,i,k,l,intFound,ok,Poi,j,m,op,oq,fo;
1008 if EC.size<> -1 then return EC.size; fi;
1009 #SetVerbose("Elliptic",1);
1010 found := FALSE;
1011 N:= 1;
1012 q := Size(EC.F);
1013 if q<11 then return Length(Points(EC)); fi;
1014 Poi:= [];
1015 j := 0;
1016 fo :=0;
1017 while found = FALSE do
1018 P := RandomPoint(EC);
1019 o := _Order(P);
1020 k := LCM(o,N);
1021 if k <> N then
1022 Poi[Length(Poi)+1]:= P;
1023 elif k=N then
1024 if q<37 and fo >4 then
1025 return Length(Points(EC)); #verhindert Endlosschleife fuer den Fall dass 2*untere Schranke <=obere Schranke, tritt nur fuer q < 37 auf
1026 fi;
1027 for i in [1..Length(Poi)] do
1028 if(_Order(Poi[i])mod o = 0 or o mod _Order(Poi[i]) = 0 ) then
1029 if _IndPoints(P,Poi[i]) then
1030 k := k*o;
1031 Append(Poi,[P]);
1032 fi;
1033 fi;
1034 od;
1035 fi;
1036 i:=0;
1037 intFound := FALSE;
1038 ok := TRUE;
1039 while ok do
1040 i:= i+1;
1041 l := i*k;
1042 if ((l >= q+1 - Floor(2*(Sqrt(q)))) and (l <= q+1 + 2*(Ceiling(Sqrt(q))))) then
1043 if intFound = FALSE then intFound := TRUE; fo:=fo+1;
1044 else ok := FALSE;N:= LCM(o,N);
1045 fi;
1046 elif (l > q+1 + 2*(Ceiling(Sqrt(q)))) then
1047 if intFound = TRUE then found := TRUE;N:= l-k;
1048 fi;
1049 ok := FALSE;
1050 N:= LCM(o,N);
1051 fi;
1052 od;
1053 j := j+1;
1054 od;
1055 EC.size := N;
1056 return N;
1057 end;
1058
1059
1060 _SizeSmall:= function(E)
1061 local sum,p,i,L,a,b;
1062 p:= Size(BaseRing(E));
1063 L:= E.L;
1064 if Length(L) <> 2 then
1065 if Length(L) = 5 and L[1]=0 and L[2]= 0 and L[3] = 0 then
1066 a:=Coerce(Integers(),L[4]);
1067 b:= Coerce(Integers(),L[5]);
1068 else return Error("missing case");
1069 fi;
1070 else
1071 a:= Coerce(Integers(),L[1]);
1072 b:= Coerce(Integers(),L[2]);
1073 fi;
1074 sum:=0;
1075 for i in [0..p-1] do
1076 sum := sum+ LegendreSymbol(i^3+a*i+b, p);
1077 od;
1078 E.size:= 1+p+sum;
1079 return E.size;
1080 end;
1081
1082 _SizeTest:= function(k)
1083 local i,j,exi,si,E;
1084 for i in [2..k] do
1085 for j in [1..k] do
1086 exi:= ExistsEllipticCurve(k,[i,j]);
1087 if exi then
1088 E:= exi.ext1;
1089 si:= Size(E);
1090 Print("i: ",i," j: ",j," Size: ",si,"\n");
1091 if Length(Points(E)) <> si then
1092 return Error("wrong calculation!");
1093 fi;
1094 fi;
1095 od;
1096 od;
1097 end;
1098
1099 QuadraticTwist:= function(E)
1100 local BR,L,a,b,E2,a1,a3,a2,a4,a6,f,g,x,T,T,found,n,p,r,i,Et;
1101 BR:=BaseRing(E);
1102 L:=Coefficients(E);
1103 p:=Size(BR);
1104 if Characteristic(BR)>3 then
1105 if Length(L) <>2 then
1106 E2:=ShortWeierstrassNormalForm(E);
1107 else E2:=E;
1108 fi;
1109 elif Characteristic(BR) =2 then
1110 if L[2]<> 0 or L[4] <>0 then
1111 a1:=L[1];
1112 a2:=L[3];
1113 a3:=L[2];
1114 a4:=L[4];
1115 a6:=L[5];
1116 f:=_GetPoly(BR,L);
1117 x:=Parent(f).1;
1118 g:=Evaluate(f,a1^3*x+(a1^2*a4+a3^2)/(a1^3));
1119 g:=Evaluate(_SwapVars(g),a1^2*x+a3/a1);
1120 g:=g/(a1^6);
1121 b:=Coerce(BR,Coefficient(g,2));
1122 a:=Coefficient(Coefficient(g,0),0);
1123 E2:=EllipticCurve(BR,[1,0,b,0,a]);
1124 L:=Coefficients(E2);
1125 else
1126 E2:=E;
1127 b:=L[3];
1128 a:=L[5];
1129 fi;
1130 if Trace(b)=0 then
1131 i:=0;
1132 T:=PrimitiveElement(BR);
1133 found:=false;
1134 while not found do
1135 if Trace(T^i) = 1 then
1136 Et:= EllipticCurve(BR,[1,0,T^i,0,a]);
1137 if E.size <> -1 then
1138 Et.size:= 2*p+2 - E.size;
1139 fi;
1140 return Et;
1141 fi;
1142 i:=i+1;
1143 od;
1144 else Et:= EllipticCurve(BR,[1,0,0,0,a]);
1145 if E.size <> -1 then
1146 Et.size:= 2*p+2 - E.size;
1147 fi;
1148 return Et;
1149
1150 fi;
1151 fi;
1152 n:=1;
1153 if Characteristic(BR) >3 then
1154 if IsPrime(p) then
1155 while not LegendreSymbol(n,p) = -1 do
1156 n:= n+1;
1157 od;
1158
1159 else n:= PrimitiveElement(BR);
1160 g:=n;
1161 r:=0;
1162 while not AllRoots(n,2) = Sequence([]) and r <> p do
1163 r:= r+1;
1164 n:= g*n;
1165 od;
1166 fi;
1167 L:= E2.L;
1168 if Length(L) <> 2 then
1169 return Error ("missing case");
1170 fi;
1171 a:= L[1];
1172 b:= L[2];
1173 Et:= EllipticCurve(BR,[(n^2)*a,(n^3)*b]);
1174 if E.size <> -1 then
1175 Et.size:= 2*p+2 - E.size;
1176 fi;
1177 return Et;
1178
1179 fi;
1180 end;
1181
1182 EllipticCurveFromJInvariant:= function(j)
1183 local BR,b,a,b,p,i,s;
1184 BR:=Parent(j);
1185 if Characteristic(BR) =2 then
1186 return EllipticCurve(BR,[1,0,0,0,1/j]);
1187 fi;
1188 p:= PrimitiveElement(BR);
1189 i:=0;
1190 s:=Size(BR);
1191 while i<s do
1192 i:=i+1;
1193 a:= p^i;
1194 b:= AllRoots((a^3*(48^3-64*j)/(27*16*j)),2);
1195 if not Length(b)=0 then
1196 return EllipticCurve(BR,[a,b[1]]);
1197 fi;
1198 od;
1199 return Error("no curve with j-invariant ",j," over ",BR," found");
1200 end;
1201
1202
1203
1204 _Size3:= function(EC)
1205 local sum,i,p,EC2,n,found,P,N2,P2,N,q,o,o2,k2,intFound2,l2,ok2,i,k,l,intFound,ok,Poi1,j,m,op,oq,fo,ords1,time,upBo,loBo,Poi2,ords2,BR,g,r,j,k,pri,EC,tri,_wei;
1206 #time:= Realtime();
1207 if EC.size<> -1 then return EC.size; fi;
1208 #SetVerbose("Elliptic",1);
1209 if Length(Coefficients(EC))= 5 and Characteristic(BaseRing(EC))>3 then
1210 _wei:=_WeierstrassNormalForm(EC.F,EC.L);
1211 EC.size:= _Size3(EllipticCurve(_wei[1],_wei[2]));
1212 return EC.size;
1213 fi;
1214 BR:= BaseRing(EC);
1215 tri:=Characteristic(BR)=3;
1216 p:= Size(BR);
1217 if IsPrime(p) then
1218 if p <230 then
1219 return _SizeSmall(EC);
1220 fi;
1221 fi;
1222 #else
1223 EC2:=QuadraticTwist(EC);
1224 found := 0;
1225 N:= 1;
1226 if not tri then
1227 N2:=1;
1228 fi;
1229 q := Size(EC.F);
1230 Poi1:= [];
1231 ords1:=[];
1232 if not tri then
1233 Poi2:=[];
1234 ords2:=[];
1235 fi;
1236 j := 0;
1237 upBo:=q+1 + Ceiling(2*(Sqrt(q)));
1238 loBo:= q+1 - Floor(2*(Sqrt(q)));
1239 while found = 0 do
1240 #Print(Length(Poi),"\n");
1241 P := RandomPoint(EC);
1242 if not tri then
1243 P2:= RandomPoint(EC2);
1244 fi;
1245 while P in Poi1 do
1246 P := RandomPoint(EC);
1247 od;
1248 if not tri then
1249 while P2 in Poi2 do
1250 P2 := RandomPoint(EC2);
1251 od;
1252 fi;
1253 o := _Order(P);
1254 k := LCM(o,N);
1255 if not tri then
1256 o2 := _Order(P2);
1257 k2 := LCM(o2,N2);
1258 fi;
1259 #Print("lcm, k: ",k,"\n");
1260 if k <> N then
1261 Poi1[Length(Poi1)+1]:= P;
1262 #Print(P,"\n");
1263 ords1:=Append(ords1,[o]);
1264 elif k=N then
1265 i:=0;
1266 while i <Length(Poi1) do
1267 i:=i+1;
1268 if(ords1[i] mod o = 0 or o mod ords1[i] = 0 ) then
1269 if not _IndPoints(P,Poi1[i]) then
1270 i:= Length(Poi1)+2;
1271 fi;
1272 fi;
1273 od;
1274 if i = Length(Poi1) then
1275 k := k*o;
1276 #Print("indpoi, k: ",k,"\n");
1277 #Print(P,"\n");
1278 Poi1:=Append(Poi1,[P]);
1279 ords1:= Append(ords1,[o]);
1280 fi;
1281 fi;
1282 if not tri then
1283 if k2 <> N2 then
1284 Poi2[Length(Poi2)+1]:= P2;
1285 #Print(P,"\n");
1286 ords2:=Append(ords2,[o2]);
1287 elif k2=N2 then
1288 i:=0;
1289 while i <Length(Poi2) do
1290 i:=i+1;
1291 if(ords2[i] mod o2 = 0 or o2 mod ords2[i] = 0 ) then
1292 if not _IndPoints(P2,Poi2[i]) then
1293 i:= Length(Poi2)+2;
1294 fi;
1295 fi;
1296 od;
1297 if i = Length(Poi2) then
1298 k2 := k2*o2;
1299 #Print("indpoi, k: ",k,"\n");
1300 #Print(P,"\n");
1301 Poi2:=Append(Poi2,[P2]);
1302 ords2:= Append(ords2,[o2]);
1303 fi;
1304 fi;
1305 fi;
1306 i:=0;
1307 intFound := FALSE;
1308 intFound2:= FALSE;
1309 ok := TRUE;
1310 ok2:=true;
1311 while ok and ok2 do
1312 i:= i+1;
1313 l := i*k;
1314 if not tri then
1315 l2:= i*k2;
1316 fi;
1317 if ((l >= loBo) and (l <= upBo)) then
1318 if intFound = FALSE then intFound := TRUE;
1319 else ok := FALSE;N:= LCM(o,N);
1320 fi;
1321 elif (l > upBo) then
1322 if intFound = TRUE then found := 1;N:= l-k;
1323 fi;
1324 ok := FALSE;
1325 N:= LCM(o,N);
1326 fi;
1327 if not tri then
1328 if ((l2 >= loBo) and (l2 <= upBo)) then
1329 if intFound2 = FALSE then intFound2 := TRUE;
1330 else ok2 := FALSE;N2:= LCM(o2,N2);
1331 fi;
1332 elif (l2 > upBo) then
1333 if intFound2 = TRUE then found := 2;N2:= l2-k2;
1334 fi;
1335 ok2 := FALSE;
1336 N2:= LCM(o2,N2);
1337 fi;
1338 fi;
1339 od;
1340 j := j+1;
1341 od;
1342 if found = 1 then
1343 EC.size := N;
1344 #Print(Poi1,"\n");
1345 #Print(Poi2,"\n");
1346 #Print("Time: ",Realtime()-time);
1347 return N;
1348 elif found = 2 then
1349 N:=2*(p+1)-N2;
1350 EC.size:= N;
1351 #Print(Poi1,"\n");
1352 #Print(Poi2,"\n");
1353 #Print("Time: ",Realtime()-time);
1354 return N;
1355 fi;
1356 #fi;
1357 #else return Error("missing case");
1358 #fi;
1359 end;
1360
1361
1362 _Size:=function(EC)
1363 local found,P,N,q,o,i,k,l,intFound,ok,Poi,j,m,op,oq,fo,ords,time,upBo,loBo;
1364 time:=Realtime();
1365 if EC.size<> -1 then return EC.size; fi;
1366 #SetVerbose("Elliptic",1);
1367 found := FALSE;
1368 N:= 1;
1369 q := Size(EC.F);
1370 if q<11 then return Length(Points(EC)); fi;
1371 Poi:= [];
1372 ords:=[];
1373 j := 0;
1374 fo :=0;
1375 upBo:=q+1 + Ceiling(2*(Sqrt(q)));
1376 loBo:= q+1 - Floor(2*(Sqrt(q)));
1377 while found = FALSE do
1378 #Print(Length(Poi),"\n");
1379 P := RandomPoint(EC);
1380 while P in Poi do
1381 P := RandomPoint(EC);
1382 od;
1383 o := _Order(P);
1384 k := LCM(o,N);
1385 #Print("lcm, k: ",k,"\n");
1386 if k <> N then
1387 Poi[Length(Poi)+1]:= P;
1388 Print(P,"\n");
1389 ords:=Append(ords,[o]);
1390 elif k=N then
1391 if q<37 and fo >4 then
1392 return Length(Points(EC)); #verhindert Endlosschleife fuer den Fall dass 2*untere Schranke <=obere Schranke, tritt nur fuer q < 37 auf
1393 fi;
1394 i:=0;
1395 while i <Length(Poi) do
1396 i:=i+1;
1397 if(ords[i] mod o = 0 or o mod ords[i] = 0 ) then
1398 if not _IndPoints(P,Poi[i]) then
1399 i:= Length(Poi)+2;
1400 fi;
1401 fi;
1402 od;
1403 if i = Length(Poi) then
1404 k := k*o;
1405 #Print("indpoi, k: ",k,"\n");
1406 Print(P,"\n");
1407 Poi:=Append(Poi,[P]);
1408 ords:= Append(ords,[o]);
1409 fi;
1410 fi;
1411 i:=0;
1412 intFound := FALSE;
1413 ok := TRUE;
1414 while ok do
1415 i:= i+1;
1416 l := i*k;
1417 if ((l >= loBo) and (l <= upBo)) then
1418 if intFound = FALSE then intFound := TRUE; fo:=fo+1;
1419 else ok := FALSE;N:= LCM(o,N);
1420 fi;
1421 elif (l > upBo) then
1422 if intFound = TRUE then found := TRUE;N:= l-k;
1423 fi;
1424 ok := FALSE;
1425 N:= LCM(o,N);
1426 fi;
1427 od;
1428 j := j+1;
1429 od;
1430 EC.size := N;
1431 Print("time: ",Realtime()-time,"\n");
1432 return N;
1433 end;
1434
1435 #find k with k*A = B
1436 DiscreteLog:= function(A,B)
1437 local N,m,L,i,j,R,found;
1438
1439 if A.parent <>B.parent then
1440 Print("Points not defined over the same curve.\n");
1441 return;
1442 fi;
1443 N := _Size(A.parent);
1444 m := Ceiling(Sqrt(N));
1445 L := [];
1446 for i in [1..m+1] do
1447 L[i] := _PointMult(i-1,A);
1448 od;
1449 for j in [1..m+1] do
1450 R := B - _PointMult((j-1)*m,A);
1451 for i in [1..m+1] do
1452 if R = L[i] then
1453 return ((i-1) + (j-1)*m) mod _Order(A);
1454 fi;
1455 od;
1456 od;
1457 Print("There is no k with k*",A,"=",B,".\n");
1458 end;
1459
1460
1461
1462
1463
1464 _ECSubEq:= function(L1,L2)
1465 local i,j,found;
1466 if Length(L1)<>Length(L2) then
1467 return FALSE;
1468 fi;
1469 for i in [1..Length(L1)] do
1470 found := FALSE;
1471 j := 0;
1472 while (found =FALSE and j<Length(L2)) do
1473 j:= j+1;
1474 if L2[j]= L1[i] then
1475 found := TRUE;
1476 fi;
1477 od;
1478 if found = FALSE then
1479 Print(L1[i],"\n");
1480 return FALSE;
1481 fi;
1482 od;
1483 return TRUE;
1484 end;
1485
1486
1487 _DisSub2:= function(L1,L2)
1488 local i,j;
1489 for i in [1..Length(L1)-1] do
1490 for j in [1..Length(L2)-1] do
1491 if L1[i] = L2[j] then return FALSE;
1492 fi;
1493 od;
1494 od;
1495 return TRUE;
1496 end;
1497
1498 _Generator:= function(EC)
1499 local N,Points,P,o1,o,sub,R,op;
1500 N:= _Size(EC);
1501 P:= RandomPoint(EC);
1502 Points:= [P];
1503 o1:= _Order(P);
1504 if o1=N then return Points;
1505 else sub:= Subgroup(Points[1]);
1506 while o1<N do
1507 P:= RandomPoint(EC);
1508 op:= _Order(P);
1509 if op = N/o1 then
1510 if _DisSub2(sub,Subgroup(P)) then
1511 Points[2] := P;
1512 o1 := N;
1513 fi;
1514 elif op = N then Points[1]:= P; return Points;
1515 elif op>o1 then Points[1]:= P; o1:=op;sub:= Subgroup(P);
1516 # else P:=Points[1]+P; op:= _Order(P);
1517 # if op>o1 then Points[1]:= P; o1:=op;sub:= Subgroup(P);fi;
1518 fi;
1519 od;
1520 if Length(Points)=1 then return Points;
1521 else
1522 R:=Points[1]+Points[2];
1523 if _Order(R) = N then
1524 return [R];
1525 else return Points;
1526 fi;
1527 fi;
1528 fi;
1529 end;
1530
1531
1532 _HasSingularPoints := function(F,L)
1533 local x,y,f,f1,f2,g1,g2,L1,L2,i,j,Li,fx,fy,kx;
1534 if Coerce(F,_Discriminant(F,L))<>0 then
1535 return FALSE;
1536 elif Characteristic(F) <>2 then
1537 f:= _GetPoly(F,L);
1538 f1:= Evaluate(f,0);
1539 Li:= Roots(f1);
1540 if Length(Li)=3 then
1541 return FALSE;
1542 else
1543 for i in [1..Length(Li)] do
1544 if Li[i][2] >1
1545 then Print("Singular point found: [",Li[i][1],",0]\n");
1546 return TRUE;
1547 fi;
1548 od;
1549 fi;
1550 return Error("disc =0 but found no singular point\n");
1551 else
1552 kx:= PolynomialAlgebra(F);
1553 f:= _GetPoly(F,L);
1554 fy:= Derivative(f);
1555 fx:= Derivative(_SwapVars(f));
1556 L:= ElementToSequence(fy);
1557 fy:= Coerce(kx,L[1]);
1558 L1:= Roots(fy);
1559 fx:= Evaluate(fx,L1[1][1]);
1560 L2:= Roots(fx);
1561 y:= L2[1][1];
1562 x:= L1[1][1];
1563 if(Evaluate(Evaluate(f,y),x)) = 0 then
1564 Print("Singular point found: [",x,",",y,"]\n");
1565 return TRUE;
1566 fi;
1567
1568 fi;
1569 return;
1570 end;
1571
1572 _Divisor :=function(P)
1573 local Pl;
1574 Pl:=_Place(P);
1575 return Divisor(Pl);
1576 end;
1577
1578 _GroundPlace:= function(P)
1579 local base;
1580 base:= Parent(P);
1581 while not Is(Type(base),fld^fin) do
1582 base:= BaseRing(base);
1583 od;
1584 return base;
1585 end;
1586
1587 MakePoint:= function(EC,L)
1588 local f,p,x,y,P,d,Pl;
1589 if L=INFTY or L[1] = INFTY then
1590 return rec(x := "INFTY POINT",parent:= EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
1591 fi;
1592 x:=Coerce(EC.F,L[1]);
1593 y:=Coerce(EC.F,L[2]);
1594 p:= _GetPoly(EC.F,EC.L);
1595 f:= Evaluate(p,y);
1596 if(Evaluate(f,x) <>0 ) then
1597 return Error([x,y],": Given coordinates do not define a point on given curve\n");
1598 fi;
1599 P:=rec(x:=x,y:=y,parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
1600 if Is(Type(EC.F),fld^fin) then
1601 Pl:= _Place(P);
1602 d:=Degree(Pl);
1603 if (d <>1) then
1604 return Error("Given coordinates do not define a point on given curve\n");
1605 fi;
1606 fi;
1607 return P;
1608 end;
1609
1610
1611 _CoercePoint:= function(E,P)
1612 if IsList(P) then return MakePoint(E,P);fi;
1613 return MakePoint(E,[PX(P),PY(P)]);
1614 end;
1615
1616
1617 MakePoints:= function(EC,x)
1618 local p,f,R,P,Pl,d,i;
1619 if x=INFTY then
1620 return rec(x := "INFTY POINT",parent:= EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
1621 fi;
1622 x:= Coerce(EC.F,x);
1623 p:= _SwapVars(_GetPoly(EC.F,EC.L));
1624 f:= Evaluate(p,x);
1625 R:= Roots(f);
1626 P:= [];
1627 for i in [1..Length(R)] do
1628 P[i]:=rec(x:=x,y:=R[i][1],parent := EC,type := elt-grp^ell,operations := _elt_grp_ell_ops);
1629 Pl:= _Place(P[i]);
1630 d:=Degree(Pl);
1631 if (d <>1) then
1632 Remove_(P,i);
1633 fi;
1634 od;
1635 return P;
1636 end;
1637
1638 _Coeffs:=function(E)
1639 return Sequence(E.L);
1640 end;
1641
1642
1643 _BaseRing:= function(EC)
1644 return EC.F;
1645 end;
1646
1647 InstallDocumentation(
1648 rec(
1649 kind := "FUNCTION",
1650 name := "EllipticCurve",
1651 sin := [[elt-ord^rat,"q"],[list,"L"]],
1652 sou := [[grp^ell]],
1653 short :=
1654 "Constructs an elliptic curve over the finite field with q elements. If L has length 2, the curve is "+
1655 "defined by the equation y^2 = x^3 + L[1]*x + L[2], if L has length 5, the curve is "+
1656 "defined by the equation y^2 + L[1]*x*y + L[2]*y = x^3 + L[3]*x^2 + L[4]*x + L[5]. ",
1657 ex:=["EC:= EllipticCurve(5,[1,2]);\n"]
1658 ));
1659
1660 InstallDocumentation(
1661 rec(
1662 kind := "FUNCTION",
1663 name := "Points",
1664 sin := [[grp^ell,"E"]],
1665 sou := [[list,"L"]],
1666 short :="Returns all Points of E as a list.",
1667 ex:=["EC:= EllipticCurve(5,[1,2]);;\nPoints(EC);\n"]
1668 ));
1669
1670
1671 InstallDocumentation(
1672 rec(
1673 kind := "FUNCTION",
1674 name := "MakePoint",
1675 sin := [[grp^ell,"EC"],[list,"L"]],
1676 sou := [[elt-grp^ell,"P"]],
1677 short :="Returns,if possible, a point on EC, defined by the coordinates given in L.",
1678 ex:=["EC:= EllipticCurve(5,[1,2]);;\nMakePoint(EC,[4,0]);\n"]
1679 ));
1680
1681 InstallDocumentation(
1682 rec(
1683 kind := "FUNCTION",
1684 name := "MakePoints",
1685 sin := [[grp^ell,"EC"],[elt-fld^fin,"x"]],
1686 sou := [[elt-grp^ell,"P"]],
1687 short :="Returns the points on EC with x-coordinate x.",
1688 ex:=["EC:= EllipticCurve(5,[1,2]);;\nMakePoints(EC,1);\n"]
1689 ));
1690
1691 InstallDocumentation(
1692 rec(
1693 kind := "FUNCTION",
1694 name := "RandomPoint",
1695 sin := [[grp^ell,"EC"]],
1696 sou := [[elt-grp^ell,"P"]],
1697 short :="Returns a random point of EC, but never the point at infinity.",
1698 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nRandomPoint(EC);\n"]
1699 ));
1700
1701 InstallDocumentation(
1702 rec(
1703 kind := "FUNCTION",
1704 name := "RandomPointWithInf",
1705 sin := [[grp^ell,"EC"]],
1706 sou := [[elt-grp^ell,"P"]],
1707 short :="Returns a random point of EC.",
1708 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nRandomPointWithInf(EC);\n"]
1709 ));
1710
1711
1712 InstallDocumentation(
1713 rec(
1714 kind := "FUNCTION",
1715 name := "Point",
1716 sin := [[elt-pls/fld^fun,"p"],[elt-grp^ell,"EC"]],
1717 sou := [[elt-grp^ell,"P"]],
1718 short :="Returns the point on EC corresponding to p.",
1719 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nRandomPoint(EC);\n"]
1720 ));
1721
1722 InstallDocumentation(
1723 rec(
1724 kind := "FUNCTION",
1725 name := "EllipticFunctionField",
1726 sin := [[grp^ell,"EC"]],
1727 sou := [[fld^fun,"F"]],
1728 short :="Returns the function field of genus one corresponding to EC.",
1729 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nRandomPoint(EC);\n"]
1730 ));
1731
1732 InstallDocumentation(
1733 rec(
1734 kind := "FUNCTION",
1735 name := "ShortWeierstrassNormalForm",
1736 sin := [[grp^ell,"EC"]],
1737 sou := [[grp^ell,"EC2"]],
1738 short :="Returns the isomorphic curve EC2 defined by a polynomial of the form y^2=x^3+a*x+b.",
1739 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nRandomPoint(EC);\n"]
1740 ));
1741
1742 InstallDocumentation(
1743 rec(
1744 kind := "FUNCTION",
1745 name := "jInvariant",
1746 sin := [[elt-grp^ell,"EC"]],
1747 sou := [[elt-fld^fin,"j"]],
1748 short :="Returns the j-invariant of EC.",
1749 ex:= ["EC:= EllipticCurve(5,[1,2]);;\njInvariant(EC);\n"]
1750 ));
1751
1752 InstallDocumentation(
1753 rec(
1754 kind := "FUNCTION",
1755 name := "Subgroup",
1756 sin := [[elt-grp^ell,"EC"]],
1757 sou := [[list,"L"]],
1758 short :="Returns the subgroup generated by EC.",
1759 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nP:= RandomPoint(EC);;\n"+"Subgroup(P);\n"]
1760 ));
1761
1762 InstallDocumentation(
1763 rec(
1764 kind := "FUNCTION",
1765 name := "DiscreteLog",
1766 sin := [[elt-grp^ell,"P"],[elt-grp^ell,"Q"]],
1767 sou := [[elt-ord^rat,"n"]],
1768 short :="Returns, if possible, n satisfying n*P = Q.",
1769 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nP:= Generator(EC)[1];;\n"+"Q:= RandomPoint(EC);;\nDiscreteLog(P,Q);\n"]
1770 ));
1771
1772 InstallMethod(
1773 rec(
1774 kind := "FUNCTION",
1775 name := "Discriminant",
1776 sin := [[grp^ell,"EC"]],
1777 sou := [[elt-fld^fin,"d"]],
1778 short :="The discriminant of EC.",
1779 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nDiscriminant(EC);\n"]
1780 ),_Discriminant2);
1781
1782
1783 InstallMethod(
1784 rec(
1785 kind := "FUNCTION",
1786 name := "Divisor",
1787 sin := [[elt-grp^ell,"P"]],
1788 sou := [[elt-dvs/fld^fun,"D"]],
1789 short :="The divisor corresponding to P.",
1790 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nP:= RandomPoint(EC);;\nDivisor(P);"]
1791 ),_Divisor);
1792
1793 InstallMethod(
1794 rec(
1795 kind := "FUNCTION",
1796 name := "BaseRing",
1797 sin := [[grp^ell,"EC"]],
1798 sou := [[fld^fin,"F"]],
1799 short :="The field over which EC is defined.",
1800 ex:= ["EC:= EllipticCurve(5,[1,2]);;\nF:= BaseRing(EC);;\n"]
1801 ),_BaseRing);
1802
1803
1804 InstallMethod(
1805 rec(
1806 kind:="FUNCTION",
1807 name:="Order",
1808 sin:=[[elt-grp^ell,"P"]],
1809 sou:=[[elt-ord^rat,"n"]],
1810 short:=
1811 "The order of the point P, only works for \"small\" curves. ",
1812 ex:=["EC:= EllipticCurve(5,[1,2]);;\nP:= RandomPoint(EC);;\n"+"Order(P);\n"
1813 ],
1814 see:=[]), _Order);
1815
1816
1817 InstallMethod(
1818 rec(
1819 kind:="FUNCTION",
1820 name:="Generator",
1821 sin:=[[grp^ell,"EC"]],
1822 sou:=[[list,"L"]],
1823 short:=
1824 "Returns a list of generator(s) of EC, only works for \"small\" curves. ",
1825 ex:=["EC:= EllipticCurve(5,[1,2]);;\nL:= Generator(EC);\n"
1826 ],
1827 see:=[]), _Generator);
1828
1829 InstallMethod(
1830 rec(
1831 kind:="FUNCTION",
1832 name:="Place",
1833 sin:=[[elt-grp^ell,"P"]],
1834 sou:=[[elt-pls/fld^fun,"p"]],
1835 short:=
1836 "The place of degree one corresponding to P ",
1837 ex:=["EC:= EllipticCurve(5,[1,2]);;\nP:= RandomPoint(EC);;\n"+"Place(P);\n"
1838 ],
1839 see:=[]), _Place);
1840
1841 InstallMethod(
1842 rec(
1843 kind:="FUNCTION",
1844 name:="Size",
1845 sin:=[[grp^ell,"P"]],
1846 sou:=[[elt-ord^rat,"n"]],
1847 short:=
1848 "The number of points belonging to EC, only works for \"small\" curves. ",
1849 ex:=["EC:= EllipticCurve(5,[1,2]);;\nSize(EC);\n"
1850 ],
1851 see:=[]), _Size3);
1852 InstallMethod(
1853 rec(
1854 kind:="FUNCTION",
1855 name:="Coerce",
1856 sin:=[[grp^ell,"EC"],[elt-grp^ell,"P"]],
1857 sou:=[[elt-grp^ell,"P2"]],
1858 short:=
1859 "Coerces P into EC. ",
1860 ex:=["EC:= EllipticCurve(5,[1,2]);;\nP:= RandomPoint(EC);\n"+
1861 "Coerce(EllipticCurve(5^2,[1,2]),P);" ],
1862 see:=[]), _CoercePoint);
1863
1864 InstallMethod(
1865 rec(
1866 kind:="FUNCTION",
1867 name:="Coerce",
1868 sin:=[[grp^ell,"EC"],[list,"L"]],
1869 sou:=[[elt-grp^ell,"P"]],
1870 short:=
1871 "Coerces L into EC. ",
1872 ex:=["EC:= EllipticCurve(5,[1,2]);;\nP:= RandomPoint(EC);\n"+
1873 "Coerce(EllipticCurve(5^2,[1,2]),P);" ],
1874 see:=[]), _CoercePoint);
1875
1876 InstallMethod(
1877 rec(
1878 kind:="FUNCTION",
1879 name:="Zero",
1880 sin:=[[grp^ell,"EC"]],
1881 sou:=[[elt-grp^ell,"P"]],
1882 short:=
1883 "Returns the point at infinity. ",
1884 ex:=["EC:= EllipticCurve(5,[1,2]);;\nP:= Zero(EC);\n"
1885 ],
1886 see:=[]), _Zero);
1887
1888 InstallMethod(
1889 rec(
1890 kind:="FUNCTION",
1891 name:="Coefficients",
1892 sin:=[[grp^ell,"EC"]],
1893 sou:=[[seq(),"L"]],
1894 short:=
1895 "The coefficients of EC as a sequence (of Length 2 resp. 5). ",
1896 ex:=["EC:= EllipticCurve(5,[1,2]);;\nL:= Coefficients(EC);\n"
1897 ],
1898 see:=[]), _Coeffs);
1899
1900 InstallMethod(
1901 rec(
1902 kind:="FUNCTION",
1903 name:="Parent",
1904 sin:=[[elt-grp^ell,"P"]],
1905 sou:=[[grp^ell,"EC"]],
1906 short:=
1907 "The default parent of P ",
1908 ex:=["EC:= EllipticCurve(5,[1,2]);;\nP:= RandomPoint(EC);;\n"+"Parent(P)=EC;\n"
1909 ],
1910 see:=[]), _Parent);