Software for finding prime constellations

[ATH]

My programs tend to be very messy, so I'll just explain in pseudocode how I trial factor in a test program I made.

For example this 19-tuple pattern:
p=[0 4 6 10 16 18 24 28 30 34 40 46 48 54 58 60 66 70 76] where the first p is: p=217153 (modulo 510510).

I choose a starting point and find the nearest n=217153+76 (modulo 510510), because I use the last prime in the 19-tuple pattern, and then I want to check n+k*510510 for k>=0.

So I create a bitset array of k-values to sieve with a size around 1M-2M which seemed most efficient. If a sieving prime divides any one of the 19 prime pattern positions I can exclude this k-value from the array as well as every p'th k-value from this first value: k +c*p for c>=0.

Since 510510 = 17# I starting sieving at p=19.

for p=19 to sievelimit; p=nextprime
{
  a:=n%p
  b:=510510%p

  for k=0 to k<p					I check (n + k*510510)%p for 0 <= k <= p-1
    if (a<=76)						Since my n is last prime position, if a = n%p <=76 then p divides inside the pattern somewhere.
      for i=(76-a)%p to i<=76 and i=i+p			As long as p<=76 there can be several values to check inside the 19-tuple pattern. (a is the last prime position mod p)
	  if i is equal to a 19-tuple pattern position	then I know current k value and all k+c*p for c>=0 can be reset in the array, because p divides this pattern position in each of them.
	    for j=k to j<=bitsetsize and j=j+p
              reset bitset[j]
	    break					No need to check rest of i values once I found one.
    a=(a+b)%p
}

Since I trial factor all 19 different positions in the pattern, the k-values are quickly eliminated, and I only need to trial factor up to like around p=59-83 until there are very few valid k-values left to PRP.
Then I PRP the rest of the k-candidates starting from the front or back of the 19-tuple pattern for each k-value.
As soon as a composite is found, the current k-value can be skipped, so most k-values requires only 1 or a few PRP tests.

[Puzzle-Peter]

Robert Gerbizc's Polysieve was not originally written for this kind of stuff, but it can handle these problems using the right input and it can be customized to do it more easily.
Not sure about performance at k=19 compared with the code the record holders used but it was good enough for a record at k=12.

[Cybertronic]

This unknown prime 19-tuplet is very hard.
I found with my 2 PCs this smallest prime 17-tuplet

1341829940444122313597407 + d, d = 0, 4, 10, 12, 16, 22, 24, 30, 36, 40, 42, 46, 52, 54, 60, 64, 66
after 3 weeks or so. The offset is 3.4e23.

I believe, to reach this offset for 19-tuplet is 10x faster...~2 days, or 1.5e24 per month (or an offset with 6e17+ per sec).

My estimation for the 1st one is 5*10^26 or 330 months (27 years).

100 PCs should be found it in 180 days

[User140242]

With this code fixed the number of elements and the delta max finds evential k-tuples less than n.

To find any k-tuples with elements in two consecutive search blocks it is necessary to increase the wheel from 210 to 2310 in the input values.

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <stdint.h>
#include <time.h> 

const int64_t PrimesBase[5]={2,3,5,7,11};
const int64_t n_PB_max = 5;

const int64_t del_bit[8] =
{
  ~(1 << 0),~(1 << 1),~(1 << 2),~(1 << 3),
  ~(1 << 4),~(1 << 5),~(1 << 6),~(1 << 7)
};

const int64_t bit_count[256] =
{
  0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
};

int64_t Euclidean_Diophantine( int64_t coeff_a, int64_t  coeff_b)
{
    // return y in  Diophantine equation  coeff_a x + coeff_b y  = 1
    int64_t k=1;
    std::vector<int64_t> div_t; 
    std::vector<int64_t> rem_t;
    std::vector<int64_t> coeff_t;
    div_t.push_back(coeff_a);
    rem_t.push_back(coeff_b);
    coeff_t.push_back((int64_t)0);
    div_t.push_back((int64_t)div_t[0]/rem_t[0]);
    rem_t.push_back((int64_t)div_t[0]%rem_t[0]);
    coeff_t.push_back((int64_t)0);
    while (rem_t[k]>1)
    {
        k=k+1;
        div_t.push_back((int64_t)rem_t[k-2]/rem_t[k-1]);
        rem_t.push_back((int64_t)rem_t[k-2]%rem_t[k-1]);
        coeff_t.push_back((int64_t)0);
    }
    k=k-1;
    coeff_t[k]=-div_t[k+1];
    if (k>0)
        coeff_t[k-1]=(int64_t)1;
    while (k > 1)
    {
        k=k-1;
        coeff_t[k-1]=coeff_t[k+1];
        coeff_t[k]+=(int64_t)(coeff_t[k+1]*(-div_t[k+1]));
    }
    if (k==1)
        return (int64_t)(coeff_t[k-1]+coeff_t[k]*(-div_t[k]));
    else
        return (int64_t)(coeff_t[0]);
}

void segmented_bit_sieve_wheel(uint64_t n,int64_t max_bW,int64_t num_c,int64_t delta_max)
{

    int64_t sqrt_n = (int64_t) std::sqrt(n);


    int64_t n_PB=(int64_t)3;    
    int64_t bW=(int64_t)30;
    //get bW base wheel equal to p1*p2*...*pn <=max_bW  with n=n_PB
    while(n_PB<n_PB_max&&(bW*PrimesBase[n_PB]<=std::min(max_bW,sqrt_n)))
    {
        bW*=PrimesBase[n_PB];
        n_PB++;
    }


    if (n>1+PrimesBase[n_PB-1]){

        int64_t k_end = (n < bW) ? (int64_t)2 :(int64_t) (n/(uint64_t)bW+1);
        int64_t k_sqrt = (int64_t) std::sqrt(k_end/bW)+1;

        //find possible remainder of base module
        std::vector<char> Remainder_i_t(bW+1,true); 
        for (int64_t i=0; i< n_PB;i++)
            for (int64_t j=PrimesBase[i]*PrimesBase[i];j< bW+1;j+=PrimesBase[i])
                Remainder_i_t[j]=false;
        std::vector<int64_t> RW;
        for (int64_t j=PrimesBase[n_PB-1]+1;j< bW+1;j++)
            if (Remainder_i_t[j]==true)
                RW.push_back(-bW+j);
        RW.push_back(1);
        int64_t  nR=RW.size();
        std::vector<int64_t> P_c(nR,0);
        std::vector<int64_t> C1(nR*nR);
        std::vector<int64_t> C2(nR*nR);
        for (int64_t j=0; j<nR-2; j++)
        {
            int64_t rW_t,rW_t1;
            rW_t1=Euclidean_Diophantine(bW,-RW[j]);
            for (int64_t i=0; i<nR; i++)
            {
                if (i==j)
                {
                    C2[nR*i+j]=0;
                    C1[nR*i+j]=RW[j]+1;
                }
                else if(i==nR-3-j )
                {
                    C2[nR*i+j]=1;
                    C1[nR*i+j]=RW[j]-1;
                }
                else
                {
                    rW_t=(int64_t)(rW_t1*(-RW[i]))%bW;
                    if (rW_t>1)
                        rW_t-=bW;
                    C1[nR*i+j]=rW_t+RW[j];
                    C2[nR*i+j]=(int64_t)(rW_t*RW[j])/bW+1;
                    if (i==nR-1)
                        C2[nR*i+j]-=1;
                }
            }
            C2[nR*j+nR-2]=(int64_t)1;
            C1[nR*j+nR-2]=-(bW+RW[j])-1;
            C1[nR*j+nR-1]=RW[j]+1;
            C2[nR*j+nR-1]=(int64_t)0;
        }
        for (int64_t i=nR-2; i<nR; i++)
        {
            C2[nR*i+nR-2]=(int64_t)0;
            C1[nR*i+nR-2]=-RW[i]-1;
            C1[nR*i+nR-1]=RW[i]+1;
            C2[nR*i+nR-1]=(int64_t)0;
        }


        int64_t  nB=nR/8;
        int64_t segment_size=1; 
        int64_t p_mask_i=(int64_t)4;
        for (int64_t i=0; i<p_mask_i;i++)
            segment_size*=(bW+RW[i]);  //if bW=30  =7*11*13*17
        while (segment_size<k_sqrt && p_mask_i<7)
        {
            segment_size*=(bW+RW[p_mask_i]);  //if bW=30 max value =7*11*13*17*19*23*29
            p_mask_i++;
        }
        
        int64_t segment_size_b=nB*segment_size;
        std::vector<uint8_t> Primes(nB+segment_size_b, 0xff);
        std::vector<uint8_t> Segment_i(nB+segment_size_b, 0xff);
 
        int64_t  pb,mb,mmin,ib,i,jb,j,k,kb;
        int64_t kmax = (int64_t) std::sqrt(segment_size/bW)+(int64_t)1;
        for (k =(int64_t)1; k  <= kmax; k++)
        {
            kb=k*nB;            
            for (jb = 0; jb<nB; jb++)
            {
                for (j = 0; j<8; j++)
                {
                    if(Primes[kb+jb] & (1 << j))
                    {
                        for (ib = 0; ib<nB; ib++)
                        {
                            for (i = 0; i<8; i++)
                            {
                                pb=nB*(bW*k+RW[j+jb*8]);
                                mmin=nB*(bW*k*k + k*C1[(i+ib*8)*nR+j+jb*8] + C2[(i+ib*8)*nR+j+jb*8]);
                                for (mb =mmin; mb <= segment_size_b && mb>=(int64_t)0; mb +=pb )
                                    Primes[mb+ib] &= del_bit[i];
                                if (pb<nB*(bW+RW[p_mask_i]) && k_end>segment_size)
                                {
                                    mb-=segment_size_b;
                                    while (mb<(int8_t)0)
                                        mb+=pb;
                                    for (; mb <= segment_size_b; mb +=pb )
                                        Segment_i[mb+ib] &= del_bit[i];
                                }
                            }
                        }
                    }
                }
            }
        }
        for (int64_t kc = 1; kc < std::min (segment_size+1,k_end); kc++)
        {
            int64_t count_p=0;
            for (int64_t ic = 0; ic<nB; ic++)
                count_p+=bit_count[Primes[nB*kc+ic]]; 
            if (count_p>=num_c)
            {
                int64_t i_P_c=0;
                for (int64_t ic = 0; ic<nB; ic++)
                    for (int64_t ic1 = 0; ic1  < 8; ic1++)
                        if(Primes[nB*kc+ic]& (1 << ic1))
                        {
                            P_c[i_P_c]=RW[ic1+ic*8]+kc*bW;
                            i_P_c++;
                        }
                for (int64_t ic = 0; ic<=i_P_c-num_c; ic++)
                {
                    if (P_c[ic+num_c-1]-P_c[ic]<=delta_max)
                    {
                        if (P_c[ic+num_c-1]-P_c[ic]<delta_max)
                            std::cout<< "delta max: "<<P_c[ic+num_c-1]-P_c[ic]<<std::endl;
                        std::cout<< "prime: "<<P_c[ic]<<" + d"<<std::endl;
                        std::cout<<std::endl;
                        std::cout<< "d: 0";
                        for (int64_t ic1 = 1; ic1<num_c; ic1++)
                                std::cout<<" , "<<P_c[ic+ic1]-P_c[ic];
                        std::cout<<std::endl;
                        std::cout<<std::endl;
                    }
                }
            }
        }

 
        if (k_end>segment_size) 
        {
            int64_t k_low,kb_low; 
            std::vector<uint8_t> Segment_t(nB+segment_size_b);
            for (int64_t k_low = segment_size; k_low < k_end; k_low += segment_size)
            {
                kb_low=k_low*nB;
                for (kb = (int64_t)0; kb <(nB+segment_size_b); kb++)
                    Segment_t[kb]=Segment_i[kb]; 
                kmax=(std::min(segment_size,(int64_t)std::sqrt((k_low+segment_size)/bW)+2));
                j=p_mask_i;
                for(k=(int64_t)1; k<=kmax;k++)
                {
                    kb=k*nB;
                    for (jb = 0; jb<nB; jb++)
                    {
                        for (; j < 8; j++)
                        {
                            if (Primes[kb+jb]& (1 << j))
                            {
                                for (ib = 0; ib<nB; ib++)
                                {
                                    for (i = 0; i < 8; i++)
                                    {
                                        pb=bW*k+RW[j+jb*8];
                                        mmin=-k_low+bW*k*k+ k*C1[(i+ib*8)*nR+j+jb*8] + C2[(i+ib*8)*nR+j+jb*8];
                                        if (mmin<0)
                                            mmin=(mmin%pb+pb)%pb;
                                        mmin*=nB;
                                        pb*=nB;
                                        for (mb =mmin; mb <= segment_size_b; mb += pb)
                                            Segment_t[mb+ib] &= del_bit[i];
                                    }
                                }
                            }
                        }
                        j=(int64_t)0;
                    }
                }
                for (int64_t kc = 1; kc <std::min (segment_size+1,k_end-k_low); kc++)
                {
                    int64_t count_p=0;
                    for (int64_t ic = 0; ic<nB; ic++)
                        count_p+=bit_count[Segment_t[nB*kc+ic]]; 
                    if (count_p>=num_c)
                    {
                        int64_t i_P_c=0;
                        for (int64_t ic = 0; ic<nB; ic++)
                            for (int64_t ic1 = 0; ic1  < 8; ic1++)
                                if(Segment_t[nB*kc+ic]& (1 << ic1))
                                {
                                    P_c[i_P_c]=RW[ic1+ic*8]+(kc+k_low)*bW;
                                    i_P_c++;
                                }
                        for (int64_t ic = 0; ic<=i_P_c-num_c; ic++)
                        {
                            if (P_c[ic+num_c-1]-P_c[ic]<=delta_max)
                            {
                                if (P_c[ic+num_c-1]-P_c[ic]<delta_max)
                                    std::cout<< "delta max: "<<P_c[ic+num_c-1]-P_c[ic]<<std::endl;
                                std::cout<< "prime: "<<P_c[ic]<<" + d"<<std::endl;
                                std::cout<<std::endl;
                                std::cout<< "d: 0";
                                for (int64_t ic1 = 1; ic1<num_c; ic1++)
                                    std::cout<<" , "<<P_c[ic+ic1]-P_c[ic];
                                std::cout<<std::endl;
                                std::cout<<std::endl;
                            }

                        }
                    }
                }
            }

        }
    }

}

int main()
{
    //segmented_bit_sieve_wheel(n,max_bW, num tuple , delta max) with max base wheel max_bW= 30 , 210 , 2310
    segmented_bit_sieve_wheel(100000,210, 8, 26);   

    return 0;
}

I corrected a error of course numberof elements must be <= 8 if using wheel 30 and <= 48 if using wheel 210

[User140242]

In reference to the code in post #15 I want to specify that I have only modified a sieve that I had made to count the prime numbers.

In theory it works for n<2^63 and maybe for n<2^64 if you change std::vector<int64_t>P_c(nR, 0); in std::vector<uint64_t>P_c (nR, 0); but for such large n is slow, you have to use multithreading.

I am not an expert and I do not have the necessary experience but I think you can make m different blocks for multithreading in the second part of the code, using m different Segment_t vectors and increasing for each block k_low by m*segment_size.

It is also possible to increase the size of the wheel by adding prime numbers to PrimesBase the problem is that for p>13 the size (number of classes)^2 of the arrays C1 and C2 becomes too large, you have to modify the code a little to use two vectors of size (number of classes) to be calculated from time to time for each single class used and this decreases the speed.

Although by increasing the size of the wheel it is possible to reach higher values of n because the algorithm works with numbers smaller than max(sqrt(n),n/bW) and by increasing bW it is possible to increase n always modifying the code and entering sqrt(n) and n/bW as input.

[ATH]

Originally Posted by Cybertronic

I believe, to reach this offset for 19-tuplet is 10x faster...~2 days, or 1.5e24 per month (or an offset with 6e17+ per sec).

How would you get 6e17+ per sec? On many cores or on GPU or do you have a highly optimized assembly program for this?

On 1 core on a very old 5960X corei7 I'm getting 6.3e13 per sec. I know my program is pretty basic without assembly, but I was expecting a very good program to be 10x to 100x faster at most, not 10000x faster.

[Cybertronic]

Originally Posted by ATH

How would you get 6e17+ per sec? On many cores or on GPU or do you have a highly optimized assembly program for this?

On 1 core on a very old 5960X corei7 I'm getting 6.3e13 per sec. I know my program is pretty basic without assembly, but I was expecting a very good program to be 10x to 100x faster at most, not 10000x faster.

Okay, forgot to say 6e17 (or more) is on 32 threads (two Ryzen 7 1700 @3GHz, round 3e16 per core and sec,without SMT.).
This target runs with 53# blocksize under my sparse freebasic code....note, I'm not a programming expert.
The searching is systematic with all 32 threads simultaneous...for this task I use an offset-range of 53#*2000=6.5e22...tooks maybe a day or so.


BTW, here is my german k-tuplet thread on "Matroids Matheplanet", 12 pages, there are a lot of running informations.

https://matheplanet.de/matheplanet/n...232720&start=0




I can try to recode the 17-tuplet for this special 19-tuplet ..so I have the exact rate/s. I have not try this yet.


After I found the smallest 7k triplet 2nd kind, I will look for the smallest 55-digit prime 12-tuplet. This will take 1 week to reach an offset 1e21.

[Cybertronic]

Testsieve is done.

Pattern: d=00,04,06,10,12,16,24,30,34,40,42,46,52,54,60,66,70,72,76


I get per core 3.36e16 / s [ 1e18/s ! with 32 threads] . The program works in cycles and one is done in 11.2 days with complete offset up to 32589158477190044730000
Deep sieving max. p<48000

My estimation is one 19-tuplet ~3000 cycles. So ~10 days with 3000 cores. RAM in using 50MB per task.



Even I found:
14 / 19 conditions are prime 28026688763307365621377
12 / 19 conditions are prime 15642812740222210642987
9 / 19 conditions are prime 20107537050380847377197
11 / 19 conditions are prime 3128586964390514346607

[Cybertronic]

I wrote on my page:
k=19 s=76 B={0 4 6 10 12 16 24 30 34 40 42 46 52 54 60 66 70 72 76}
smallest is unknown
But one is known !

622803914376064301858782434517 + d, d = 0, 4, 6, 10, 12, 16, 24, 30, 34, 40, 42, 46, 52, 54, 60, 66, 70, 72, 76
I must rewrite the code for the smallest nontrivial , this is unknown .... d = 0, 4, 6, 10, 16, 18, 24, 28, 30, 34, 40, 46, 48, 54, 58, 60, 66, 70, 76 [trivial: 13,17,...,89] ; 217153 (modulo 510510). For this pattern, the speed is 3x faster... ~ 1e17 /s and core.

[Cybertronic]

Edit:
Sieve for d = 0, 4, 6, 10, 16, 18, 24, 28, 30, 34, 40, 46, 48, 54, 58, 60, 66, 70, 76 is also ready.

here an example with 14/19 are prime , first 10 conditions are prime

4008486315568214747203+d,d=0, 4, 6, 10, 16, 18, 24, 28, 30, 34, (40), 46, (48), (54), 58, 60, (66), 70, (76)

Ratio 1:2,7


Running time 5 days per cycle and one core up to 32589158477190044730000

[Cybertronic]

I have run now 32 threads...will see, how many conditions in order are true ... up to 32589158477190044730000
Tooks me 4h, so the offset progress (averaged) is 7e16/s and core


After 30min: first with 13/19 in order.
3718239377799223934593+d, d=00,04,06,10,16,18,24,28,30,34,40,46,48 ; [like primes 13,..,61]