(* On a Power Mac, this Mathematica program takes about 40 minutes *)
nn = PrimePi[8191];
t1 = {2};
t5 = {31, 1801};
t7 = {127};
t13 = {8191}; found = Union[t1, t5, t7, t13]; t2 = {};
Do[ps = Position[(2^Range[k/2] + 1)/k, _?(IntegerQ[#] &), 1, 1];
If[ps != {}, AppendTo[found, k]; AppendTo[t2, {k, ps[[1, 1]]}]], {k,
Prime[Range[2, nn]]}]; Length[found]
t3 = {};
Do[Do[If[! MemberQ[found, k],
ps = Position[(2^j + 2^Range[j - 1] + 1)/k, _?(IntegerQ[#] &), 1, 1];
If[ps != {}, AppendTo[found, k];
AppendTo[t3, {k, j, ps[[1, 1]]}]; Break[]]], {j, 2, k/3}], {k,
Prime[Range[2, nn]]}]; Length[found]
t4 = {};
Do[Do[If[! MemberQ[found, k],
ps = Position[(2^j + 2^i + 2^Range[i - 1] + 1)/k, _?(IntegerQ[#] &), 1, 1];
If[ps != {}, AppendTo[found, k];
AppendTo[t4, {k, j, i, ps[[1, 1]]}]; Break[]]], {j, 3, k/10}, {i, 2, j}],
{k, Prime[Range[2, nn]]}]; Length[found]
s1 = {{2, 1}};
s2 = Table[{i[[1]], (2^i[[2]] + 1)/i[[1]]}, {i, t2}];
s3 = Table[{i[[1]], (2^i[[2]] + 2^i[[3]] + 1)/i[[1]]}, {i, t3}];
s4 = Table[{i[[1]], (2^i[[2]] + 2^i[[3]] + 2^i[[4]] + 1)/i[[1]]}, {i,
t4}];
s5 = Transpose[{t5, Table[1, {Length[t5]}]}];
s7 = {Append[t7, 1]};
s13 = {Append[t13, 1]};
seq2 = Sort[Join[s1, s2, s3, s4, s5, s7, s13]];
tAll = Transpose[seq2][[2]];