mersenneforum.org How many ways can you code an LL test
 Register FAQ Search Today's Posts Mark Forums Read

 2015-12-31, 13:05 #1 science_man_88     "Forget I exist" Jul 2009 Dumbassville 203008 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<

1,s%P,s))-1)%P ); s==0 } where d= digits(241) =[2,4,1] and Code: LL2(p,d)={ my(s=6,P=1<1,s%P,s))%P ); s==0 } Code: lucaslehmer(p) ={ s=Mod(4,2^p-1); for(x=1,p-2, s=s^2-2 ); if(s==0,1,0) } etc. how many ways can you code it ? Last fiddled with by science_man_88 on 2015-12-31 at 13:27

 2015-12-31, 13:52 #2 R. Gerbicz     "Robert Gerbicz" Oct 2005 Hungary 31·47 Posts 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)}
2015-12-31, 14:22   #3
science_man_88

"Forget I exist"
Jul 2009
Dumbassville

26·131 Posts

Quote:
 Originally Posted by R. Gerbicz 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)}
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.

 2015-12-31, 19:42 #4 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 16148 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.
2015-12-31, 20:31   #5
science_man_88

"Forget I exist"
Jul 2009
Dumbassville

26×131 Posts

Quote:
 Originally Posted by danaj 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.
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)))
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.

Last fiddled with by science_man_88 on 2015-12-31 at 20:36

2016-01-01, 03:34   #6
LaurV
Romulan Interpreter

Jun 2011
Thailand

33×347 Posts

Quote:
 Originally Posted by R. Gerbicz A short and faster example, still in PARI-GP: Code: mylucas(p)={mp=1<=mp,x-=mp);if(x<0,x+=mp)); return(x==0)}
Fixed that for you, for few microseconds gain

2016-01-01, 03:45   #7
science_man_88

"Forget I exist"
Jul 2009
Dumbassville

838410 Posts

Quote:
 Originally Posted by LaurV Fixed that for you, for few microseconds gain
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:
(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.
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.

Last fiddled with by science_man_88 on 2016-01-01 at 03:51

2016-01-01, 04:11   #8
axn

Jun 2003

132F16 Posts

Quote:
 Originally Posted by R. Gerbicz 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)}
Quote:
 Originally Posted by LaurV Fixed that for you, for few microseconds gain
Another micro optimization is to get rid of
Code:
if(x<0,x+=mp)

2016-01-01, 09:22   #9
xilman
Bamboozled!

"πΊππ·π·π­"
May 2003
Down not across

22×3×887 Posts

Quote:
 Originally Posted by science_man_88 etc. how many ways can you code it ?
Countably infinite.

2016-01-01, 09:49   #10
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

26618 Posts

Quote:
 Originally Posted by LaurV Fixed that for you, for few microseconds gain
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:
 Originally Posted by axn Another micro optimization is to get rid of Code: if(x<0,x+=mp)
Yes, it is true, but you need a proof that you can skip that.

Last fiddled with by R. Gerbicz on 2016-01-01 at 09:50

2016-01-01, 12:55   #11
science_man_88

"Forget I exist"
Jul 2009
Dumbassville

26×131 Posts

Quote:
 Originally Posted by R. Gerbicz Yes, it is true, but you need a proof that you can skip that.
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 $mp^2 = 2^{2p}-(2^{p+1})+1$ 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.

Last fiddled with by science_man_88 on 2016-01-01 at 13:19

 Similar Threads Thread Thread Starter Forum Replies Last Post Dubslow Soap Box 17 2012-05-14 08:51 R.D. Silverman Soap Box 137 2012-01-07 07:52 science_man_88 Puzzles 9 2010-07-30 21:22 geoff Programming 7 2006-08-18 12:16 GP2 Hardware 11 2004-01-21 03:01

All times are UTC. The time now is 04:39.

Thu Apr 15 04:39:54 UTC 2021 up 6 days, 23:20, 0 users, load averages: 1.53, 1.70, 1.76