How to resolve the algorithm Lucas-Lehmer test step by step in the PARI/GP programming language

Problem Statement
Step by Step Solution
Sourcecode

Problem Statement

Lucas-Lehmer Test: for

{\displaystyle p}

an odd prime, the Mersenne number

− 1

{\displaystyle 2^{p}-1}

is prime if and only if

− 1

{\displaystyle 2^{p}-1}

divides

S ( p − 1 )

{\displaystyle S(p-1)}

where

S ( n + 1 )

( S ( n )

)

− 2

{\displaystyle S(n+1)=(S(n))^{2}-2}

, and

S ( 1 )

{\displaystyle S(1)=4}

Calculate all Mersenne primes up to the implementation's maximum precision, or the 47th Mersenne prime (whichever comes first).

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Lucas-Lehmer test step by step in the PARI/GP programming language

Source code in the pari/gp programming language

LL(p)={
  my(m=Mod(4,1<
  for(i=3,p,m=m^2-2);
  m==0
};

search()={
  print("2^2-1");
  forprime(p=3,43112609,
    if(LL(p), print("2^"p"-1"))
  )
};

/* ll(p): input odd prime 'p'. */
/* returns '1' if 2^p-1 is a Mersenne prime. */
ll(p) = { 
    /* trial division up to a reasonable depth (time ratio tdiv/llt approx. 0.2) */
    my(l=log(p), ld=log(l));       
    forprimestep(q = 1, sqr(ld)^(l/log(2))\4, p+p,
        if(Mod(2,q)^p == 1, return)
    );
    /* Lucas-Lehmer test with fast modular reduction. */ 
    my(s=4, m=2^p-1, n=m+2); 
    for(i = 3, p,
        s = sqr(s);
        s = bitand(s,m)+ s>>p;
        if(s >= m, s -= n, s -= 2)
    );
    !s
};      /* end ll */


/* get Mersenne primes in range [a,b] */
llrun(a, b) = {
    my(t=0, c=0, p=2, thr=default(nbthreads));
    if(a <= 2,
        c++;
        printf("#%d\tM%d\t%3dh, %2dmin, %2d,%03d ms\n", c, p, t\3600000, t\60000%60, t\1000%60, t%1000);
        a = 3;
    );
    gettime();
    parforprime(p= a, b, ll(p), d,       /* ll(p) -> d  copy from parallel world into real world. */
        if(d,
            t += gettime()\thr;
            c++;    
            printf("#%d\tM%d\t%3dh, %2dmin, %2d,%03d ms\n", c, p, t\3600000, t\60000%60, t\1000%60, t%1000)
        )
    )
};   /* end llrun */


\\  export(ll);     /* if running ll as script */


llrun(2, 132049)