The OEIS is supported by the many generous donors to the OEIS Foundation.
%I #35 Feb 16 2025 08:32:41
%S 105,1419,1729,1885,4505,5719,15387,24211,25085,27559,31929,54205,
%T 59081,114985,207177,208681,233569,287979,294409,336611,353977,448585,
%U 507579,982513,1012121,1073305,1242709,1485609,2089257,2263811,2953711,3077705,3506371,3655861,3812599
%N Zeisel numbers.
%C Pick any integers A and B and consider the linear recurrence relation given by p(0) = 1, p(i + 1) = A*p(i) + B. If for some n > 2, p(1), p(2), ..., p(n) are distinct primes, then the product of these primes is called a Zeisel number.
%H M. F. Hasler and Lars Blomberg, <a href="/A051015/b051015.txt">Table of n, a(n) for n = 1..9607</a> (first 70 terms from _M. F. Hasler_)
%H Kevin S. Brown, <a href="https://www.mathpages.com/home/kmath015.htm">Zeisel Numbers</a>, MathPages website.
%H OEIS Wiki, <a href="https://oeis.org/wiki/Zeisel_numbers">Zeisel numbers</a>.
%H Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/ZeiselNumber.html">Zeisel Number</a>.
%H Wikipedia, <a href="http://en.wikipedia.org/wiki/Zeisel_number">Zeisel number</a>.
%H Helmut Zeisel, <a href="https://groups.google.com/g/sci.math/c/m3iN-ygRjr0/m/ATnWPFeV01QJ">Primes of the form 2^(k-1)+k</a>, sci.math newsgroup, February 24, 1994.
%t maxTerm = 3*10^7; ZeiselQ[n_] := Module[{a, b, pp, eq, r}, If[PrimeQ[n] || ! SquareFreeQ[n], False, pp = Join[{1}, FactorInteger[n][[All, 1]]]; If[Length[pp] <= 3, False, eq = Thread[Rest[pp] == b + a*Most[pp]]; r = Reduce[eq, {a, b}, Integers]; r =!= False]]]; p = 3; A051015 = Reap[While[p^3 < maxTerm, q = NextPrime[p]; While[p*q^2 < maxTerm, If[ ! IntegerQ[a = (q - p)/(p - 1)] || !IntegerQ[b = (p^2 - q)/(p - 1)], q = NextPrime[q]; Continue[]]; r = b + a*q; n = r*p*q; While[PrimeQ[r] && n < maxTerm, Sow[n]; r = b + a*r; n *= r]; q = NextPrime[q]]; p = NextPrime[p]]][[2, 1]]; A051015 = Select[Sort[A051015], ZeiselQ] (* _Jean-François Alcover_, Oct 31 2012, with much help from _Giovanni Resta_ *)
%o (PARI) is_A051015(n)={#(n=factor(n)~)>2 & vecmax(n[2,])==1 & denominator(n[2,1]=(n[1,3]-n[1,2])/(n[1,2]-n[1,1]))==1 & #Set(n[1,]-n[2,1]*concat(1,vecextract(n[1,],"^-1")))==1} \\ - _M. F. Hasler_, Oct 31 2012
%o (Haskell)
%o a051015 n = a051015_list !! (n-1)
%o a051015_list = filter zeisel [3, 5 ..] where
%o zeisel x = 0 `notElem` ds && length ds > 2 &&
%o all (== 0) (zipWith mod (tail ds) ds) && all (== q) qs
%o where q:qs = (zipWith div (tail ds) ds)
%o ds = zipWith (-) (tail ps) ps
%o ps = 1 : a027746_row x
%o -- _Reinhard Zumkeller_, Dec 15 2014
%Y Cf. A027746, A061422, A252094 (A values), A252095 (B values).
%K nonn
%O 1,1
%A _Eric W. Weisstein_
%E More terms from _David Wasserman_, Feb 19 2002
%E Extended by _Karsten Meyer_, Jun 08 2006, but values were incorrect. _M. F. Hasler_, Oct 31 2012
%E Values up to a(70) computed by _Jean-François Alcover_ and double-checked by _M. F. Hasler_, Oct 31 2012
%E Values < 10^15 by _Lars Blomberg_, Nov 02 2012