I wrote a code to produce twin primes, however it is pretty strenuous for my computer and I need some high numbers to check a conjecture. Is there a better way to get twin primes. Note: these twin primes are represented by 6m's where 6m+-1 are prime.
2 次查看(过去 30 天)
显示 更早的评论
n = 1:120000;
m61 = n*6-1;
Twins = [];
for i=1:numel(m61);
z = isprime(m61(i));
y = isprime(m61(i)+2);
if z == 1 && y == 1;
Twins(:,(numel(Twins)+1)) = m61(i)+1;
else
end
end
采纳的回答
John D'Errico
2018-12-8
编辑:John D'Errico
2018-12-9
Ok, it looks like you are still interested in this, so here is an approach targetted more at what I hear from your last comment.
Simple would be, if you just want to find all twin primes less than some value, just use primes to get all primes, then find all primes that are exactly 2 larger than the one before.
N = 1e7;
P = primes(N);
ind = find(diff(P) == 2);
TP = [P(ind)',P(ind+1)'];
size(TP)
ans =
58980 2
TP(1:10,:)
ans =
3 5
5 7
11 13
17 19
29 31
41 43
59 61
71 73
101 103
107 109
But there is no formula to find all such twin primes, nor to predict when a twin prime pair will occur, except that it must be a pair of the form 6*m+[-1 1].
Yes, you could use isprime to test each pair of numbers of that form, but that would be extremely inefficient, far less good than just calling primes directly.
So, is there something slightly better? For example, could we write a variation of a sieve, targetted specifically at twin primes? For example, suppose I wanted to find all twin primes less than 1e9? primes(1e9) is probably too big for me to do on my computer. But I know that twin primes must be of the form 6*m +/-1.
tic
% find all twin prime pairs where both are less than N
N = 1e9;
Plist = primes(sqrt(N));
% Keep the small twin primes from Plist
twin_primes = Plist(diff(Plist) == 2);
twin_primes = [twin_primes;twin_primes+2]';
% We don't need to check 6*m-1 for divisibility by 2 or 3.
Plist(1:2) = [];
% m will correspond to the index of a potential twin prime, 6*m-1.
% At the end of this comutation, if m(k) is true, then
% 6*m(k) +/- 1 will be a twin prime. Note that all of this will miss the
% very first twin prime pair, [3 5]. But all other twin prime pairs are of
% the form 6*m +/- 1.
m = true(1,floor((N-2)/6));
for p_i = Plist
% which elements of m correspond to 6*m-1 that are divisible by p_i?
% Find them using the modular inverse of 6, modulo p_i.
p_i_inv = minv(6,p_i);
m1 = p_i_inv;
if p_i == (6*m1 - 1)
m1 = m1 + p_i;
end
m(m1:p_i:end) = false;
m2 = mod((p_i - 1)*p_i_inv,p_i);
if p_i == (6*m2 - 1)
m2 = m2 + p_i;
end
m(m2:p_i:end) = false;
end
m = find(m);
twin_primes =[twin_primes;[6*m-1;6*m+1]'];
toc
So reasonably fast.
Elapsed time is 9.773307 seconds.
This ran pretty quickly on my computer, since the only RAM it required was a vector of size N/6 bytes. So for N=1e9, m only required 0.17 gigabytes. And the loop was also pretty efficient.
Note that the above code uses my minv tool, which is found in my vpi toolbox, although when used on doubles, it returns double precision output.
whos m twin_primes
Name Size Bytes Class Attributes
m 1x3424019 27392152 double
twin_primes 3424506x2 54792096 double
So it found every one of the 3.4 million twin prime pairs for which neither of them exceed 1e9. But since it did not need to explicitly test any of them explicitly for primality, nor did it need to generate the entire list of primes below 1e9, the computation was actually pretty fast. While my computer is moderately old, I'd guess I could easily push N out further out than 1e10. (I just tried it for N=1e10, and it took 110 seconds to find all 27412679 twin prime pairs less than 1e10. Not bad for a 10 year old computer.)
I'm not sure if this will be of any real value to you, but it will indeed generate all such twin primes less than N.
0 个评论
更多回答(1 个)
John D'Errico
2016-9-18
编辑:John D'Errico
2016-9-18
I recall writing a toy that worked as a partial Euclidean prime sieve to find twin primes. The point being that you can then search for relatively large primes. The domain I was playing in was searching for twin primes with thousands of decimal digits, far too large for the tools in MATLAB. I was using my VPIJ tool to test it out as a VPI replacement, but the ideas are applicable for smaller numbers too. I can dig up the code I was using if you are interested.
Edit, let me see if I can remember the basic ideas in that exploration. For example, here is a number with roughly 100 decimal digits.
N = prod(vpij(primes(250)))
N =
256041159035492609053110100510385311995538591998443060216114576417920917800321526504084465112487730
log10(N)
ans =
98.4083097844727
The point is that N is the product of LOTS of small primes.
Now, lets get a list of twin primes that are between 250 and 10000.
plist = primes(10000)';
plist(plist < 250) = [];
k = find(diff(plist) == 2);
tp = plist(k);
tp'
ans =
Columns 1 through 12
269 281 311 347 419 431 461 521 569 599 617 641
Columns 13 through 24
659 809 821 827 857 881 1019 1031 1049 1061 1091 1151
Columns 25 through 36
1229 1277 1289 1301 1319 1427 1451 1481 1487 1607 1619 1667
Columns 37 through 48
1697 1721 1787 1871 1877 1931 1949 1997 2027 2081 2087 2111
Columns 49 through 60
2129 2141 2237 2267 2309 2339 2381 2549 2591 2657 2687 2711
Columns 61 through 72
2729 2789 2801 2969 2999 3119 3167 3251 3257 3299 3329 3359
Columns 73 through 84
3371 3389 3461 3467 3527 3539 3557 3581 3671 3767 3821 3851
Columns 85 through 96
3917 3929 4001 4019 4049 4091 4127 4157 4217 4229 4241 4259
Columns 97 through 108
4271 4337 4421 4481 4517 4547 4637 4649 4721 4787 4799 4931
Columns 109 through 120
4967 5009 5021 5099 5231 5279 5417 5441 5477 5501 5519 5639
Columns 121 through 132
5651 5657 5741 5849 5867 5879 6089 6131 6197 6269 6299 6359
Columns 133 through 144
6449 6551 6569 6659 6689 6701 6761 6779 6791 6827 6869 6947
Columns 145 through 156
6959 7127 7211 7307 7331 7349 7457 7487 7547 7559 7589 7757
Columns 157 through 168
7877 7949 8009 8087 8219 8231 8291 8387 8429 8537 8597 8627
Columns 169 through 180
8819 8837 8861 8969 8999 9011 9041 9239 9281 9341 9419 9431
Columns 181 through 188
9437 9461 9629 9677 9719 9767 9857 9929
The vector tp contains the first member of each pair of twin primes in that interval. I can convince you they are all twin primes with this pair or tests:
all(isprime(tp))
ans =
1
all(isprime(tp+2))
ans =
1
So, now consider the pair of numbers
[N+tp(i), N+tp(i)+2]
For every element of tp. Are those numbers both prime? The idea is both of those numbers are not divisible by any small prime less than 250, since we know that N is divisible by EVERY prime less than 250. And since tp(i) is prime, as is tp(i)+2, we know the sum cannot have any small prime factors, nor is the sum divisible by the twin prime offset.
Using the above approach, it took only a moderate amount of time to generate the twin prime pair below, which are primes with roughly 1000 decimal digits.
N = prod(vpij(primes(2360)));
isprime(N+[2414411, 2414413])
ans =
1 1
N
N =
2372745398987024026246661730569400071940852598530913706750961313093613019661136660043304214953229588549115824235879381842077106653008440157763482515826755896463281987961979543338530673344287018463219016115004917088496306089189981437869746067778819705086659226608214905773265406521005821190696164966769062789270607392466810411611364927251425724094853199855468930458472451779485739850642137721909433918071274371636467786383202340657583081021199102890278550492075053970663402100153191893555813143629199328494500822016926378226575514315855644658990142857007620141512167568953260848559217323755974661494616621625481469726141488042628605904007828553835113010863602272246388549781675554738720289081288962182158908888074895912630253490848795442738869861959492164311090600752598942985608896520433889214529842303378943565378331077082821046504392936708043242708752032968977683539225351550077540671477173072122185514120758317186048432316935013230715793661245186916334985219768354365429487437645709510479749387570
Again, the computations were done using my VPIJ tool. VPI is on the file exchange, and VPIJ is fully working, though still in beta as far as I am concerned.
Or, you can use a similar scheme to find smaller twin primes, but still pretty large.
N = prod(primes(41))
N =
304250263527210
log2(N)
ans =
48.1122518409569
So N is moderately large.
plist = primes(1000)';
plist(plist < 42) = [];
k = find(diff(plist) == 2);
tp = plist(k);
Can we find any twin primes here?
find(all(isprime(N + [tp,tp+2]),2))
ans =
12
YES! What are they?
N + [tp(12),tp(12)+2]
ans =
304250263527479 304250263527481
isprime(ans)
ans =
1 1
So the pair shown above is indeed a 15 digit twin prime pair, and far too large to have been efficiently found using a brute force search through lists of primes. (You don't want to generate primes(1e15) for a search like that.) And the nice thing is the entire search took far more time to type in the code than it did to find the twin prime pair itself. Again, the basic idea here was to employ an implicit sieve that essentially finds blocks of numbers that are known not to have any small prime factors. That increases the odds that we might find a twin prime pair in the vicinity.
4 个评论
John D'Errico
2016-9-19
编辑:John D'Errico
2016-9-19
You have made no conjecture. Is your conjecture that all twin prime pairs beyond (3,5) are of the form 6*m +/- 1?
If so, that is trivially provable. No large twin primes need be generated. In fact, you cannot prove a conjecture by finding cases where it holds true.
To prove a twin prime pair is of the form 6m+/-1 is easy.
NO members of a twin primes pair can be of the form of 6*m, 6*m+2, 6*m+4. That seems trivial, because they are all divisible by 2, and 2 is not a twin prime. (Twin primes are a pair of primes separated by 2.)
Therefore all twin prime members must be of the form 6*m+1, 6*m+3, or 6*m+5. Note the latter case can be written in the form 6*m-1, for some m.
Since 6*m+3 is always divisible by 3 (and is composite) unless m=0, then no twin prime pair besides (3,5) has a member of the form 6*m+3.
Therefore all twin prime pairs beyond (3,5) can be written in the form (6*m-1,6*m+1), for some m.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Discrete Math 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!