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.

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

[Jesse Crotts]

Nevermind, I found a much faster way.

if true
  % code
end
P = primes(10000000);
P2 = P+2;
Twins = intersect(P,P2)-1;
Twins(1) = [];

Answer 1

[John D'Errico]

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.

Answer 2

[John D'Errico]

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]

A fun thing about primes is they are easy to play with, to find primes with interesting properties. You can learn a bit about number theory along the way.

My personal favorite prime is:

2 + 1723*17^1717

But then I'm a fan of the number 17. My other favorite prime is the SECOND prime number in the family:

(n+1)*17^n-1

It is sort of a modified Cullen prime .

https://en.wikipedia.org/wiki/Cullen_number

In fact, I don't know what that number is, since I've never looked past 23000 or so. For n = 23000, one is working with numbers that have roughly 28000 decimal digits. (Quiz: what value of n generates the FIRST prime number in that family? It is still moderately large.)

[John D'Errico]

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]

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]

I believe that the conjecture I was playing with was something like: you can get every even number by adding two non twin primes or adding 2 twin primes. Its been a long time since I did that stuff though. I can't really remember exactly what my conjecture was but I think it was a stronger case of the goldbach conjecture.