cp's OEIS Frontend

This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.

A051015 Zeisel numbers.

Original entry on oeis.org

105, 1419, 1729, 1885, 4505, 5719, 15387, 24211, 25085, 27559, 31929, 54205, 59081, 114985, 207177, 208681, 233569, 287979, 294409, 336611, 353977, 448585, 507579, 982513, 1012121, 1073305, 1242709, 1485609, 2089257, 2263811, 2953711, 3077705, 3506371, 3655861, 3812599
Offset: 1

Views

Author

Keywords

Comments

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.

Crossrefs

Cf. A027746, A061422, A252094 (A values), A252095 (B values).

Programs

  • Haskell
    a051015 n = a051015_list !! (n-1)
    a051015_list = filter zeisel [3, 5 ..] where
       zeisel x = 0 `notElem` ds && length ds > 2 &&
             all (== 0) (zipWith mod (tail ds) ds) && all (== q) qs
             where q:qs = (zipWith div (tail ds) ds)
                   ds = zipWith (-) (tail ps) ps
                   ps = 1 : a027746_row x
    -- Reinhard Zumkeller, Dec 15 2014
  • Mathematica
    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 *)
  • 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
    

Extensions

More terms from David Wasserman, Feb 19 2002
Extended by Karsten Meyer, Jun 08 2006, but values were incorrect. M. F. Hasler, Oct 31 2012
Values up to a(70) computed by Jean-François Alcover and double-checked by M. F. Hasler, Oct 31 2012
Values < 10^15 by Lars Blomberg, Nov 02 2012