20151231, 13:05  #1 
"Forget I exist"
Jul 2009
Dumbassville
20300_{8} Posts 
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:
lucas1(p) = { my(u = 4, q = 1<<p  1,k=2); while(k<p, u = (sqr(u)2) % q;k++); u == 0; } Code:
lucas2(p) = { my(u = 4, q = 1<<p  1,k=2); until(k==p, u = (sqr(u)2) % q;k++); u == 0; } Code:
lucas3(p) = { my(u = 4, q = 1<<p  1,k=2); until(k==p, u = (norm(u)2) % q;k++); u == 0; } Code:
LL1(p,d)={ my(s=6,P=1<<p1); for(x=1,p1, s=(fromdigits(d,if(s%P>1,s%P,s))1)%P ); s==0 } and Code:
LL2(p,d)={ my(s=6,P=1<<p1); for(x=1,p1, s=subst(Pol([2,4,0],y),y,if(s%P>1,s%P,s))%P ); s==0 } Code:
lucaslehmer(p) ={ s=Mod(4,2^p1); for(x=1,p2, s=s^22 ); if(s==0,1,0) } Last fiddled with by science_man_88 on 20151231 at 13:27 
20151231, 13:52  #2 
"Robert Gerbicz"
Oct 2005
Hungary
31·47 Posts 
A short and faster example, still in PARIGP:
Code:
mylucas(p)={mp=2^p1; x=4;for(i=1,p2,x=sqr(x); hi=shift(x,p);lo=bitand(x,mp); x=lo+hi2; if(x>=mp,x=mp);if(x<0,x+=mp)); return(x==0)} 
20151231, 14:22  #3 
"Forget I exist"
Jul 2009
Dumbassville
2^{6}·131 Posts 
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.

20151231, 19:42  #4 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
1614_{8} Posts 
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.

20151231, 20:31  #5  
"Forget I exist"
Jul 2009
Dumbassville
2^{6}×131 Posts 
Quote:
Code:
parforprime(x=1,9941,c=mylucas(x);if(c,print(x))) Last fiddled with by science_man_88 on 20151231 at 20:36 

20160101, 03:34  #6 
Romulan Interpreter
Jun 2011
Thailand
3^{3}×347 Posts 

20160101, 03:45  #7 
"Forget I exist"
Jul 2009
Dumbassville
8384_{10} Posts 
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<<p1; x=4;for(i=1,p2,x=sqr(x); hi=shift(x,p);lo=bitand(x,mp); x=lo+hi2; if(x>=mp,x=mp);if(x<0,x+=mp)); return(x==0)} %136 = (p)>mp=1<<p1;x=4;for(i=1,p2,x=sqr(x);hi=shift(x,p);lo=bitand(x,mp);x=lo+hi2;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<<x1,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:
(23:39) gp > mylucas(p)=mp=2^p1;x=4;for(i=1,p2,x=sqr(x);hi=shift(x,p);lo=bitand(x,mp);x=lo+hi2;if(x>=mp,x=mp);if(x<0,x+=mp));return(x==0) %138 = (p)>mp=2^p1;x=4;for(i=1,p2,x=sqr(x);hi=shift(x,p);lo=bitand(x,mp);x=lo+hi2;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<<x1,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. Last fiddled with by science_man_88 on 20160101 at 03:51 
20160101, 04:11  #8  
Jun 2003
132F_{16} Posts 
Quote:
Code:
if(x<0,x+=mp) 

20160101, 09:22  #9 
Bamboozled!
"πΊππ·π·π"
May 2003
Down not across
2^{2}×3×887 Posts 

20160101, 09:49  #10 
"Robert Gerbicz"
Oct 2005
Hungary
2661_{8} Posts 
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<<1000000001; and my version of y=2^1000000001;
it is quite likely that pari is using shifts for the calculation of 2^n and no multiplication. Yes, it is true, but you need a proof that you can skip that. Last fiddled with by R. Gerbicz on 20160101 at 09:50 
20160101, 12:55  #11 
"Forget I exist"
Jul 2009
Dumbassville
2^{6}×131 Posts 
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 x2 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 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^p1) are possible x values. of course other than 1 the next step will take care of this.
Last fiddled with by science_man_88 on 20160101 at 13:19 
Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
The Joys of Cracked.com: 5 Ways We Ruined the Occupy Wall Street Generation  Dubslow  Soap Box  17  20120514 08:51 
RDS's unique pedagogic ways  R.D. Silverman  Soap Box  137  20120107 07:52 
ways to get rid of oil spills  science_man_88  Puzzles  9  20100730 21:22 
Help needed to test Athlon64 code  geoff  Programming  7  20060818 12:16 
Creative ways to achieve Athlon 64 / Opteron optimization  GP2  Hardware  11  20040121 03:01 