[seqfan] Re: A051015 Zeisel numbers

[Giovanni Resta]

On 10/31/2012 08:39 AM, Jean-François Alcover wrote:
> 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 ?

I think that your program is correct.
I prefer to use local variables, to avoid possible errors.
Clearly the sequence in OEIS has errors.

If I may make a comment about your program, is that
it is quite elegant but very slow.

A naive version in which you actually compute a,b from
the first 2 terms and then test the others:

zeiQ[n_] := Catch[
   Block[{a, b, f = FactorInteger[n]},
    If[Length[f] < 3 || Max[Last /@ f] > 1, Throw[False]];
    f = First /@ f;
    a = (-f[[1]] + f[[2]])/(f[[1]] - 1);
    b = (f[[1]]^2 - f[[2]])/(f[[1]] - 1);
    If[! IntegerQ[a] || ! IntegerQ[b], Throw[False]];
    Do[If[a f[[i - 1]] + b != f[[i]], Throw[False]], {i, 3,
      Length[f]}]; True]];

is at least 10 times faster.

Clearly, if you want to build a large list of such numbers, it
may be better to generate them directly, instead of testing them.
Something like this (not tested, just written down...)

L = {}; p = 3; Lim = 10^9;
While[p^3 < Lim,
   q = NextPrime[p];
   While[p q^2 < Lim,
    If[! IntegerQ[a = (q - p)/(p - 1)] ||
       ! IntegerQ[b = (p^2 - q)/(p - 1)], q = NextPrime[q]; Continue[]];
    r = a q + b; n = r p q;
    While[PrimeQ[r] && n < Lim, AppendTo[L, n]; r = a r + b; n *= r];
    q = NextPrime[q]];
   p = NextPrime[p]];
(* sorting and testing just to be sure *)
L = Select[Sort[L], zeiQ]

Giovanni Resta