a(n) = {
	
	maxp = 11;
	r = remaindersuptop(n, maxp);
	d = diff(r);
	s = vecprod(primes(primepi(maxp)))/3;
	d = concat(d, r[1] + s - r[#r]);
	inc = 10^12;
	printat = inc;
	forstep(i = r[1], oo, d,
		if(i > printat,
			print1(i \ inc", ");
			printat += inc;
		);
		if(is(i, n),
			return(i)
		)
	)
}

is(c, n) = {
	for(k = 1, n, 
		if(!ispseudoprime(9*k*c + 1),
			return(0)
		)
	);
	1
}

remaindersuptop(n, {maxp=31}) = {
	my(res = List(Mod(0, 1)), pr = primes(primepi(maxp)), q);
	pr = pr[^2];
	for(i = 1, #pr,
		q = #res;
		rems = remaindersmodp(pr[i], n); 
		nres = List(); 	
		for(j = 1, q, 
			for(k = 1, #rems,
				c = chinese(res[j], Mod(rems[k], pr[i]));
				listput(nres, c)
			)
		);
		res = nres;
	);
	for(i = 1, #res, 
		res[i] = lift(res[i])
	);
	Set(res)
	
}

remaindersmodp(p, n) = {
	my(res = List(), i);
	for(i = 0, p-1,
		if(isrem(i, p, n),
			listput(res, i)
		)
	);
	res
	
}

isrem(r, p, n) = {
	for(k = 1, n, 
		c = lift(Mod(9*k*r+1, p));
		if(c == 0,
			return(0)
		)
	);
	return(1)
}

diff(v) = vector(#v-1, i, v[i+1]-v[i]);