%I #35 Oct 25 2020 22:51:06
%S 2,18,250,6174,3660250,1542294,2839714,41154,117793122328750,
%T 7978057537338,2898701538750,33734898,29688151506250,21107677374,
%U 69834458642125879757481250,3999523458421521342
%N a(n) is the least natural k which is a multiple of prime(n) such that for some m >= 0, phi(k) = rad(k)^m, where phi(k) = A000010(k) and rad(k) = A007947(k).
%C The number m mentioned above is usually referred to as the order of the corresponding number a(n). The sequence of these orders is in A337776.
%C The algorithm suggested here for the calculation of a(n) starts its work from prime(n).
%C Numbers k such that phi(k) = rad(k)^m with m >= 1 are given in A211413. - _Andrew Howroyd_, Sep 21 2020
%D J.-M. De Koninck, Ces nombres qui nous fascinent, Entry 108, p. 38, Ellipses, Paris 2008.
%D J.-M. De Koninck & A. Mercier, 1001 Problèmes en Théorie Classique Des Nombres, Problème 745 ; pp 95; 317-8, Ellipses Paris 2004.
%H J.-M. De Koninck, <a href="https://www.jstor.org/stable/4145084">When the Totient Is the Product of the Squared Prime Divisors: Problem 10966</a>, Amer. Math. Monthly, 111 (2004), p. 536.
%e For n=12 the initial prime is prime(12) = 37 and a(12) = 33734898 because phi(33734898) = 10941048, rad(33734898) = 222 and 222^3 = 10941048 and there is no smaller number satisfying the requirements. The order of a(12) is 3.
%t nn = 16;
%t Sar = Table[0, {nn}]; Sar[[1]] = 2;
%t (*It is a list oh the sequence A337775*)
%t OrdSar = Table[0, {nn}]; OrdSar[[1]] = 0;
%t (*It is a sequence A337776 - the orders of members in sequence A337775*) For[Index = 2, Index <= nn, Index++,
%t InitialPrime = Prime[Index];
%t InitialInteger = InitialPrime - 1;
%t InitialArray = FactorInteger[InitialInteger];
%t For[i = 1, i <= Length[InitialArray], i++,
%t CurrentArray =
%t FactorInteger[InitialArray[[-i, 1]] - 1] ~Join~ InitialArray;
%t InitialInterger =
%t Product[CurrentArray[[k, 1]] ^ CurrentArray[[k, 2]], {k, 1,
%t Length[CurrentArray]}];
%t InitialArray = FactorInteger[InitialInterger];
%t ];
%t InitialArray = InitialArray ~Join~ {{InitialPrime, 0}};
%t Ord = Max[InitialArray[[All, 2]]];
%t Lint = Product[
%t Power[InitialArray[[k, 1]], Ord - InitialArray[[k, 2]] + 1], {k,
%t 1, Length[InitialArray]}];
%t radn = Product[InitialArray[[k, 1]], {k, 1, Length[InitialArray]}];
%t Sar[[Index]] = Lint;
%t OrdSar[[Index]] = Ord;
%t ];
%t Print["Sar= ", Sar]
%t Print["OrdSar= ", OrdSar]
%o (PARI) rad(n) = factorback(factorint(n)[, 1]);
%o isok(k) = my(phik=eulerphi(k), radk=rad(k), x=logint(phik, radk)); radk^x == phik;
%o a(n) = {my(p=prime(n), k=p); while (!isok(k), k+=p); k;} \\ _Michel Marcus_, Sep 23 2020
%Y Cf. A000010 (phi), A000040 (prime), A007947 (rad), A023503, A024619, A105261, A211413.
%K nonn
%O 1,1
%A _Vladislav Shubin_, Sep 20 2020