![]() |
|
|
#1 |
|
Sep 2004
10000101012 Posts |
I was intrigued by http://primes.utm.edu/glossary/page....ortunateNumber and http://www-lipn.univ-paris13.fr/~ban...factorial.html and I thought I would try a little fortunate computing myself. I made this pari script myself and checked up to k=300. What do you guys think? I thought I might see if I could do it with gmp in 64-bit which should be way faster.
{ p=psum=2; q=nextprime(p+2); print(q-p); for (k=2, 3000, p=nextprime(p+1); psum=psum*p; q=nextprime(psum+2); print(k); if(!isprime(q-psum), print(q-psum); ) ) } |
|
|
|
|
|
#2 |
|
Sep 2004
10258 Posts |
Running overnight it got to 550. I bet if I increased its memory above 4 megs it would go significantly faster...
To summarize for the people that don't want to read the links, the fortunate series is conjectured to be all prime numbers, and it has been checked to about n=2000, but someone said there might have been a few errors in that result. |
|
|
|
|
|
#3 |
|
Sep 2004
10000101012 Posts |
I made some optimizations and it seems to be around ~10% faster. Didn't really do official benchmarks.
allocatemem allocatemem allocatemem allocatemem { arr=[2, 3, 5, 7, 11, 13, ..., 9973]; psum=1; for (k=1, 3000, p=arr[k]; psum=psum*p; q=nextprime(psum+3); print(psum+3); f=q-psum; print(k " " f); if(!isprime(f), print(k); k=3000;) ); } |
|
|
|
|
|
#4 |
|
Oct 2007
Manchester, UK
23·59 Posts |
How about just using a forprime(2,10000,...) loop? You can work out in advance which prime you want to stop at, or you could include a counter in the loop and break() it when it reaches a specified point (though this would slow the code down more).
Also, I find that using only 1 letter variables, and in general as fewer characters as possible, helps to speed up the code in Pari/GP. |
|
|
|
|
|
#5 | ||
|
Sep 2004
13×41 Posts |
Quote:
Quote:
It seems most of the time (99%) is spent on the nextprime() method. The primes are ~1400 digits around n=400. Would anyone care to help me get started with GMP 64-bit? I believe that would significantly increase speed. I downloaded the source. I have visual studio 2008 and vista x64. Last fiddled with by Joshua2 on 2009-01-27 at 04:19 |
||
|
|
|
|
|
#6 |
|
Oct 2007
Manchester, UK
23×59 Posts |
Of the two alternatives I gave, only the second would introduce any slow down. Though it would only be very slight, and certainly nowhere near as much as is introduced by printing out a couple of lines worth of progress every iteration.
However, since the loop only gets through very few iterations, you're only losing milliseconds for the print()s, and microseconds for the extra variable being iterated. A forprime() loop would eliminate the need for the prime array entirely. It is a built in loop in Pari: http://pari.math.u-bordeaux.fr/docht....html#forprime Here is some sample code that will run through all primes below 1000 (168 iterations): Code:
P=1;
k=0;
forprime(a=2,1000,
k++;
P*=a;
q=nextprime(P+2);
if(isprime(q-P)==0,
print(k);
break()
)
);
print("Done ",k," iterations.");
Code:
P=1;
k=0;
forprime(a=2,1000,
k++;
P*=a;
q=nextprime(P+2);
while(isprime(q)==0,
q=nextprime(q+2)
);
if(isprime(q-P)==0,
print(k);
break()
)
);
print("Done ",k," iterations.");
Code:
l=1000; \\ The number of iterations to run.
P=1;
k=0;
forprime(a=2,500000,
k++;
P*=a;
q=nextprime(P+2);
while(isprime(q)==0,
q=nextprime(q+2)
);
r=q-P;
if(isprime(r)==0,
print("F_k NOT PRIME\nk =\t",k,"\np =\t",a,"\nF_k =\t",r,"\nF_k NOT PRIME\n");
break()
,
print("k =\t",k,"\np =\t",a,"\nF_k =\t",r,"\n")
);
if(k>=l,
break()
)
);
print("Done ",k," iterations.");
|
|
|
|
|
|
#7 |
|
Sep 2004
10000101012 Posts |
Cool. I read that this was verified up to n=2000 with mathematica, and n=1300 with maple. Those programs should have been even slower than pari or at least not faster. Maybe they just left it running for a month or two. I was hoping to raise the bounds of how far it had been proved. It seems 64-bit would be my best speedup.
I am also considering various other conjectures. Instead of starting with P=1 as in yours, start with 2, 3, or whatever. I have tried up to 40 for the first 100 or so iterations and they are all prime as well. This isn't obvious is it? Edit: I discovered a couple P starting values (both odd and even) that result in the first number of the series being 9, after that the series appear to be all primes. Edit2: P=221 2nd number is 9, then fine. Up to P=2000 ignoring cases where the first or 2nd is 9. Edit3: P=2018 3rd is 9. Maybe these things can't produce anything but prime or 9? Randomly testing way out 150064 4th is 9. Last fiddled with by Joshua2 on 2009-01-27 at 05:24 |
|
|
|
|
|
#8 |
|
Oct 2007
Manchester, UK
135710 Posts |
The code I posted actually starts with P=2, it's just that I initialise it to 1, then on the first iteration the value is set to P=P*2, so P=2 for the first iteration. The next it is set to P=P*3, so P=6, and so on.
If you want to start with a different value, the simply change the first number in the forprime loop, eg: Code:
forprime(a=7,500000, As far as compiling the code, GP2C would probably be worth a look. The site claims that code runs 3 to 4 times faster that way. http://pari.math.u-bordeaux.fr/pub/p...gp2c/gp2c.html I don't think that GMP contains any deterministic primality tests, so if you were to go the GMP route, you would have to either write your own implementation, or else find an open source one from somewhere (such as Pari for instance ). However, I could be mistaken about GMP there.467 is the first iteration to cross into 1400+ digits by the way. |
|
|
|
|
|
#9 |
|
Sep 2004
13×41 Posts |
That code you gave seems a lot slower. 15s to get to k=75, whereas mine gets to k=100 in 4.5s. It appears to be all that isprime to check the nextprime.
Unfortunately, gp2c seems to require linux, at least that tutorial is for linux. I can't just type make or do links with ln on vista :). What I'm talking about with the P thing, called i currently on my version, is a multiplier on the primorial. Instead of 1*2*3*5*7*11... have 2*2*3*5*7 or P*2*3*5*7. Where P is 100 or 100000 or whatever. I tried having P being restricted to prime numbers, but I am still coming up with the first term being 9 or 2nd, 3rd, extra as P gets large. After it gets the 9 out of its system it seems to be fine. I might call something like this the Jordan generalized Fortune conjecture :) Let P be the product of the first n primes. If q is the smallest prime greater than kP+1, then q-kP is prime or 9 (or instead prime after n is sufficiently large). Last fiddled with by Joshua2 on 2009-01-27 at 06:28 |
|
|
|
|
|
#10 | |
|
Sep 2004
13×41 Posts |
Quote:
"However, all of them are probabilistic tests (Mathematica uses multiple Rabin-Miller and Lucas, Pari GP uses Baillie-Pomerance-Selfridge-Wagstaff and Lucas) and the best known deterministic test is about 1000 times slower." I'm not sure if the slowdown is worth it or not. I suppose I'll do some of each. I attached the pari of the your program modified to test the generalized conjecture. Last fiddled with by Joshua2 on 2009-01-27 at 06:50 |
|
|
|
|
|
|
#11 |
|
(loop (#_fork))
Feb 2006
Cambridge, England
23×11×73 Posts |
I think you have very much underestimated the sloth of guaranteed primality tests for large general numbers. The product of the first 2000 primes is about 10^7483; checking such a number with a competent parallel implementation of ECPP (NB I don't know if such a thing exists) would take about a week, using gp would take to a fair approximation forever.
for(c=10,20,print(nextprime(c*10^A)-c*10^A)) (eg eleven sets of {some probably-primality tests + one true-primality}) takes on 2.4GHz core2 Code:
A= 999 222 seconds A=1499 880 seconds A=1999 2290 seconds A=2499 5404 seconds A=2999 8870 seconds |
|
|
|
![]() |
| Thread Tools | |
Similar Threads
|
||||
| Thread | Thread Starter | Forum | Replies | Last Post |
| Let's attack the Bayesian-ECM-bounds again | fivemack | Math | 34 | 2021-08-05 10:55 |
| Nuke attack - run or hide? | MooMoo2 | Soap Box | 40 | 2018-01-19 23:48 |
| TF: A job half done? | davieddy | Lounge | 35 | 2010-10-01 20:18 |
| Attack of the Cosmic Rays | S485122 | Hardware | 3 | 2010-08-24 01:19 |
| Attack of the Killer Zombies | ewmayer | Lounge | 12 | 2007-01-30 05:56 |