rd[n_] := rd[n] = RealDigits[ N[1 + Sum[1 - 1/Zeta[j], {j, 2, 2^n}], 400]][[1]]; rd[n = 4]; While[rd[n] =!= rd[n-1], n++]; Niven = FromDigits[{rd[n], 1}];
A033151 = ContinuedFraction[Niven]; a[n_] := Position[
A033151, n][[1, 1]]; Table[a[n], {n, 1, 18}] (*
Jean-François Alcover, Oct 31 2012 *)