[seqfan] A051015 Zeisel numbers

[Jean-François Alcover]

Dear SeqFans,

I'm trying to code a Mma program for A051015 (Zeisel numbers)
Unfortunately, output from Mma shows a few discrepances with data.
Can someone help me to find what's wrong with my program ?

http://oeis.org/A051015

[assuming all terms are odd]
ZeiselQ[n_] :=
 If[PrimeQ[n] || ! SquareFreeQ[n], False,
 pp = Join[{1}, FactorInteger[n][[All, 1]]];
 If[Length[pp] <= 3, False,
 eq = Thread[Rest[pp] == Most[pp] a + b];
 r = Reduce[eq, {a, b}, Integers]; r =!= False]];
Reap[For[n = 1, n <= 3*10^7, n = n + 2,
 If[ZeiselQ[n], Print[n]; Sow[n]]]][[2, 1]]

Absent from data :
{3812599, 4293793, 4813879, 6622609, 8835799, 10235499, 13534129,
14662681, 16596371, 18651417, 20138379, 22549639, 25711399}

ZeiselQ[3812599]
True

pp
{1, 7, 151, 3607}

eq
{7 == a + b, 151 == 7 a + b, 3607 == 151 a + b}

r
a == 24 && b == -17


Absent from Reap output :
{721907}

ZeiselQ[721907]
False

pp
{1, 37, 109, 179}

eq
{37 == a + b, 109 == 37 a + b, 179 == 109 a + b}

r
False