mersenneforum.org  

Go Back   mersenneforum.org > Math Stuff > Computer Science & Computational Number Theory

Reply
 
Thread Tools
Old 2011-04-21, 11:36   #23
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

22·5·373 Posts
Default

Quote:
Originally Posted by CRGreathouse View Post
Of course it matters for this thread, since that's the whole topic. But I would say that it matters to me personally: I frequently find myself in the position of needing to test the primality of many (trillions) of small numbers, usually below 64 bits. At the moment I use PARI which takes about 3 microseconds. (Well, if you abuse it enough you can get to that speed anyway; the current version doesn't yet know about Jan's result and so if you ask it to prove primality it takes more like a millisecond.)

At those speeds even just a billion numbers take an hour. Reducing that to 15 or 20 minutes would be great.
Agreed.
R.D. Silverman is offline   Reply With Quote
Old 2011-04-22, 11:01   #24
Xitami
 
Apr 2010

168 Posts
Default

:) PARI/GP can write programs in C
Code:
ifelse(n,p,k)={ my(c);
	if(p<k,
		if(n>0,print);
	 	for(i=0,n,print1(" "));
		c=(p+k)>>1;
		print1("if(a1262["c"]<n)");  
		ifelse(n+1, c+1, k);
		for(i=0,n,print1(" "));
		print1("else");
		ifelse(n+1, p, c)
	,
		print(" return a1262["p"]==n;"); 
	);
}
fun(k)={
	print("int test(unsigned long long n) {");
	lt=eq=0;
	ifelse(0,1,k);
	print("} ");
}
fun(6)
Output (table is not necesary)
Code:
int test(unsigned long long n) {
 if(a1262[3]<n)
  if(a1262[5]<n) return a1262[6]==n;
  else
   if(a1262[4]<n) return a1262[5]==n;
   else return a1262[4]==n;
 else
  if(a1262[2]<n) return a1262[3]==n;
  else
   if(a1262[1]<n) return a1262[2]==n;
   else return a1262[1]==n;
}
only 2314 sprp to base 2 less then 2^32
only log2(2314)+1=12 tests in 2314*2 rows of programm

310 base 2 and 101 pseudoprimes below 145046817709

Last fiddled with by Xitami on 2011-04-22 at 11:16
Xitami is offline   Reply With Quote
Old 2011-04-22, 14:22   #25
Xitami
 
Apr 2010

168 Posts
Default

Quote:
Originally Posted by ATH View Post
[...]But I just took the list of the 31,894,014 SPRP base2 < 2^64 [...]
1. where i can find this list?
2. how many pseudoprimes survive test to two bases ?
Xitami is offline   Reply With Quote
Old 2011-04-22, 18:27   #26
ATH
Einyen
 
ATH's Avatar
 
Dec 2003
Denmark

2×5×172 Posts
Default

The full list is here:
http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html
but this include all 118,968,378 Fermat base2 PRP + Miller Rabin SPRP base2.

The best constant second base I found so far is 9375 where 1,147,253 (0.9643%) survives.

More exotic choices does a bit better. For number n if you use base: floor(n/5)+49 then 1,092,753 survives. Base: floor(Sqrt(n))+1: 1,023,767 survives.
ATH is offline   Reply With Quote
Old 2011-04-26, 13:14   #27
Xitami
 
Apr 2010

2×7 Posts
Default

Quote:
Originally Posted by zerothbase View Post
[...]This is definitely the weakest part of the entire code.[...]
Code:
Code:
typedef unsigned long long int u64;

inline u64 mulmod(u64 i, u64 j, u64 k){
        u64 p, r = 0;
        while(j > 0){
                if(j&1) //  r = (i + r) % k
                        if(0==r) r=i;
                        else {  r=k-r;
                                if(i>=r) r=i-r;
                                else r=k-r+i;}
                // i = (i + i) % k
                if(i>k-i) i=(i-k)+i;
                else i+=i; 
                j>>=1;}
        return r;}
no MOD, no overflow
Xitami is offline   Reply With Quote
Old 2011-04-27, 18:33   #28
Jeff Gilchrist
 
Jeff Gilchrist's Avatar
 
Jun 2003
Ottawa, Canada

7×167 Posts
Default

Quote:
Originally Posted by ATH View Post
This has already been done: http://www.trnicely.net/misc/bpsw.html
"Furthermore, preliminary analysis by Gilchrist of Feitsma's database, further extended, indicates that no Baillie-PSW pseudoprime (standard or strong) exists below 2^64 (24 October 2009; double checking of this result is in progress), approximately 1.8446744e19."

But I just took the list of the 31,894,014 SPRP base2 < 2^64 and did a Lucas-Selfridge test on them, it took around 15min, and none of them are lucas PRP. So assuming doublechecking proves that the SPRP list is complete there is BPSW pseudoprime below 2^64.
Jan Feitsma has posted the results as "preliminary" because nobody has verified the math of his algorithm to ensure that all the base 2 psp up to 2^64 have been enumerated. David Cleaver (WraithX on these boards) and I verified the pseudoprime database with Jan's algorithm using independently written code so we have double checked that his algorithm generated the same set of psp using different code. But Jan has not proven that his algorithm in fact enumerates 100% of the psp up to 2^64. Jan was originally in contact with Wil Galway (the one hosting the 2^64 database now) to review his math but for whatever reason his algorithm was never evaluated.

Back in 2009 when David and I finished double checking the psp database, I was able to confirm there were no BPSW psp up to 2^64 with his database (assuming we did in fact enumerate all of them) and posted that here: http://gilchrist.ca/jeff/factoring/pseudoprimes.html

Charles Greathouse has also confirmed no BPSW psp exist using the BPSW implementation in PARI with Jan Feitsma's database so two different implementations of BPSW have also been tried.

Last fiddled with by Jeff Gilchrist on 2011-04-27 at 18:39
Jeff Gilchrist is offline   Reply With Quote
Old 2011-04-28, 04:03   #29
CRGreathouse
 
CRGreathouse's Avatar
 
Aug 2006

22·32·163 Posts
Default

Quote:
Originally Posted by Jeff Gilchrist View Post
Charles Greathouse has also confirmed no BPSW psp exist using the BPSW implementation in PARI with Jan Feitsma's database so two different implementations of BPSW have also been tried.
I know that many people have checked this, but I don't know which implementations of BPSW were used (though as you point out, it's at least two).
CRGreathouse is online now   Reply With Quote
Old 2011-05-24, 16:57   #30
JohnFullspeed
 
May 2011
France

7×23 Posts
Default 2^32 primes

Why are you not using an array with all the primes smaller tahnn2^32?
You have only 200 000 000 of primes so using a simple RLE the file
is 200 000 000 of byte

Two hours to compute the llst (i can give the file) and you can factorized
20 digits number in less than one minute.

For the classic m^k mod you can use a fast algo

Code:
function E_M(m, K, N: int64): int64;
var b, a: int64;
begin
 b := 1;
 a := m;
 if odd(k) then
  b := A;
 k := k shr 1;
 repeat
  A := (A * A) mod N;
  if odd(k) then
   b := (A * b) mod N;
  k := k shr 1;
 until k = 0;
 result := b;
end;
It's also good for RSA encryption

John

Last fiddled with by JohnFullspeed on 2011-05-24 at 17:00
JohnFullspeed is offline   Reply With Quote
Old 2011-05-24, 18:38   #31
CRGreathouse
 
CRGreathouse's Avatar
 
Aug 2006

586810 Posts
Default

Quote:
Originally Posted by JohnFullspeed View Post
Two hours to compute the llst (i can give the file)
Two hours?!? It takes PARI 11 seconds on this machine. (Technically that's "just" to 4.26e9 since this computer is running 32-bit gp.)

Quote:
Originally Posted by JohnFullspeed View Post
and you can factorized 20 digits number in less than one minute.
Sure, as long as it's not a semiprime where both factors are between 2^32 and 10^10. But you can do much better -- using rho and SQUFOF, PARI can factor worst-case numbers of that size in about 200 milliseconds, and ten times faster for typical numbers.

But primality testing is much faster than factoring! PARI only takes about 2 milliseconds to prove that a number is prime; typical composites are extremely easy to detect, taking only 20 *microseconds* by my calculations (I had to test 100,000 to get good timing!).

So trial division works at that size, but we're looking for something really fast here.
CRGreathouse is online now   Reply With Quote
Old 2011-05-24, 20:51   #32
Xitami
 
Apr 2010

2×7 Posts
Default

A := (A * A) mod N;
John, how long is A*A?
Xitami is offline   Reply With Quote
Old 2011-05-25, 07:32   #33
JohnFullspeed
 
May 2011
France

7×23 Posts
Default a:=a*a mod N

10000 bits if you want....

In fact you can use this mrethod for othezr problems as divisioins..

The method is:
Let A and B be big positive integerss
Let ? an sllower opoeraareationn
A ? B = C

you work on the binary representation
each bit of A
If te bit is 0 you make a specific code
if the bit is 1 you make and other code
you pass to, the Bit after
unil A have bit


Code:
function E_M(m, K, N: int64): int64; var b, a: int64; begin  b := 1;  a := m;  if odd(k) then   b := A;  k := k shr 1; <= first  bit of tha operande  repeat   A := (A * A) mod N;  <= in all casae Bit=0 or Bit=1   if odd(k) then       <= the bit is 1    b := (A * b) mod N;   k := k shr 1;        <+ nexT bit  until k = 0;  <= no new bit  (end)  result := b; end;
if you remove the mod you have a speed exopanciation
Is't easy to do a faster div....

A := (A * A); <= in all casae Bit=0 or Bit=1 if odd(k) then <= the bit is 1 b := (A * b); k := k shr 1; <+ nexT bit until k = 0; <= no new bit (end)


You understood that A can have 1000 bios you make 'only'
1000+5000 divmod

John
JohnFullspeed is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Do normal adults give themselves an allowance? (...to fast or not to fast - there is no question!) jasong jasong 35 2016-12-11 00:57
How to Check if Non-Mersenne Number isPrime? FloatingPoint Operation Billion Digits 39 2015-10-21 02:15
How fast is the dog? Andi47 Puzzles 20 2009-04-01 02:35
I wonder how fast this is... ixfd64 Hardware 1 2005-11-21 21:28
Fast way to square??? maheshexp Math 2 2004-05-29 01:54

All times are UTC. The time now is 07:52.

Tue Aug 11 07:52:22 UTC 2020 up 25 days, 3:39, 1 user, load averages: 1.44, 1.59, 1.70

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, 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.