OFFSET
1,1
COMMENTS
For every k >= 0, the sequence includes 4^k + 2^(k+1), with m = 2^k + 1. - David Wasserman, Jan 29 2004
So a(13) <= 4198400. - Michel Marcus, Aug 10 2014
Are there other terms like 4374 that are not of this form? - Michel Marcus, Aug 10 2014
FORMULA
G.f.: Conjecture: Q(0)/x - 1/x where Q(k)= 1 + 2^k*x/(1 - 2*x/(2*x + 2^k*x/Q(k+1) )); (continued fraction ). - Sergei N. Gladkovskii, Apr 10 2013
EXAMPLE
With n=3 and m=2, rad(3) = rad(3) and rad(2) = rad(4), so 3 is in the sequence.
MAPLE
rad:= n -> convert(numtheory:-factorset(n), `*`):
count:= 0: lastr:= rad(1):
for n from 2 to 10^7 do
newr:= rad(n);
P[lastr, newr]:= n-1;
if assigned(P[newr, lastr]) then
count:= count+1; A[count]:= n-1; M[count]:= P[newr, lastr];
fi;
lastr:= newr;
od:
seq(A[n], n=1..count); # Robert Israel, Aug 10 2014
MATHEMATICA
(* Recomputation up to a(13), assuming m of the form 2^k+1 *)
rad[n_] := rad[n] = Select[Divisors[n], SquareFreeQ][[-1]];
okQ[n_] := Module[{r = rad[n], r1 = rad[n+1], k, m}, For[k = 0, k < Log[2, n-1], k++, m = 2^k+1; If[r == rad[m+1] && rad[m] == r1, Return[True]]]; False];
Reap[For[n = 1, n <= 5*10^6, n++, If[okQ[n], Print[n]; Sow[n]]]][[2, 1]] (* Jean-François Alcover, Apr 11 2019 *)
PROG
(PARI) lista(nn) = {v = vector(nn, i, rad(i)); for (n=1, nn-1, ok = 0; if (n % 2, ma = 2, ma = 1); forstep (m = ma, n-1, 2, if ((v[n] == v[m+1]) && (v[m] == v[n+1]), ok = 1; break); ); if (ok, print1(n, ", ")); ); } \\ Michel Marcus, Aug 10 2014
CROSSREFS
KEYWORD
nonn,more
AUTHOR
Naohiro Nomoto, Oct 29 2003
EXTENSIONS
More terms from David Wasserman, Jan 29 2004
a(13) confirmed by Robert Israel, Aug 10 2014
STATUS
approved