1maple_mode(1);
2 M:=t->(cos(t)/sin(t)^3,sin(t)/sin(t)^3);
3 f:=(x,y,z)->x^3-3*x*y*z-z^3;
4 C2:=implicitplot(f(x,y,1),x,y);
5/*  Pour ajuster le pas d'un trace implicit on utilise xstep et ystep�pour abscisse et ordonnee, quels que soient les noms de variables.*/
6 C2Z:=implicitplot(f(x,1,z),x,z,xstep=0.01,ystep=0.01);
7/* ---------------------------------------------------------------------------------� Mini texte Courbes elliptiques�---------------------------------------------------------------------------------*/
8/* Etudier la documentation de factorint dans l'aide de pari-gp: Il combine�plusieurs methodes, dont celle de Lenstra: ECM*/
9/*  egalite de 2 vecteurs a facteur pres:*/
10 egal:=proc(P,Q)
11 a:=expand(P[1]*Q[2]-P[2]*Q[1]);
12 b:=expand(P[3]*Q[2]-P[2]*Q[3]);
13 c:=expand(P[1]*Q[3]-P[3]*Q[1]);
14 if [a,b,c]=[0,0,0] then true; else false; fi;
15end proc;
16/* On cree une procedure simplif que l'on pourra modifier ensuite:*/
17simplif:=proc(R)
18d:=igcd(R[1],R[2],R[3]);
19normal(R/d);end proc;
20/* On prend un courbe  C d'equation: y^2=x^3+ax^2+bx+c*/
21 pointresiduel:=proc(a,b,c,P,Q)
22local d,C,QQ,M,sol,R;
23 C:=-y^2*z+x^3+a*x^2*z+b*x*z^2+c*z^3;
24if egal(P,Q) then
25//la tangente est P+tQQ sauf si P est omega. Autrement dit, QQ est le point
26//d'intersection de la tangente avec la droite de l'infini, ie c'est le vecteur
27//directeur. Sauf lorsque la tangente est la droite de l'infini ie QQ=omega
28QQ:=[subs(x=P[1],y=P[2],z=P[3],diff(C,y)),subs(x=P[1],y=P[2],z=P[3],-diff(C,x)),0];
29M:=expand(P+t*QQ);d:=subs(x=M[1],y=M[2],z=M[3],C);
30sol:=normal(d/t^2);
31//le vecteur QQ doit etre non nul sinon c'est que la courbe est singuliere.
32if egal(QQ,[0,1,0]) then R:=[0,1,0]; else R:=expand(-coeff(sol,t,1)*P+coeff(sol,t,0)*QQ);fi;
33else
34//On parametre (PQ) par P+tQ, 0 et l'infini sont solutions.
35M:=expand(P+t*Q);d:=subs(x=M[1],y=M[2],z=M[3],C);
36//sol doit etre un poly de degre 1 puisque 0 et l'infini sont solutions.
37sol:=normal(d/t);
38//ce vecteur doit etre non nul.
39R:=expand(-coeff(sol,t,1)*P+coeff(sol,t,0)*Q);
40fi;
41end proc;
42/* exemple:*/
43P:=[0,0,1];Q:=[1,1,1];omega:=[0,1,0];
44/* verification egal:*/
45 if egal(P,Q) then print(eg) else  print(noneg) fi;
46 if egal(P,P) then print(eg) else  print(noneg) fi;
47 pointresiduel(1,-1,0,P,Q);pointresiduel(1,-1,0,P,omega);pointresiduel(1,-1,0,P,omega);
48 pointresiduel(1,-1,0,P,P);pointresiduel(1,-1,0,omega,omega);
49/*  addition:*/
50plus:=proc(a,b,c,P,Q)
51local R;
52R:=pointresiduel(a,b,c,P,Q);
53R[2]:=-R[2];
54R;
55end proc;
56plus(1,-1,0,omega,omega);plus(1,-1,0,P,omega);plus(1,-1,0,P,P);R:=plus(1,-1,0,P,Q);
57/* On programme la mult par n en O(log n) additions.*/
58nplus:=proc(a,b,c,P,n)
59local Y,m,X;
60Y:=[0,1,0];X:=P;m:=n;
61while m>0 do
62 if odd(m) then Y:=plus(a,b,c,X,Y);X:=plus(a,b,c,X,X);m:=(m-1)/2;
63  else X:=plus(a,b,c,X,X);m:=m/2;
64  fi;
65  end do;
66  Y;
67end proc;
68/*  Ex: P est d'ordre 2, Q d'ordre 3, P+Q d'ordre 6 */
69 nplus(1,-1,0,P,2),nplus(1,-1,0,Q,2),nplus(1,-1,0,Q,3),nplus(1,-1,0,R,3),nplus(1,-1,0,R,2),nplus(1,-1,0,R,6);
70 pari()
71 [x1,x2]:=pari_ellpow([0,1,0,2,-15],[2,1],5);
72 [y1,y2,y3]:=nplus(1,2,-15,[2,1,1],5);
73 normal([x1,x2]-[y1/y3,y2/y3]);#doit etre nul
74/* Application \`a la  factorisation.*/
75egalomega:=proc(P,N)
76 Q:=[0,1,0];
77 a:=expand(P[1]*Q[2]-P[2]*Q[1]);
78 b:=expand(P[3]*Q[2]-P[2]*Q[3]);
79 c:=expand(P[1]*Q[3]-P[3]*Q[1]);
80 g:=igcd(a,b,c,N);
81 if g<>1 and g<>N  then print("diviseur",g):true; else false; fi;
82end proc;
83//On redefinit simplifier pour travailler modulo N
84simplif:=proc(R)
85d:=igcd(R[1],R[2],R[3],N);
86 if d<>1 and d<>N then print("diviseur",d); break;
87else  R mod N;fi;end proc;
88/* On redefinit plus pour qu'il s'arrete si l'on obtient�omega modulo un diviseur strict de N.*/
89plus:=proc(a,b,c,P,Q)
90R:=pointresiduel(a,b,c,P,Q);
91R[2]:=-R[2];
92if egalomega(R,N) then break; fi;
93R;
94end proc;
95/* factorisation d'un entier via les courbes elliptiques*/
96n:=lcm(seq(i,i=1..19));
97/* On cherche c tel que (2,1) soit sur y^2=x^3+bx+c*/
98P:=[2,1,1];N:=nextprime(20000)*nextprime(40000);
99N:=nextprime(1234567)*nextprime(7654321);n:=lcm(seq(i,i=1..50));
100