The cubic Frobenius primality test as a Pari/gp script

I've written the cubic Frobenius primality test as a Pari/gp script. It can be loaded with the command:
bash$  gp script.gp

and compiled with the  command:
bash$ gp2c-run -g script.gp

I tested the compiled script with four Mersenne primes. An extract of the gp session as well as 
the script used are copied below.




david@HPLaptop:~/mar31/dec04$ gp2c-run  -g my_script2900a.gp


                                                        GP/PARI CALCULATOR Version 2.15.3 (released)
                                        

parisizemax = 16000000000, primelimit = 500000
? isfrob3select(2^521-1)    \\ Mersenne prime
%1 = 1
? isfrob3select(2^11213-1)  \\ Mersenne prime
%2 = 1
? ##
  ***   last result computed in 7,642 ms.
? isfrob3select(2^19937-1)   \\ Mersenne prime
%3 = 1
? ##
  ***   last result computed in 33,256 ms.
? isfrob3select(2^86243-1)   \\ Mersenne prime
%4 = 1
? ##
  ***   last result computed in 18min, 26,963 ms.






Copy of the script my_script2900.gp:


isfrob3select(n) =
{
    if(((n%2)==0)&&(n>2), return(0));
    my(conductors = [7, 9, 13, 19, 37, 61, 67, 79, 97, 103, 139, 163, 199, 241, 271, 337, 349, 379, 409, 421, 463, 523, 1087], 
       r_values = [7, 3, 13, 19, 37, 61, 67, 79, 97, 103, 139, 163, 199, 241, 271, 337, 349, 379, 409, 421, 463, 523, 1087],  
       s_values = [7, 1, 13, 19, 37, 183, 201, 79, 97, 309, 139, 163, 199, 1205, 813, 2359, 349, 1895, 2045, 2947, 3241, 1569, 7609]);

    for (i = 1, #conductors,
        my(N = conductors[i], r = r_values[i], s = s_values[i], rootD=sqrtint(4*r^3-27*s^2));  
        if((gcd(n,N)>1)&& (n>r), return(0));    
        my(witness = n % N);
        while (!isprime(witness), witness += N);
        
        if (polisirreducible(Mod(1, witness)*X^3 - Mod(r, witness)*X - Mod(s, witness)),            
            my(M = matrix(3, 3), v3 = matrix(3, 1));
            M[1, 1] = Mod(0, n); M[1, 2] = Mod(r, n); M[1, 3] = Mod(s, n);
            M[2, 1] = Mod(1, n); M[2, 2] = M[1, 1]; M[2, 3] = M[1, 1];
            M[3, 1] = M[1, 1]; M[3, 2] = Mod(1, n); M[3, 3] = M[1, 1];
            v3[1, 1] = Mod(2*r, n); v3[2, 1] = Mod(0, n); v3[3, 1] = Mod(3, n);          
            my(v = digits(n, 2), l = length(v), Mpp = M, ap2);
            for (j = 2, l,
                ap2 = Mpp * Mpp;
                if (v[j], ap2 = M * ap2);
                Mpp = ap2;
            );
            my(output = Mpp * v3);           
            my(test1 = (output[2, 1] == Mod(-r, n)),
               test2 = (output[3, 1] == Mod(0, n)),
               a=((-3*s + rootD)/2)%n, b=((-3*s - rootD)/2)%n, test3=((output[1,1]==Mod(a,n))||(output[1,1]==Mod(b,n))));
            if (test1 && test2 && test3, return(1)); 
        )
    );
    return(0); 
}