Asked by Jesse Crotts
on 18 Sep 2016

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

Answer by John D'Errico
on 8 Dec 2018 at 23:23

Edited by John D'Errico
on 9 Dec 2018 at 0:03

Accepted Answer

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.

Sign in to comment.

Answer by John D'Errico
on 18 Sep 2016

Edited by John D'Errico
on 18 Sep 2016

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.

John D'Errico
on 19 Sep 2016

Please don't add answers every time you want to make a comment. Answer (from Jesse) moved to a comment:

"Thank You for posting how to find large twin primes. However my conjecture involves knowing all twin primes up to a particular point. Thanks for the idea though."

John D'Errico
on 19 Sep 2016

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.

Jesse Crotts
on 8 Dec 2018 at 21:29

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Jesse Crotts (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/303515-i-wrote-a-code-to-produce-twin-primes-however-it-is-pretty-strenuous-for-my-computer-and-i-need-som#comment_391879

Sign in to comment.