EliptDva[m_,a4_,a6_,P_,Q_]:=
Module[{a=P[[1]],b=P[[2]],c=Q[[1]],d=Q[[2]],s,t,e,f},
If[P=={m,m},Print["P=0"];Return[Q]];
If[Q=={m,m},(*Print["Q=0"];*)Return[P]];
If[Mod[b^2,m]≠Mod[a^3+a4*a+a6,m] ,Print["?P"];Return[{m,0}]];
If[Mod[d^2,m]≠Mod[c^3+a4*c+a6,m] ,Print["?Q"];Return[{m,0}]];
If[Mod[a-c,m]≠0,
s=Mod[(b-d)*PowerMod[(a-c),-1,m],m];
t=Mod[(b*c-d*a)*PowerMod[(c-a),-1,m],m];
e=Mod[s^2-a-c,m];f=Mod[-s*e-t,m];
Return[{e,f}],
If[Mod[b-d,m]≠0,
Return[{m,m}],
s=Mod[(3*a^2+a4)*PowerMod[(2*b),-1,m],m];
t=Mod[(-a^3+a4*a+2*a6)*PowerMod[(2*b),-1,m],m]
]; (* b≠d *)
e=Mod[s^2-a-c,m];f=Mod[-s*e-t,m];
Return[{e,f}];
]; (* a≠c *)
] (* Module *)

Print["E:y^2=x^3+x+1, m=5, generator (2,3):"]
R1={2,3};
R2=R1;
m=5;a4=0;a6=1;
For[i=2,i≤8,i++,
R2=EliptDva[m,a4,a6,R1,R2];
Print[i,"(2,3)=",R2];
];
Print["E:y^2=x^3+x+1, m=109, generator (0,1):"]
R1={0,1}
R2=R1;
m=109;a4=1;a6=1;
For[i=2,i≤125,i++,
R2=EliptDva[m,a4,a6,R1,R2];
If[R2=={39,45},Print[i,R1,"=",R2]];
];