How many ways can you code an LL test
 

including examples in PARI/gp 2.8.0 we have:

[CODE]lucas(p) =
{ my(u = 4, q = 1<<p - 1);
for(k=3, p, u = (sqr(u)-2) % q);
u == 0;
}[/CODE]

[CODE]lucas1(p) =
{ my(u = 4, q = 1<<p - 1,k=2);
while(k<p, u = (sqr(u)-2) % q;k++);
u == 0;
}[/CODE]

[CODE]lucas2(p) =
{ my(u = 4, q = 1<<p - 1,k=2);
until(k==p, u = (sqr(u)-2) % q;k++);
u == 0;
}[/CODE]

[CODE]lucas3(p) =
{ my(u = 4, q = 1<<p - 1,k=2);
until(k==p, u = (norm(u)-2) % q;k++);
u == 0;
}[/CODE]

then without just changing the example given we can do:
[CODE]
LL1(p,d)={
my(s=6,P=1<<p-1);
for(x=1,p-1,
s=(fromdigits(d,if(s%P>1,s%P,s))-1)%P
);
s==0
}[/CODE]

where d= digits(241) =[2,4,1]
and
[CODE]
LL2(p,d)={
my(s=6,P=1<<p-1);
for(x=1,p-1,
s=subst(Pol([2,4,0],y),y,if(s%P>1,s%P,s))%P
);
s==0
}[/CODE]

[CODE]lucaslehmer(p) ={
s=Mod(4,2^p-1);
for(x=1,p-2,
s=s^2-2
);
if(s==0,1,0)
}[/CODE]

etc. how many ways can you code it ?

A short and faster example, still in PARI-GP:
[CODE]
mylucas(p)={mp=2^p-1;
x=4;for(i=1,p-2,x=sqr(x);
hi=shift(x,-p);lo=bitand(x,mp);
x=lo+hi-2;
if(x>=mp,x-=mp);if(x<0,x+=mp));
return(x==0)}
[/CODE]

[QUOTE=R. Gerbicz;420686]A short and faster example, still in PARI-GP:
[CODE]
mylucas(p)={mp=2^p-1;
x=4;for(i=1,p-2,x=sqr(x);
hi=shift(x,-p);lo=bitand(x,mp);
x=lo+hi-2;
if(x>=mp,x-=mp);if(x<0,x+=mp));
return(x==0)}
[/CODE][/QUOTE]

wow that is fast ( like 3 times how fast lucas runs on my system) and potentially parallelizable although I don't see how bitand(x,mp) will be different than x unless x>mp maybe that could be coded in as a speedup. also your two separate if statements about x can use the switch case way of writing it.

The C+GMP example on RosettaCode shows quite a few basic pretests, and also speeding up the modulo (since it is a special form), which it looks like the admirably short example from Gerbicz does.

[QUOTE=danaj;420719]The C+GMP example on RosettaCode shows quite a few basic pretests, and also speeding up the modulo (since it is a special form), which it looks like the admirably short example from Gerbicz does.[/QUOTE]

I know a few things that speed it up still though, at least for my machine. Changing sqr() to norm(), changing shift(x,-p) to x>>p seems to speed it up ( my thought is the doesn't make it figure out -p in gp. same with i's range it can be done with i=3,p so that gp can just recall the value of p. and lo=bitand(mp,x) and checking mp<=x in the first if got me a slight if not imaginary speedup ( seems to speed the whole test up by a few hundred milliseconds though in theory so does not checking mylucas for sophie germain primes. either from within mylucas itself or before calling it in the test should speed things up) as well in my testing with:

[CODE]parforprime(x=1,9941,c=mylucas(x);if(c,print(x)))[/CODE]

edit:and of course I wish I knew enough about FFT to code it. or get the -2 for free. I can't find the convolution stuff on wikipedia/mersennewiki I read at one point.

[QUOTE=R. Gerbicz;420686]A short and faster example, still in PARI-GP:
[CODE]
mylucas(p)={mp=[COLOR=Red]1<<p[/COLOR]-1;
x=4;for(i=1,p-2,x=sqr(x);
hi=shift(x,-p);lo=bitand(x,mp);
x=lo+hi-2;
if(x>=mp,x-=mp);if(x<0,x+=mp));
return(x==0)}
[/CODE][/QUOTE]
Fixed that for you, for few microseconds gain :razz:

[QUOTE=LaurV;420779]Fixed that for you, for few microseconds gain :razz:[/QUOTE]

I think you'll find that it depends on what version of PARI/gp and architecture example:

[CODE](23:25) gp > mylucas(p)={mp=1<<p-1;
x=4;for(i=1,p-2,x=sqr(x);
hi=shift(x,-p);lo=bitand(x,mp);
x=lo+hi-2;
if(x>=mp,x-=mp);if(x<0,x+=mp));
return(x==0)}
%136 = (p)->mp=1<<p-1;x=4;for(i=1,p-2,x=sqr(x);hi=shift(x,-p);lo=bitand(x,mp);x=lo+hi-2;if(x>=mp,x-=mp);if(x<0,x+=mp));return(x==0)
(23:38) gp > parforprime(x=1,9941,if(x%4==3 && isprime(x<<1+1),0,factor(1<<x-1,1<<min(x\2,12))[1,1]!=[;],0,mylucas(x)),c,if(c,print(x)))
5
7
13
17
19
31
61
89
107
127
521
607
1279
2203
2281
3217
4253
4423
9689
9941
(23:39) gp > ##
*** last result computed in 30,161 ms.
[/CODE]
[CODE]
(23:39) gp > mylucas(p)=mp=2^p-1;x=4;for(i=1,p-2,x=sqr(x);hi=shift(x,-p);lo=bitand(x,mp);x=lo+hi-2;if(x>=mp,x-=mp);if(x<0,x+=mp));return(x==0)
%138 = (p)->mp=2^p-1;x=4;for(i=1,p-2,x=sqr(x);hi=shift(x,-p);lo=bitand(x,mp);x=lo+hi-2;if(x>=mp,x-=mp);if(x<0,x+=mp));return(x==0)
(23:39) gp > parforprime(x=1,9941,if(x%4==3 && isprime(x<<1+1),0,factor(1<<x-1,1<<min(x\2,12))[1,1]!=[;],0,mylucas(x)),c,if(c,print(x)))
5
7
13
17
19
31
61
89
107
127
521
607
1279
2203
2281
3217
4253
4423
9689
9941
(23:40) gp > ##
*** last result computed in 30,023 ms.[/CODE]

so technically that change on mine ( 64 bit OS, i3 2130 processor, and 64 bit 2.8.0 (development 18332-50485cf) gp) slows it down if doing runs of more than one number.

[QUOTE=R. Gerbicz;420686]A short and faster example, still in PARI-GP:
[CODE]
mylucas(p)={mp=2^p-1;
x=4;for(i=1,p-2,x=sqr(x);
hi=shift(x,-p);lo=bitand(x,mp);
x=lo+hi-2;
if(x>=mp,x-=mp);if(x<0,x+=mp));
return(x==0)}
[/CODE][/QUOTE]

[QUOTE=LaurV;420779]Fixed that for you, for few microseconds gain :razz:[/QUOTE]

Another micro optimization is to get rid of [CODE]if(x<0,x+=mp)[/CODE]

[QUOTE=science_man_88;420677]etc. how many ways can you code it ?[/QUOTE]
Countably infinite.

[QUOTE=LaurV;420779]Fixed that for you, for few microseconds gain :razz:[/QUOTE]

I'd like to give a shorter and cleaner code. And actually I don't see any speedup, tried 6 runs to get x=1<<100000000-1; and my version of y=2^100000000-1;
it is quite likely that pari is using shifts for the calculation of 2^n and no multiplication.

[QUOTE=axn;420785]Another micro optimization is to get rid of [CODE]if(x<0,x+=mp)[/CODE][/QUOTE]

Yes, it is true, but you need a proof that you can skip that.

[QUOTE=R. Gerbicz;420816]
Yes, it is true, but you need a proof that you can skip that.[/QUOTE]

my basic thought is it would be better to have another condition to test because the only time I see that x<0 is when x<2 if the original squaring leaves it at 2 then x-2 is 0 only when x=0 ( the termination condition) or x=1 can it go negative (x=1 leading to -1,-1,-1,-1,-1 and x=0 leading to -2,2,2,2,2,2,2) x will never grow to more than [TEX]mp^2 = 2^{2p}-(2^{p+1})+1[/TEX] if the value is kept in check (edit: or since hi and lo are used 2^(p+1)-2 because that's what it brings it less than). the first check will only produce a value that needs the second check if 1 or 2^p mod (2^p-1) are possible x values. of course other than -1 the next step will take care of this.

[QUOTE=xilman;420808]Countably infinite.[/QUOTE]

okay then what's the absolute fastest way you can code it.

[QUOTE=science_man_88;420843]okay then what's the absolute fastest way you can code it.[/QUOTE]

If by 'fast' you are referring to execution speed of the tests, several of us have been working on that for many years, and continue to do so. You are welcome to compare timings of your efforts to ours.

If OTOH by 'fast' you mean 'speed of coding it up', then you got us beat, hands-down. (But only if you neglect the coding time others put into the bignum package you are using.)

[QUOTE=ewmayer;420926]If by 'fast' you are referring to execution speed of the tests, several of us have been working on that for many years, and continue to do so. You are welcome to compare timings of your efforts to ours.

If OTOH by 'fast' you mean 'speed of coding it up', then you got us beat, hands-down. (But only if you neglect the coding time others put into the bignum package you are using.)[/QUOTE]

mine are all PARI/gp and suck even compared to the lucas example.

I've thought about the fact that in the squaring version if[TEX] s_n=a{M_n}[/TEX] then [TEX]s_{n+1}=a^2{M_n}^2-2[/TEX] which mod the next mersenne number is [TEX]a^22^{n-1}-2[/TEX] but that's not all that helpful when a's value is the larger part of that typically and in the [TEX]s_n=2{s_{n-1}}^2-1[/TEX] version(s) of the test you can find the relation that the next value is also [TEX]2{s_{n-1}-1}^2+4{s_{n-1}-1}+1[/TEX] which is like the way mersenne numbers that have exponents in a cunningham chain of the first kind follow but I have a hard time finding a useful relationship to use to mod that one.

Here's another crappy version:
[CODE]LL(p)=my(s=Mod(2,2^p-1));for(x=1,p-2,s=vecprod([2,s-1,s+1])+1);s[/CODE]

Crappy-performing but functional bc code:
[code]define lltest(p) {
auto i,m,r; i = p-2; m = 2^p-1; r = 4;
while(i--) {
r = (r^2-2)%m;
}
return(r==0);
}[/code]

In Ruby....

[CODE]def LL(p)
s, n = 4, 2**p-1
3.upto(p) { s=(s**2-2)%n }
s==0
end

print "LL exponent : "
p = gets.to_i
puts LL(p)
[/CODE]

Earlier in the thread someone mentioned the Rosetta Code page, with dozens of languages. None particularly efficient, though.

[url]http://rosettacode.org/wiki/Lucas-Lehmer_test[/url]

[QUOTE=GP2;494490]Earlier in the thread someone mentioned the Rosetta Code page, with dozens of languages. None particularly efficient, though.

[url]http://rosettacode.org/wiki/Lucas-Lehmer_test[/url][/QUOTE]

Ah! I did not handle exponent 2. Also [c](p-2).times[/c] maybe neater then [c]3.upto(p)[/c].

I guess the linked page is too small for George's source :wink:

[QUOTE=paulunderwood;494491]Ah! I did not handle exponent 2.[/QUOTE]

No need - LL assumes p an odd prime.