Thanks to the people who provided the references for gmp.
My own (very simple) gmp-based Lucas-Lehmer program follows below:
I wrote it mostly in order to make sure I understood what mprime was doing,
especially the Res64. As I mentionned, it is *WAY* slower than mprime,
even for primes far below 1,000,000. This is not at all surprising,
of course, given the algorithm that gmp uses for multiplication.
Even if gmp had a faster algorithm, mprime would *still* be faster
since it was optimized in machine-language.
Oh, btw, this was done in Linux. It probably works anywhere you have
gmp and a C compiler, though.
==========================================================================
/*************************************************************************
lucaslehmer.c
Compute this program thus:
gcc lucaslehmer.c -lm -lgmp -O3 -Wall -o lucaslehmer
************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
unsigned long int p, i, j;
mpz_t pp, ll, aux;
int
main(int argc, char *argv[])
{
extern unsigned long int p, i, j;
extern mpz_t pp, ll, aux;
int k;
mpz_init(pp);
mpz_init(ll);
mpz_init(aux);
if ( argc == 1 )
{
fprintf(stderr,"lucaslehmer: usage: lucaslehmer p1 p2 ... pn\n");
}
else for ( k = 1 ; k < argc ; k++ )
{
p = atol(argv[k]);
mpz_set_ui(ll,(long)(4));
mpz_ui_pow_ui(aux,(long)(2),p);
mpz_sub_ui(pp,aux,(long)(1));
for ( j = 1 ; j + 2 <= p ; j++ )
{
mpz_powm_ui(aux,ll,(long)(2),pp);
mpz_sub_ui(ll,aux,(long)(2));
}
if ( mpz_sgn(ll) == 0 )
{
fprintf(stdout,"M%ld is a prime!\n",p);
}
else
{
mpz_fdiv_r_2exp(aux,ll,(long)(64));
fprintf(stdout,"M%ld is not a prime. Res64: ",p);
mpz_out_str(stdout,(int)(16),aux);
fprintf(stdout,"\n");
}
}
return(0);
}
