"Fossies" - the Fresh Open Source Software Archive 
"Fossies" - the Fresh Open Source Software Archive A hint: This file contains one or more very long lines, so maybe it is better readable using the pure text view mode that shows the contents as wrapped lines within the browser window.
1 #debugging function 2 Assert:= function(b,s) 3 if b then 4 return; 5 else return Error("Assumption: ",s,"failed!"); 6 fi; 7 end; 8 9 #some functions that are forwarded 10 MultMap:=function(a) 11 return a; 12 end; 13 14 IsogenyFromSubgroup:=function(a) 15 return a; 16 end; 17 18 extFrob:=function(a) 19 return a; 20 end; 21 22 IsogenyFromKernel:=function(a) 23 return a; 24 end; 25 26 IdealFactorization:=function(a) 27 return a; 28 end; 29 30 MakeSmooth2:=function(a) 31 return a; 32 end; 33 34 PicardNumber:=function(a) 35 return a; 36 end; 37 38 PseudoRandom:=function(a) 39 return a; 40 end; 41 42 FindIso:=function(a) 43 return a; 44 end; 45 46 #true if the isogeny f is the zero-map 47 IsZero:=function(f) 48 return "zero" in RecFields(f); 49 end; 50 51 52 #isogeny-operations 53 _iso_ops:=rec(); 54 _iso_ops.Print:= function(f) 55 if IsZero(f) then 56 if f.domain = f.codomain then 57 Print(" ZeroMap on ",f.domain,"\n"); 58 else Print("ZeroMap from ",f.domain, "to ",f.codomain,"\n"); 59 fi; 60 else 61 if "fro" in RecFields(f) then 62 Print("Isogeny from ",f.domain," to ",f.codomain,"(x,y)|->(x^",Degree(f),", y^",Degree(f),")\n"); 63 return; 64 fi; 65 AssignNames_(CoefficientRing(Parent(f.xpol)),["x"]); 66 AssignNames_(Parent(f.ypol),["y"]); 67 Print("Isogeny from ",f.domain," to ",f.codomain,"(x,y)|->(",f.xpol,","); 68 Print(f.ypol,")\n"); 69 fi; 70 end; 71 72 73 #Frobenius-isogeny of degree q 74 FrobeniusToOtherCurve:=function(E,q) 75 local E2,j,L,BR,L2,i,_frob,frob,x,y,j2,Funf; 76 j:=jInvariant(E); 77 L:=Coefficients(E); 78 BR:=BaseRing(E); 79 x:=PolynomialAlgebra(BR).1; 80 Funf:=EllipticFunctionField(E); 81 y:=Funf.1; 82 if GCD(q,Characteristic(BR)) = 1 then 83 return Error; 84 fi; 85 L2:=[]; 86 for i in [1..Length(L)] do 87 L2[i]:=L[i]^q; 88 od; 89 E2:=EllipticCurve(BR,L2); 90 _frob:=function(P) 91 if P = Zero(Parent(P)) then 92 return Zero(E2); 93 else return MakePoint(E2,[PX(P)^q, PY(P)^q]); 94 fi; 95 end; 96 frob:=Map(E,E2,_frob); 97 frob.type:= iso; 98 frob.fro:=true; 99 frob.operations:= _iso_ops; 100 frob.degree:=q; 101 frob.xpol:=Coerce(Funf,x^q); 102 frob.ypol:=y^q; 103 return frob; 104 end; 105 106 107 # computes invariants a_i and b_i 108 ECInvariants:=function(E) 109 local L,Inv; 110 L:=Coefficients(E); 111 Inv:=rec(); 112 if Length(L) = 5 then 113 Inv.a1:=L[1]; 114 Inv.a3:=L[2]; 115 Inv.a2:=L[3]; 116 Inv.a4:=L[4]; 117 Inv.a6:=L[5]; 118 else 119 Inv.a1:=0; 120 Inv.a3:=0; 121 Inv.a2:=0; 122 Inv.a4:=L[1]; 123 Inv.a6:=L[2]; 124 fi; 125 Inv.b2:= Inv.a1^2+4*Inv.a2; 126 Inv.b4:=2*Inv.a4+Inv.a1*Inv.a3; 127 Inv.b6:=Inv.a3^2+4*Inv.a6; 128 Inv.b8:=Inv.a1^2*Inv.a6+4*Inv.a2*Inv.a6-Inv.a1*Inv.a3*Inv.a4+Inv.a2*Inv.a3^2-Inv.a4^2; 129 return Inv; 130 end; 131 132 133 # the generic point (x,y) of a curve 134 GenericPoint:= function(E) 135 local i,Funf,y0,x0,EE,x,y,coLi; 136 Funf:=_OldFunf(E); 137 coLi:=[]; 138 for i in [1..Length(E.L)] do 139 coLi[i]:=E.L[i]; 140 od; 141 EE:= EllipticCurve(Funf,coLi); 142 y:= BaseRing(EE).1; 143 x:= BaseField(BaseRing(EE)).1; 144 return MakePoint(EE,[x,y]); 145 end; 146 147 148 # computes an isomorphic curve in short WeierstrassNormalForm 149 ShortWNF:=function(E) 150 local L,Inv,E2,Iso,xmap,ymap,P,PA,PA2,x,y,f,ff,co; 151 L:=Coefficients(E); 152 if Length(L) = 2 then 153 return rec(base:=E,ext1:=Isomorphism(E,E,1)); 154 else Inv:=ECInvariants(E); 155 E2:=_WeierstrassNormalForm(BaseRing(E),L); 156 E2:=EllipticCurve(BaseRing(E),E2[2]); 157 xmap:=function(P) 158 return 36*PX(P)+3*Inv.b2; 159 end; 160 ymap:=function(P) 161 return (2*PY(P)+Inv.a1*PX(P)+Inv.a3)*108; 162 end; 163 Iso:=Map(E,E2, 164 function(P) 165 if P = Zero(E) then 166 return Zero(E2); 167 fi; 168 return Coerce(E2,[xmap(P),ymap(P)]); 169 end 170 ); 171 P:=GenericPoint(E); 172 Iso.type:= iso; 173 Iso.operations:= _iso_ops; 174 Iso.xpol:= xmap(P); 175 Iso.ypol:= ymap(P); 176 Iso.degree:= 1; 177 return rec(base:=E2,ext1:=Iso); 178 fi; 179 end; 180 181 182 # returns -f if f is an isogeny 183 Negate:=function(f) 184 local fmap,Iso,E2,a1,a3,L; 185 fmap:=function(P) 186 return -f(P); 187 end; 188 E2:=Codomain(f); 189 L:=Coefficients(E2); 190 if Length(L) = 2 then 191 a1:=0; 192 a3:=0; 193 else 194 a1:=L[1]; 195 a3:=L[2]; 196 fi; 197 Iso:=Map(Domain(f),E2,fmap); 198 Iso.type:=iso; 199 Iso.operations:=_iso_ops; 200 Iso.xpol:=f.xpol; 201 Iso.ypol:=f.ypol-a1*f.xpol-a3; 202 Iso.degree:=f.degree; 203 return Iso; 204 end; 205 206 207 208 # returns true if E is of the form y^2+xy=x^3+bx^2+a over GF(2,n) and trace(b) =1 or b =0 209 210 InNormalFormChar2:=function(E) 211 local BR,L,b1,b2,b; 212 BR:=BaseRing(E); 213 L:=Coefficients(E); 214 if not Characteristic(BR) =2 then 215 return Error("function only applicable for curves over fields of genus 2"); 216 fi; 217 b1:= (L[1]=1 and L[2] = 0 and L[4] =0); 218 b:=L[3]; 219 if Is(Type(Parent(b)),fld^fun) then 220 b2:= (Trace(Coerce(CoefficientRing(CoefficientRing(Parent(b))),b)) <> 0 ) or (L[3] = 0); 221 else 222 b2:=(Trace(b) <>0) or (b =0); 223 fi; 224 return b1 and b2; 225 end; 226 227 228 #isomorphic curve of the form y^2+xy=x^3+bx^2+a over GF(2,r), n defines direction of isomorphism 229 230 NormalFormChar2:=function(E,n) 231 local L,f,BR,x,a1,a2,a3,a4,a6,g,b,a,s,E2,Iso,_Iso,xmap,ymap,P,EE,b2; 232 if InNormalFormChar2(E) then 233 return rec(base:=E,ext1:=Isomorphism(E,E,1)); 234 fi; 235 if n <> 1 and n <> 2 then 236 return Error(" n must be in [1,2] "); 237 fi; 238 L:=Coefficients(E); 239 BR:=BaseRing(E); 240 f:=_GetPoly(BR,L); 241 x:=Parent(f).1; 242 a1:=L[1]; 243 a3:=L[2]; 244 a2:=L[3]; 245 a4:=L[4]; 246 a6:=L[5]; 247 b:=((a1^3)*a3+(a1^4)*a2)/(a1^6); 248 s:=0; 249 if Trace(b) = 0 then 250 s:=Roots(PolynomialAlgebra(BR).1^2+PolynomialAlgebra(BR).1+b)[1][1]; 251 b:=0; 252 fi; 253 a:=(a1^4*a4^2+a3^4+a1^3*a3^3+a1^4*a2*a3^2+a1^5*a3*a4+a1^6*a6)/(a1^12); 254 E2:=EllipticCurve(BR,[1,0,b,0,a]); 255 E2.size:=E.size; 256 if n = 1 then 257 P:=GenericPoint(E); 258 xmap:= function(P) 259 if PX(P) = INFTY then 260 return INFTY; 261 fi; 262 return (PX(P)+a3/a1)/a1^2; 263 end; 264 ymap:= function(P) 265 if PY(P) = INFTY then 266 return INFTY; 267 fi; 268 return (PY(P)+(a1^2*a4+a3^2)/a1^3+s*xmap(P))/a1^3; 269 end; 270 _Iso:= function(P) 271 if P = Zero(E) then 272 return Zero(E2); 273 else return Coerce(E2,[xmap(P),ymap(P)]); 274 fi; 275 end; 276 Iso:=Map(E,E2,_Iso); 277 else 278 P:=GenericPoint(E2); 279 xmap:= function(P) 280 if PX(P) = INFTY then 281 return INFTY; 282 fi; 283 return a1^2*PX(P)+a3/a1; 284 end; 285 ymap:= function(P) 286 if PY(P) = INFTY then 287 return INFTY; 288 fi; 289 return a1^3*PY(P)+(a1^2*a4+a3^2)/a1^3+s*xmap(P); 290 end; 291 _Iso:= function(P) 292 if P = Zero(E) then 293 return Zero(E2); 294 else return Coerce(E,[xmap(P),ymap(P)]); 295 fi; 296 end; 297 Iso:=Map(E2,E,_Iso); 298 fi; 299 Iso.type:= iso; 300 Iso.operations:= _iso_ops; 301 Iso.xpol:= xmap(P); 302 Iso.ypol:= ymap(P); 303 Iso.degree:= Degree(Iso.xpol)/2; 304 AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]); 305 AssignNames_(Parent(Iso.ypol),["y"]); 306 return rec(base:=E2,ext1:=Iso); 307 end; 308 309 310 #recursion for division polynomials 311 _ff := function(EC,m) 312 local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n; 313 L:= Coefficients(EC); 314 if Length(L)=2 then 315 a1:=0; 316 a3:=0; 317 a2:=0; 318 a4:=L[1]; 319 a6:=L[2]; 320 else 321 a1:=L[1]; 322 a3:=L[2]; 323 a2:=L[3]; 324 a4:=L[4]; 325 a6:=L[5]; 326 fi; 327 b2:=a1^2+4*a2; 328 b4:=a1*a3+2*a4; 329 b6:=a3^2+4*a6; 330 b8:=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2; 331 kx:= PolynomialAlgebra(BaseRing(EC)); 332 x:= kx.1; 333 F:= EllipticFunctionField(EC); 334 y:=F.1; 335 if m<=0 then return 0; 336 elif m =1 then return 1; 337 elif m=2 then return 2*y+a1*x+a3; 338 elif m=3 then return 3*x^4+b2*x^3+3*b4*x^2+3*b6*x+b8; 339 elif m=4 then return (2*x^6+b2*x^5+5*b4*x^4+10*b6*x^3+10*b8*x^2+(b2*b8-b4*b6)*x+b4*b8-b6^2)*(2*y+a1*x+a3); 340 elif IsOdd(Coerce(Integers(),m)) then 341 n:=(m-1)/2; 342 return _ff(EC,n+2)*_ff(EC,n)^3-_ff(EC,n-1)*_ff(EC,n+1)^3; 343 else 344 n:=m/2; 345 return (_ff(EC,n+2)*_ff(EC,n-1)^2-_ff(EC,n-2)*_ff(EC,n+1)^2)*_ff(EC,n)/(2*y+a1*x+a3); 346 fi; 347 end; 348 349 350 #recursion for division polys over char 2 351 _ff2 := function(EC,m) 352 local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n,E,Iso,li; 353 L:= Coefficients(EC); 354 a2:=L[3]; 355 a6:=L[5]; 356 x:=PolynomialAlgebra(BaseRing(EC)).1; 357 F:= EllipticFunctionField(EC); 358 y:=F.1; 359 if m<=0 then return 0; 360 elif m =1 then return 1; 361 elif m=2 then return x; 362 elif m=3 then return 3*x^4+x^3+a6; 363 elif m=4 then return x^6+a6*x^2; 364 elif IsOdd(Coerce(Integers(),m)) then 365 n:=(m-1)/2; 366 return _ff2(EC,n+2)*_ff2(EC,n)^3-_ff2(EC,n-1)*_ff2(EC,n+1)^3; 367 else 368 n:=m/2; 369 return (_ff2(EC,n+2)*_ff2(EC,n-1)^2-_ff2(EC,n-2)*_ff2(EC,n+1)^2)*_ff2(EC,n)/x; 370 fi; 371 end; 372 373 #Division polynomials for characteristic 2 374 _DivisionPolynomial2:=function(EC,m) 375 local PA,L,a2,a6,Phi,Phi2,els,ro,kx,x,F,y,E,Iso; 376 L:= Coefficients(EC); 377 PA:= PolynomialAlgebra(BaseRing(EC)); 378 x:= PA.1; 379 return Coerce(PA,_ff2(EC,m)); 380 end; 381 382 #returns Division polynomial of degree m 383 DivisionPolynomial:=function(EC,m) 384 local PA,L,a1,a3,Phi,Phi2,els,ro,kx,x,F,y; 385 if Characteristic(BaseRing(EC)) =2 then 386 return _DivisionPolynomial2(EC,m); 387 fi; 388 L:= Coefficients(EC); 389 if Length(L) <>2 then 390 return Error("division polynomial: expecting elliptic curve in short WNF"); 391 fi; 392 a1:=L[1]; 393 a3:=L[2]; 394 kx:= PolynomialAlgebra(BaseRing(EC)); 395 x:= kx.1; 396 F:= EllipticFunctionField(EC); 397 y:=F.1; 398 PA:= PolynomialAlgebra(BaseRing(EC)); 399 if m mod 2 = 0 then 400 Phi:= Coerce(PA,_ff(EC,m)/_ff(EC,2)); 401 if Characteristic(BaseRing(EC))>3 then 402 if Length(L) = 2 then 403 return Phi*(x^3+a1*x+a3); 404 else return Error("case not implemented yet"); 405 fi; 406 elif Characteristic(BaseRing(EC))=3 then return Error("case not implemented yet"); 407 else Phi2:= -a1*x-a3; 408 fi; 409 return Phi*Phi2; 410 elif IsOdd(m) then 411 return Coerce(PA,_ff(EC,m)); 412 else 413 return Coerce(PA,_ff(EC,m)/_ff(EC,2)); 414 fi; 415 end; 416 417 418 419 420 #the characteristic polynomial of the frobenius endomorphism 421 CharacteristicFrobeniusPolynomial:= function(E) 422 local t,size,p,Zx; 423 size:= Size(E); 424 p:= Size(BaseRing(E)); 425 Zx:= PolynomialAlgebra(Integers()); 426 t:= p+1-size; 427 return Coerce(Zx,[p,-t,1]); 428 end; 429 430 431 #returns true if l is an Elkies prime of EC 432 IsElkiesPrime:=function(EC,l) 433 local mo,F,PA,ro,nr,i,ch,O; 434 ch:= CharacteristicFrobeniusPolynomial(EC); 435 O:=EquationOrder(ch); 436 if Conductor(O) mod l = 0 then 437 return rec(base:=false); 438 fi; 439 PA:= PolynomialAlgebra(GF(l)); 440 ch:= Coerce(PA,ch); 441 ro:=Roots(ch); 442 return rec(base:=ro <>Sequence([]),ext1:=ro); 443 end; 444 445 446 #returns Elkies primes <=37 of E 447 ElkiesPrimes:= function(E) 448 local i,el,bool; 449 el:= []; 450 for i in [2..Minimum(Size(BaseRing(E))-1,37)] do 451 if IsPrime(i) then 452 bool:=IsElkiesPrime(E,i); 453 if bool.base then 454 el:= Append(el,[i]); 455 fi; 456 fi; 457 od; 458 return el; 459 end; 460 461 #returns Elkies primes <=59 of E 462 ElkiesPrimes59:= function(E) 463 local i,el,elli,bool; 464 el:= []; 465 elli:=[]; 466 for i in [2..Minimum(Size(BaseRing(E))-1,59)] do 467 if IsPrime(i) then 468 bool:=IsElkiesPrime(E,i); 469 if bool.base then 470 el:= Append(el,[i]); 471 fi; 472 fi; 473 od; 474 return el; 475 end; 476 477 478 479 480 #the torsion subgroup EC[l]/K 481 TorsionSubgroup:=function(EC,l) 482 local di,i,ro,li,litmp; 483 di:= DivisionPolynomial(EC,l); 484 ro:= Roots(di); 485 if ro = [] then 486 return [Zero(EC)]; 487 else li:=[]; 488 for i in [1..Length(ro)] do 489 litmp:= MakePoints(EC,ro[i][1]); 490 Append_(li,litmp); 491 od; 492 Append_(li,[Zero(EC)]); 493 return li; 494 fi; 495 end; 496 497 #recursion for evaluation of Division polynomial at certain point in char 2 498 _DivisionPolynomialEval2 := function(EC,m,P) 499 local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n,E,Iso1,nor,Iso2,re,P2; 500 nor:=InNormalFormChar2(EC); 501 if not nor then 502 E:=NormalFormChar2(EC,1); 503 Iso1:=E.ext1; 504 Iso2:=NormalFormChar2(EC,2).ext1; 505 E:=Base(E); 506 else 507 Iso1:=false; 508 E:=EC; 509 fi; 510 L:= Coefficients(E); 511 a2:=L[3]; 512 a6:=L[5]; 513 if Iso1 = false then 514 x:=PX(P); 515 y:=PY(P); 516 else 517 x:=PX(Iso1(P)); 518 y:=PY(Iso1(P)); 519 fi; 520 if m<=0 then return 0; 521 elif m =1 then 522 re:=1; 523 if Iso1 = false then 524 return re; 525 else 526 P2:=MakePoints(E,re)[1]; 527 return PX(Iso2(P2)); 528 fi; 529 elif m=2 then 530 re:=x; 531 if Iso1 = false then 532 return re; 533 else 534 P2:=MakePoints(E,re)[1]; 535 return PX(Iso2(P2)); 536 fi; 537 elif m=3 then re:= x^4+x^3+a6; 538 if Iso1 = false then 539 return re; 540 else 541 P2:=MakePoints(E,re)[1]; 542 return PX(Iso2(P2)); 543 fi; 544 elif m=4 then re:= x^6+a6*x^2; 545 if Iso1 = false then 546 return re; 547 else 548 P2:=MakePoints(E,re)[1]; 549 return PX(Iso2(P2)); 550 fi; 551 elif IsOdd(Coerce(Integers(),m)) then 552 n:=(m-1)/2; 553 re:= _DivisionPolynomialEval2(EC,n+2,P)*_DivisionPolynomialEval2(EC,n,P)^3+_DivisionPolynomialEval2(EC,n-1,P)*_DivisionPolynomialEval2(EC,n+1,P)^3; 554 if Iso1 = false then 555 return re; 556 else 557 P2:=MakePoints(E,re)[1]; 558 return PX(Iso2(P2)); 559 fi; 560 else 561 n:=m/2; 562 re:= (_DivisionPolynomialEval2(EC,n+2,P)*_DivisionPolynomialEval2(EC,n-1,P)^2+_DivisionPolynomialEval2(EC,n-2,P)*_DivisionPolynomialEval2(EC,n+1,P)^2)*_DivisionPolynomialEval2(EC,n,P)/x; 563 if Iso1 = false then 564 return re; 565 else 566 P2:=MakePoints(E,re)[1]; 567 return PX(Iso2(P2)); 568 fi; 569 fi; 570 end; 571 572 573 574 575 #recursion for evaluation of Division polynomial at certain point 576 _DivisionPolynomialEval := function(EC,m,P) 577 local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n; 578 if Characteristic(BaseRing(EC))=2 then 579 return _DivisionPolynomialEval2(EC,m,P); 580 fi; 581 L:= Coefficients(EC); 582 if Length(L)=2 then 583 a1:=0; 584 a3:=0; 585 a2:=0; 586 a4:=L[1]; 587 a6:=L[2]; 588 else 589 a1:=L[1]; 590 a3:=L[2]; 591 a2:=L[3]; 592 a4:=L[4]; 593 a6:=L[5]; 594 fi; 595 b2:=a1^2+4*a2; 596 b4:=a1*a3+2*a4; 597 b6:=a3^2+4*a6; 598 b8:=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2; 599 kx:= PolynomialAlgebra(BaseRing(EC)); 600 x:= PX(P); 601 F:= EllipticFunctionField(EC); 602 y:= PY(P); 603 if m<=0 then return 0; 604 elif m =1 then return 1; 605 elif m=2 then return 2*y+a1*x+a3; 606 elif m=3 then return 3*x^4+b2*x^3+3*b4*x^2+3*b6*x+b8; 607 elif m=4 then return (2*x^6+b2*x^5+5*b4*x^4+10*b6*x^3+10*b8*x^2+(b2*b8-b4*b6)*x+b4*b8-b6^2)*(2*y+a1*x+a3); 608 elif IsOdd(Coerce(Integers(),m)) then 609 n:=(m-1)/2; 610 return _DivisionPolynomialEval(EC,n+2,P)*_DivisionPolynomialEval(EC,n,P)^3-_DivisionPolynomialEval(EC,n-1,P)*_DivisionPolynomialEval(EC,n+1,P)^3; 611 else 612 n:=m/2; 613 return (_DivisionPolynomialEval(EC,n+2,P)*_DivisionPolynomialEval(EC,n-1,P)^2-_DivisionPolynomialEval(EC,n-2,P)*_DivisionPolynomialEval(EC,n+1,P)^2)*_DivisionPolynomialEval(EC,n,P)/(2*y+a1*x+a3); 614 fi; 615 end; 616 617 618 619 620 621 622 #Conductor of an order 623 _IntCond:= function(O) 624 local Omax; 625 Omax:= MaximalOrder(O); 626 return Index(Omax,O); 627 end; 628 629 630 #the i-th power of the frobenius endomorphism of E 631 Frobenius:= function(E,i) 632 local ma,p,frob,EE,P,q; 633 P:= GenericPoint(E); 634 EE:= Parent(P); 635 p:= Size(BaseRing(E)); 636 q:= p^i; 637 frob:= function(P1) 638 local P2,E2; 639 E2:= Parent(P1); 640 return MakePoint(E2,[PX(P1)^q, PY(P1)^q]); 641 end; 642 ma:= Map(E,E,frob); 643 ma.type:= iso; 644 ma.fro:=true; 645 ma.operations:= _iso_ops; 646 ma.degree:=q; 647 return ma; 648 end; 649 650 651 #isogeny of degree 0, i.e. phi(P) = O for all P in E1(K) 652 NuMap:=function(E1,E2) 653 local ma,ze; 654 ze:=function(P) 655 return Zero(E2); 656 end; 657 ma:=Map(E1,E2,ze); 658 ma.type:=iso; 659 ma.operations:= _iso_ops; 660 ma.degree:=0; 661 ma.zero:=true; 662 return ma; 663 end; 664 665 #converts integer to binary list 666 _IntToBinList:=function(n) 667 local L; 668 L:=[]; 669 while n>0 do 670 Append_(L,[n mod 2]); 671 n:= Div(n,2); 672 od; 673 return L; 674 end; 675 #converts integer to finite field element, T generator of the field, only for characteristic 2 676 _IntToFldElt:=function(n,T) 677 local T,K,L,i,el; 678 L:=_IntToBinList(n); 679 el:=0; 680 for i in [1..Length(L)] do 681 el:=el+L[i]*T^(i-1); 682 od; 683 return el; 684 end; 685 686 687 688 #returns the dual isogeny of f 689 DualIsogeny:= function(f) 690 local VR,ba,li,M,apoly,bpoly,polys,i,d,isodu,xpoly,xmap,E2,Fu,tmp,x,fu,mux, xpol,ypol,ymap,ypoly,y,yoben,yunten,BR,E,ymupol,Funf,xpol,muma,xmuma,xoben,xunten,PA,xmupol,Iso,f,dy,fact,P,j,dego1,dego2,degu1,degu2,GenP,_Iso,d2,xpart,ypart,xdenom,ydenom,yy,fact,ro,E3,psi,isotest; 691 E:= f.domain; 692 E2:= f.codomain; 693 BR:= BaseRing(E); 694 PA:= PolynomialAlgebra(BR); 695 polys:= []; 696 xpol:= f.xpol; 697 d:= Coerce(Integers(),Degree(f)); 698 if not IsPrime(d) then 699 return Error("assuming prime degree for dual isogeny-computation"); 700 fi; 701 muma:=MultMap(E,d); 702 for i in [1..d] do 703 polys[i]:= xpol^i; 704 od; 705 M:= [xpol^0]; 706 for i in [2..d+1] do 707 M[i]:= polys[i-1]; 708 od; 709 mux:=-Coerce(Parent(xpol),muma.xpol); 710 M:= Append(M,[mux]); 711 for i in [2..d] do 712 M:= Append(M,[mux*polys[i-1]]); 713 od; 714 VR:= Relations(M,BR,1); 715 ba:= Basis(VR); 716 Assert(ba<> Sequence([])," Vectorspace of dimension > 0 "); 717 li:= List(ba[1]); 718 apoly:= []; 719 bpoly:= []; 720 for i in [1..d+1] do 721 apoly:= Append(apoly,[li[i]]); 722 od; 723 for i in [d+2..Length(li)] do 724 bpoly:= Append(bpoly,[li[i]]); 725 od; 726 xunten:= Coerce(PA,bpoly); 727 ro:=Roots(xunten,SplittingField(xunten))[1][1]; 728 P:=MakePoints(ChangeBaseRing(E2,Extension(Parent(ro),2)),ro)[1]; 729 isodu:=IsogenyFromSubgroup(P,E2,d); 730 E3:=Codomain(isodu); 731 if not E3 = E then 732 isodu:= Composition(Isomorphism(E3,E,1),isodu); 733 fi; 734 isotest:=Composition(isodu,f); 735 P:=RandomPoint(E); 736 while 2*P = Zero(E) do 737 P:=RandomPoint(E); 738 od; 739 if d*P = isotest(P) then 740 return isodu; 741 else Assert(d*P = -isotest(P), "dual assumption"); 742 fi; 743 return Negate(isodu); 744 end; 745 746 747 748 #multiplication-by-m-map 749 MultMap:= function(E1,m) 750 local muma,mu,mult,P,EE,xmuma,ymuma,xmu,ymu,ymuma2,cha,E,Iso1,Iso2,E2,di,Psi,Psip,Psim,Phi,Psim2,omega,x,y; 751 if m = 0 then return NuMap(E1,E1); fi; 752 if m = 1 then return Isomorphism(E1,E1,1);fi; 753 cha:= Characteristic(BaseRing(E1)); 754 if cha = 2 then 755 E:=NormalFormChar2(E1,1); 756 Iso1:=E.ext1; 757 E:=Base(E); 758 Iso2:=NormalFormChar2(E1,2).ext1; 759 else 760 if Length(Coefficients(E1)) = 2 then 761 Iso1:=false; 762 E:=E1; 763 else 764 E:=ShortWNF(E1); 765 Iso1:=E.ext1; 766 Iso2:=DualIsogeny(Iso1); 767 E:=Base(E); 768 fi; 769 fi; 770 P:= GenericPoint(E); 771 EE:= Parent(P); 772 if cha = 2 then 773 Psi:= _ff(E,m); 774 Psip:= _ff(E,m+1); 775 Psim:= _ff(E,m-1); 776 Phi:= Psip*Psim; 777 Psim2:=_ff(E,m-2); 778 xmuma:= function(P2) 779 local psi,phi,BR,PA; 780 if Is(Type(PX(P2)),elt-fld^fin) then 781 BR:=BaseRing(Parent(P2)); 782 PA:=PolynomialAlgebra(BR); 783 psi:=Coerce(PA,Psi); 784 if Evaluate(psi,PX(P2)) = 0 then 785 return "INFTY POINT"; 786 fi; 787 phi:=Coerce(PA,Phi); 788 return PX(P2)+Evaluate(phi/psi^2,PX(P2)); 789 fi; 790 if Evaluate(Psi,PX(P2)) = 0 then 791 return "INFTY POINT"; 792 fi; 793 return PX(P2)+Evaluate((Phi/Psi^2),PX(P2)); 794 end; 795 ymuma:= function(P2) 796 local Phi2,phi2,psi,psim,psip,psim2,BR,PA; 797 if Is(Type(PX(P2)),elt-fld^fin) then 798 BR:=BaseRing(Parent(P2)); 799 PA:=PolynomialAlgebra(BR); 800 psi:=Coerce(PA,Psi); 801 psim:=Coerce(PA,Psim); 802 psip:=Coerce(PA,Psip); 803 psim2:=Coerce(PA,Psim2); 804 phi2:= (PX(P2)^2+PX(P2)+PY(P2))*Evaluate(psim*psi*psip,PX(P2))+Evaluate(psim2*psip^2,PX(P2)); 805 return PX(P2)+PY(P2)+ phi2/(PX(P2)*Evaluate(psi^3,PX(P2))); 806 fi; 807 Phi2:= (PX(P2)^2+PX(P2)+PY(P2))*Evaluate(Psim*Psi*Psip,PX(P2))+Evaluate(Psim2*Psip^2,PX(P2)); 808 return PX(P2)+PY(P2)+ Phi2/(PX(P2)*Evaluate(Psi^3,PX(P2))); 809 end; 810 else 811 Psi:= _ff(E,m); 812 Psip:= _ff(E,m+1); 813 Psim:= _ff(E,m-1); 814 if Is(Type(Psi),elt-fld^fun) then 815 x:=CoefficientRing(Parent(Psi)).1; 816 else 817 x:=PolynomialAlgebra(BaseRing(E)).1; 818 fi; 819 Phi:= x*Psi^2-Psip*Psim; 820 omega:=_ff(E,m+2)*Psim^2-_ff(E,m-2)*Psip^2; 821 y:=EllipticFunctionField(E).1; 822 xmuma:= function(P2) 823 local BR,PA,phi,psi; 824 if Is(Type(PX(P2)),elt-fld^fin) then 825 BR:=BaseRing(Parent(P2)); 826 PA:=PolynomialAlgebra(BR); 827 psi:=Coerce(PA,Psi); 828 if Evaluate(psi,PX(P2)) = 0 then 829 return "INFTY POINT"; 830 fi; 831 phi:=Coerce(PA,Phi); 832 return Evaluate(phi/psi^2,PX(P2)); 833 fi; 834 if Evaluate(Psi,PX(P2)) = 0 then 835 return "INFTY POINT"; 836 fi; 837 return Evaluate(Phi/Psi^2,PX(P2)); 838 end; 839 ymuma:= function(P2) 840 local omega2,psi,BR,PA; 841 if Is(Type(PX(P2)),elt-fld^fin) then 842 BR:=BaseRing(Parent(P2)); 843 PA:=PolynomialAlgebra(BR); 844 omega2:=Coerce(PA,omega); 845 omega2:=Evaluate(omega2,PX(P2))/(4*PY(P2)); 846 psi:=Coerce(PA,Psi); 847 return omega2/(Evaluate(psi,PX(P2))^3); 848 fi; 849 omega2:= Evaluate(omega,PX(P2))/(4*PY(P2)); 850 return omega2/(Evaluate(Psi,PX(P2)))^3; 851 end; 852 fi; 853 muma:=function(P1) 854 local P2,EE,P; 855 if P1 = Zero(E) then 856 return P1; 857 fi; 858 P:= GenericPoint(E); 859 EE:= Parent(P); 860 P2:= Coerce(EE,P1); 861 xmu:= xmuma(P1); 862 ymu:=ymuma(P1); 863 if xmu = "INFTY POINT" then 864 return Zero(E); 865 fi; 866 return MakePoint(E,[xmu,ymu]); 867 end; 868 mu:= Map(E,E,muma); 869 mu.type:=iso; 870 mu.operations:= _iso_ops; 871 if Characteristic(BaseRing(E))>3 then 872 mu.xpol:= Phi/Psi^2; 873 mu.ypol:= omega/(4*y*Psi^3); 874 else 875 x:=PolynomialAlgebra(BaseRing(E)).1; 876 y:=EllipticFunctionField(E).1; 877 mu.xpol:=x+Phi/Psi^2; 878 mu.ypol:=x+y+(((x^2+x+y)*Psim*Psi*Psip)+Psim2*Psip^2)/(x*Psi^3); 879 fi; 880 if GCD(cha,Coerce(Integers(),m))<>1 then 881 mu.degree:= Degree(mu.xpol)/2; 882 else mu.degree:=m^2; 883 fi; 884 if Iso1 = false then 885 return mu; 886 else return 887 Composition(Iso2,Composition(mu,Iso1)); 888 fi; 889 end; 890 891 #returns true if Z[pi] maximal order 892 EndRingIsMaximal:= function(EC) 893 local q,eq,nop,t,eq,orde,z; 894 z:= PolynomialAlgebra(Integers()).1; 895 q:= Size(EC.F); 896 nop:= Size(EC); 897 t:= q+1-nop; 898 eq:= EquationOrder(z^2-t*z+q); 899 orde:=rec(base:=IsMaximal(eq)); 900 if orde.base then 901 orde.ext1:= eq; 902 fi; 903 return orde; 904 end; 905 906 907 908 909 910 911 912 #returns an isomorphism between E1 and E2 over an extension of ground field of degree n 913 Isomorphism:= function(E1,E2,n) 914 local _Iso,L1,L2,xmap,ymap,P,a,b,A,B,j,delta,u1,u2,Iso,u,i,K,E12,E22,s,a1,a2,x,Iso1,Iso2; 915 if not jInvariant(E1) = jInvariant(E2) then 916 return Error(E1," and ", E2," are not isomorphic."); 917 fi; 918 if E1 = E2 then 919 Iso:=Map(E1,E2,function(P) return P; end); 920 Iso.type:= iso; 921 Iso.operations:= _iso_ops; 922 P:= GenericPoint(E1); 923 Iso.xpol:= PX(P); 924 Iso.ypol:= PY(P); 925 Iso.degree:= Degree(Iso.xpol)/2; 926 AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]); 927 AssignNames_(Parent(Iso.ypol),["y"]); 928 return Iso; 929 fi; 930 K:=Base(Extension(BaseRing(E1),n)); 931 if Characteristic(K)=2 then 932 E12:=NormalFormChar2(ChangeBaseRing(E1,K),1); 933 Iso1:=E12.ext1; 934 E12:=ChangeBaseRing(Base(E12),K); 935 E22:=NormalFormChar2(ChangeBaseRing(E2,K),2); 936 Iso2:=E22.ext1; 937 E22:=ChangeBaseRing(Base(E22),K); 938 else 939 E12:=ChangeBaseRing(E1,K); 940 E22:=ChangeBaseRing(E2,K); 941 fi; 942 L1:= Coefficients(E12); 943 L2:= Coefficients(E22); 944 P:= GenericPoint(E12); 945 if Length(L1) =2 then 946 a:= L1[1]; 947 b:= L1[2]; 948 A:= L2[1]; 949 B:= L2[2]; 950 if a = 0 then 951 u:=AllRoots(b/B,6); 952 if u = Sequence([]) then return Error("no isomorphism defined over BaseRing"); 953 else 954 if u[1] =1 then 955 u := u[2]; 956 else u:= u[1]; 957 fi; 958 fi; 959 else 960 u1:=AllRoots(A/a,4); 961 if u1 = Sequence([]) then return Error("no isomorphism defined over BaseRing"); 962 else 963 if b=0 then 964 if u1[1] =1 then 965 u := u1[2]; 966 else u:= u1[1]; 967 fi; 968 else u2 := AllRoots(B/b,6); 969 if u2 =Sequence([]) then return Error("no isomorphism defined over BaseRing"); 970 fi; 971 i:=0; 972 j:=0; 973 while i <Length(u1) do 974 i:=i+1; 975 if not u1[i] =1 then 976 j:=0; 977 while j <Length(u2) do 978 j:=j+1; 979 if u2[j] = u1[i] then 980 u:= u1[i]; 981 j:= Length(u2); 982 i:= Length(u1); 983 fi; 984 od; 985 fi; 986 od; 987 fi; 988 fi; 989 fi; 990 xmap:= function(P) 991 if PX(P) = INFTY then return INFTY; fi; 992 return PX(P)/(u^2); 993 end; 994 ymap:= function(P) 995 if PY(P) = INFTY then return INFTY; fi; 996 return PY(P)/(u^3); 997 end; 998 _Iso:= function(P) 999 if P = Zero(E12) then 1000 return Zero(E22); 1001 else return Coerce(E22,[xmap(P),ymap(P)]); 1002 fi; 1003 end; 1004 Iso:=Map(E12,E22,_Iso); 1005 Iso.type:= iso; 1006 Iso.operations:= _iso_ops; 1007 Iso.xpol:= xmap(P); 1008 Iso.ypol:= ymap(P); 1009 Iso.degree:= Degree(Iso.xpol)/2; 1010 AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]); 1011 AssignNames_(Parent(Iso.ypol),["y"]); 1012 return Iso; 1013 elif Characteristic(K) =2 then 1014 Assert(L1[2]=0 and L1[4]=0," curve in characteristic 2 form "); 1015 a1:=L1[3]; 1016 a2:=L2[3]; 1017 x:=PolynomialAlgebra(K).1; 1018 s:=Roots(x^2+x+a1+a2); 1019 if s = Sequence([]) then 1020 return Error("no isomorphism defined over BaseRing"); 1021 else s:=s[1][1]; 1022 fi; 1023 xmap:= function(P) 1024 if PX(P) = INFTY then return INFTY; fi; 1025 return PX(P); 1026 end; 1027 ymap:= function(P) 1028 if PY(P) = INFTY then return INFTY; fi; 1029 return PY(P)-s*PX(P); 1030 end; 1031 _Iso:= function(P) 1032 if P = Zero(E12) then 1033 return Zero(E22); 1034 else return Coerce(E22,[xmap(P),ymap(P)]); 1035 fi; 1036 end; 1037 Iso:=Map(E12,E22,_Iso); 1038 Iso.type:= iso; 1039 Iso.operations:= _iso_ops; 1040 Iso.xpol:= xmap(P); 1041 Iso.ypol:= ymap(P); 1042 Iso.degree:= Degree(Iso.xpol)/2; 1043 AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]); 1044 AssignNames_(Parent(Iso.ypol),["y"]); 1045 return Composition(Iso2,Composition(Iso,Iso1)); 1046 else return Error("missing case"); 1047 fi; 1048 end; 1049 1050 1051 #new base ring is extension of degree n 1052 ChangeBaseRingN:=function(E,n) 1053 local K; 1054 K:=Extension(BaseRing(E),n); 1055 return ChangeBaseRing(E,K); 1056 end; 1057 1058 #G subgroup of E(K), computes E/G 1059 GetCodomain:= function(E,G) 1060 local A1,A2,A3,A4,A6,E2,F2,xx,M,a4,a6,a1,a3,a2,b2,b4,b6,nu,i,tq,uq,S,t,w,tQ,E2; 1061 E2:=Parent(G[1]); 1062 F2:= []; 1063 xx:=[]; 1064 M:=[]; 1065 if Length(E.L) = 2 then 1066 a4:= E.L[1]; 1067 a6:= E.L[2]; 1068 a1:= 0; 1069 a3:= 0; 1070 a2:= 0; 1071 else a1:= E.L[1]; 1072 a3:= E.L[2]; 1073 a2:= E.L[3]; 1074 a4:= E.L[4]; 1075 a6:= E.L[5]; 1076 fi; 1077 b2:= a1^2+4*a2; 1078 b4:=a1*a3+2*a4; 1079 b6:= a3^2+4*a6; 1080 nu:= Zero(E2); 1081 for i in [1..Length(G)] do 1082 if G[i]<>nu then 1083 if 2*G[i] = nu then 1084 F2:=Append(F2,[G[i]]); 1085 elif not G[i].x in xx then 1086 xx:=Append(xx,[G[i].x]); 1087 M:=Append(M,[MakePoints(E2,[G[i].x])[1]]); 1088 fi; 1089 fi; 1090 od; 1091 tq:= function(Q) 1092 if 2*Q = Zero(Parent(Q)) then 1093 return (3*Q.x^2+2*a2*Q.x+a4-a1*Q.y); 1094 else return 6*Q.x^2+b2*Q.x+b4; 1095 fi; 1096 end; 1097 uq:= function(Q) 1098 return 4*Q.x^3+b2*Q.x^2+2*b4*Q.x+b6; 1099 end; 1100 S:= Append(F2,M); 1101 t:=0; 1102 w:=0; 1103 for i in [1..Length(S)] do 1104 Q:= S[i]; 1105 tQ:= tq(Q); 1106 t:= t+tQ; 1107 w:= w+(uq(Q)+PX(Q)*tQ); 1108 od; 1109 A1:= a1; 1110 A2:=a2; 1111 A3:= a3; 1112 A4:= a4- 5*t; 1113 A6:= a6-b2*t-7*w; 1114 if A1 = 0 and A3 =0 and A2 = 0 then 1115 E2:= EllipticCurve(BaseRing(E),[A4,A6]); 1116 else 1117 E2:= EllipticCurve(BaseRing(E),[A1,A3,A2,A4,A6]); 1118 fi; 1119 return E2; 1120 end; 1121 1122 #the degree of an isogeny 1123 _IsoDegree:= function(f) 1124 return f.degree; 1125 end; 1126 1127 #returns E[l] 1128 CompleteTorsionSubgroup:=function(E,l) 1129 local poi,E2,di,ro,nu,i,Sp,PA,F,din,ex,k,li,l,BR; 1130 di:= DivisionPolynomial(E,l); 1131 BR:= BaseRing(E); 1132 Sp:= Extension(SplittingField(di),2); 1133 di:=Coerce(PolynomialAlgebra(Sp),di); 1134 ro:= Roots(di); 1135 E2:= ChangeBaseRing(E,Sp); 1136 poi:= []; 1137 for k in [1..Length(ro)] do 1138 li:= MakePoints(E2,ro[k][1]); 1139 if (Length(li) <> 2 and not IsCoercible(Integers(),l/2)) then return Error("something went wrong"); 1140 fi; 1141 poi:= Append(poi,li); 1142 od; 1143 Append_(poi,[Zero(E2)]); 1144 return poi; 1145 end; 1146 1147 #return isogeny phi:E ->E/P if P has order 2 1148 IsoFromSub2:=function(P,E,E2) 1149 local coeffs,a1,a2,a3,a4,a6,b2,b4,b6,xmap,ymap,_Iso,PP,Zx,EE,polx,poly,Iso,FFchange,f,PA2,xdenom,xnom,xpol,els,f1,f2,ydenom,ynom,y,ypol,PA,P2; 1150 coeffs:=Coefficients(E); 1151 if Length(coeffs) = 2 then 1152 a1:=0; 1153 a2:=0; 1154 a3:=0; 1155 a4:=coeffs[1]; 1156 a6:=coeffs[2]; 1157 else 1158 a1:=coeffs[1]; 1159 a2:=coeffs[3]; 1160 a3:=coeffs[2]; 1161 a4:=coeffs[4]; 1162 a6:=coeffs[5]; 1163 fi; 1164 b2:=a1^2+4*a2; 1165 b4:=a1*a3+2*a4; 1166 b6:=a3^2+4*a6; 1167 PA:=PolynomialAlgebra(BaseRing(E)); 1168 xmap:=function(Q) 1169 local PR; 1170 PR:=Coerce(Parent(Q),P); 1171 return PX(Q)+PX(PR+Q) - PX(PR); 1172 end; 1173 ymap:=function(Q) 1174 local PR; 1175 PR:=Coerce(Parent(Q),P); 1176 return PY(Q)+PY(PR+Q) - PY(PR); 1177 end; 1178 _Iso:=function(Q) 1179 if Q = P or Q = Zero(E) then 1180 return Zero(E2); 1181 else 1182 return MakePoint(E2,[xmap(Q),ymap(Q)]); 1183 fi; 1184 end; 1185 PP:= GenericPoint(Parent(P)); 1186 Zx:= PolynomialAlgebra(Integers()); 1187 EE:= Parent(PP); 1188 P2:=Coerce(EE,PP); 1189 polx:= xmap(PP); 1190 poly:=ymap(PP); 1191 Iso:= Map(E,E2,_Iso); 1192 Iso.operations:= _iso_ops; 1193 Iso.type:= iso; 1194 FFchange:=true; 1195 if Parent(polx)= EllipticFunctionField(E) then 1196 Iso.xpol:= polx; 1197 Iso.ypol:= poly; 1198 FFchange:=false; 1199 else 1200 f:=Eltseq(polx)[1]; 1201 PA2:=PolynomialAlgebra(CoefficientRing(Parent(f))); 1202 xdenom:=Denominator(f); 1203 xnom:=Coerce(PA,Coerce(PA2,xdenom*f)); 1204 xdenom:=Coerce(PA,xdenom); 1205 Iso.xnom:=xnom; 1206 Iso.xdenom:=xdenom; 1207 xpol:=Coerce(EllipticFunctionField(E),xnom/xdenom); 1208 els:=Eltseq(poly); 1209 f1:= els[1]; 1210 f2:= els[2]; 1211 xdenom:=Denominator(f1); 1212 xnom:=Coerce(PA,Coerce(PA2,xdenom*f1)); 1213 xdenom:=Coerce(PA,xdenom); 1214 ydenom:=Denominator(f2); 1215 ynom:=Coerce(PA,Coerce(PA2,ydenom*f2)); 1216 ydenom:=Coerce(PA,ydenom); 1217 y:=EllipticFunctionField(E).1; 1218 ypol:=xnom/xdenom+(ynom/ydenom)*y; 1219 Iso.xpol:=xpol; 1220 Iso.ypol:=ypol; 1221 fi; 1222 Iso.degree:= Degree(Iso.xpol)/2; 1223 Iso.sub:= [P]; 1224 AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]); 1225 AssignNames_(Parent(Iso.ypol),["y"]); 1226 if not FFchange then 1227 els:= Eltseq(Iso.xpol); 1228 xdenom:= Denominator(els[1]); 1229 Iso.xnom:=Coerce(PolynomialAlgebra(BaseRing(Parent(P2))), els[1]*xdenom); 1230 Iso.xdenom:= xdenom; 1231 fi; 1232 return Iso; 1233 end; 1234 1235 1236 #returns phi:E->E/<P> of degree l 1237 IsogenyFromSubgroup:= function(P,E,l) 1238 local PA,P2,EE,i,Iso,xmap,ymap,E2,_Iso,E,G,L,P,Zx,polx,xdenom,dseq,nseq,xnom,ydenom,ynom,els,BR,poly,f,y,ydenom,ynom,els,f1,f2,FFchange,xpol,ypol,PA2,x,y,mytime, coeffs,a1,a2,a3,a4,a6,b2,b4,b6,R,F2,xcoords,ze,t,w,S,tqtmp,tq,uq,gqy,gqx,xxunten,d,fro; 1239 fro:=extFrob(E,Degree(BaseRing(Parent(P)),BaseRing(E))); 1240 if not fro(P) in Subgroup(P) then 1241 return Error("Subgroup is not galois-invariant"); 1242 fi; 1243 if P.x = "INFTY POINT" then 1244 return Isomorphism(E,E,1); 1245 fi; 1246 if l=2 then 1247 S:=[P]; 1248 else 1249 d:=Coerce(Integers(),(l-1)/2); 1250 S:=[P]; 1251 for i in [2..d] do 1252 Append_(S,[i*P]); 1253 od; 1254 fi; 1255 BR:=BaseRing(E); 1256 PA:=PolynomialAlgebra(BR); 1257 coeffs:=Coefficients(E); 1258 if Length(coeffs) = 2 then 1259 a1:=0; 1260 a2:=0; 1261 a3:=0; 1262 a4:=coeffs[1]; 1263 a6:=coeffs[2]; 1264 else 1265 a1:=coeffs[1]; 1266 a2:=coeffs[3]; 1267 a3:=coeffs[2]; 1268 a4:=coeffs[4]; 1269 a6:=coeffs[5]; 1270 fi; 1271 b2:=a1^2+4*a2; 1272 b4:=a1*a3+2*a4; 1273 b6:=a3^2+4*a6; 1274 gqx:=function(P) 1275 return 3*PX(P)^2+2*a2*PX(P)+a4-a1*PY(P); 1276 end; 1277 gqy:=function(P) 1278 return -2*PY(P) - a1*PX(P) -a3; 1279 end; 1280 uq:=function(P) 1281 return 4*PX(P)^3+b2*PX(P)^2+2*b4*PX(P)+b6; 1282 end; 1283 tq:=function(P) 1284 if 2*P = Zero(Parent(P)) then 1285 return gqx(P); 1286 else 1287 return 6*PX(P)^2+b2*PX(P)+b4; 1288 fi; 1289 end; 1290 t:=0; 1291 w:=0; 1292 for i in [1..Length(S)] do 1293 tqtmp:=tq(S[i]); 1294 t:=t+tqtmp; 1295 w:=w+uq(S[i])+tqtmp*PX(S[i]); 1296 od; 1297 if Length(E.L) =2 then 1298 E2:= EllipticCurve(BaseRing(E),[a4-5*t,a6-b2*t-7*w]); 1299 else 1300 E2:= EllipticCurve(BaseRing(E),[a1,a3,a2,a4-5*t,a6-b2*t-7*w]); 1301 fi; 1302 E2.size:=E.size; 1303 if 2*P = Zero(Parent(P)) then 1304 return IsoFromSub2(P,E,E2); 1305 fi; 1306 P2:=P; 1307 xmap:= function(P1) 1308 local sum,i,xx,xxunten; 1309 xx := PX(P1); 1310 sum:=xx; 1311 for i in [1..Length(S)] do 1312 if sum = INFTY then 1313 return sum; 1314 fi; 1315 xxunten:=xx-PX(S[i]); 1316 if xxunten = 0 then 1317 return INFTY; 1318 fi; 1319 sum := sum +tq(S[i])/xxunten+uq(S[i])/(xxunten^2); 1320 od; 1321 return sum; 1322 end; 1323 ymap:= function(P1) 1324 local sum,i,yy,Q,xx,xxunten; 1325 xx:=PX(P1); 1326 yy := PY(P1); 1327 sum:=yy; 1328 for i in [1..Length(S)] do 1329 if sum = INFTY then 1330 return sum; 1331 fi; 1332 Q:=S[i]; 1333 xxunten:=xx-PX(Q); 1334 if xxunten = 0 then 1335 return INFTY; 1336 fi; 1337 sum := sum -(uq(Q)*((2*yy+a1*xx+a3)/(xxunten)^3) + tq(Q)*((a1*xxunten+yy-PY(Q))/(xxunten)^2)+(a1*uq(Q)-gqx(Q)*gqy(Q))/(xxunten)^2); 1338 od; 1339 return sum; 1340 end; 1341 _Iso:=function(P1) 1342 local P2,m; 1343 if P1 = Zero(Parent(P1)) then 1344 return Zero(E2); 1345 else 1346 if xmap(P1) = INFTY then 1347 return Zero(E2); 1348 fi; 1349 fi; 1350 return MakePoint(E2,[xmap(P1),ymap(P1)]); 1351 end; 1352 P:= GenericPoint(Parent(P)); 1353 Zx:= PolynomialAlgebra(Integers()); 1354 EE:= Parent(P); 1355 L:= []; 1356 for i in [1..Length(S)] do 1357 S[i]:=Coerce(EE,S[i]); 1358 od; 1359 polx:= xmap(P); 1360 poly:=ymap(P); 1361 E2.size:=Size(E); 1362 Iso:= Map(E,E2,_Iso); 1363 Iso.operations:= _iso_ops; 1364 Iso.type:= iso; 1365 FFchange:=true; 1366 if Parent(polx)= EllipticFunctionField(E) then 1367 Iso.xpol:= polx; 1368 Iso.ypol:= poly; 1369 FFchange:=false; 1370 else 1371 f:=Eltseq(polx)[1]; 1372 PA2:=PolynomialAlgebra(CoefficientRing(Parent(f))); 1373 xdenom:=Denominator(f); 1374 xnom:=Coerce(PA,Coerce(PA2,xdenom*f)); 1375 xdenom:=Coerce(PA,xdenom); 1376 Iso.xnom:=xnom; 1377 Iso.xdenom:=xdenom; 1378 xpol:=Coerce(EllipticFunctionField(E),xnom/xdenom); 1379 els:=Eltseq(poly); 1380 f1:= els[1]; 1381 f2:= els[2]; 1382 xdenom:=Denominator(f1); 1383 xnom:=Coerce(PA,Coerce(PA2,xdenom*f1)); 1384 xdenom:=Coerce(PA,xdenom); 1385 ydenom:=Denominator(f2); 1386 ynom:=Coerce(PA,Coerce(PA2,ydenom*f2)); 1387 ydenom:=Coerce(PA,ydenom); 1388 y:=EllipticFunctionField(E).1; 1389 ypol:=xnom/xdenom+(ynom/ydenom)*y; 1390 Iso.xpol:=xpol; 1391 Iso.ypol:=ypol; 1392 fi; 1393 Iso.degree:= Degree(Iso.xpol)/2; 1394 AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]); 1395 AssignNames_(Parent(Iso.ypol),["y"]); 1396 if not FFchange then 1397 els:= Eltseq(Iso.xpol); 1398 xdenom:= Denominator(els[1]); 1399 Iso.xnom:=Coerce(PolynomialAlgebra(BaseRing(Parent(P2))), els[1]*xdenom); 1400 Iso.xdenom:= xdenom; 1401 fi; 1402 Codomain(Iso).size:=E.size; 1403 return Iso; 1404 end; 1405 1406 #frobenius applicable to points of E(GF(q^i)) 1407 extFrob:= function(E,i) 1408 local E2,BR,F,fro,re; 1409 BR:= BaseRing(E); 1410 F:= Extension(BR,i); 1411 E2:= ChangeBaseRing(E,F); 1412 fro:= Frobenius(E,1); 1413 re:= Map(E2,E2,fro.base); 1414 re.domain:=E2; 1415 re.codomain:=E2; 1416 re.type:=iso; 1417 re.operations:=fro.operations; 1418 re.fro:=fro.fro; 1419 re.degree:=fro.degree; 1420 return re; 1421 end; 1422 1423 1424 #subfunction for computing isogenies between elliptic curves over fields of characteristic 2 1425 B:=function(i,j) 1426 return (Binomial(i,j) mod 2); 1427 end; 1428 1429 1430 #returns elliptic curve with j-invariant j and size n 1431 ECFromJGivenNumber:=function(j,n) 1432 local E,P,i,O,q,t,l1,l2,E2; 1433 E:=EllipticCurveFromJInvariant(j); 1434 O:=Zero(E); 1435 q:=Size(Parent(j)); 1436 t:=q+1-n; 1437 P:=RandomPoint(E); 1438 l1:=n; 1439 l2:=q+1+t; 1440 while true do 1441 if not l1*P = O then 1442 if not l2*P = O then 1443 return Error("there is no elliptic curve with j-invariant ",j," and size ",n," over ",Parent(j)); 1444 fi; 1445 E2:= QuadraticTwist(E); 1446 E2.size:=n; 1447 return E2; 1448 elif not l2*P = O then 1449 E.size:=n; 1450 return E; 1451 else P:=RandomPoint(E); 1452 fi; 1453 od; 1454 end; 1455 1456 1457 #subfunction for computing isogenies between elliptic curves over fields of characteristic 2 1458 Fillp:=function(p,alpha,beta,PAMult,d,x) 1459 local ok,k,k1,k2,elimvars,sum,sum1,sum2,i,ck,bk,f,coeffs,vars,T,elim,co,varelim,i2; 1460 ok:=false; 1461 k:=1; 1462 k1:=0; 1463 k2:=1; 1464 elimvars:=[]; 1465 while k<d-2*k do 1466 sum:=0; 1467 sum2:=0; 1468 for i in [1..k] do 1469 sum:=sum+(p[i]^2)*(p[d-k+i]^2)*alpha^(2*(i-1)); 1470 od; 1471 for i in [1..Floor(k/2)] do 1472 sum2:=sum2+p[k-2*i+1]*B(d-k+2*i,i); 1473 od; 1474 sum:=AllRoots(alpha,4)[1]*sum+(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k)*sum2); 1475 ck:=sum/(alpha^(2*k)*AllRoots(alpha,4)[1]); 1476 bk:=(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k))/(alpha^(2*k)*AllRoots(alpha,4)[1]); 1477 f:=ck/(bk^2); 1478 if Is(Type(f),elt-fld^fin) then 1479 f:=MultPol(PAMult,[0],f); 1480 fi; 1481 if Is(Type(f),elt-alg^mult) then 1482 coeffs:=Coefficients(f); 1483 vars:=Variables(f); 1484 T:=1; 1485 while T<> 0 do 1486 T:=0; 1487 for i in [1..Length(coeffs)] do 1488 if coeffs[i]=0 then 1489 coeffs[i]:=Coerce(GF(2),0); 1490 fi; 1491 if Trace(coeffs[i],GF(2)) = 1 then 1492 T:=T+MultPol(PAMult,vars[i],Coerce(GF(2),1)); 1493 fi; 1494 od; 1495 if T<>0 then 1496 elim :=false; 1497 coeffs:=Coefficients(T); 1498 vars:=Variables(T); 1499 i:=0; 1500 while i <Length(vars) do 1501 i:=i+1; 1502 if Length(vars[i]) = 1 and vars[i][1]<>0 then 1503 elim :=true; 1504 for i2 in [1..Length(vars)] do 1505 if i2 <> i then 1506 if vars[i][1] in vars[i2] then 1507 elim:=false; 1508 fi; 1509 fi; 1510 od; 1511 if elim then 1512 elimvars:=Append(elimvars,[vars[i][1]]); 1513 co:=coeffs[i]; 1514 varelim:=vars[i][1]; 1515 Remove_(vars,i); 1516 Remove_(coeffs,i); 1517 f:=MultPol(PAMult,vars,coeffs)/co; 1518 for i2 in [1..Length(p)] do 1519 if (Type(GetEntry(p,i2))=elt-fld^fin or Type(GetEntry(p,i2))=elt-alg^mult) then 1520 if Is(Type(p[i2]),elt-alg^mult) then 1521 p[i2]:=Elimination_(p[i2],varelim,f); 1522 fi; 1523 fi; 1524 od; 1525 i:=Length(vars)+1; 1526 fi; 1527 fi; 1528 od; 1529 if not elim then 1530 return Error("not yet implemented"); 1531 fi; 1532 sum:=0; 1533 sum2:=0; 1534 for i in [1..k] do 1535 sum:=sum+(p[i]^2)*(p[d-k+i]^2)*alpha^(2*(i-1)); 1536 od; 1537 for i in [1..Floor(k/2)] do 1538 sum2:=sum2+p[k-2*i+1]*B(d-k+2*i,i); 1539 od; 1540 sum:=AllRoots(alpha,4)[1]*sum+(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k)*sum2); 1541 ck:=sum/(alpha^(2*k)*AllRoots(alpha,4)[1]); 1542 bk:=(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k))/(alpha^(2*k)*AllRoots(alpha,4)[1]); 1543 f:=ck/(bk^2); 1544 coeffs:=Coefficients(f); 1545 vars:=Variables(f); 1546 fi; 1547 od; 1548 sum:=0; 1549 for i in [1..Length(coeffs)] do 1550 sum:=sum+MultPol(PAMult,vars[i],Roots(x^2+x+coeffs[i])[1][1]); 1551 od; 1552 sum:=sum*bk; 1553 p[k+1]:=sum+MultPol(PAMult,[k],bk); 1554 k:=k+1; 1555 if not k>d-2*k+1 then 1556 k1:=k1+1; 1557 sum1:=0; 1558 sum2:=0; 1559 for i in [1..k1] do 1560 sum1:=sum1+p[d-2*k1+2*i]*B(d-2*k1-1+2*i,i)*alpha^(2*i); 1561 sum2:=sum2+p[d-2*k1+2*(i-1)+1]*B(d-2*k1+2*(i-1),i-1)*alpha^(2*(i-1)); 1562 od; 1563 sum2:=sum2+p[d+1]*B(d,k1)*alpha^(2*k1); 1564 p[d-2*k1]:= alpha^(4*k1+1-2*d)*p[k1+1]^4+sum1+alpha*sum2; 1565 if not k>d-2*k then 1566 k2:=k2+1; 1567 sum1:=0; 1568 sum2:=0; 1569 for i in [1..k2] do 1570 sum1:=sum1+p[d-2*k2+2*i]*B(d+1-2*k2+2*(i-1),(i-1))*alpha^(2*(i-1)); 1571 sum2:=sum2+p[d-2*k2+2*i+1]*B(d-2*k2+2*i,i)*alpha^(2*i); 1572 od; 1573 p[d-2*k2+1]:=p[d-k2+1]^4+alpha*sum1+sum2; 1574 fi; 1575 fi; 1576 fi; 1577 od; 1578 return [p,elimvars,k,k1,k2]; 1579 end; 1580 1581 #subfunction for computing isogenies between elliptic curves over fields of characteristic 2 1582 _Solve:=function(p,m,elimvars,M,BR,PAMult) 1583 local i,elim,j,le,f,L,c,easy,ellist,k,g,summ1,sum2,Ls,cs,pos,sols,sum1,sum2; 1584 i:=m; 1585 elim:=[]; 1586 while i >1 do 1587 while i in elimvars do 1588 i:=i-1; 1589 m:=m-1; 1590 od; 1591 j:=0; 1592 le:=m-i; 1593 if i >1 then 1594 while j <Length(M) and (Length(elim) = le) do 1595 j:=j+1; 1596 f:=M[j]; 1597 L:=Variables(f); 1598 if i in Last(L) then 1599 c:=Coefficients(f); 1600 easy:=true; 1601 ellist:=[]; 1602 for k in [1..Length(L)-1] do 1603 if i in L[k] then 1604 easy:=false; 1605 Append_(ellist,[k]); 1606 fi; 1607 od; 1608 if easy then 1609 if Length(Last(L)) = 1 then 1610 g:= f/Last(c); 1611 g:=MultPol(PAMult,Remove(Variables(g),Length(L)),Remove(Coefficients(g),Length(L))); 1612 else return Error("case not yet implemented"); 1613 fi; 1614 for k in [1..Length(M)] do 1615 M[k]:=Elimination_(M[k],Last(L)[1],g); 1616 od; 1617 Append_(elim,[[Last(L)[1],g]]); 1618 Remove_(M,j); 1619 i:=i-1; 1620 else 1621 Append_(ellist,[k+1]); 1622 sum1:=0; 1623 sum2:=0; 1624 for k in [1..Length(ellist)] do 1625 sum1:= sum1+MultPol(PAMult,L[ellist[k]],c[ellist[k]]); 1626 od; 1627 sum2:=MultPol(PAMult,Variables(f),Coefficients(f))-sum1; 1628 Ls:=Variables(sum1); 1629 cs:=Coefficients(sum1); 1630 sum1:=0; 1631 for k in [1..Length(Ls)] do 1632 pos:=Position(Ls[k],i); 1633 sum1:=sum1+MultPol(PAMult,Remove(Ls[k],pos),cs[k]); 1634 od; 1635 g:=MultPol(PAMult,Variables(sum2),Coefficients(sum2),sum1); 1636 for k in [1..Length(M)] do 1637 M[k]:=Elimination_(M[k],Last(L)[1],g); 1638 od; 1639 Append_(elim,[[Last(L)[1],g]]); 1640 Remove_(M,j); 1641 i:=i-1; 1642 fi; 1643 else 1644 return Error("not yet implemented"); 1645 fi; 1646 od; 1647 fi; 1648 od; 1649 Assert(i=1," all variables except for one eliminated "); 1650 sols:=[0,1]; 1651 for i in [1..Length(M)] do 1652 if sols = [] then 1653 return Error(" no solution found "); 1654 fi; 1655 j:=0; 1656 while j <Length(sols) do 1657 j:=j+1; 1658 if "denom" in RecFields(M[i]) then 1659 if Coerce(BR,Elimination_(M[i],1,sols[j]))/Coerce(BR,Elimination_(M[i].denom,1,sols[j])) <> 0 then 1660 Remove_(sols,j); 1661 j:=j-1; 1662 fi; 1663 else 1664 if Coerce(BR,Elimination_(M[i],1,sols[j])) <> 0 then 1665 Remove_(sols,j); 1666 j:=j-1; 1667 fi; 1668 fi; 1669 od; 1670 od; 1671 Assert(Length(sols) =1 , " only one solution found "); 1672 sols[1]:=Tuple([1,sols[1]]); 1673 i:=Length(elim); 1674 for j in [1..Length(elimvars)] do 1675 sols[elimvars[j]]:= Tuple([elimvars[j],"nothing"]); 1676 od; 1677 while i >=1 do 1678 k:=elim[i][1]; 1679 f:=elim[i][2]; 1680 if "denom" in RecFields(f) then 1681 j:=k-1; 1682 while j in elimvars do 1683 j:=j-1; 1684 od; 1685 while j >1 do 1686 while j in elimvars do 1687 j:=j-1; 1688 od; 1689 g:=Elimination_(f,j,sols[j][2]); 1690 f:=MultPol(PAMult,Variables(g),Coefficients(g),Elimination_(f.denom,j,sols[j][2])); 1691 j:=j-1; 1692 od; 1693 f:=Coerce(BR,Elimination_(f,1,sols[1][2])); 1694 Assert(f = 1 or f = 0 , " right solution"); 1695 sols[k]:=Tuple([k,f]); 1696 else 1697 j:=k-1; 1698 while j in elimvars do 1699 j:=j-1; 1700 od; 1701 while j >1 do 1702 while j in elimvars do 1703 j:=j-1; 1704 od; 1705 f:=Elimination_(f,j,sols[j][2]); 1706 j:=j-1; 1707 od; 1708 f:=Coerce(BR,Elimination_(f,1,sols[1][2])); 1709 Assert(f = 1 or f = 0 , " right solution"); 1710 sols[k]:=Tuple([k,f]); 1711 fi; 1712 i:=i-1; 1713 od; 1714 for i in [1..Length(p)] do 1715 while not IsConstant(p[i]) do 1716 k:=Last(Variables(p[i]))[1]; 1717 p[i]:=Elimination_(p[i],k,sols[k][2]); 1718 od; 1719 od; 1720 return p; 1721 end; 1722 1723 1724 #computes isogeny of degree l, starting from EE, to curve whose j-invariant is the m-th root of the l-th modular polynomial, characteristic 2, Tr(L[3])=0 1725 IsoFromJChar2new:=function(EE,l,m) 1726 local d,L,size,E,Iso1,E4,BR,s,PA,x,PA2,mo,j,ro,PAMult,E2,b,a,alpha,beta,j2,mytime, p,elimvars,k,k1,k2,M,sum1,sum2,i,q,Q,P,spliq,testP,G,K,R,H,xmap,ymap,_Iso,re,PP,fr,lambda; 1727 if not IsOdd(l) or not IsPrime(l) then 1728 return Error("Algorithm implemented only for odd primes"); 1729 fi; 1730 d:=Coerce(Integers(),(l-1)/2); 1731 L:=Coefficients(EE); 1732 size:=EE.size; 1733 E:=NormalFormChar2(EE,1); 1734 Iso1:=E.ext1; 1735 E4:=Base(E); 1736 L:= Coefficients(E4); 1737 BR:=BaseRing(EE); 1738 s:=Size(BR); 1739 E:=E4; 1740 E.size:=size; 1741 PA:=PolynomialAlgebra(BR); 1742 x:=PA.1; 1743 PA2:=PolynomialAlgebra(PA); 1744 mo:=Coerce(PA2,ModularPolynomial(l)); 1745 j:=jInvariant(E); 1746 ro:=Roots(Evaluate(mo,j)); 1747 if Length(ro) = 0 then 1748 return Error("There is no ",l,"-isogenous curve"); 1749 fi; 1750 if m =1 then 1751 j2:=ro[1][1]; 1752 elif Length(ro) <m then 1753 return Error("The modular polynomial does not have k roots"); 1754 else j2:= ro[m][1]; 1755 fi; 1756 PAMult:=MultPolAlg(BR,Floor(d/3)); 1757 E2:=EllipticCurveFromJInvariant(j2); 1758 b:=j2^(-1); 1759 L:= Coefficients(E); 1760 Assert(L[1]=1 and L[2]=0 and L[3]=0 and L[4] = 0 and L[5] <>0, " curve in normal form"); 1761 a:= L[5]; 1762 alpha:=AllRoots(a,4)[1]; 1763 beta:=AllRoots(b,4)[1]; 1764 p:=[]; 1765 p[d+1]:=1; 1766 p[d]:=alpha+beta; 1767 if d >1 then 1768 if IsOdd(d) then 1769 p[d-1]:= p[d]^4+alpha*p[d]+alpha^2; 1770 else 1771 p[d-1]:= p[d]^4+alpha*p[d]; 1772 fi; 1773 fi; 1774 p[1]:=AllRoots(alpha^(2*d)+alpha^(2*d-1)*p[d],4)[1]; 1775 p:=Fillp(p,alpha,beta,PAMult,d,x); 1776 elimvars:=p[2]; 1777 k:=p[3]; 1778 k1:=p[4]; 1779 k2:=p[5]; 1780 p:=p[1]; 1781 M:=[]; 1782 m:=k-1; 1783 if m > 0 then 1784 for k in [k1+1..Floor((d-1)/2)] do 1785 sum1:=0; 1786 sum2:=0; 1787 for i in [1..k+1] do 1788 sum1:=sum1+p[d-2*k+2*(i-1)]*B(d-2*k-1+2*(i-1),i-1)*alpha^(2*(i-1)); 1789 sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 1790 od; 1791 Append_(M,[p[k+1]^4+alpha^(2*d-4*(k)-1)*sum1+alpha^(2*d-4*(k))*sum2]); 1792 od; 1793 for k in [k2+1..Floor(d/2)] do 1794 sum1:=0; 1795 sum2:=0; 1796 for i in [1..k] do 1797 sum1:=sum1+p[d-2*k+2*i]*B(d+1-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 1798 sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 1799 od; 1800 sum2:= sum2+p[d+1]*B(d,k)*alpha^(2*k); 1801 Append_(M,[p[d-k+1]^4+alpha*sum1+sum2]); 1802 od; 1803 Assert(IsDense(p)," all coefficients of p computed "); 1804 p:=_Solve(p,m,elimvars,M,BR,PAMult); 1805 fi; 1806 q:=[]; 1807 q[d+1]:=1; 1808 for i in [1..d+1] do 1809 q[i]:=Coerce(BR,(AllRoots(alpha,4)[1]/AllRoots(beta,4)[1])*Sqrt(alpha)^(d-2*(i-1))*p[d-(i-1)+1]); 1810 od; 1811 Q:=q[1]^2; 1812 for i in [2..d+1] do 1813 Q:=Q+q[i]^2*x^(i-1); 1814 od; 1815 P:=Coerce(BR,p[1])^2; 1816 for i in [2..d+1] do 1817 P:=P+Coerce(BR,p[i])^2*x^(i-1); 1818 od; 1819 spliq:=SplittingField(Q); 1820 testP:=MakePoints(ChangeBaseRing(E,Extension(spliq,2)),HasRoot(Q,spliq).ext1)[1]; 1821 G:=x*P^2; 1822 K:=(P^2)*Q; 1823 R:=x*Derivative(P*Q); 1824 H:=x*R*P+Sqrt(b)*Q^3+Sqrt(a)*P^2*Q; 1825 xmap:=function(P) 1826 local xx,q; 1827 xx:=PX(P); 1828 q:=Evaluate(Q^2,xx); 1829 if q = 0 then 1830 return "INFTY"; 1831 else 1832 return Evaluate(G,xx)/q; 1833 fi; 1834 end; 1835 ymap:= function(P) 1836 local xx,yy; 1837 xx:=PX(P); 1838 yy:=PY(P); 1839 return (Evaluate(H,xx)+yy*Evaluate(K,xx))/Evaluate(Q^3,xx); 1840 end; 1841 E2:=ChangeBaseRing(E2,BR); 1842 E2.size:=E.size; 1843 _Iso:=function(P) 1844 local x; 1845 if P = Zero(E) then 1846 return Zero(E2); 1847 else 1848 x:=xmap(P); 1849 if x = "INFTY" then 1850 return Zero(E2); 1851 else 1852 return MakePoint(E2,[x,ymap(P)]); 1853 fi; 1854 fi; 1855 end; 1856 re:= Map(E,E2,_Iso); 1857 re.type:= iso; 1858 re.operations:= _iso_ops; 1859 PP:=GenericPoint(E); 1860 re.xpol:=xmap(PP); 1861 re.ypol:=ymap(PP); 1862 re.degree:= l; 1863 i:= Degree(BaseRing(Parent(testP)),BaseRing(E)); 1864 fr:= MakePoint(Parent(testP),[PX(testP)^s,PY(testP)^s]); 1865 ro:= Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E))); 1866 lambda:= Coerce(Integers(),ro[1][1]); 1867 if not (lambda*testP=fr) then 1868 lambda:=Coerce(Integers(),ro[2][1]); 1869 Assert(lambda*testP = fr, " eigenvalue found "); 1870 fi; 1871 re.ei:= lambda; 1872 re:=Composition(re,Iso1); 1873 Codomain(re).size:=EE.size; 1874 return re; 1875 end; 1876 1877 1878 #computes isogeny of degree l, starting from EE, to curve whose j-invariant is the m-th root of the l-th modular polynomial, characteristic 2 1879 IsoFromJChar2:=function(EE,l,m) 1880 local d,L,E,Iso1,E3,BR,s,PA,x,PA2,mo,j,ro,j2,PAMult,E2,b,alpha,beta,p,a,k,k1,k2,elimvars,M,m,sum1,sum2,i,q,Q,K,P,Iso11,Iso2,Iso,re,Kxy,Alf,defi,Alf2,fr,lambda; 1881 if not IsOdd(l) or not IsPrime(l) then 1882 return Error("Algorithm implemented only for odd primes"); 1883 fi; 1884 d:=Coerce(Integers(),(l-1)/2); 1885 L:=Coefficients(EE); 1886 if Trace(L[3],GF(2)) = 0 then 1887 return IsoFromJChar2new(EE,l,m); 1888 else 1889 E:=NormalFormChar2(EE,2); 1890 fi; 1891 Iso1:=E.ext1; 1892 E3:=Base(E); 1893 L:= Coefficients(E3); 1894 E2:=QuadraticTwist(E3); 1895 E:=E2; 1896 BR:=BaseRing(E); 1897 s:=Size(BR); 1898 PA:=PolynomialAlgebra(BR); 1899 x:=PA.1; 1900 PA2:=PolynomialAlgebra(PA); 1901 mo:=Coerce(PA2,ModularPolynomial(l)); 1902 j:=jInvariant(E); 1903 ro:=Roots(Evaluate(mo,j)); 1904 if Length(ro) = 0 then 1905 return Error("There is no ",l,"-isogenous curve"); 1906 fi; 1907 if m =1 then 1908 j2:=ro[1][1]; 1909 elif Length(ro) <m then 1910 return Error("The modular polynomial does not have k roots"); 1911 else 1912 j2:= ro[m][1]; 1913 fi; 1914 PAMult:=MultPolAlg(BR,Floor(d/3)); 1915 E2:=EllipticCurveFromJInvariant(j2); 1916 b:=j2^(-1); 1917 L:= Coefficients(E); 1918 Assert(L[1]=1 and L[2]=0 and L[3]=0 and L[4] = 0 and L[5] <>0, " Form of the curve "); 1919 a:= L[5]; 1920 alpha:=AllRoots(a,4)[1]; 1921 beta:=AllRoots(b,4)[1]; 1922 p:=[]; 1923 p[d+1]:=1; 1924 p[d]:=alpha+beta; 1925 if d >1 then 1926 if IsOdd(d) then 1927 p[d-1]:= p[d]^4+alpha*p[d]+alpha^2; 1928 else 1929 p[d-1]:= p[d]^4+alpha*p[d]; 1930 fi; 1931 fi; 1932 p[1]:=AllRoots(alpha^(2*d)+alpha^(2*d-1)*p[d],4)[1]; 1933 p:=Fillp(p,alpha,beta,PAMult,d,x); 1934 elimvars:=p[2]; 1935 k:=p[3]; 1936 k1:=p[4]; 1937 k2:=p[5]; 1938 p:=p[1]; 1939 M:=[]; 1940 m:=k-1; 1941 if m > 0 then 1942 for k in [k1+1..Floor((d-1)/2)] do 1943 sum1:=0; 1944 sum2:=0; 1945 for i in [1..k+1] do 1946 sum1:=sum1+p[d-2*k+2*(i-1)]*B(d-2*k-1+2*(i-1),i-1)*alpha^(2*(i-1)); 1947 sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 1948 od; 1949 Append_(M,[p[k+1]^4+alpha^(2*d-4*(k)-1)*sum1+alpha^(2*d-4*(k))*sum2]); 1950 od; 1951 for k in [k2+1..Floor(d/2)] do 1952 sum1:=0; 1953 sum2:=0; 1954 for i in [1..k] do 1955 sum1:=sum1+p[d-2*k+2*i]*B(d+1-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 1956 sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 1957 od; 1958 sum2:= sum2+p[d+1]*B(d,k)*alpha^(2*k); 1959 Append_(M,[p[d-k+1]^4+alpha*sum1+sum2]); 1960 od; 1961 Assert(IsDense(p)," all coefficients of p computed "); 1962 p:=_Solve(p,m,elimvars,M,BR,PAMult); 1963 fi; 1964 q:=[]; 1965 q[d+1]:=1; 1966 for i in [1..d+1] do 1967 q[i]:=Coerce(BR,(AllRoots(alpha,4)[1]/AllRoots(beta,4)[1])*Sqrt(alpha)^(d-2*(i-1))*p[d-(i-1)+1]); 1968 od; 1969 Q:=q[1]^2; 1970 for i in [2..d+1] do 1971 Q:=Q+q[i]^2*x^(i-1); 1972 od; 1973 K:=Extension(SplittingField(Q),2); 1974 ro:=HasRoot(Q,K); 1975 P:=MakePoints(ChangeBaseRing(E,K),ro.ext1)[1]; 1976 Iso11:=NormalFormChar2(ChangeBaseRing(EE,K),2); 1977 Iso2:=Composition(Iso11.ext1,Isomorphism(ChangeBaseRing(E,K),Base(Iso11),1)); 1978 P:=Iso2(P); 1979 E:= EE; 1980 Assert(l*P=Zero(Parent(P))," successful computation of Q "); 1981 Iso:=IsogenyFromSubgroup(P,E,l); 1982 E3:=Codomain(Iso); 1983 re:= Map(ChangeBaseRing(E,BR),ChangeBaseRing(Codomain(Iso),BR),Base(Iso)); 1984 re.type:= iso; 1985 re.operations:= Iso.operations; 1986 Kxy:= PolynomialAlgebra(PolynomialAlgebra(BR)); 1987 Alf:= Parent(Iso.xpol); 1988 defi:= DefiningPolynomial(Alf); 1989 Alf2:= FunctionField(Coerce(Kxy,defi)); 1990 re.xpol:= Coerce(Alf2,Eltseq(Iso.xpol)); 1991 re.ypol:= Coerce(Alf2,Eltseq(Iso.ypol)); 1992 re.degree:= Iso.degree; 1993 i:= Degree(BaseRing(Parent(P)),BaseRing(E)); 1994 fr:= MakePoint(Parent(P),[PX(P)^s,PY(P)^s]); 1995 ro:= Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E))); 1996 lambda:= Coerce(Integers(),ro[1][1]); 1997 if not (lambda*P=fr) then 1998 lambda:=Coerce(Integers(),ro[2][1]); 1999 Assert(lambda*P = fr, " eigenvalue found "); 2000 fi; 2001 re.ei:= lambda; 2002 if not InNormalFormChar2(E3) then 2003 re:=Composition(NormalFormChar2(E3,1).ext1,re); 2004 Codomain(re).size:=EE.size; 2005 return re; 2006 fi; 2007 re:=Composition(re,Iso1); 2008 Codomain(re).size:=EE.size; 2009 return re; 2010 end; 2011 2012 2013 #computes isogeny of degree 2, starting from E to curve whose j-invariant is the k-th root of the second modular polynomial 2014 IsoFromJ2:= function(E,k) #berechnet Isogenie vom Grad 2 2015 local to,Iso,re,Kxy,Alf,defi,Alf2,l,P,i,R,fr,fro,lambda,eifound,BR,f,E2,j,found,m,mo,fo,ro,fro; 2016 to:= CompleteTorsionSubgroup(E,2); 2017 mo:= ModularPolynomial(2); 2018 j:=Roots(Evaluate(Coerce(PolynomialAlgebra(PolynomialAlgebra(BaseRing(E))),mo),jInvariant(E)))[k][1]; 2019 for i in [1..Length(to)] do 2020 P:=to[i]; 2021 fro:=extFrob(E,Degree(BaseRing(Parent(P)),BaseRing(E))); 2022 if fro(P) = P then 2023 Iso:= IsogenyFromSubgroup(P,E,2); 2024 if jInvariant(Codomain(Iso)) = j then 2025 Iso.ei:=1; 2026 return Iso; 2027 fi; 2028 fi; 2029 od; 2030 Iso:= FrobeniusToOtherCurve(E,2); 2031 Iso.ei:=0; 2032 Codomain(Iso).size:=E.size; 2033 return Iso; 2034 end; 2035 2036 2037 #subfunction for endomorphism ring computation, computes conductor of End(E) 2038 FindEndRingIndexFromEqOrder:= function(E) 2039 local jE,f,O,i,l,fact,a,m,ind,d,divpol,x,BR,divap,divam,divas,h,las,p,j,r,t,atmp,modpoly,PA,g,nr,ro,n,co,n1,n2,mi,n3,prev,p,n2set,n,to,fro,k,atmp,P,done; 2040 f:= CharacteristicFrobeniusPolynomial(E); 2041 t:= -Coefficient(f,1); 2042 O:= EquationOrder(f); 2043 m:= Conductor(O); 2044 if m = 1 then 2045 return 1; 2046 fi; 2047 d:= Discriminant(MaximalOrder(O)); 2048 if d mod 4 = 0 then 2049 a:= t/2; 2050 elif d mod 4 = 1 then 2051 a:= (t+m)/2; 2052 else return Error("something went wrong"); 2053 fi; 2054 a:= Coerce(Integers(),a); 2055 ind:=1; 2056 fact:= Factorization(m); 2057 BR:= BaseRing(E); 2058 x:= PolynomialAlgebra(BR).1; 2059 p:= Size(BR); 2060 PA:= PolynomialRing(PolynomialRing(BR)); 2061 jE:= jInvariant(E); 2062 prev:=[]; 2063 n2set:=false; 2064 for i in [1..Length(fact)] do 2065 l:= fact[i][1]; 2066 r:= fact[i][2]; 2067 n:=0; 2068 j:=jE; 2069 modpoly:= Coerce(PA,ModularPolynomial(l)); 2070 las:=false; 2071 while n <=r and not las do 2072 prev:=Append(prev,[j]); 2073 n:=n+1; 2074 g:= Evaluate(modpoly,j); 2075 ro:= Roots(g); 2076 nr:= 0; 2077 for co in [1..Length(ro)] do 2078 nr:= nr+ ro[co][2]; 2079 od; 2080 if nr <> l+1 then 2081 las:= true; 2082 n:=n-1; 2083 else 2084 j:= ro[1][1]; 2085 p:=1; 2086 while (j in prev or j =0) and not las do 2087 p:= p+1; 2088 if p>Length(ro) then 2089 las:=true; 2090 n:=n-1; 2091 else 2092 j:= ro[p][1]; 2093 fi; 2094 od; 2095 fi; 2096 od; 2097 n1:=n; 2098 n:=0; 2099 j:=jE; 2100 las:=false; 2101 prev:=[]; 2102 while n <=r and not las do 2103 prev:=Append(prev,[j]) ; 2104 n:=n+1; 2105 g:= Evaluate(modpoly,j); 2106 ro:= Roots(g); 2107 nr:= 0; 2108 for co in [1..Length(ro)] do 2109 nr:= nr+ ro[co][2]; 2110 od; 2111 if nr <> l+1 then 2112 las:= true; 2113 n:=n-1; 2114 else 2115 j:= ro[2][1]; 2116 p:= 1; 2117 while (j in prev or j = 0) and n2set=false do 2118 p:=p+1; 2119 if p>Length(ro) then 2120 n2:=n1; 2121 n2set:=true; 2122 else 2123 j:= ro[p][1]; 2124 fi; 2125 od; 2126 fi; 2127 od; 2128 if not n2set then 2129 n2:= n; 2130 fi; 2131 if not (n1>r and n2>r) then 2132 mi:= Minimum([n1,n2]); 2133 if not mi = 0 then las:= true; 2134 ind:=ind*l^mi; 2135 fi; 2136 else ind:= ind*(l^r); 2137 fi; 2138 od; 2139 return Coerce(Integers(),m/ind); 2140 end; 2141 2142 2143 #computes End(E) 2144 EndRing:= function(E) 2145 local c,Omax,f,g,O,OE,alpha,beta,fact,D,x; 2146 if IsSupersingular(E) then 2147 return Error("Cannot compute Endomorphism Ring for supersingular curve"); 2148 fi; 2149 f:= CharacteristicFrobeniusPolynomial(E); 2150 O:= EquationOrder(f); 2151 c:=Conductor(O); 2152 if not c = 1 then 2153 fact:=Factorization(c); 2154 if Last(fact)[1] > 59 then 2155 return FAILURE; 2156 fi; 2157 fi; 2158 c:= FindEndRingIndexFromEqOrder(E); 2159 Omax:= MaximalOrder(O); 2160 OE:=EquationOrder(MinimalPolynomial(c*Omax.2)); 2161 Assert(Conductor(OE) = c, "Order has expected conductor"); 2162 return OE; 2163 end; 2164 2165 2166 2167 #computes isogeny of degree l, starting from E, to curve whose j-invariant is k-th root of l-th modular polynomial 2168 IsogenyFromJInvariant:= function(E,l,k) 2169 local sum,jTilDer,p1,coeffs,k,PA2,m,C,A,K,CFun,ro,AFun,unf,gen,bo,f, xval,sub,BR,jDer,a,b,ModPoly,PA,j2,j,E2,a,b,a2,d,b2,Phix,Phiy,Phixx,Phiyy,Phixy,J,E4,E6,E4Til,E6Til,c,cTil,i,j,sum1,sum2,divpolyfact,insum,Kxy,Alf,Alf2,defi,Iso,re,R,fro,fr,P,lambda,eifound,E2,i,ro,r, mytime,esize,xx,ff,kk,sumA,jj,f2; 2170 esize:=Size(E); 2171 if l = 2 then 2172 return IsoFromJ2(E,k); 2173 fi; 2174 if Characteristic(BaseRing(E)) = 2 then 2175 return IsoFromJChar2(E,l,k); 2176 fi; 2177 if IsSupersingular(E) then 2178 return Error("cannot compute isogeny for supersingular curve"); 2179 fi; 2180 Assert(Length(Coefficients(E)) =2, "curve in short WeierstrassNormalForm"); 2181 a:= E.L[1]; 2182 b:= E.L[2]; 2183 c:=[]; 2184 cTil:=[]; 2185 BR:= BaseRing(E); 2186 j:=jInvariant(E); 2187 if j = Coerce(BR,0) or j = Coerce(BR,1728) then 2188 return Error("cannot compute isogeny for j-Invariant ",j); 2189 fi; 2190 PA:=PolynomialAlgebra(BR); 2191 PA2:= PolynomialAlgebra(PA); 2192 ModPoly:= Coerce(PA2,ModularPolynomial(l)); 2193 ro:=Roots(Evaluate(ModPoly,j));#Print(ro, "\n"); 2194 if Length(ro) = 0 then 2195 return Error("There is no ",l,"-isogenous curve"); 2196 fi; 2197 E:= ChangeBaseRing(E,BR); 2198 E.size:=esize; 2199 PA:= PolynomialAlgebra(BR); 2200 if k =1 then 2201 j2:=ro[1][1]; 2202 elif Length(ro) <k then 2203 return Error("The modular polynomial does not have k roots"); 2204 else 2205 j2:= ro[k][1]; 2206 fi; 2207 if j2 = 0 or j2 = 1728 then 2208 return FAILURE; 2209 fi; 2210 Phix:= Derivative(ModPoly); 2211 if Evaluate(Evaluate(Phix,j),j2) = 0 then 2212 return FAILURE; 2213 fi; 2214 Phiy:=(Derivative(_SwapVars(ModPoly))); 2215 Phiyy:= _SwapVars(Derivative(Phiy)); 2216 Phiy:= _SwapVars(Phiy); 2217 Phixx:= Derivative(Phix); 2218 Phixy:= _SwapVars(Derivative(_SwapVars(Phix))); 2219 E4:=-48*a; 2220 E6:= 864*b; 2221 jDer:= -j*E6/E4; 2222 jTilDer:= -((jDer*Evaluate(Evaluate(Phix,j),j2))/(l*Evaluate(Evaluate(Phiy,j),j2))); 2223 jTilDer:= Coerce(BR,jTilDer); 2224 a2:= -(jTilDer^2)/(48*j2*(j2-1728)); 2225 b2:= -(jTilDer^3)/(864*j2^2*(j2-1728)); 2226 E4Til:= -48*a2; 2227 E6Til:= 864*b2; 2228 J:=(jDer^2*Evaluate(Evaluate(Phixx,j),j2)); 2229 J:= J+2*l*jDer*jTilDer*Evaluate(Evaluate(Phixy,j),j2); 2230 J:= J+l^2*jTilDer^2*Evaluate(Evaluate(Phiyy,j),j2); 2231 J:= -J/(jDer*Evaluate(Evaluate(Phix,j),j2)); 2232 p1:= l/2*J+l/4*(E4^2/E6-l*E4Til^2/E6Til)+l/3*(E6/E4-l*E6Til/E4Til); 2233 E2:=EllipticCurve(BR,[l^4*a2,l^6*b2]); 2234 c[1]:=-a/5; 2235 cTil[1]:=-(l^4*a2)/5; 2236 c[2]:= -b/7; 2237 cTil[2]:= -(l^6*b2)/7; 2238 d:= Coerce(Integers(),(l-1)/2); 2239 for i in [3..d+1] do 2240 sum1 :=0; 2241 sum2:=0; 2242 for j in [1..i-2] do 2243 sum1:= sum1 + c[j]*c[i-1-j]; 2244 sum2:= sum2 + cTil[j]*cTil[i-1-j]; 2245 od; 2246 c[i] := sum1*3/((i-2)*(2*i+3)); 2247 cTil[i] := sum2*3/((i-2)*(2*i+3)); 2248 od; 2249 f:= Coerce(PA,Append([0],c)); 2250 C:= function(k,j) 2251 local f2,i,elseq; 2252 f2:=f^(k-j); 2253 elseq:= Eltseq(f2); 2254 if elseq=Sequence([]) then 2255 return 0; 2256 elif Length(elseq)<=j then 2257 return 0; 2258 else 2259 return elseq[j+1]; 2260 fi; 2261 end; 2262 A:=[]; 2263 xx:= PA.1; 2264 ff:=-1/2*p1*xx; 2265 for kk in [1..d] do 2266 ff:= ff-(cTil[kk]-l*c[kk])/((2*kk+1)*(2*kk+2))*xx^(kk+1); 2267 od; 2268 sumA:=0; 2269 for jj in [0..d] do 2270 sumA:= sumA+ff^jj/(Factorial(jj)); 2271 od; 2272 A:= Eltseq(sumA); 2273 coeffs:=[]; 2274 coeffs[d+1]:=1; 2275 for i in [1..d] do 2276 m:= d-i+1; 2277 sum:=0; 2278 for k in [1..i] do 2279 insum:=0; 2280 for j in [0..k] do 2281 insum:= insum+Binomial(d-i+k,k-j)*C(k,j); 2282 od; 2283 sum:=sum+insum*coeffs[d-i+k+1]; 2284 od; 2285 coeffs[m]:=A[i+1]-sum; 2286 od; 2287 divpolyfact := Coerce(PA,coeffs); 2288 if Characteristic(BaseRing(E))>6*d+1 then 2289 return IsogenyFromKernel(E,E2,divpolyfact); 2290 fi; 2291 K := SplittingField( divpolyfact ); 2292 PA := PolynomialAlgebra(K); 2293 xval:= HasRoot(Coerce(PA,divpolyfact)).ext1; 2294 E2:= ChangeBaseRing(E,Extension(K,2)); 2295 E:= ChangeBaseRing(E,BR); 2296 E.size:=esize; 2297 P:=MakePoints(E2,xval)[1]; 2298 Assert( l * P = Zero(E2), "l * P = O"); 2299 Iso:=IsogenyFromSubgroup(P,E,l); 2300 re:= Map(E,ChangeBaseRing(Codomain(Iso),BR),Base(Iso)); 2301 Codomain(re).size:=esize; 2302 re.type:= iso; 2303 re.operations:= Iso.operations; 2304 Kxy:= PolynomialAlgebra(PolynomialAlgebra(BR)); 2305 Alf:= Parent(Iso.xpol); 2306 defi:= DefiningPolynomial(Alf); 2307 Alf2:= FunctionField(Coerce(Kxy,defi)); 2308 re.xpol:= Coerce(Alf2,Eltseq(Iso.xpol)); 2309 re.ypol:= Coerce(Alf2,Eltseq(Iso.ypol)); 2310 re.degree:= Iso.degree; 2311 i:= Degree(BaseRing(Parent(P)),BaseRing(E)); 2312 if i =1 then 2313 re.ei:=1; 2314 else 2315 fro:= extFrob(E,i); 2316 fr:= fro(P); 2317 ro:= Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E))); 2318 lambda:= Coerce(Integers(),ro[1][1]); 2319 if not (lambda*P=fr) then 2320 lambda:=Coerce(Integers(),ro[2][1]); 2321 Assert(lambda*P = fr, " eigenvalue found "); 2322 fi; 2323 re.ei:= lambda; 2324 fi; 2325 Codomain(re).size:=E.size; 2326 return re; 2327 end; 2328 2329 2330 #computes eigenvalue of frobenius on subgroup in kernel of isogeny of degree 2, starting from E to curve whose j-invariant is the k-th root of the second modular polynomial 2331 SubgroupFromJ2:=function(E,k) 2332 local to,j,j2,ModPoly,PA,PA2,ro,Iso,BR,i; 2333 to:= TorsionSubgroup(E,2); 2334 j:=jInvariant(E); 2335 BR:=BaseRing(E); 2336 ModPoly:= ModularPolynomial(2); 2337 PA:=PolynomialAlgebra(BR); 2338 PA2:= PolynomialAlgebra(PA); 2339 ModPoly:= Coerce(PA2,ModPoly); 2340 ro:=Roots(Evaluate(ModPoly,j)); 2341 j2:= ro[k][1]; 2342 i:=0; 2343 while i<(Length(to) - 1) do 2344 i:=i+1; 2345 Iso:=IsogenyFromSubgroup(to[i],E,2); 2346 if jInvariant(Codomain(Iso)) = j2 then 2347 return [2,1]; 2348 fi; 2349 od; 2350 if GCD(2,Characteristic(BR))<>1 then 2351 Iso:=FrobeniusToOtherCurve(E,2); 2352 if jInvariant(Codomain(Iso)) = j2 then 2353 return [2,0]; 2354 fi; 2355 fi; 2356 return Error("something went wrong"); 2357 end; 2358 2359 #computes eigenvalue of frobenius on subgroup in kernel of isogeny of degree l, starting from E to curve whose j-invariant is the m-th root of the l-th modular polynomial, characteristic of base ring is 2 2360 SubFromJChar2:=function(EE,l,m) 2361 local L,a,j,mo,PA,PA2,BR,j2,ro,K,k1,P,k2,d,E2,b,q,alpha,beta,x,f,ok,k,i,g,pi,b1,xlist,PAtmp,xtmp,q,Q,Iso,E3,p,sum,sum2,co,ck,bk,PAMult,coeffs,T,vars,fel,i2,elim,sum1,M,PXM,V,W,n,pi,L2,c,easy,ellist,sols,le,Ls,cs,pos,re,Kxy,Alf,defi,Alf2,fr,sub,fro,lambda,r,s,E,quad,Iso1,Iso2,E4,Iso11,EC,elimlist,elimvars,mytime,ro,lambda,lambdax,li,de,nu,roro,poi,xval,gen; 2362 if not IsOdd(l) or not IsPrime(l) then 2363 return Error("Algorithm implemented only for odd primes"); 2364 fi; 2365 d:=Coerce(Integers(),(l-1)/2); 2366 L:=Coefficients(EE); 2367 if Trace(L[3],GF(2)) = 0 then 2368 E:=NormalFormChar2(EE,1); 2369 else 2370 E:=NormalFormChar2(EE,2); 2371 fi; 2372 Iso1:=E.ext1; 2373 E4:=Base(E); 2374 L:= Coefficients(E4); 2375 quad:=L[3]<>0; 2376 if quad then 2377 quad:=-1; 2378 E2:=QuadraticTwist(E4); 2379 E:=E2; 2380 else 2381 E:=E4; 2382 quad:=1; 2383 fi; 2384 BR:=BaseRing(E); 2385 s:=Size(BR); 2386 PA:=PolynomialAlgebra(BR); 2387 x:=PA.1; 2388 PA2:=PolynomialAlgebra(PA); 2389 mo:=Coerce(PA2,ModularPolynomial(l)); 2390 j:=jInvariant(E); 2391 ro:=Roots(Evaluate(mo,j)); 2392 if Length(ro) = 0 then 2393 return Error("There is no ",l,"-isogenous curve"); 2394 fi; 2395 if m =1 then 2396 j2:=ro[1][1]; 2397 elif Length(ro) <m then 2398 return Error("The modular polynomial does not have k roots"); 2399 else j2:= ro[m][1]; 2400 fi; 2401 EC:=ECFromJGivenNumber(j2,Size(EE)); 2402 PAMult:=MultPolAlg(BR,Floor(d/3)); 2403 E2:=EllipticCurveFromJInvariant(j2); 2404 b:=j2^(-1); 2405 L:= Coefficients(E); 2406 if not L[1]=1 and L[2]=0 and L[3]=0 and L[4] = 0 and L[5] <>0 then 2407 return Error("Algorithm implemented for curve of the form Y^2+XY=X^3+a"); 2408 fi; 2409 a:= L[5]; 2410 alpha:=AllRoots(a,4)[1]; 2411 beta:=AllRoots(b,4)[1]; 2412 p:=[]; 2413 p[d+1]:=1; 2414 p[d]:=alpha+beta; 2415 if d>1 then 2416 if IsOdd(d) then 2417 p[d-1]:= p[d]^4+alpha*p[d]+alpha^2; 2418 else 2419 p[d-1]:= p[d]^4+alpha*p[d]; 2420 fi; 2421 fi; 2422 p[1]:=AllRoots(alpha^(2*d)+alpha^(2*d-1)*p[d],4)[1]; 2423 p:=Fillp(p,alpha,beta,PAMult,d,x); 2424 elimvars:=p[2]; 2425 k:=p[3]; 2426 k1:=p[4]; 2427 k2:=p[5]; 2428 p:=p[1]; 2429 M:=[]; 2430 m:=k-1; 2431 if m > 0 then 2432 for k in [k1+1..Floor((d-1)/2)] do 2433 sum1:=0; 2434 sum2:=0; 2435 for i in [1..k+1] do 2436 sum1:=sum1+p[d-2*k+2*(i-1)]*B(d-2*k-1+2*(i-1),i-1)*alpha^(2*(i-1)); 2437 sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 2438 od; 2439 Append_(M,[p[k+1]^4+alpha^(2*d-4*(k)-1)*sum1+alpha^(2*d-4*(k))*sum2]); 2440 od; 2441 for k in [k2+1..Floor(d/2)] do 2442 sum1:=0; 2443 sum2:=0; 2444 for i in [1..k] do 2445 sum1:=sum1+p[d-2*k+2*i]*B(d+1-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 2446 sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1)); 2447 od; 2448 sum2:= sum2+p[d+1]*B(d,k)*alpha^(2*k); 2449 Append_(M,[p[d-k+1]^4+alpha*sum1+sum2]); 2450 od; 2451 Assert(IsDense(p)," all coefficients of p computed "); 2452 p:=_Solve(p,m,elimvars,M,BR,PAMult); 2453 fi; 2454 q:=[]; 2455 q[d+1]:=1; 2456 for i in [1..d+1] do 2457 q[i]:=Coerce(BR,(AllRoots(alpha,4)[1]/AllRoots(beta,4)[1])*Sqrt(alpha)^(d-2*(i-1))*p[d-(i-1)+1]); 2458 od; 2459 Q:=q[1]^2; 2460 for i in [2..d+1] do 2461 Q:=Q+q[i]^2*x^(i-1); 2462 od; 2463 ro:=Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E))); 2464 if Length(ro) = 2 and ro[1][1] = - ro[2][1] then 2465 K := SplittingField( Q ); 2466 PA := PolynomialAlgebra(K); 2467 xval:= HasRoot(Coerce(PA,Q)).ext1; 2468 E2:= ChangeBaseRing(E,Extension(K,2)); 2469 P:=MakePoints(E2,xval)[1]; 2470 fro:=extFrob(E,Degree(BaseRing(E2),BR)); 2471 lambda:=Coerce(Integers(),ro[1][1]); 2472 if lambda*P = fro(P) then 2473 return [l,quad*lambda]; 2474 else 2475 return [l,quad*(Coerce(Integers(),ro[2][1]))]; 2476 fi; 2477 fi; 2478 lambda:=ro[1][1]; 2479 if lambda = 1 then 2480 roro:= HasRoot(Q); 2481 if roro then 2482 poi:=MakePoints(E,roro.ext1); 2483 if poi <> [] then 2484 return [l,quad]; 2485 else 2486 return [l,quad*Coerce(Integers(),ro[2][1])]; 2487 fi; 2488 else return [l,quad*Coerce(Integers(),ro[2][1])]; 2489 fi; 2490 fi; 2491 lambdax:=MultMap(E,lambda).xpol; 2492 li:=List(lambdax); 2493 de:=Denominator(li[1]); 2494 nu:=de*li[1]; 2495 de:=Coerce(Parent(Q),de); 2496 nu:=Coerce(Parent(Q),nu); 2497 if not nu mod Q = de*(Modexp(x,Size(BR),Q)) mod Q then 2498 lambda:=ro[2][1]; 2499 fi; 2500 return [l,quad*Coerce(Integers(),lambda)]; 2501 end; 2502 2503 #computes eigenvalue of frobenius on subgroup in kernel of isogeny of degree l, starting from E to curve whose j-invariant is the k-th root of the second modular polynomial 2504 SubgroupFromJInvariant:= function(E,l,k) 2505 local sum,jTilDer,p1,coeffs,k,PA2,m,C,A,K,CFun,ro,AFun,unf,gen,bo,f,ENeu,xx,ff,kk,sumA,jj, xval,sub,BR,jDer,a,b,ModPoly,PA,j2,j,E2,a,b,a2,d,b2,Phix,Phiy,Phixx,Phiyy,Phixy,J,E4,E6,E4Til,E6Til,c,cTil,i,j,sum1,sum2,divpolyfact,insum,Kxy,Alf,Alf2,defi,Iso,re,R,fro,fr,P,lambda,eifound,E2,i,gen,ro,lambda,lambdax,li,de,nu,x,roro,poi; 2506 Assert(IsPrime(l)," computing prime order subgroup"); 2507 if IsSupersingular(E) then 2508 return Error("cannot compute isogeny for supersingular curve"); 2509 fi; 2510 if l = 2 then 2511 return SubgroupFromJ2(E,k); 2512 fi; 2513 if Characteristic(BaseRing(E))=2 then 2514 return SubFromJChar2(E,l,k); 2515 fi; 2516 Assert(Length(E.L)=2,"curve in short WNF"); 2517 a:= E.L[1]; 2518 b:= E.L[2]; 2519 c:=[]; 2520 cTil:=[]; 2521 BR:= BaseRing(E); 2522 j:=jInvariant(E); 2523 if j = Coerce(BR,0) or j = Coerce(BR,1728) then 2524 return Error("cannot compute isogeny for j-Invariant ",j); 2525 fi; 2526 ModPoly:= ModularPolynomial(l); 2527 PA:=PolynomialAlgebra(BR); 2528 PA2:= PolynomialAlgebra(PA); 2529 ModPoly:= Coerce(PA2,ModPoly); 2530 ro:=Roots(Evaluate(ModPoly,j));#Print(ro, "\n"); 2531 if Length(ro) = 0 then 2532 return Error("There is no ",l,"-isogenous curve"); 2533 fi; 2534 PA:= PolynomialAlgebra(BR); 2535 if k =1 then 2536 j2:=ro[1][1]; 2537 elif Length(ro) =1 then 2538 return Error("The modular polynomial has only one root, try k=1"); 2539 else j2:= ro[2][1]; 2540 fi; 2541 Phix:= Derivative(ModPoly); 2542 if Evaluate(Evaluate(Phix,j),j2) = 0 then 2543 return FAILURE; 2544 fi; 2545 Phiy:=(Derivative(_SwapVars(ModPoly))); 2546 Phiyy:= _SwapVars(Derivative(Phiy)); 2547 Phiy:= _SwapVars(Phiy); 2548 Phixx:= Derivative(Phix); 2549 Phixy:= _SwapVars(Derivative(_SwapVars(Phix))); 2550 E4:=-48*a; 2551 E6:= 864*b; 2552 jDer:= -j*E6/E4; 2553 jTilDer:= -((jDer*Evaluate(Evaluate(Phix,j),j2))/(l*Evaluate(Evaluate(Phiy,j),j2))); 2554 jTilDer:= Coerce(BR,jTilDer); 2555 a2:= -(jTilDer^2)/(48*j2*(j2-1728)); 2556 b2:= -(jTilDer^3)/(864*j2^2*(j2-1728)); 2557 E4Til:= -48*a2; 2558 E6Til:= 864*b2; 2559 J:=(jDer^2*Evaluate(Evaluate(Phixx,j),j2)); 2560 J:= J+2*l*jDer*jTilDer*Evaluate(Evaluate(Phixy,j),j2); 2561 J:= J+l^2*jTilDer^2*Evaluate(Evaluate(Phiyy,j),j2); 2562 J:= -J/(jDer*Evaluate(Evaluate(Phix,j),j2)); 2563 p1:= l/2*J+l/4*(E4^2/E6-l*E4Til^2/E6Til)+l/3*(E6/E4-l*E6Til/E4Til); 2564 ENeu:=EllipticCurve(BR,[l^4*a2,l^6*b2]); 2565 c[1]:=-a/5; 2566 cTil[1]:=-(l^4*a2)/5; 2567 c[2]:= -b/7; 2568 cTil[2]:= -(l^6*b2)/7; 2569 d:= Coerce(Integers(),(l-1)/2); 2570 for i in [3..d+1] do 2571 sum1 :=0; 2572 sum2:=0; 2573 for j in [1..i-2] do 2574 sum1:= sum1 + c[j]*c[i-1-j]; 2575 sum2:= sum2 + cTil[j]*cTil[i-1-j]; 2576 od; 2577 c[i] := sum1*3/((i-2)*(2*i+3)); 2578 cTil[i] := sum2*3/((i-2)*(2*i+3)); 2579 od; 2580 C:= function(k,j) 2581 local f,i,elseq; 2582 f:= Coerce(PA,Append([0],Sublist(c,[1..j]))); 2583 f:=f^(k-j); 2584 elseq:= Eltseq(f); 2585 if elseq=Sequence([]) then 2586 return 0; 2587 elif Length(elseq)<=j then 2588 return 0; 2589 else 2590 return elseq[j+1]; 2591 fi; 2592 end; 2593 A:=[]; 2594 xx:= PA.1; 2595 ff:=-1/2*p1*xx; 2596 for kk in [1..d] do 2597 ff:= ff-(cTil[kk]-l*c[kk])/((2*kk+1)*(2*kk+2))*xx^(kk+1); 2598 od; 2599 sumA:=0; 2600 for jj in [0..d] do 2601 sumA:= sumA+ff^jj/(Factorial(jj)); 2602 od; 2603 A:= Eltseq(sumA); 2604 coeffs:=[]; 2605 coeffs[d+1]:=1; 2606 for i in [1..d] do 2607 m:= d-i+1; 2608 sum:=0; 2609 for k in [1..i] do 2610 insum:=0; 2611 for j in [0..k] do 2612 insum:= insum+Binomial(d-i+k,k-j)*C(k,j); 2613 od; 2614 sum:=sum+insum*coeffs[d-i+k+1]; 2615 od; 2616 coeffs[m]:=A[i+1]-sum; 2617 od; 2618 divpolyfact := Coerce(PA,coeffs); 2619 ro:=Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E))); 2620 if Length(ro) = 2 and ro[1][1] = - ro[2][1] then 2621 K := SplittingField( divpolyfact ); 2622 PA := PolynomialAlgebra(K); 2623 xval:= HasRoot(Coerce(PA,divpolyfact)).ext1; 2624 E2:= ChangeBaseRing(E,Extension(K,2)); 2625 P:=MakePoints(E2,xval)[1]; 2626 fro:=extFrob(E,Degree(BaseRing(E2),BR)); 2627 lambda:=Coerce(Integers(),ro[1][1]); 2628 if lambda*P = fro(P) then 2629 return [l,lambda]; 2630 else 2631 return [l,Coerce(Integers(),ro[2][1])]; 2632 fi; 2633 fi; 2634 lambda:=Coerce(Integers(),ro[1][1]); 2635 if lambda = 1 then 2636 roro:= HasRoot(divpolyfact); 2637 if roro then 2638 poi:=MakePoints(E,roro.ext1); 2639 if poi <> [] then 2640 return [l,1]; 2641 else 2642 return [l,Coerce(Integers(),ro[2][1])]; 2643 fi; 2644 else 2645 return [l,Coerce(Integers(),ro[2][1])]; 2646 fi; 2647 fi; 2648 lambdax:=MultMap(E,lambda).xpol; 2649 li:=List(lambdax); 2650 de:=Denominator(li[1]); 2651 nu:=de*li[1]; 2652 de:=Coerce(Parent(divpolyfact),de); 2653 nu:=Coerce(Parent(divpolyfact),nu); 2654 x:=PolynomialAlgebra(BR).1; 2655 if not nu mod divpolyfact = de*(Modexp(x,Size(BR),divpolyfact)) mod divpolyfact then 2656 lambda:=ro[2][1]; 2657 fi; 2658 return [l,Coerce(Integers(),lambda)]; 2659 end; 2660 2661 2662 2663 2664 2665 2666 #computes isogeny of degree 2 from an ideal 2667 IsoIdeals2:= function(fact,E1,E2,isoli) 2668 local E,n,mo,ma,I,i,lev,k,jtmp,j,j2,ro,path,Iso,c,O; 2669 if Length(isoli)=0 then 2670 E:=E1; 2671 else 2672 E:= Codomain(isoli[Length(isoli)]); 2673 fi; 2674 O:=EndRing(E); 2675 c:= Conductor(O); 2676 n:= fact[2]; 2677 I:=fact[1]; 2678 mo:= Coerce(PolynomialAlgebra(PolynomialAlgebra(BaseRing(E))),ModularPolynomial(2)); 2679 j:=jInvariant(E); 2680 j2:= jInvariant(E2); 2681 Extend(j); 2682 j.prev:=false; 2683 ro:=Roots(Evaluate(mo,j)); 2684 j.next:=[]; 2685 for i in [1..Length(ro)] do 2686 if Conductor(EndRing(EllipticCurveFromJInvariant(ro[i][1])))=c then 2687 Append_(j.next,[ro[i][1]]); 2688 k:=Length(j.next); 2689 Extend(j.next[k]); 2690 j.next[k].level:=1; 2691 j.next[k].prev:=j; 2692 j.next[k].ro:=i; 2693 j.next[k].checked:=0; 2694 j.checked:=0; 2695 fi; 2696 od; 2697 j.level:=0; 2698 while not j.base=j2 do 2699 i:= Length(j.next); 2700 if j.checked=i then 2701 j:=j.prev; 2702 else 2703 lev:=j.level; 2704 j.checked:=j.checked+1; 2705 j:= j.next[j.checked]; 2706 j.level:=lev+1; 2707 if not j.level =n then 2708 ro:=Roots(Evaluate(mo,j)); 2709 j.next:=[]; 2710 for i in [1..Length(ro)] do 2711 if Conductor(EndRing(EllipticCurveFromJInvariant(ro[i][1])))=c then 2712 Append_(j.next,[ro[i][1]]); 2713 k:= Length(j.next); 2714 Extend(j.next[k]); 2715 jtmp:=j.next[k]; 2716 jtmp.level:=lev+2; 2717 jtmp.prev:=j; 2718 jtmp.ro:=i; 2719 jtmp.checked:=0; 2720 fi; 2721 od; 2722 elif j.base = j2 then 2723 path:=[]; 2724 path[1]:=Tuple([j.base,j.ro]); 2725 while not j.prev = false do 2726 j:=j.prev; 2727 if not j.prev = false then 2728 Append_(path,[Tuple([j.base,j.ro])]); 2729 fi; 2730 od; 2731 j.base:=j2; 2732 else 2733 j:=j.prev; 2734 fi; 2735 fi; 2736 od; 2737 if path=[] then 2738 path[1]:=Tuple([j.base,j.ro]); 2739 while not j.prev = false do 2740 j:=j.prev; 2741 if not j.prev=false then 2742 Append_(path,[Tuple([j.base,j.ro])]); 2743 fi; 2744 od; 2745 fi; 2746 path:= Reversed(path); 2747 for i in [1..Length(path)] do 2748 Append_(isoli,[IsoFromJ2(E,path[i][2])]); 2749 E:= Codomain(isoli[Length(isoli)]); 2750 od; 2751 if not E = E2 then 2752 Append_(isoli,[Isomorphism(E,E2,1)]); 2753 fi; 2754 Iso:=isoli[1]; 2755 for i in [2..Length(isoli)] do 2756 Iso:=Composition(isoli[i],Iso); 2757 od; 2758 return Iso; 2759 end; 2760 2761 #returns eigenvalue of frobenius on kernel subgroup of f, if it has been assigned 2762 FrobeniusEigenvalue:=function(f) 2763 if "ei" in RecFields(f) then 2764 return f.ei; 2765 else 2766 return Error("eigenvalue of frobenius has not been assigned to ", f); 2767 fi; 2768 end; 2769 2770 2771 #computes isogeny between E1 and E2 from a factored product of ideals 2772 IsoFromIdeals:= function(fact,E1,E2) 2773 local bo,nr2,O,PA,ch,pi,I,r,no, elk,Iso,j,f,BR,l,i,P,isoli,n,lambda,Iso,fro,to,E,li,iso1,iso2,j,ro,k,m,tmp,c,PA2,O,c,pi1,ff,Omax1,OE,Omax,pi,rlist,to2,ll; 2774 nr2:=false; 2775 E:= E1; 2776 isoli:=[]; 2777 j:= jInvariant(E2); 2778 ff:=CharacteristicFrobeniusPolynomial(E1); 2779 O:=EndRing(E); 2780 c:=Conductor(O); 2781 if c > 1 then 2782 pi1:=EquationOrder(ff).2; 2783 Omax1:=MaximalOrder(Parent(pi1)); 2784 OE:=EquationOrder(MinimalPolynomial(c*Omax1.2)); 2785 Omax:=MaximalOrder(OE); 2786 ff:=Coerce(PolynomialAlgebra(Omax),ff); 2787 pi:=Roots(ff); 2788 rlist:=List(pi[1][1]); 2789 if rlist[2]>0 then 2790 pi:=Coerce(OE,pi[1][1]); 2791 else 2792 pi:=Coerce(OE,pi[2][1]); 2793 fi; 2794 else 2795 OE:=EquationOrder(ff); 2796 Omax:=MaximalOrder(OE); 2797 pi:=Coerce(Omax,OE.2); 2798 OE:=Omax; 2799 fi; 2800 i:=1; 2801 while i < Length(fact) do 2802 if Norm(fact[i][1]) = Norm(fact[i+1][1]) then 2803 if fact[i][2] = fact[i+1][2] then 2804 Remove_(fact,i); 2805 Remove_(fact,i); 2806 elif fact[i][2]>fact[i+1][2] then 2807 fact[i][2]:=fact[i][2]-fact[i+1][2]; 2808 Remove_(fact,i+1); 2809 i:=i+1; 2810 else 2811 fact[i+1][2]:=fact[i+1][2]-fact[i][2]; 2812 Remove_(fact,i); 2813 i:=i+1; 2814 fi; 2815 else i:=i+1; 2816 fi; 2817 od; 2818 i:=0; 2819 PA2:= PolynomialAlgebra(PolynomialAlgebra(BaseRing(E1))); 2820 while i <Length(fact) do 2821 i:=i+1; 2822 I:=(List(Basis(fact[i][1])[1])[1])*OE+(List(Basis(fact[i][1])[2])[1]+List(Basis(fact[i][1])[2])[2]*OE.2)*OE; 2823 l:= Norm(I); 2824 r:=fact[i][2]; 2825 no:=0; 2826 while no<r do 2827 no:=no+1; 2828 if l > 2 then 2829 ll:=List(pi-Basis(I)[2]); 2830 if ll[2] = 0 then 2831 lambda:=Coerce(Integers(),ll[1]); 2832 else 2833 ro:=Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E1))); 2834 lambda:=Coerce(Integers(),ro[1][1]); 2835 if not pi -lambda in I then 2836 lambda:=Coerce(Integers(),ro[2][1]); 2837 fi; 2838 fi; 2839 iso1:= SubgroupFromJInvariant(E,l,1); 2840 if iso1[2] mod l = lambda mod l then 2841 isoli:= Append(isoli,[IsogenyFromJInvariant(E,l,1)]); 2842 else 2843 iso2:= IsogenyFromJInvariant(E,l,2); 2844 isoli:= Append(isoli,[iso2]); 2845 fi; 2846 E:= Codomain(isoli[Length(isoli)]); 2847 if jInvariant(E) = j then 2848 i := Length(fact); 2849 no:=r+1; 2850 fi; 2851 else 2852 to2:=TorsionSubgroup(E,2); 2853 if Length(to2)<>2 then 2854 return IsoIdeals2(Tuple([fact[i][1],r-(no-1)]),E1,E2,isoli); 2855 fi; 2856 lambda:=Coerce(Integers(),List(pi-Basis(I)[2])[1]); 2857 if lambda mod 2 = 0 then 2858 iso1:=FrobeniusToOtherCurve(E,2); 2859 else 2860 iso1:=IsogenyFromSubgroup(to2[1],E,2); 2861 fi; 2862 isoli:= Append(isoli,[iso1]); 2863 E:= Codomain(isoli[Length(isoli)]); 2864 if jInvariant(E) = j then 2865 i := Length(fact); 2866 fi; 2867 fi; 2868 od; 2869 od; 2870 if not E = E2 then 2871 if not jInvariant(E) = j then 2872 nr2:=true; 2873 else 2874 isoli:= Append(isoli,[Isomorphism(E,E2,1)]); 2875 fi; 2876 fi; 2877 Iso:= isoli[1]; 2878 for i in [2..Length(isoli)] do 2879 Iso:= Composition(isoli[i],Iso); 2880 od; 2881 return Iso; 2882 end; 2883 2884 2885 2886 #The composition f1(f2()) 2887 IsoComp:= function(f2,f1) 2888 local i,xelseq2,xdenom2,xnom2,xpol,re,yelseq2,ynom2,ydenom2; 2889 re:=[]; 2890 xelseq2:= Eltseq(f2.xpol); 2891 xdenom2:= Denominator(xelseq2[1]); 2892 xnom2:= xdenom2*xelseq2[1]; 2893 re[1]:= Evaluate(xnom2,f1.xpol)/Evaluate(xdenom2,f1.xpol); 2894 yelseq2:= Eltseq(f2.ypol); 2895 re[2]:= 0; 2896 for i in [1..Length(yelseq2)] do 2897 ydenom2:= Denominator(yelseq2[i]); 2898 ynom2:= ydenom2*yelseq2[i]; 2899 re[2]:= (re[2]+(Evaluate(ynom2,f1.xpol)/Evaluate(ydenom2,f1.xpol)))*(f1.ypol^(i-1)); 2900 od; 2901 return re; 2902 end; 2903 2904 2905 #subfunction for composition 2906 _Decrease:=function(Iso,K) 2907 local E1,E2,E11,E22,Iso2,yp,el,K2,el1,el2; 2908 yp:=Iso.ypol; 2909 el:=Eltseq(yp); 2910 K2:=BaseRing(Domain(Iso)); 2911 el1:=Coerce(PolynomialAlgebra(K2),el[1]*Denominator(el[1])); 2912 el2:=Coerce(PolynomialAlgebra(K2),el[2]*Denominator(el[2])); 2913 if not Base(IsCoercible(PolynomialAlgebra(K),el2)) or not Base(IsCoercible(PolynomialAlgebra(K),el1)) then 2914 return false; 2915 fi; 2916 E1:=Domain(Iso); 2917 E2:=Codomain(Iso); 2918 E11:=ChangeBaseRing(E1,K); 2919 E22:=ChangeBaseRing(E2,K); 2920 Iso2:=Copy(Iso); 2921 Iso2.domain:=ChangeBaseRing(Domain(Iso),K); 2922 Iso2.codomain:=ChangeBaseRing(Codomain(Iso),K); 2923 return Iso2; 2924 end; 2925 2926 #subfunction for composition 2927 _Increase:=function(Iso,K) 2928 local E1,E2,E11,E22,Iso2; 2929 E1:=Domain(Iso); 2930 E2:=Codomain(Iso); 2931 E11:=ChangeBaseRing(E1,K); 2932 E22:=ChangeBaseRing(E2,K); 2933 Iso2:=Copy(Iso); 2934 Iso2.domain:=ChangeBaseRing(Domain(Iso),K); 2935 Iso2.codomain:=ChangeBaseRing(Codomain(Iso),K); 2936 return Iso2; 2937 end; 2938 2939 #composition Iso1(Iso2) 2940 _Composition:= function(Iso2,Iso1) 2941 local Iso,_Iso,P,EE,comp,Iso22,Iso11; 2942 if not Codomain(Iso1) = Domain(Iso2) then 2943 if IsSubset(BaseRing(Codomain(Iso1)),BaseRing(Domain(Iso2))) then 2944 Iso11:=Iso1; 2945 Iso22:=_Decrease(Iso2,BaseRing(Codomain(Iso1))); 2946 if Iso22 = false then 2947 Iso11:=_Increase(Iso1,BaseRing(Codomain(Iso2))); 2948 Iso22:=Iso2; 2949 fi; 2950 elif IsSubset(BaseRing(Domain(Iso2)),BaseRing(Codomain(Iso1))) then 2951 Iso11:=_Decrease(Iso1,BaseRing(Domain(Iso2))); 2952 Iso22:=Iso2; 2953 if Iso11 = false then 2954 Iso22:=_Increase(Iso2,BaseRing(Codomain(Iso1))); 2955 Iso11:=Iso1; 2956 fi; 2957 else 2958 return Error("composition: domain of ", Iso2," is not codomain of ",Iso1); 2959 fi; 2960 if not Codomain(Iso11) = Domain(Iso22) then 2961 return Error("composition: domain of ", Iso2," is not codomain of ",Iso1); 2962 fi; 2963 else 2964 Iso11:=Iso1; 2965 Iso22:=Iso2; 2966 fi; 2967 if IsZero(Iso11) or IsZero(Iso22) then 2968 return NuMap(Iso11.domain,Iso22.codomain); 2969 fi; 2970 P:= GenericPoint(Iso11.domain); 2971 EE:=Parent(P); 2972 _Iso:=function(P1) 2973 return Iso22(Iso11(P1)); 2974 end; 2975 Iso:= Map(Iso11.domain,Iso22.codomain,_Iso); 2976 Iso.type:=iso; 2977 Iso.degree:=Iso11.degree*Iso22.degree; 2978 comp:= IsoComp(Iso22,Iso11); 2979 Iso.xpol:= comp[1]; 2980 Iso.ypol:= comp[2]; 2981 Iso.operations:= _iso_ops; 2982 if Degree(Iso1) = 1 then 2983 if "ei" in RecFields(Iso2) then 2984 Iso.ei:=Iso2.ei; 2985 fi; 2986 elif Degree(Iso2) = 1 then 2987 if "ei" in RecFields(Iso1) then 2988 Iso.ei:=Iso1.ei; 2989 fi; 2990 fi; 2991 return Iso; 2992 end; 2993 2994 2995 #returns true if I is L-smooth, i.e. factors have norms given in L 2996 IsSmooth:= function(I,O,L,E) 2997 local i,m,smooth,fact,Idoben,ba,Omax; 2998 fact:=IdealFactorization(I); 2999 if fact = FAILURE then 3000 return false; 3001 fi; 3002 smooth:= true; 3003 for i in [1..Length(fact)] do 3004 if not Norm (fact[i][1]) in L then 3005 smooth:= false; 3006 i:= Length(fact)+1; 3007 fi; 3008 od; 3009 return smooth; 3010 end; 3011 3012 3013 MakeSmooth:=function(I,O,L,E) 3014 local Omax,Cl,F,G,b,L1,bool,c,L2,M,fact,r,j,Ma,H,L3,e,II,ten,primeIds,i,l,p,ba,L3,k,pp,fact,gens,hnf,a,l,MM,m; 3015 if IsSmooth(I,O,L,E) then 3016 return I; 3017 fi; 3018 if Conductor(O)<>1 then return Error("smoothing ideals only implemented for ideals of maximal order"); 3019 fi; 3020 Omax:=MaximalOrder(O); 3021 Cl:=ClassGroup(Omax); 3022 G:=Cl.ext1; 3023 F:=Inverse(G); 3024 ba:=Basis(I); 3025 a:=Coerce(Omax,ba[1])*Omax+(Coerce(Omax,ba[2]))*Omax; 3026 L1:=List(F(a)); 3027 bool:=true; 3028 r:=Length(L1); 3029 gens:=[]; 3030 for i in [1..r] do 3031 fact:=Factorization(G(Cl.(i))); 3032 gens[i]:=1*Omax; 3033 for j in [1..Length(fact)] do 3034 gens[i]:=gens[i]*fact[j][1]; 3035 od; 3036 od; 3037 M:=[]; 3038 primeIds:=[]; 3039 c:= Conductor(EquationOrder(CharacteristicFrobeniusPolynomial(E))); 3040 for i in [1..Length(L)] do 3041 l:=L[i]; 3042 fact:= IdealFactorization(l*Omax); 3043 if GCD(Norm(fact[1][1]),c)=1 then 3044 Append_(primeIds,[fact[1][1]]); 3045 if Length(fact) = 2 then 3046 Append_(primeIds,[fact[2][1]]); 3047 fi; 3048 fi; 3049 od; 3050 MM:=[]; 3051 for i in [1..r] do 3052 L2:=List(F(primeIds[i])); 3053 for j in [1..r] do 3054 if j=i then 3055 Append_(L2,[1]); 3056 else 3057 Append_(L2,[0]); 3058 fi; 3059 od; 3060 MM[i]:=L2; 3061 od; 3062 M:=[]; 3063 for i in [1..r] do 3064 Append_(M,MM[i]); 3065 od; 3066 Ma:=Matrix(Z,r,2*r,M); 3067 H:=Base(HermiteForm(Ma)); 3068 bool:=true; 3069 k:=r; 3070 hnf:=false; 3071 while k<Length(primeIds) and not hnf do 3072 for i in [1..r] do 3073 if H[i][i]<>1 then 3074 bool:=false; 3075 fi; 3076 od; 3077 if not bool then 3078 k:=k+1; 3079 L2:=List(F(primeIds[k])); 3080 for j in [1..k-1] do 3081 Append_(L2,[0]); 3082 Append_(MM[j],[0]); 3083 od; 3084 Append_(L2,[1]); 3085 MM[k]:=L2; 3086 M:=[]; 3087 for i in [1..k] do 3088 Append_(M,MM[i]); 3089 od; 3090 Ma:=Matrix(Z,k,k+r,M); 3091 H:=Base(HermiteForm(Ma)); 3092 bool:=true; 3093 else hnf:=true; 3094 fi; 3095 od; 3096 L3:=[]; 3097 for i in [1..r] do 3098 m:=r-i+1; 3099 e:=H[m][m]; 3100 pp:=1*Omax; 3101 for j in [r+1..r+k] do 3102 pp:=pp*primeIds[j-r]^H[m][j]; 3103 od; 3104 for l in [1..i-1] do 3105 pp:=pp*L3[l]^(-H[m][i-l]); 3106 od; 3107 pp:=pp^(-e); 3108 L3[i]:=Denominator(pp)*pp; 3109 od; 3110 p:=1*Omax; 3111 Reverse_(L3); 3112 for i in [1..r] do 3113 e:=L1[i]; 3114 if e <>0 then 3115 p:=p*L3[i]^e; 3116 fi; 3117 od; 3118 return p; 3119 end; 3120 3121 #old MakeSmooth, tries to make I L-smooth 3122 MakeSmooth2:= function(I,O,L,E) 3123 local Omax, fact,a,b,l,ro,f,phi,primeIds,h,ra,found,factored,pos,i,c,bli,posl,j,tr,fact,c,L2,Idsoben,ba,Iinv,m,p,n,r,k,p1,t,e,p2,e1,e2,Omax,Cl,G,F,primes; 3124 if IsSmooth(I,O,L,E) then 3125 return I; 3126 fi; 3127 c:= Conductor(EquationOrder(CharacteristicFrobeniusPolynomial(E))); 3128 tr:= 0; 3129 primeIds:=[]; 3130 Idsoben:=[]; 3131 Iinv:=IdInvert(I); 3132 Omax:= MaximalOrder(O); 3133 h:= PicardNumber(O); 3134 phi:= Basis(Omax)[2]; 3135 found:=false; 3136 primes:=PrimesUpTo(1000); 3137 for i in [1..Length(primes)] do 3138 if IsElkiesPrime(E,primes[i]) then 3139 fact:= IdealFactorization(primes[i]*O); 3140 if fact<> FAILURE then 3141 Append_(primeIds,[fact[1][1]]); 3142 if Length(fact) = 2 then 3143 Append_(primeIds,[fact[2][1]]); 3144 fi; 3145 fi; 3146 fi; 3147 od; 3148 m:=Length(primeIds); 3149 t:=Ceiling(Root(Root(Size(BaseRing(E)),4),m)); 3150 a:=[]; 3151 i:=1; 3152 while i <=Length(primeIds) do 3153 n:=Norm(primeIds[i]); 3154 Append_(a,[n]); 3155 Append_(a,[[primeIds[i],0]]); 3156 i:=i+1; 3157 if Norm(primeIds[i]) = n then 3158 Append_(a[Length(a)],[primeIds[i],0]); 3159 i:=i+1; 3160 fi; 3161 od; 3162 b:= I; 3163 i:=0; 3164 while not found do 3165 p:=Random(primeIds); 3166 n:=Norm(p); 3167 r:=Random(t); 3168 b:= IdReduction(b*(p^r)); 3169 pos:=Position(a,n); 3170 if a[pos+1][1] = p then 3171 a[pos+1][2]:=a[pos+1][2]+ r; 3172 else a[pos+1][4]:=a[pos+1][4]+ r; 3173 fi; 3174 fact:=IdealFactorization(b); 3175 i:=i+1; 3176 if fact <> FAILURE then 3177 if IsSmooth(b,O,L,E) then 3178 found:=true; 3179 fi; 3180 fi; 3181 od; 3182 p:=1*O; 3183 k:=1; 3184 i:=1; 3185 while i <=Length(a) do 3186 p1:=a[i+1][1]; 3187 e1:=a[i+1][2]; 3188 if IsPrincipal(p1^2) then 3189 p2:=p1; 3190 e2:=e1; 3191 else 3192 p2:=a[i+1][3]; 3193 e2:=a[i+1][4]; 3194 fi; 3195 if k<=Length(fact) then 3196 if p1 = fact[k][1] then 3197 e:=fact[k][2]-e1; 3198 k:=k+1; 3199 if e>=0 then 3200 p:=p*(p1^e); 3201 else p:=p*p2^e; 3202 fi; 3203 else p:=p*p2^e1; 3204 fi; 3205 else p:=p*p2^e1; 3206 fi; 3207 if k<=Length(fact) then 3208 if p2 = fact[k][1] then 3209 e:=fact[k][2]-e2; 3210 k:=k+1; 3211 if e>=0 then 3212 p:=p*(p2^e); 3213 else p:=p*p1^e; 3214 fi; 3215 else p:=p*p1^e2; 3216 fi; 3217 else p:=p*p1^e2; 3218 fi; 3219 i:=i+2; 3220 od; 3221 return p; 3222 end; 3223 3224 3225 #returns the domain of a map 3226 _Domain:= function(f) 3227 return f.domain; 3228 end; 3229 3230 #returns the codomain of a map 3231 _Codomain:= function(f) 3232 return f.codomain; 3233 end; 3234 3235 3236 3237 3238 #true if there is more than one root of Psi_l(j,x) 3239 ExistsSecondRoot:= function(l,j) 3240 local ro,f,F,PA; 3241 F:= Parent(j); 3242 PA:= PolynomialAlgebra(PolynomialAlgebra(F)); 3243 f:= Coerce(PA,ModularPolynomial(l)); 3244 f:= Evaluate(f,j); 3245 ro:= Roots(f); 3246 if Length(ro) =1 then 3247 return FALSE; 3248 elif Length(ro) >= 2 then 3249 return TRUE; 3250 else 3251 return Error("ExistsSecondRoot: l must be an Elkies Prime"); 3252 fi; 3253 end; 3254 3255 #returns the n-th root of Psi_l(j,x) 3256 ModPolyRoot:= function(l,j,n) 3257 local ro,f,F,PA; 3258 F:= Parent(j); 3259 PA:= PolynomialAlgebra(PolynomialAlgebra(F)); 3260 f:= Coerce(PA,ModularPolynomial(l)); 3261 f:= Evaluate(f,j); 3262 ro:= Roots(f); 3263 if n = 1 then 3264 return [ro[n][1],1]; 3265 elif n > Length(ro) then 3266 return [ro[1][1],1]; 3267 else 3268 return [ro[2][1],n]; 3269 fi; 3270 end; 3271 3272 3273 #tests if an isogeny of E is computable 3274 TestElkiesPrimes:= function(E) 3275 local el,i,j,c; 3276 el:= ElkiesPrimes59(E); 3277 j:= jInvariant(E); 3278 for i in [1..Length(el)] do 3279 if Type(SubgroupFromJInvariant(E,el[i],1)[1])=elt-ord^rat then 3280 return [el[i],1]; 3281 fi; 3282 if ExistsSecondRoot(el[i],j) then 3283 if Type(SubgroupFromJInvariant(E,el[i],2)[1]) = elt-ord^rat then 3284 return [el[i],2]; 3285 fi; 3286 fi; 3287 od; 3288 return FAILURE; 3289 end; 3290 3291 #computes the picard-number of O 3292 PicardNumber:=function(O) 3293 local h,p,i,fact,c,Omax,UE,UO,l,pic,Id,factid; 3294 if IsMaximal(O) then 3295 return ClassNumber(O); 3296 fi; 3297 Omax:=MaximalOrder(O); 3298 h:=ClassNumber(Omax); 3299 c:=Conductor(O); 3300 UE:=UnitGroup(O); 3301 UO:=UnitGroup(Omax); 3302 if Size(Base(UO)) = INFTY and Size(Base(UE)) = INFTY then 3303 l:=1; 3304 else 3305 l:=Size(Base(UO))/Size(Base(UE)); 3306 fi; 3307 pic:=1; 3308 fact:=Factorization(c); 3309 for i in [1..Length(fact)] do 3310 p:=fact[i][1]; 3311 Id:=p*Omax; 3312 factid:=Factorization(Id); 3313 if Length(factid) = 2 then 3314 pic:=fact[i][2]*pic*(1-1/p); 3315 elif fact[i][2] = 1 then 3316 pic:=fact[i][2]*pic*(1+1/p); 3317 fi; 3318 od; 3319 return Coerce(Integers(),pic*h*(c/l)); 3320 end; 3321 3322 3323 #tests if an isogeny of degree l between E and E2 (with j(E) = j and j(E2) = j2) is computable 3324 SubComputable:= function(j,j2,l) 3325 local ModPoly,BR,PA,PA2,Phix; 3326 ModPoly:= ModularPolynomial(l); 3327 BR:= Parent(j); 3328 PA:=PolynomialAlgebra(BR); 3329 PA2:= PolynomialAlgebra(PA); 3330 ModPoly:= Coerce(PA2,ModPoly); 3331 Phix:= Derivative(ModPoly); 3332 if Evaluate(Evaluate(Phix,j),j2) = 0 then 3333 return false; 3334 else return true; 3335 fi; 3336 end; 3337 3338 #tries to find a connecting pseudo-random function between E1 and E2 3339 PseudoRandom2:= function(E1,E2,r,L) 3340 local j, i, jr, k, n, j2,f,found,l,f,m,o,tes,ma,O,c,L2; 3341 O:=EndRing(E1); 3342 c:=Conductor(O); 3343 k:=-1; 3344 m:=Length(L); 3345 found:= false; 3346 while not found do 3347 k:= k+1; 3348 if k =m then 3349 L2:=L; 3350 for i in [Last(L)..59] do 3351 if IsPrime(i) then 3352 if IsElkiesPrime(E1,i) then 3353 L2:=Append(L,[i]); 3354 fi; 3355 fi; 3356 od; 3357 if L2 = L then 3358 return Error(" no connectiong pseudo random function found"); 3359 else 3360 return PseudoRandom(E1,E2,r,L2); 3361 fi; 3362 fi; 3363 j:= jInvariant(E1); 3364 for i in [1..r] do 3365 l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1]; 3366 n:= (Coerce(Integers(),(j+k+1)) mod 2) +1; 3367 j2:= ModPolyRoot(l,j,n)[1]; 3368 f:=0; 3369 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3370 while ((not ma) or j2 = 0 or j2 = 1728 or (not SubComputable(j,j2,l))) do 3371 f:=f+1; 3372 l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1]; 3373 n:= (Coerce(Integers(),(j+k+1)) mod 2) +1; 3374 j2:= ModPolyRoot(l,j,n)[1]; 3375 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3376 if f = 5 then 3377 tes:= TestElkiesPrimes(E); 3378 if tes = FAILURE then 3379 return Error("Isogeny cannot be computed"); 3380 else j2 := ModPolyRoot(tes[1],j,tes[2])[1]; 3381 fi; 3382 fi; 3383 od; 3384 j:= j2; 3385 od; 3386 jr:=j; 3387 j:= jInvariant(E2); 3388 o:=0; 3389 while o < r do 3390 o:= o+1; 3391 l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1]; 3392 n:= (Coerce(Integers(),(j+k+1)) mod 2) +1; 3393 j2:= ModPolyRoot(l,j,n)[1]; 3394 f:=0; 3395 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3396 while ((not ma) or j2 = 0 or j2 = 1728 or (not SubComputable(j,j2,l))) do 3397 f:=f+1; 3398 l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1]; 3399 n:= (Coerce(Integers(),(j+k+1)) mod 2) +1; 3400 j2:= ModPolyRoot(l,j,n)[1]; 3401 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3402 if f = 5 then tes:= TestElkiesPrimes(E); 3403 if tes = FAILURE then 3404 return Error("Isogeny cannot be computed"); 3405 else j2 := ModPolyRoot(tes[1],j,tes[2])[1]; 3406 fi; 3407 fi; 3408 od; 3409 j:= j2; 3410 if j = jr then 3411 return [k,2,true,jr]; 3412 fi; 3413 od; 3414 od; 3415 end; 3416 3417 #tries to find a connecting pseudo-random function between E1 and E2 3418 PseudoRandom59:= function(E1,E2,r,L) 3419 local j, i, jr, k, n, j2,f,found,l,f,m,o,tes,ma,O,c; 3420 k:=-1; 3421 O:=EndRing(E1); 3422 c:= Conductor(O); 3423 L:= ElkiesPrimes59(E1); 3424 m:=Length(L); 3425 found:= false; 3426 while not found do 3427 k:= k+1; 3428 if k =m then 3429 return PseudoRandom2(E1,E2,r,L); 3430 fi; 3431 j:= jInvariant(E1); 3432 for i in [1..r] do 3433 l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1]; 3434 n:= (Coerce(Integers(),(j+k)) mod 2) +1; 3435 j2:= ModPolyRoot(l,j,n)[1]; 3436 f:=0; 3437 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3438 while ((not ma) or j2 = 0 or j2 = 1728 or ( not SubComputable(j,j2,l))) do 3439 f:=f+1; 3440 l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1]; 3441 n:= (Coerce(Integers(),(j+k)) mod 2) +1; 3442 j2:= ModPolyRoot(l,j,n)[1]; 3443 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3444 if f = 5 then tes:= TestElkiesPrimes(E); 3445 if tes = FAILURE then 3446 return Error("Isogeny cannot be computed"); 3447 else j2 := ModPolyRoot(tes[1],j,tes[2])[1]; 3448 fi; 3449 fi; 3450 od; 3451 j:= j2; 3452 od; 3453 jr:=j; 3454 j:= jInvariant(E2); 3455 o:=0; 3456 while o < r do 3457 o:= o+1; 3458 l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1]; 3459 n:= (Coerce(Integers(),(j+k)) mod 2) +1; 3460 j2:= ModPolyRoot(l,j,n)[1]; 3461 f:=0; 3462 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3463 while ((not ma) or j2 = 0 or j2 = 1728 or (not SubComputable(j,j2,l))) do 3464 f:=f+1; 3465 l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1]; 3466 n:= (Coerce(Integers(),(j+k)) mod 2) +1; 3467 j2:= ModPolyRoot(l,j,n)[1]; 3468 ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c; 3469 if f = 5 then 3470 tes:= TestElkiesPrimes(E); 3471 if tes = FAILURE then 3472 return Error("Isogeny cannot be computed"); 3473 else j2 := ModPolyRoot(tes[1],j,tes[2])[1]; 3474 fi; 3475 fi; 3476 od; 3477 j:= j2; 3478 if j = jr then 3479 return [k,1,true,jr]; 3480 fi; 3481 od; 3482 od; 3483 end; 3484 3485 #tries to find a connecting pseudo-random function between E1 and E2 3486 PseudoRandom:= function(E1,E2,r,L,polys) 3487 local j, i, jr, k, n,loop,c,xalt, j2,f,found,l,f,m,o,tes,ma,O,c,eval,BR,q,p,L2,poss,mi,ind,PA2,r2,degs,els; 3488 m:=Length(L); 3489 BR:= BaseRing(E1); 3490 found:= false; 3491 k:=-1; 3492 poss:=[]; 3493 loop:=Last(L) > 37; 3494 eval:=function(l,j,n) 3495 local ro; 3496 ro:= Roots(Evaluate(polys[Position(L,l)],j)); 3497 if n = 1 then 3498 return [ro[n][1],1,Length(ro)]; 3499 elif (n =2 and Length(ro) =1) then 3500 return [ro[1][1],1,Length(ro)]; 3501 else return [ro[2][1],2,Length(ro)]; 3502 fi; 3503 end; 3504 q:=1; 3505 Print("von E1\n"); 3506 while not found do 3507 k:= k+1; 3508 if k =m then 3509 if poss=[] then 3510 L2:=L; 3511 PA2:=Parent(polys[1]); 3512 for i in [Last(L)+1..59] do 3513 if IsPrime(i) then 3514 if IsElkiesPrime(E1,i) then 3515 Append_(L2,[i]); 3516 Append_(polys,[Coerce(PA2,ModularPolynomial(i))]); 3517 fi; 3518 fi; 3519 od; 3520 if Length(L2) = m then 3521 return Error(" no connectiong pseudo random function found"); 3522 else return PseudoRandom(E1,E2,r,L2,polys); 3523 fi; 3524 elif Length(poss)= 1 then 3525 return poss[1][1]; 3526 else mi:=poss[1][3]; 3527 ind:=1; 3528 for i in [2..Length(poss)] do 3529 if poss[i][3] < mi then 3530 mi:=poss[i][3]; 3531 ind:=i; 3532 fi; 3533 od; 3534 return poss[ind][1]; 3535 fi; 3536 fi; 3537 j:= jInvariant(E1); 3538 i :=0; 3539 while i <r do 3540 i:=i+1; 3541 els:=Eltseq(j); 3542 q:=Sum(function(x) return Coerce(Integers(),x); end,els); 3543 l:= L[(q+k+1) mod m +1]; 3544 n:= ((q+k) mod 2) +1; 3545 j2:= eval(l,j,n); 3546 f:=0; 3547 while (j2[1] = 0 or j2[1] = 1728 or (not SubComputable(j,j2[1],l)) or j2[1] = j) do 3548 f:=f+1; 3549 l:= L[(q+k+f+1) mod m +1]; 3550 n:= ((q+k) mod 2) +1; 3551 j2:= eval(l,j,n); 3552 if f = 5 then tes:= TestElkiesPrimes(E1); 3553 if tes = FAILURE then 3554 return Error("Isogeny cannot be computed"); 3555 else j2 := eval(tes[1],j,tes[2]); 3556 fi; 3557 fi; 3558 od; 3559 j:= j2[1]; 3560 Print(j,"\n"); 3561 od; 3562 jr:=j; 3563 j:= jInvariant(E2); 3564 Print("von E2\n"); 3565 o:=0; 3566 degs:=0; 3567 while o < r do 3568 o:= o+1; 3569 els:=Eltseq(j); 3570 q:=Sum(function(x) return Coerce(Integers(),x); end,els); 3571 l:= L[(q+k+1) mod m +1]; 3572 n:= ((q+k) mod 2) +1; 3573 j2:= eval(l,j,n); 3574 f:=0; 3575 while (j2[1] = 0 or j2[1] = 1728 or not (SubComputable(j,j2[1],l))or j2[1]=j) do 3576 f:=f+1; 3577 l:= L[(q+f+k+1) mod m +1]; 3578 n:= ((q+k) mod 2) +1; 3579 j2:= eval(l,j,n); 3580 if f = 5 then tes:= TestElkiesPrimes(E1); 3581 if tes = FAILURE then 3582 return Error("Isogeny cannot be computed"); 3583 else j2 := eval(tes[1],j,tes[2]); 3584 fi; 3585 fi; 3586 od; 3587 degs:=degs+l; 3588 j:= j2[1]; 3589 Print(j,"\n"); 3590 if j = jr then 3591 return [k,1,loop,jr]; 3592 fi; 3593 od; 3594 od; 3595 end; 3596 3597 3598 3599 #returns ideal that represents an isogeny between E1 and E2 3600 GetIsoIdeal:= function(E1,E2,c) 3601 local pos,p,r,k,I1,I2,n,E,I,f,h,BR,EC,fro,P,mu,ro,frob,R,lambda,sub,nu,su,tes,nr,eq,D,a,fact,pos,ma,PA2,eval,mopos,index,co,Etmp,m,cha2,Iso,l,el,Isos,Ids,dom,Id,Cl,O,Ge,Gen,ma,Id,pi1,pi,su,j,i,Id,ord,o,elli,pos,j1,j2,jInvs,found,jr,ch,id,PA,phi,c,Omax1,Omax,polys,pri,q,jr,L,mytime,loop,OE,w,r1,r2,ba,ff,rlist,xalt,els,mytime; 3602 BR:= BaseRing(E1); 3603 pri:=PrimitiveElement(BR); 3604 PA2:=PolynomialAlgebra(PolynomialAlgebra(BR)); 3605 ff:=CharacteristicFrobeniusPolynomial(E1); 3606 if c > 1 then 3607 pi1:=EquationOrder(ff).2; 3608 Omax1:=MaximalOrder(Parent(pi1)); 3609 OE:=EquationOrder(MinimalPolynomial(c*Omax1.2)); 3610 Omax:=MaximalOrder(OE); 3611 ff:=Coerce(PolynomialAlgebra(Omax),ff); 3612 pi:=Roots(ff); 3613 rlist:=List(pi[1][1]); 3614 if rlist[2]>0 then 3615 pi:=Coerce(OE,pi[1][1]); 3616 else 3617 pi:=Coerce(OE,pi[2][1]); 3618 fi; 3619 else 3620 OE:=EquationOrder(ff); 3621 Omax:=MaximalOrder(OE); 3622 pi:=Coerce(Omax,OE.2); 3623 OE:=Omax; 3624 fi; 3625 phi:=Omax.2; 3626 PA:=PolynomialAlgebra(Omax); 3627 ch:= CharacteristicFrobeniusPolynomial(E1); 3628 el:= ElkiesPrimes(E1); 3629 j:= jInvariant(E1); 3630 m:=Length(el); 3631 h:= PicardNumber(OE); 3632 r:= Ceiling(Sqrt(h)); 3633 k:= -1; 3634 found:=false; 3635 I1:=1*OE; 3636 I2:=1*OE; 3637 E:=E1; 3638 m:=Length(el); 3639 found:= false; 3640 polys:=[]; 3641 for i in [1..Length(el)] do 3642 polys[i]:= Coerce(PA2,ModularPolynomial(el[i])); 3643 od; 3644 k:=PseudoRandom(E1,E2,r,el,polys); 3645 nr:= k[2]-1; 3646 jr:=k[4]; 3647 if k[3] then el:=ElkiesPrimes59(E1); 3648 polys:=[]; 3649 for i in [1..Length(el)] do 3650 polys[i]:= Coerce(PA2,ModularPolynomial(el[i])); 3651 od; 3652 fi; 3653 k:=k[1]; 3654 m:=Length(el); 3655 eval:=function(l,j,n) 3656 local ro; 3657 ro:= Roots(Evaluate(polys[Position(el,l)],j)); 3658 if n = 1 then 3659 return [ro[n][1],1,Length(ro)]; 3660 elif n > Length(ro) then 3661 return [ro[1][1],1,Length(ro)]; 3662 else return [ro[n][1],n,Length(ro)]; 3663 fi; 3664 end; 3665 Print("computing subgroups\n"); 3666 while not found do 3667 E:=E1; 3668 j:=jInvariant(E1); 3669 i:=0; 3670 while i <r do 3671 i:=i+1; 3672 els:=Eltseq(j); 3673 q:=Sum(function(x) return Coerce(Integers(),x); end,els); 3674 l:= el[(q+k+1) mod m +1]; 3675 n:= (Coerce(Integers(),(q+k)) mod 2) +1; 3676 j2:= eval(l,j,n); 3677 su:= (SubgroupFromJInvariant(E,l,j2[2])); 3678 f:=0; 3679 while (j2[1] = 0 or j2[1] = 1728 or not Type(su) = list or j2[1]=j) do 3680 f:=f+1; 3681 l:= el[(q+f+k+1) mod m +1]; 3682 n:= (Coerce(Integers(),(q+k)) mod 2) +1; 3683 j2:= eval(l,j,n); 3684 su:= SubgroupFromJInvariant(E,l,j2[2]); 3685 if f = 5 then tes:= TestElkiesPrimes(E); 3686 if tes = FAILURE then 3687 return Error("Isogeny cannot be computed"); 3688 else su:=SubgroupFromJInvariant(E,tes[1],tes[2]); 3689 fi; 3690 fi; 3691 od; 3692 lambda:=su[2]; 3693 Id:=l*OE+(pi-lambda)*OE; 3694 I1:= IdReduction(Id*I1); 3695 j:=j2[1]; 3696 Print("j: ",j, "l: ",l,"\n"); 3697 E:=ECFromJGivenNumber(j,Size(E)); 3698 if j = jr then 3699 i:=r; 3700 fi; 3701 od; 3702 Print("von E2\n"); 3703 jr:=j; 3704 j:= jInvariant(E2); 3705 E:=E2; 3706 nu:=0; 3707 while nu <r do 3708 nu:=nu+1; 3709 els:=Eltseq(j); 3710 q:=Sum(function(x) return Coerce(Integers(),x); end,els); 3711 l:= el[(q+k+1) mod m +1]; 3712 n:= ((q+k) mod 2) +1; 3713 j2:= eval(l,j,n); 3714 su:= (SubgroupFromJInvariant(E,l,j2[2])); 3715 f:=0; 3716 while j2[1] =0 or j2[1] = 1728 or Type(su)<>list or j2[1]=j do 3717 f:=f+1; 3718 l:= el[(q+f+k+1) mod m +1]; 3719 n:= ((q+k) mod 2) +1; 3720 j2:= eval(l,j,n); 3721 su:= SubgroupFromJInvariant(E,l,j2[2]); 3722 if f = 5 then 3723 tes:= TestElkiesPrimes(E); 3724 if tes = FAILURE then 3725 return Error("Isogeny cannot be computed"); 3726 else su:=SubgroupFromJInvariant(E,tes[1],tes[2]); 3727 fi; 3728 fi; 3729 od; 3730 Id:=l*OE+(pi-su[2])*OE; 3731 I1:=IdReduction(IdInvert(Id)*I1); 3732 j:=j2[1]; 3733 Print("j: ",j, "l: ",l,"\n"); 3734 E:=ECFromJGivenNumber(j,Size(E)); 3735 if j = jr then 3736 found:=true; 3737 nu := r+1; 3738 fi; 3739 od; 3740 od; 3741 Print(I1); 3742 return MakeSmooth(I1,OE,el,E1); 3743 end; 3744 3745 3746 3747 3748 #converts a hex-number to an integer 3749 HexToInt:=function(L) 3750 local i,sum,k,tmp; 3751 sum:=0; 3752 k:=Length(L); 3753 for i in [1..k] do 3754 tmp:=L[k-i+1]; 3755 if not (ASCII(tmp)-48) in [0..9] then 3756 tmp:=ASCII(tmp)-87; 3757 else tmp:=ASCII(tmp)-48; 3758 fi; 3759 sum:=sum+tmp*16^(i-1); 3760 od; 3761 return sum; 3762 end; 3763 3764 3765 3766 3767 3768 3769 3770 #factorization of an ideal prime to the conductor of its order 3771 IdealFactorization:=function(Id) 3772 local O,Omax,fact,Ioben,ba,i,fact2,II,ten,d,Id2; 3773 O:=Order(Id); 3774 d:=Denominator(Id); 3775 Assert(d = 1," integral ideal "); 3776 if IsMaximal(O) then 3777 return Factorization(Id); 3778 else 3779 Omax:=MaximalOrder(O); 3780 ba:=Basis(Id); 3781 Ioben:=Coerce(Omax,ba[1])*Omax+(Coerce(Omax,ba[2]))*Omax; 3782 fact:=Factorization(Ioben); 3783 fact2:=[]; 3784 for i in [1..Length(fact)] do 3785 II:=Intersection(O,fact[i][1]); 3786 if GCD(Norm(II),Conductor(O))<>1 then 3787 return FAILURE; 3788 fi; 3789 ten:=TwoElementNormal(II); 3790 II:=Ideal(O,[Base(ten),ten.ext1]); 3791 II:=Ideal(O,Basis(II)); 3792 Append_(fact2,[Tuple([II,fact[i][2]])]); 3793 od; 3794 fi; 3795 return fact2; 3796 end; 3797 3798 3799 #subfunction of FindIso 3800 _FindIso:= function(E1,E2,c) 3801 local O,Id,fact; 3802 if not BaseRing(E1) = BaseRing(E2) then 3803 return Error("FindIso: given curves have different base rings"); 3804 elif not Size(E1) = Size(E2) then 3805 return Error("FindIso:= given curves are not isogenous"); 3806 elif IsSupersingular(E1) or IsSupersingular(E2) then 3807 return Error("Algorithm does not work for supersingular curves"); 3808 elif jInvariant(E1) = jInvariant(E2) then 3809 return Isomorphism(E1,E2,1); 3810 fi; 3811 Id:= GetIsoIdeal(E1,E2,c); 3812 if GCD(Norm(Id),c)<>1 then 3813 Id:=FindIso(E1,E2,true); 3814 fi; 3815 fact:= IdealFactorization(Id); 3816 if fact = Sequence([]) then 3817 return Error("something went wrong: Ideal representation of isogeny is ",Id); 3818 fi; 3819 return IsoFromIdeals(fact,E1,E2); 3820 end; 3821 3822 3823 3824 #tries to find an isogeny between E1 and E2 3825 FindIso:= function(E1,E2) 3826 local O,EE1,EE2, PA,BR,E3,E4,c1,c2,OE1,OE2,fact1,fact2,len,i,k,j,Psi,ro,l,E,Iso,Iso1,Iso2,Iso11,Iso21,n,Otmp,r,c,s,kgV,mytime; 3827 mytime:=Realtime(); 3828 if IsSupersingular(E1) or IsSupersingular(E2) then 3829 return Error("Algorithm does not work for supersingular curves"); 3830 fi; 3831 if BaseRing(E1) <> BaseRing(E2) then 3832 return Error("Curves are defined over different Base Rings"); 3833 fi; 3834 if Characteristic(BaseRing(E1))>3 then 3835 EE1:=ShortWNF(E1); 3836 EE2:=ShortWNF(E2); 3837 elif Characteristic(BaseRing(E1)) =2 then 3838 EE1:=NormalFormChar2(E1,1); 3839 EE2:=NormalFormChar2(E2,1); 3840 else return Error("Algorithm not implemented for Characteristic 3"); 3841 fi; 3842 E1:=Base(EE1); 3843 E2:=Base(EE2); 3844 if EndRingIsMaximal(E1) and EndRingIsMaximal(E2) then 3845 Iso:= _FindIso(E1,E2,1); 3846 return Composition(EE2.ext1,Composition(Iso,EE1.ext1)); 3847 fi; 3848 OE1:= EndRing(E1); 3849 if OE1 = FAILURE then 3850 return Error("cannot compute endomorphism ring: conductor too large"); 3851 fi; 3852 OE2:= EndRing(E2); 3853 if OE2 = FAILURE then 3854 return Error("cannot compute endomorphism ring: conductor too large"); 3855 fi; 3856 c1:= Conductor(OE1); 3857 c2:= Conductor(OE2); 3858 kgV:=Lcm(c1,c2); 3859 fact1:=Factorization(Coerce(Integers(),kgV/c1)); 3860 fact2:=Factorization(Coerce(Integers(),kgV/c2)); 3861 Iso1:=[]; 3862 Iso2:=[]; 3863 BR:= BaseRing(E1); 3864 PA:=PolynomialAlgebra(PolynomialAlgebra(BR)); 3865 j:= jInvariant(E1); 3866 E:= E1; 3867 O:=OE1; 3868 s:=Size(E1); 3869 if Length(fact1)>0 then 3870 for i in [1..Length(fact1)] do 3871 l:= fact1[i][1]; 3872 Psi:= Coerce(PA,ModularPolynomial(l)); 3873 for k in [1..fact1[i][2]] do 3874 ro:= Roots(Evaluate(Psi,j)); 3875 if Length(ro) = 1 then 3876 Iso:=IsogenyFromJInvariant(E,l,1); 3877 elif Length(ro) = l+1 then 3878 n:=0; 3879 while n in [0..l] do 3880 n:=n+1; 3881 if ro[n][1] <>0 then 3882 Otmp:= EndRing(ECFromJGivenNumber(ro[n][1],s)); 3883 if Conductor(Otmp)< Conductor(O) then 3884 Iso:= IsogenyFromJInvariant(E,l,n); 3885 if (not Iso = FAILURE) then 3886 n:= l+3; 3887 fi; 3888 fi; 3889 fi; 3890 od; 3891 if n = l+1 then 3892 return Error("no suitable curve found"); 3893 fi; 3894 fi; 3895 if Iso= FAILURE then 3896 return Error("cannot compute isogenous maximal endRing curve for E1"); 3897 fi; 3898 Iso1:= Append(Iso1,[Iso]); 3899 E:= Codomain(Iso); 3900 j:= jInvariant(E); 3901 O:= EndRing(E); 3902 od; 3903 od; 3904 Iso11:= Iso1[1]; 3905 for i in [2..Length(Iso1)] do 3906 Iso11:= Composition(Iso1[i],Iso11); 3907 od; 3908 else Iso11:= Isomorphism(E1,E1,1); 3909 fi; 3910 if Length(fact2)>0 then 3911 E:= E2; 3912 j:= jInvariant(E2); 3913 O:=OE2; 3914 for i in [1..Length(fact2)] do 3915 l:= fact2[i][1]; 3916 Psi:= Coerce(PA,ModularPolynomial(l)); 3917 for k in [1..fact2[i][2]] do 3918 ro:= Roots(Evaluate(Psi,j)); 3919 if Length(ro) = 1 then 3920 Iso:=IsogenyFromJInvariant(E,l,1); 3921 elif Length(ro) = l+1 then 3922 n:=0; 3923 while n in [0..l] do 3924 n:=n+1; 3925 if ro[n][1] <>0 then 3926 Otmp:= EndRing(ECFromJGivenNumber(ro[n][1],s)); 3927 if Conductor(Otmp)< Conductor(O) then 3928 Iso:= IsogenyFromJInvariant(E,l,n); 3929 if not Iso= FAILURE then 3930 n:= l+3; 3931 fi; 3932 fi; 3933 fi; 3934 od; 3935 if n = l+1 then 3936 return Error("something went wrong 2"); 3937 fi; 3938 fi; 3939 if Iso= FAILURE then 3940 return Error("cannot compute isogenous maximal endRing curve for E2"); 3941 fi; 3942 Iso2:= Append(Iso2,[Iso]); 3943 E:= Codomain(Iso); 3944 j:= jInvariant(E); 3945 O:=EndRing(E); 3946 od; 3947 od; 3948 len:= Length(Iso2); 3949 Iso21:= DualIsogeny(Iso2[len]); 3950 for i in [1..len-1] do 3951 Iso21:= Composition(DualIsogeny(Iso2[len-i]),Iso21); 3952 od; 3953 else Iso21:= Isomorphism(E2,E2,1); 3954 fi; 3955 E3:= Codomain(Iso11); 3956 E4:= Domain(Iso21); 3957 c:=Conductor(EndRing(E3)); 3958 Print("init:",Realtime()-mytime,"\n"); 3959 Iso:= _FindIso(E3,E4,c); 3960 Iso:= Composition(Iso,Iso11); 3961 Iso:= Composition(Iso21,Iso); 3962 return Composition(EE2.ext1,(Composition(Iso,EE1.ext1))); 3963 end; 3964 3965 3966 #direkte magma-Uebersetzung 3967 3968 WeierstrassFunction:= function(E,prec) 3969 local Q,z,coeffs,a4,a6,k,c,ctmp,h,ztmp; 3970 if (2*prec+1) >= Characteristic(BaseRing(E)) then return Error( 3971 "Request for division by zero in kernel polynomial computation."); 3972 fi; 3973 Q := LaurentSeriesRing(BaseRing(E),2*prec+2); 3974 z := Q.1; 3975 coeffs:=Coefficients(E); 3976 a4:=coeffs[1]; 3977 a6:=coeffs[2]; 3978 c := [ -a4/5, -a6/7 ]; 3979 for k in [3..prec-1] do 3980 ctmp:=0; 3981 for h in [1..k-2] do 3982 ctmp:=ctmp+c[h]*c[k-1-h]; 3983 od; 3984 Append_(c, [3*ctmp/ ((k-2) * (2*k + 3))]); 3985 od; 3986 ztmp:=0; 3987 for k in [1..prec-1] do 3988 ztmp:= ztmp+c[k]*z^(2*k); 3989 od; 3990 return z^-2 +ztmp; 3991 end; 3992 3993 3994 KernelPolynomial:=function(E,l,Eb,p1) 3995 local FF,xtr,deg,coeffs,coeffs2,e4,e6,e4_tilde,e6_tilde,c,c_tilde,P,x,psi2,dpsi2,mulist,k,pp,np,ii,i,M,B,v_raw,V,p0,v,ps,tlist,g,tmp; 3996 FF := BaseRing(E); 3997 xtr := 2; 3998 deg := Div(l-1, 2); 3999 if (2*(deg+xtr)+3) >= Characteristic(FF) then 4000 return Error( 4001 "Request for division by zero in kernel polynomial computation."); 4002 fi; 4003 coeffs:=Coefficients(E); 4004 coeffs2:=Coefficients(Eb); 4005 if Length(coeffs) <> 2 then 4006 return Error; 4007 fi; 4008 e4:=coeffs[1]; 4009 e6:=coeffs[2]; 4010 e4_tilde:=coeffs2[1]; 4011 e6_tilde := coeffs2[2]; 4012 c := WeierstrassFunction(E,deg+xtr+1); 4013 c_tilde := WeierstrassFunction(Eb,deg+xtr+1); 4014 P := PolynomialRing(FF); x := P.1; 4015 psi2 := 4*x^3 + 4*e4*x + 4*e6; 4016 dpsi2 := 6*x^2 + 2*e4; 4017 mulist := [ 2*e4 + 6*x^2 ]; 4018 for k in [2..deg+xtr] do 4019 pp := Eltseq(mulist[k-1]); 4020 np := pp[2] * dpsi2; 4021 for ii in [3..Length(pp)] do 4022 i := ii - 1; 4023 np := np+ i*pp[i+1]*(dpsi2*x^(i-1) + (i-1)*psi2*x^(i-2)); 4024 od; 4025 Append_(mulist,[ np]); 4026 od; 4027 M:=[]; 4028 for k in [1..deg+xtr] do 4029 pp := Eltseq(mulist[k]); 4030 for i in [1..Length(pp)] do 4031 M[(k-1)*(deg+2+xtr)+i] := pp[i] * 2 / Factorial(2*k); 4032 od; 4033 for i in [Length(pp)+1..deg+2+xtr] do 4034 M[(k-1)*(deg+2+xtr)+i]:=0; 4035 od; 4036 od; 4037 M :=Coerce( KMatrixSpace(FF, deg+xtr, deg+2+xtr),M); 4038 B:=[]; 4039 for i in [1..deg+xtr] do 4040 B[i]:=Coefficient(c_tilde,2*i) - 4041 Coefficient(c,2*i); 4042 od; 4043 B:=Vector(FF,B); 4044 v_raw:= Solution(Transpose(M), B); 4045 V:=v_raw.ext1; 4046 if Dimension(V) <> 2 then 4047 return Error( "Trace term not uniquely determined in kernel polynomial computation."); 4048 fi; 4049 p0 := deg; 4050 v := v_raw + V.1 * (p0 - Eltseq(v_raw)[1]); 4051 v := v + V.2 * (p1 - Eltseq(v)[2]); 4052 ps := Eltseq(v); 4053 tlist:=[]; 4054 for i in [1..deg+1+xtr] do 4055 tlist[i]:=Coerce(FF,0); 4056 od; 4057 tlist[deg+1+xtr] := Coerce(FF,1); 4058 for k in [1..deg+xtr] do 4059 tmp := Coerce(FF,0); 4060 for i in [1..k-1] do 4061 tmp := tmp + tlist[deg+xtr-i+1] * ps[k-i+1]; 4062 od; 4063 tmp := -tmp - ps[k+1]; 4064 tlist[deg+xtr-k+1] := tmp / k; 4065 od; 4066 g := Coerce(P,Sublist(tlist,[xtr+1..deg+xtr+1])); 4067 return g; 4068 end; 4069 4070 4071 4072 4073 PolynomialExpression:=function(S2,S1) 4074 local SS,t,F,PP,x,vS1,vS2,degP,S1_powers_list,S1_i,i,powers_of_S1,P,S2tmp,S1coeff,S2coeff,Pcoeff; 4075 SS := Parent(S1); 4076 t := SS.1; 4077 F := BaseRing(Parent(S1)); 4078 PP := PolynomialRing(F); 4079 x := PP.1; 4080 vS1 := Valuation(S1); 4081 vS2 := Valuation(S2); 4082 degP := Div(vS2,vS1); 4083 if vS1 = 0 or vS2 = 0 or (vS2 mod vS1) <> 0 or degP < 0 then 4084 return Error( "Invalid data in polynomial reconstruction."); 4085 fi; 4086 S1_powers_list := [Coerce(Parent(S1),1)]; 4087 S1_i := Coerce(Parent(S1),1); 4088 for i in [1..Div(vS2,vS1) ] do 4089 S1_i := S1_i * S1; 4090 Append_(S1_powers_list, [S1_i]); 4091 od; 4092 powers_of_S1 := function(i) 4093 return S1_powers_list[i+1]; 4094 end; 4095 P := Coerce(PP,0); 4096 S2tmp := S2; 4097 S1coeff := Coefficient(S1, vS1); 4098 for i in [0..degP] do 4099 S2coeff := Coefficient(S2tmp, vS2 - i*vS1); 4100 Pcoeff := S2coeff / S1coeff; 4101 P := P + Pcoeff * x^(degP - i); 4102 S2tmp := S2tmp - Pcoeff * powers_of_S1(degP - i); 4103 if (i <> degP) and S2tmp = Coerce(Parent(S2tmp),0) then 4104 Print("not enough terms to determine poly"); 4105 return Coerce(PP,0); 4106 fi; 4107 od; 4108 return P; 4109 end; 4110 4111 IsogenyPhi:=function(E,F,psi,prec) 4112 local wE,wF,den,nu; 4113 wE := WeierstrassFunction(E,prec); 4114 wF := WeierstrassFunction(F,prec); 4115 den := Evaluate(psi,wE); 4116 nu := den^2*wF; 4117 return PolynomialExpression(nu,wE); 4118 end; 4119 4120 4121 IsogenyFromKernel:=function(E,F,psi) 4122 local wE,wF,den,nu,phi,wF2,wE2,omega,xmap,ymap,_Iso,re,prec,deg; 4123 deg:=Degree(psi); 4124 prec:=3*deg; 4125 wE := WeierstrassFunction(E,prec); 4126 wF := WeierstrassFunction(F,prec); 4127 den := Evaluate(psi,wE); 4128 nu := den^2*wF; 4129 phi:= PolynomialExpression(nu,wE); 4130 wF2:=Derivative(wF); 4131 wE2:=Derivative(wE); 4132 den:=Evaluate(psi,wE); 4133 nu:=den^3*wF2; 4134 omega:=PolynomialExpression(nu/wE2,wE); 4135 xmap:=function(P) 4136 local x; 4137 x:=PX(P); 4138 if x = INFTY then 4139 return INFTY; 4140 elif Evaluate(psi,x) = 0 then 4141 return INFTY; 4142 else return Evaluate(phi/psi^2,x); 4143 fi; 4144 end; 4145 ymap:=function(P) 4146 local x,y; 4147 x:=PX(P); 4148 y:=PY(P); 4149 return Evaluate(omega/psi^3,x)*y; 4150 end; 4151 _Iso:=function(P) 4152 if xmap(P) = INFTY then 4153 return Zero(F); 4154 else return MakePoint(F,[xmap(P),ymap(P)]); 4155 fi; 4156 end; 4157 re:=Map(E,F,_Iso); 4158 re.operations:= _iso_ops; 4159 re.type:= iso; 4160 re.ypol:=omega/psi^3*EllipticFunctionField(E).1; 4161 re.xpol:=Coerce(Parent(re.ypol),phi/psi^2); 4162 re.degree:=2*deg+1; 4163 F.size:=E.size; 4164 return re; 4165 end; 4166 4167 4168 4169 4170 4171 4172 InstallMethod( 4173 rec( 4174 kind:="FUNCTION", 4175 name:="Degree", 4176 sin:=[[iso,"f"]], 4177 sou:=[[elt-ord^rat,"d"]], 4178 short:= 4179 "The degree of f ", 4180 ex:=[], 4181 see:=[]), _IsoDegree); 4182 4183 InstallMethod( 4184 rec( 4185 kind:="FUNCTION", 4186 name:="Composition", 4187 sin:=[[iso,"phi"],[iso,"psi"]], 4188 sou:=[[iso]], 4189 short:= 4190 "The composition phi(psi()). ", 4191 ex:=[], 4192 see:=[]), _Composition); 4193 InstallMethod( 4194 rec( 4195 kind:="FUNCTION", 4196 name:="Domain", 4197 sin:=[[iso,"phi"]], 4198 sou:=[[elt-grp^ell]], 4199 short:= 4200 "The domain of phi. ", 4201 ex:=[], 4202 see:=[]), _Domain); 4203 4204 4205 InstallMethod( 4206 rec( 4207 kind:="FUNCTION", 4208 name:="Codomain", 4209 sin:=[[iso,"phi"]], 4210 sou:=[[elt-grp^ell]], 4211 short:= 4212 "The codomain of phi. ", 4213 ex:=[], 4214 see:=[]), _Codomain); 4215 4216 InstallMethod( 4217 rec( 4218 kind:="FUNCTION", 4219 name:="Conductor", 4220 sin:=[[ord^num,"O"]], 4221 sou:=[[elt-ord^rat]], 4222 short:= 4223 "The index of O in its maximal order. ", 4224 ex:=[], 4225 see:=[]), _IntCond); 4226