Software for finding prime constellations

[CRGreathouse]

Suppose I have an admissible tuple (0, n_1, n_2, ..., n_{k-1}) and I am interested in finding instances of this constellation: numbers N such that N, N + n_1, ..., N + n_{k-1} are all prime. Is there software out there to do this search efficiently?

In my immediate case, k = 19.

[CRGreathouse]

https://arxiv.org/abs/1807.08777 is relevant, but I don't know of any implementations.

[a1call]

1st SM and now Mr. Greathouse. Are the stars lining up somewhere?

My 2 cents:
PFGW it’s the most efficient primality tester and you can script it so that it checks the 1st then and if and only if Prime would test next and so on. It inserts a noticeable line in the output when all items match (are Prime).
It lacks an advanced sieving feature but you can (I have successfully, so should be a piece of cake for you) combine Pari-gp script and initiate pfgw when a sieved candidate is found.
Obvious not telling you anything new, why would you look for another software?
Nice to have you posting again. The forum had deteriorated without your inputs.

[MattcAnderson]

I wrote some Maple code to find the smallest few k-tuples for
k= 7,8,9,10,11,12,13. back in 2015
My Maple code and some prime k-tuples appear in OEIS.org
Dana Jacobsen and Norman Luhn have extended these lists (and others)
A link to this work is on my personal webpage
https://mattanderson.fun/f/prime-constellations
Also, Norman Luhn has a web page with lists of prime k-tuples, and
the form that these numbers have.
see - https://pzktupel.de/ktpatt_hl.php
for example, form Norman's webpage, about 19-tuples, if the 19-tuple is
{0 4 6 10 16 22 24 30 34 36 42 46 52 60 64 66 70 72 76} then the
smallest prime in the constellation will have the form
37 (modulo 30030)
I read this off the chart at pzktupel.de/ktpatt_hl.php
There are 4 patterns for 19-tuples that are also prime constellations.
They are shown on the chart.
Currently, there is no known prime 19-tuple with pattern
{0 4 6 10 12 16 24 30 34 40 42 46 52 54 60 66 70 72 76}
Initial primes in this 19-tuple will have the form
29917 (modulo 30030)
I pulled this information from this webpage
https://pzktupel.de/SMArchiv/smadditions.html

Maple code is slow, but easy for me to write code. I sped up the
Maple isprime() function by first testing if a number has any factors
less than 100, then, if all candidates of the k-tuple pass that test,
run the full isprime() function on all 'k' of the possible primes.
See my code for 10-tuples here
http://oeis.org/A027569

Good luck on being the first to find this 19-tuple.
I don't know of any downloadable software for finding
admissable k-tuples. Norman Luhn has his software but it is not
ready to be shared.

Also, if someone can find a certain 447-tuple, they can show a counterexample to a conjecture of
Hardy and Littlewood.
see - http://www.opertech.com/primes/k-tuples.html

Best of luck, and
have a nice day.

[Dr Sardonicus]

Originally Posted by CRGreathouse

Suppose I have an admissible tuple (0, n_1, n_2, ..., n_{k-1}) and I am interested in finding instances of this constellation: numbers N such that N, N + n_1, ..., N + n_{k-1} are all prime. Is there software out there to do this search efficiently?

In my immediate case, k = 19.

If you go to cybertronic's pzktupel site, you'll find some 19-tuples (and some longer ones) listed, along with who found them and when. Perhaps he could point you in the right direction.

[mart_r]

AFAIK, Norman Luhn, Jens Kruse Andersen and others all use highly customized code. My code for the 13-tuples some years ago also was highly customized, written in Basic, and too slow for 19-tuplets.
I haven't yet seen any general-purpose efficient software specifically for large k-tuples, but if there is one, I would also be interested

[mart_r]

I've written a Pari program that captures the gist of the most important steps in a k-tuple finding algorithm. Problem is, it is quite slow. Finding 12- or 13-tuples can take several minutes, so it's most likely not suitable for 19-tuples.

print("ktuple(pattern as vector): tries to find tuples of primes with given pattern");
ktuple(p)=
{
h=2048; \\ lower bound in number of pre-wheel-sieved values - can be changed manually if so desired
z=2^16; \\ trial division limit; tuples larger than z^2 get TD'd up to z, after that, numbers are checked sequentially if prime - can also be changed if necessary

p=vecsort(p); \\ just in case pattern is not strictly increasing
\\ check if pattern is admissible (a=1:yes, a=0:no)
a=1;
forprime(q=2,#p,
	v=vector(q);
	for(i=1,#p,
		v[q-p[i]%q]=1
	);
	v=vecsort(v);
	if(v[1]==1,
		print("pattern is not admissible / all residues mod "q" are covered");
		a=0;
		break()
	)
);
\\ find starting point s and difference pattern d; wheel of length q# for first q such that #d > h
if(a,
	s=p[1]%2+1;
	d=[2];
	q=1;
	m=2;
	while(#d<h,
		q=nextprime(q+2);
		m*=q;
		t=s;
		f=0;
		b=1;
		while(b,
			for(i=1,#p,
				b*=sign((t+p[i])%q)
			);
			if(b==0,
				f=f+1-#d*(f==#d);
				t+=d[f];
				b=1
			,
				b=0
			)
		);
		s=t;
		u=t;
		e=[];
		while(t-s<m,
			f=f+1-#d*(f==#d);
			u+=d[f];
			b=1;
			while(b,
				for(i=1,#p,
					b*=sign((u+p[i])%q)
				);
				if(b==0,
					f=f+1-#d*(f==#d);
					u+=d[f];
					b=1
				,
					e=concat(e,u-t);
					t=u;
					b=0
				)
			)
		);
		d=e;
		print1(#d" admissible pattern");
		if(#d>1,
			print1("s")
		);
		print(" in wheel of length "m)
	);
	print();
\\ residue patterns in vectors
	w=p[#p]-p[1];
	o=vector(w);
	for(i=1,#p-1,
		o[p[#p]-p[i]]=1
	);
	l=primepi(w)-primepi(q);
	q=nextprime(q+2);
	if(l>0,
		g=vector(w);
		n=vector(l);
		k=0;
		forprime(c=q,w,
			k++;
			g[c]=k;
			n[k]=vector(c);
			for(i=1,#p-1,
				j=(p[#p]-p[i])%c;
				if(j,
					n[k][j]=1
				)
			)
		)
	);
	s+=p[#p];
	f=0;
	print("searching for primes with {pattern} = "p);
\\ main search
	while(s<z^2,
		b=1;
		forprime(c=q,w,
			r=s%c;
			if(r==0,
				b=0;
				break()
			,
				if(n[g[c]][r],
					b=0;
					break()
				)
			)
		);
		if(b,
			forprime(c=w+1,sqrt(s),
				r=s%c;
				if(r==0,
					b=0;
					break()
				,
					if(r<=w,
						if(o[r],
							b=0;
							break()
						)
					)
				)
			);
			if(b,
				print(s-p[#p]" + {pattern}")
			)
		);
		f=f+1-#d*(f==#d);
		s+=d[f]
	);
	while(1,
		b=1;
		forprime(c=q,w,
			r=s%c;
			if(r==0,
				b=0;
				break()
			,
				if(n[g[c]][r],
					b=0;
					break()
				)
			)
		);
		if(b,
			forprime(c=w+1,z,
				r=s%c;
				if(r==0,
					b=0;
					break()
				,
					if(r<=w,
						if(o[r],
							b=0;
							break()
						)
					)
				)
			);
			if(b,
				for(i=1,#p,
					if(!isprime(s-p[#p]+p[i]),
						b=0;
						break()
					)
				);
				if(b,
					print(s-p[#p]" + {pattern}")
				)
			)
		);
		f=f+1-#d*(f==#d);
		s+=d[f]
	)
)
}

It can however serve as a general-purpose program for finding smaller tuples with only a given pattern:

(19:53) gp > p=[0,2,6,12,20,30,42,56,72,90,110]
%50 = [0, 2, 6, 12, 20, 30, 42, 56, 72, 90, 110]
(19:55) gp > ktuple(p)
1 admissible pattern in wheel of length 6
2 admissible patterns in wheel of length 30
6 admissible patterns in wheel of length 210
30 admissible patterns in wheel of length 2310
180 admissible patterns in wheel of length 30030
1440 admissible patterns in wheel of length 510510
12960 admissible patterns in wheel of length 9699690

searching for primes with {pattern} = [0, 2, 6, 12, 20, 30, 42, 56, 72, 90, 110]
180078317 + {pattern}
1278189947 + {pattern}
1761702947 + {pattern}
1829187287 + {pattern}
5862143447 + {pattern}

[LaurV]

Originally Posted by CRGreathouse

Suppose I have ..

Welcome back!

[science_man_88]

Originally Posted by mart_r

print("ktuple(pattern as vector): tries to find tuples of primes with given pattern");
ktuple(p)=
{
h=2048; \\ lower bound in number of pre-wheel-sieved values - can be changed manually if so desired
z=2^16; \\ trial division limit; tuples larger than z^2 get TD'd up to z, after that, numbers are checked sequentially if prime - can also be changed if necessary

p=vecsort(p); \\ just in case pattern is not strictly increasing
\\ check if pattern is admissible (a=1:yes, a=0:no)
a=1;
forprime(q=2,#p,
	v=vector(q);
	for(i=1,#p,
		v[q-p[i]%q]=1
	);
	v=vecsort(v);
	if(v[1]==1,
		print("pattern is not admissible / all residues mod "q" are covered");
		a=0;
		break()
	)
);
\\ find starting point s and difference pattern d; wheel of length q# for first q such that #d > h
if(a,
	s=p[1]%2+1;
	d=[2];
	q=1;
	m=2;
	while(#d<h,
		q=nextprime(q+2);
		m*=q;
		t=s;
		f=0;
		b=1;
		while(b,
			for(i=1,#p,
				b*=sign((t+p[i])%q)
			);
			if(b==0,
				f=f+1-#d*(f==#d);
				t+=d[f];
				b=1
			,
				b=0
			)
		);
		s=t;
		u=t;
		e=[];
		while(t-s<m,
			f=f+1-#d*(f==#d);
			u+=d[f];
			b=1;
			while(b,
				for(i=1,#p,
					b*=sign((u+p[i])%q)
				);
				if(b==0,
					f=f+1-#d*(f==#d);
					u+=d[f];
					b=1
				,
					e=concat(e,u-t);
					t=u;
					b=0
				)
			)
		);
		d=e;
		print1(#d" admissible pattern");
		if(#d>1,
			print1("s")
		);
		print(" in wheel of length "m)
	);
	print();
\\ residue patterns in vectors
	w=p[#p]-p[1];
	o=vector(w);
	for(i=1,#p-1,
		o[p[#p]-p[i]]=1
	);
	l=primepi(w)-primepi(q);
	q=nextprime(q+2);
	if(l>0,
		g=vector(w);
		n=vector(l);
		k=0;
		forprime(c=q,w,
			k++;
			g[c]=k;
			n[k]=vector(c);
			for(i=1,#p-1,
				j=(p[#p]-p[i])%c;
				if(j,
					n[k][j]=1
				)
			)
		)
	);
	s+=p[#p];
	f=0;
	print("searching for primes with {pattern} = "p);
\\ main search
	while(s<z^2,
		b=1;
		forprime(c=q,w,
			r=s%c;
			if(r==0,
				b=0;
				break()
			,
				if(n[g[c]][r],
					b=0;
					break()
				)
			)
		);
		if(b,
			forprime(c=w+1,sqrt(s),
				r=s%c;
				if(r==0,
					b=0;
					break()
				,
					if(r<=w,
						if(o[r],
							b=0;
							break()
						)
					)
				)
			);
			if(b,
				print(s-p[#p]" + {pattern}")
			)
		);
		f=f+1-#d*(f==#d);
		s+=d[f]
	);
	while(1,
		b=1;
		forprime(c=q,w,
			r=s%c;
			if(r==0,
				b=0;
				break()
			,
				if(n[g[c]][r],
					b=0;
					break()
				)
			)
		);
		if(b,
			forprime(c=w+1,z,
				r=s%c;
				if(r==0,
					b=0;
					break()
				,
					if(r<=w,
						if(o[r],
							b=0;
							break()
						)
					)
				)
			);
			if(b,
				for(i=1,#p,
					if(!isprime(s-p[#p]+p[i]),
						b=0;
						break()
					)
				);
				if(b,
					print(s-p[#p]" + {pattern}")
				)
			)
		);
		f=f+1-#d*(f==#d);
		s+=d[f]
	)
)
}

It can however serve as a general-purpose program for finding smaller tuples with only a given pattern:

(19:53) gp > p=[0,2,6,12,20,30,42,56,72,90,110]
%50 = [0, 2, 6, 12, 20, 30, 42, 56, 72, 90, 110]
(19:55) gp > ktuple(p)
1 admissible pattern in wheel of length 6
2 admissible patterns in wheel of length 30
6 admissible patterns in wheel of length 210
30 admissible patterns in wheel of length 2310
180 admissible patterns in wheel of length 30030
1440 admissible patterns in wheel of length 510510
12960 admissible patterns in wheel of length 9699690

searching for primes with {pattern} = [0, 2, 6, 12, 20, 30, 42, 56, 72, 90, 110]
180078317 + {pattern}
1278189947 + {pattern}
1761702947 + {pattern}
1829187287 + {pattern}
5862143447 + {pattern}

Looks like it should be highly paralleled.

[mart_r]

Originally Posted by science_man_88

Looks like it should be highly paralleled.

Thanks for the reply! You make me happy

Different patterns can be looked at in parallel, but manually setting a starting point would also be nice. Should be no problem for me to implement.
I have a couple more ideas for improvements. But I also have to earn some money, so give me a day or three...

[mart_r]

It would still take thousands of years to find a prime novemdecuplet with my program, but it might be useful for small searches and doublechecks.

- Working with "concat" was slowing the wheel sieve down. Fixed that, it's now way faster.
- The search itself is now also about 20% faster (for medium-sized tuples). And since the wheel can now be made larger, another significant speedup is possible for large tuples.
- A starting point can be set. Example: "ktuple([1,3,7,9],10^49)" finds the first quadruplets with 50 digits.
- (Did anyone notice that the first member of the input pattern need not be zero? You can also search for quintuplets with pattern [-4,0,2,6,8], for instance.)
- Optional write to file and optional output tracking. I've finally learned how to use the input command properly.
- There's a lot of room for optimisation by experimenting with different wheel sieve sizes and trial division depths. This is left as an exercise for the user ;)
- Any errors or bugs? Feel free to report.

ktuple(p,x)=
{
h=8192; \\ lower bound in number of pre-wheel-sieved values - can be changed manually if so desired
z=1024*#p; \\ trial division limit - can also be changed if necessary
print1("Write output file "Strchr(34)"ktuple_output.txt"Strchr(34)"? 1:yes / 0:no ");l=input();
print1("Output search status every ___ seconds? (0/default: don't output search status, only tuples) ");j=input();

p=vecsort(p); \\ just in case pattern is not strictly increasing
\\ check if pattern is admissible (a=1:yes, a=0:no)
a=1;
y=vector(#p);
forprime(q=2,#p,
	v=vector(q);
	for(i=1,#p,
		y[q]+=1-v[q-p[i]%q];
		v[q-p[i]%q]=1
	);
	v=vecsort(v);
	if(v[1]==1,
		print("pattern is not admissible / all residues mod "q" are covered");
		a=0;
		break()
	)
);
\\ find admissible initial members s[f]; generate wheel of size q# := m for first q such that d := #s[f] >= h
if(a,
	s=[p[1]%2+1];
	d=1;
	q=1;
	m=2;
	while(d<h,
		q=nextprime(q+2);
		if(q>#p,
			c=0;
			v=vector(q);
			for(i=1,#p,
				c+=1-v[q-p[i]%q];
				v[q-p[i]%q]=1
			)
		,
			c=y[q]
		);
		d*=(q-c);
		t=vector(d);
		f=0;
		u=0;
		e=0;
		b=1;
		while(e<d,
			f++;
			if(f>#s,
				f=1;
				u+=m
			);
			for(i=1,#p,
				b*=sign((u+s[f]+p[i])%q)
			);
			if(b,
				e++;
				t[e]=u+s[f]
			,
				b=1
			)
		);
		s=t;
		m*=q;
		print1(d" admissible pattern");
		if(d>1,
			print1("s")
		);
		print(" in wheel of size "m)
	);
	print();
\\ residue patterns in vectors
	w=p[#p]-p[1];
	o=vector(w);
	for(i=1,#p-1,
		o[p[#p]-p[i]]=1
	);
	f=primepi(w)-primepi(q);
	q=nextprime(q+2);
	if(f>0,
		g=vector(w);
		n=vector(f);
		k=0;
		forprime(c=q,w,
			k++;
			g[c]=k;
			n[k]=vector(c);
			for(i=1,#p-1,
				e=(p[#p]-p[i])%c;
				if(e,
					n[k][e]=1
				)
			)
		)
	);
	u=(x\m)*m+p[#p];
	if(u==p[#p],
		print("Patterns with initial member <= "precprime(q-2)" are not shown in the output!");
		print();
		if(l,
			write("ktuple_output.txt",Str("Patterns with initial member <= "precprime(q-2)" are not shown in the output!"))
		)
	);
	print("Searching for primes with {pattern} = "p);
	if(l,
		write("ktuple_output.txt",Str("Searching for primes with {pattern} = "p))
	);
\\ main search (stage 1: complete trial division)
	j=(j*1000\256)*256;
	t=0;
	gettime();
	while(u<z^2,
		for(f=1,d,
			x=u+s[f];
			b=1;
			forprime(c=q,w,
				r=x%c;
				if(r==0,
					b=0;
					break()
				,
					if(n[g[c]][r],
						b=0;
						break()
					)
				)
			);
			if(b,
				forprime(c=w+1,sqrt(x),
					r=x%c;
					if(r==0,
						b=0;
						break()
					,
						if(r<=w,
							if(o[r],
								b=0;
								break()
							)
						)
					)
				);
				if(b,
					print(x-p[#p]" + {pattern}");
					if(l,
						write("ktuple_output.txt",Str(x-p[#p]" + {pattern}"))
					)
				)
			)
		);
		u+=m;
		if(j,
			t+=gettime();
			if(t>j,
				t=0;
				print("searching ... "u-p[#p])
			)
		)
	);
\\ main search (stage 2: trial division up to z, then PRP test)
	while(1,
		for(f=1,d,
			x=u+s[f];
			b=1;
			forprime(c=q,w,
				r=x%c;
				if(r==0,
					b=0;
					break()
				,
					if(n[g[c]][r],
						b=0;
						break()
					)
				)
			);
			if(b,
				forprime(c=w+1,z,
					r=x%c;
					if(r==0,
						b=0;
						break()
					,
						if(r<=w,
							if(o[r],
								b=0;
								break()
							)
						)
					)
				);
				if(b,
					for(i=1,#p,
						if(!ispseudoprime(x-p[#p]+p[i]),
							b=0;
							break()
						)
					);
					if(b,
						print(x-p[#p]" + {pattern}");
						if(l,
							write("ktuple_output.txt",Str(x-p[#p]" + {pattern}"))
						)
					)
				)
			)
		);
		u+=m;
		if(j,
			t+=gettime();
			if(t>j,
				t=0;
				print("searching ... "u-p[#p])
			)
		)
	)
)
};
print("ktuple(pattern as vector [, starting point]): tries to find tuples of primes with given pattern [optionally start at x]")