Re: Use Pari/gp for learning Elliptic Curves and ECPP

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Thu, Jul 13, 2023 at 11:03:35PM +0200, tony.reix@laposte.net wrote:
> Hi,
> 
> I'd like to learn Elliptic Curves and ECPP by using Pari/gp.
> I'm now reading L. C. Washington's book.
> 
> I've already found:
>    http://pari.math.u-bordeaux.fr/Events/PARI2018b/talks/elliptic.pdf‌;
>    http://pari.math.u-bordeaux.fr/dochtml/html/Elliptic_curves.html
>    http://pari.math.u-bordeaux.fr/Events/PARI2017c/talks/ecc_en.pdf‌‌;
>    http://pari.math.u-bordeaux.fr/Events/PARIday2021/talks/ellrank.pdf‌‌;
>    http://pari.math.u-bordeaux.fr/Events/PARI2018/talks/ecpp.pdf‌;
> that I'll have to read and experiment with.
> Are there more Pari/gp documents I should use ?
> 
> About the last paper (by Jared Asuncion), it talks about a  "ecpp()" function, which is not provided by Pari/gp AFAIK.
> Is there some Pari/gp code implementing ECPP somewhere ?
> 
> My goal is to understand and experiment ECPP with Pari/gp with a special kind
> of numbers (Wagstaff numbers) and see if I can find "constants" for different
> instances of this kind of numbers.

Try this script that generates a prime certificate without specifying the prime number.
This works this way.

Start with a ~128bit prime, here 2^127-1
do
C=primecert(2^127-1);
then do several time in a row
C=concat(backecpp(C[1][1]),C)
to enlarge the certificate.
The prime number it certifies is C[1][1].
You can do
? primecertisvalid(C)
%2 = 1
to prove that the number is prime.

Each time you call backecpp the certified prime has nearly the double number of bits
of the previous one and you get an extremly compact certificate.

Cheers,
Bill.

check(N,D,m,q)=
{
  my(T,j,E,P,P2,P1,g);
  T=polclass(D)*Mod(1,N);
  j=polrootsmod(T,N)[1];
  E=ellinit(ellfromj(j));
  P=random(E); P2=ellmul(E,P,(m\q)); P1=ellmul(E,P2,q);
  if (#P2==2 && #P1==1, return(P));
  g=Mod(1,N);
  until(!issquare(g),g++);
  E=ellinit([g^2*E.a4,g^3*E.a6]);
  P=random(E); P2=ellmul(E,P,(m\q)); P1=ellmul(E,P2,q);
  if (#P2==2 && #P1==1, return(P));
  0;
}

backecpp(p)=
{
  my(check=check);
  my(e=logint(p,2));
  parfor(i=1,oo,
    my(c=p*(2*i+1)*2^(e-logint(2*i+1,2)-1));
    my(V=qfbsolve(Qfb(1,0,3),4*c,2));
    if(#V,
      forstep(e=-1,1,2,
        my(t=e*V[1]+2,q=t+c-1);
        if(ispseudoprime(q),return([[q,t,c/p,0,lift(check(q,-3,c,p))]])))),
    E,if(E,return(E)));
}