mersenneforum.org  

Go Back   mersenneforum.org > Fun Stuff > Lounge

Reply
 
Thread Tools
Old 2015-12-31, 13:05   #1
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

203008 Posts
Default 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;
}
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
}
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:
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
science_man_88 is offline   Reply With Quote
Old 2015-12-31, 13:52   #2
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

31·47 Posts
Default

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)}
R. Gerbicz is offline   Reply With Quote
Old 2015-12-31, 14:22   #3
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

26·131 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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.
science_man_88 is offline   Reply With Quote
Old 2015-12-31, 19:42   #4
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

16148 Posts
Default

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.
danaj is offline   Reply With Quote
Old 2015-12-31, 20:31   #5
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

26×131 Posts
Default

Quote:
Originally Posted by danaj View Post
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
science_man_88 is offline   Reply With Quote
Old 2016-01-01, 03:34   #6
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

33×347 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
A short and faster example, still in PARI-GP:
Code:
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)}
Fixed that for you, for few microseconds gain
LaurV is offline   Reply With Quote
Old 2016-01-01, 03:45   #7
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

838410 Posts
Default

Quote:
Originally Posted by LaurV View Post
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
science_man_88 is offline   Reply With Quote
Old 2016-01-01, 04:11   #8
axn
 
axn's Avatar
 
Jun 2003

132F16 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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 View Post
Fixed that for you, for few microseconds gain
Another micro optimization is to get rid of
Code:
if(x<0,x+=mp)
axn is online now   Reply With Quote
Old 2016-01-01, 09:22   #9
xilman
Bamboozled!
 
xilman's Avatar
 
"π’‰Ίπ’ŒŒπ’‡·π’†·π’€­"
May 2003
Down not across

22×3×887 Posts
Default

Quote:
Originally Posted by science_man_88 View Post
etc. how many ways can you code it ?
Countably infinite.
xilman is offline   Reply With Quote
Old 2016-01-01, 09:49   #10
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

26618 Posts
Default

Quote:
Originally Posted by LaurV View Post
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 View Post
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
R. Gerbicz is offline   Reply With Quote
Old 2016-01-01, 12:55   #11
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

26×131 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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
science_man_88 is offline   Reply With Quote
Reply

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 2012-05-14 08:51
RDS's unique pedagogic ways R.D. Silverman Soap Box 137 2012-01-07 07:52
ways to get rid of oil spills science_man_88 Puzzles 9 2010-07-30 21:22
Help needed to test Athlon64 code geoff Programming 7 2006-08-18 12:16
Creative ways to achieve Athlon 64 / Opteron optimization 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

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.