How many ways can you code an LL test

science_man_88 · 2015-12-31, 13:05

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 ?

R. Gerbicz · 2015-12-31, 13:52

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)}

science_man_88 · 2015-12-31, 14:22

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.

danaj · 2015-12-31, 19:42

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.

science_man_88 · 2015-12-31, 20:31

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.

LaurV · 2016-01-01, 03:34

Quote:

Originally Posted by R. Gerbicz

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

science_man_88 · 2016-01-01, 03:45

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.

axn · 2016-01-01, 04:11

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)

xilman · 2016-01-01, 09:22

Quote:

Originally Posted by science_man_88

etc. how many ways can you code it ?

Countably infinite.

R. Gerbicz · 2016-01-01, 09:49

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.

science_man_88 · 2016-01-01, 12:55

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.

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)}

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

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