Here we will implement an algorithm for computing the 'inverse' of Euler totient function ϕ(n). ϕ(n) is defined as the number of positive integers less than n and relatively prime to n. The inverse function for a given m will return the list of all the numbers n such that ϕ(n) = m.

How can all the solutions to the equation ϕ(n) = m be found? Naturally, let's start by looking at the decomposition of m and n into primes. Let

                                                                                n = (∏_(k = 1)^r) (Subscript[p, k]^((Subscript[a, k]))),  
      m = (∏_(k = 1)^s) (Subscript[q, k]^((Subscript[b, k]))), ϕ(n) = m

Euler function can be expressed explicitly via the prime factors of n:

                                                                              \(TraditionalForm\`ϕ(
          n) = \(∏\+\(k = 1\)\%r\( p\_k\%\(a\_k - 1\)\)(
              p\_k - 1) = \(m = ∏\+\(k = 1\)\%s q\_k\%\(b\_k\)\)\)\)                                                                   ( 1 )

It immediately follows from (1) that (Subscript[p, k]^((Subscript[a, k]) - 1)) and ((Subscript[p, k]) - 1) cannot contain any prime factors that are not in the set Q = {(Subscript[q, k])}_(k = 1)^s. Hence we can rewrite their decompositions thus:

                                                                                        \(TraditionalForm\`∀ k = \(\(1, r\)\&_\ \ a\_k = \(1 ∨ ∃ q\_i = p\_k\)\)\)
and
                                                                              \(TraditionalForm\`∀ k = \(\(1, r\)\&_\ \ p\_k = \(2 ∨ 
              p\_k = ∏\+\(i = 1\)\%\(t\_k\)q\_\(k, i\)\%\(b\_\(k, i\)\) + 1\)\)\)

where each (Subscript[q, k, i]) belongs to Q. Allowing the powers (Subscript[b, k, i]) to be zero, we can rewrite the last formula including all the (Subscript[q, k]) in the product:

                                                                                   \(TraditionalForm\`∀ k = \(\(1, r\)\&_\ \ p\_k = ∏\+\(i = 1\)\%s q\_i\%\(b\_\(k, i\)\) + 1\)\)

Another fact that follows from (1) is that (∏_1^r)((Subscript[p, k]) - 1) divides m. Therefore, for each i, k we have (Subscript[b, k, i]) ≤ (Subscript[b, i]). This means that we can find all the (Subscript[p, k]) that can possibly occur in n by simply trying all possible products of (Subscript[q, i]) raised to any power not greater than (Subscript[b, i]). And by transforming (1) we get

                                                                       \(TraditionalForm\`ϕ(
          n) = \(∏\+1\%r\(\( p\_k\%\(a\_k\)\)(p\_k - 1)\)\/p\_k = \(n \(∏\+1\%r\((1 - 
                    1\/p\_k)\)\) = \(\(m\)\(\ \)\(=>\)\)\)\)\)
                                                                    
                                                                                        => n = 
        m (∏_1^r)(((1 - (1/(Subscript[p, k]))))^(-1))                                                                                              ( 2 )

so we can get all the solutions n by substituting all possible combinations of (Subscript[p, k]) into (2). Not all n that we get this way will satisfy the equation ϕ(n) = m, but right now it is important to us that no solutions will be lost.

Another simple way to find (ϕ^(-1))(m) would be to estimate the range containing all the solutions. First, (((1 - (1/(Subscript[p, k]))))^(-1)) > 1, and second, if (Subscript[p, i]) < (Subscript[p, j]), then \(TraditionalForm\`\((1 - 1\/p\_i)\)\^\(-1\) > \((1 - 1\/p\_j)\)\^\(-1\)\). If the total number of (Subscript[p, k]) is R and they are sorted from smallest to largest: (Subscript[p, 1]) < (Subscript[p, 2]) <  ...  < (Subscript[p, R]), then 

                                                                                 ((m(1 - (1/(Subscript[p, R]))))^(-1)) ≤ 
        n ≤ 
        m (∏_(k = 1)^R)(((1 - (1/(Subscript[p, k]))))^(-1))

Now we need to estimate the values of (Subscript[p, k]). The largest is not greater than (m + 1) because each ((Subscript[p, k]) - 1) divides m (see above). And if we denote the kth prime in the row of natural numbers as (Subscript[π, k]), then (Subscript[p, k]) ≥ (Subscript[π, k]). Additionally, for all k this inequality holds: (Subscript[π, k]) > k ln(k). For k equal to one the right side is zero, so we use this inequality for k ≥ 2 (and explicitly substitute (Subscript[π, 1]) = 2), and get

                                                   \(TraditionalForm\`m \(∏\+\(k = 1\)\%R\((1 - 1\/p\_k)\)\^\(-1\)\) ≤ 
        m \(∏\+\(k = 1\)\%R\((1 - 1\/π\_k)\)\^\(-1\)\) ≤ 
        2  m \(∏\+\(k = 2\)\%R\((1 - 1\/\(k\ \(ln(k)\)\))\)\^\(-1\)\)\)  

The constant R is not greater than χ(m), the number of divisors of m (again, because m is divisible by each ((Subscript[p, k]) - 1)), and finally we get

                                                                                \(TraditionalForm\`m + 1 ≤ n ≤ 
        2  m \(∏\+\(k = 2\)\%\(χ(m)\)\((1 - 1\/\(k\ \(ln(k)\)\))\)\^\(-1\)\)\)

invphidemo1[m_] := 
      Module[
{Lp, Lq, Lb, Lpow, Lfact, 
          n}, 
{Lq, Lb} = 
          Transpose[FactorInteger[m]]; 
Lpow = 
          MapThread[
Table[#1^i, {i, 
                  0, #2}] &, 
{Lq, 
              Lb}]; 
Lp = 
          Outer[Times, Sequence @@ Lpow]; 
Lp = 
          Select[Flatten[Lp] + 1, PrimeQ]; 
Lfact = 
          Table[
(Lp[[i]]/(Lp[[i]] - 1))^
              j, 
{i, Length[Lp]}, {j, 0, 
              1}]; 
n = 
          m Outer[Times, Sequence @@ Lfact]; 
Select[
          Flatten[n], 
EulerPhi[#] == 
              m &]
]

invphidemo2[m_] := 
      Select[
Range[m + 1, 
          2  m (∏_(k = 2)^(DivisorSigma[0, m]))(((1 - (1/(k Log[k]))))^(-1))], 
EulerPhi[#] == 
            m &]

invphidemo1[12]
invphidemo2[12]

Of course, these are bad solutions. They are very slow, and the first one requires prodigious amounts of memory. The reason is that if the powers (Subscript[b, k]) are large, the set of (Subscript[p, k]) can be many times larger than the set Q, and creating the table of all the possible products of (Subscript[p, k]) is nearly impossible. The function invphi still will involve the search for all combinations of (Subscript[p, k]) that yield solutions but our goal is to make this search more efficient.

Let's denote the set of all the (Subscript[p, k]) that can possibly occur in a solution as P = {(Subscript[p, k])}_(k = 1)^R. We want to create list (Subscript[L, r]) = {(Subscript[p, (Subscript[i, 1])]) , (Subscript[p, (Subscript[i, 2])]) ,  ... , 
          (Subscript[p, (Subscript[i, r])])}, with elements from P, corresponding to one of the solutions n. This means that we'll get a solution by applying formula (2) to (Subscript[L, r]). Note that (Subscript[L, r]) denotes the list with r elements and not the rth element. This notation will be handy for us because r is different for each n and is not known in advance. Let on some step in the algorithm we already have t elements (Subscript[L, t]). Then the problem is reduced to deciding adding which (Subscript[p, k]) to (Subscript[L, t]) can yield a valid solution, and which (Subscript[p, k]) needn't be tried, being 'incompatible' with those (Subscript[p, k]) that are in the list already.

First, as was already shown, the product (∏_(k = 1)^r)((Subscript[p, (Subscript[i, k])]) - 1) divides m, so their quotient

                                                                                              (Subscript[d, r]) = 
        m (∏_(k = 1)^r)((((Subscript[p, (Subscript[i, k])]) - 1))^(-1))

is integer. Therefore, before adding (Subscript[p, k]) to (Subscript[L, t]) we should check that the resulting quotient (Subscript[d, t + 1]) is integer, and if it's not, (Subscript[p, k]) is not added.

Second, (Subscript[d, r]) contains only prime factors from Q (because (Subscript[d, r]) is m divided by something). But by formula (1), 

                                                                                                    \(TraditionalForm\`d\_r = ∏\+\(k = 1\)\%r p\_\(i\_k\)\%\(a\_\(i\_k\) - 1\)\)

It follows that only those primes can occur in (Subscript[d, r]) which belong to P ⋂ Q and are also members of (Subscript[L, r]). So on the step where we have t elements we should check if the 'current' quotient (Subscript[d, t]) contains a factor (Subscript[q, k]) which is not a member of (Subscript[L, t]). If such (Subscript[q, k]) exists, then it cannot appear in the 'final' quotient (Subscript[d, r]) unless it belongs to (Subscript[L, r]). Therefore, we can add (Subscript[q, k]) to (Subscript[L, t]) (possible only if (Subscript[q, k]) belongs to P) or eliminate (Subscript[q, k]) from the quotient on the next steps, if we want to arrive to a solution. The elimination can be done by selecting (Subscript[p, (Subscript[i, t + 1])]) such that ((Subscript[p, (Subscript[i, t + 1])]) - 1) is divisible by (Subscript[q, k]), so in the next quotient (Subscript[d, t + 1]) = 
        (Subscript[d, t])/((Subscript[p, (Subscript[i, t + 1])]) - 1) the power to which (Subscript[q, k]) appears in it will be lowered (though (Subscript[q, k]) may still be present in (Subscript[d, t + 1])). How this all benefits the algorithm is that instead of trying to add all the possible (Subscript[p, k]) that are not in the list yet, we can find out that there is a much smaller subset of P such that at least one of its elements should necessarily be appended to (Subscript[L, t]), in order to eliminate a (Subscript[q, k]). Then we may try building new lists from (Subscript[L, t]) by adding only (Subscript[p, k]) from that subset first (adding them separately, thus creating several new 'candidates' from (Subscript[L, t])). This way we can significantly reduce the amount of search without losing any of the solutions.

Another important consequence is that if want to know whether (Subscript[L, t]) corresponds to a valid solution or not, it can be checked without computing EulerPhi. The criterion of its validity is exactly the condition discussed above, namely that (Subscript[d, t]) can contain only those (Subscript[q, k]) that appear in (Subscript[L, t]).

The program invphi first finds P. P and Q are rearranged so that all their common elements, the number of which is denoted (Subscript[r, 0]), come in the beginning and have the same positions in both lists. The main structure wrk is a list of solution candidates. Each element of wrk has in turn two elements. The first one defines (Subscript[L, t]). It contains an indicator for each (Subscript[p, k]) from P showing if that (Subscript[p, k]) is already in (Subscript[L, t]), can, or can not be added to it (1, 0, -1 respectively). The second element is (Subscript[d, t]). On each step we check if there are any factors (Subscript[q, k]) in the quotient (Subscript[d, t]) that should be eliminated. If there are, then we see how many (Subscript[p, i]) can reduce the power to which (Subscript[q, k]) occurs in (Subscript[d, t]), ie how many (Subscript[p, i]) have the property that ((Subscript[p, i]) - 1) is divisible by (Subscript[q, k]). (As discussed above, another possibility is (Subscript[p, i]) = (Subscript[q, k])). We always select the (Subscript[q, k]) for which the number of those (Subscript[p, i]) is the smallest, thus trying to minimize the number of new candidates that will be added. If no (Subscript[q, k]) to be eliminated are found, then all members of P that are allowed for the current (Subscript[L, t]) are added to it. For each new element of wrk thus created we find the new quotient (Subscript[d, t + 1]), and mark as 'incompatible' those (Subscript[p, k]) that would make (Subscript[d, t + 2]) fractional. Also, if we are adding several (Subscript[p, k]) to (Subscript[L, t]), say (Subscript[p, 1]) and (Subscript[p, 2]), creating two new lists, (Subscript[L, t + 1]^((1))) = Append[(Subscript[L, t]), (Subscript[p, 1])] and (Subscript[L, t + 1]^((2))) = Append[(Subscript[L, t]), (Subscript[p, 2])], then in (Subscript[L, t + 1]^((2))) (Subscript[p, 1]) is marked as 'incompatible' too in order to avoid possible repetitions of the same sets of (Subscript[p, k]) later. When nothing further can be added to wrk, the final list of solutions is generated by selecting those elements from wrk that yield correct values of n. This is done in a similar way; we consider the corresponding quotient, and if no (Subscript[q, k]) need to be eliminated, then we have a solution.

There are some other ways to improve the performance. It can be easily noticed that if we add (Subscript[p, k]) = 2 to the list, then the quotient doesn't change on the next step. We can speed up the search by never considering (Subscript[p, k]) = 2, and only taking into consideration the possible presence of this factor later, when generating solutions from the lists of (Subscript[p, k]). Also, for selecting the (Subscript[p, i]) that can eliminate a given (Subscript[q, k]), the structure Mdiv is created in advance, in which Mdiv[[k]] is the list of indices i such that (Subscript[q, k]) divides ((Subscript[p, i]) - 1). Finally, after we have created new lists from (Subscript[L, t]), we'll never return to (Subscript[L, t]) again. When the structure wrk grows large, we can convert those elements to the final result n and remove them from wrk. This not only reduces the memory usage, but makes large computations significantly faster, because operations will be performed with smaller lists.

In Mathematica, ϕ(0) = 0, ϕ(1) = 1. invphi complies with this conventions. Also note that only the list of nonnegative solutions is returned, while in Mathematica ϕ(-n) = ϕ(n).

invphi[m_Integer] := 
      Module[
{main, init, genp, bestcand, gencand, 
          addcand, genans, 
Lp, Lq, r, s, r0, 
          Mdiv}, 

main[] := 
          Module[
{ans = {}, wrk, 
              threshold = 100, 
Lstate, Ladd, quo, indx, 
              i}, 
wrk = {{Table[0, {r}], 
                  m}}; 
wrk[[1, 1, 
                  1]] = -1; 
For[i = 1, 
              i ≤ Length[wrk], i++, 
If[
                i == 
                  threshold + 1, 
ans = {ans, 
                    genans[Take[wrk, 
                        threshold] ] }; 
wrk = 
                  Drop[wrk, threshold]; 
i = 
                  1
]; 
{Lstate, quo} = 
                wrk[[i]]; 
indx = 
                bestcand[Lstate, quo]; 
Ladd = 
                gencand[Lstate, indx]; 
wrk = 
                Join[wrk, 
                  addcand[Lstate, quo, 
                    Ladd] ]
]; 
Flatten[{ans, 
                genans[wrk]}]
]; 

init[] := 
          Module[
{Lb, Lpq}, 
{Lq, Lb} = 
              Transpose[FactorInteger[m]]; 
genp[
              Lb]; 
Lpq = 
              Intersection[Lp, Lq]; 
{Lp, 
                Lq} = Join[
                    Lpq, 
Complement[#, 
                      Lpq] ] & /@ 
{Lp, 
                  Lq}; 
{r, s, r0} = 
              Length /@ {Lp, Lq, 
                  Lpq}; 
Mdiv = Cases[
Range[r], 
x_ /; 
                      Mod[Lp[[x]] - 1, #] == 
                        0
] & /@ 
                Lq;
]; 

genp[Lb_] := 
          Module[
{Lpow, tmp}, 
Lpow = 
              MapThread[
Table[#1^i, {i, 
                      0, #2}] &, 
{Lq, 
                  Lb}]; 
Lp = {}; 
Outer[
If[
                  PrimeQ[tmp = Times[##] + 1], 
Lp = {Lp, 
                      
                      tmp}
] &, 
Sequence @@ Lpow]; 
Lp = 
              Flatten[Lp];
]; 

bestcand[Lstate_, quo_] := 
          Module[
{len = Infinity, indx = 0, cur, 
              i}, 
For[i = 1, 
              i ≤ 
                s, i++, 
If[((i ≤ 
                            r0 && 
Lstate[[i]] ≠ 1) || 
i > r0) && 
Mod[
                      quo, Lq[[i]] ] == 
                    0, 
cur = 
                  Length[Mdiv[[i]] ]; 
If[
                  cur < len, 
len = 
                    cur; 
indx = 
                    i
]
]
]; 
indx
]; 

gencand[Lstate_, indx_] := 
          Module[
{Ladd}, 
Ladd = 
If[
                indx ≠ 0, 
If[
                  indx ≤ r0, 
Prepend[
                    Mdiv[[indx]], 
                    indx], 
Mdiv[[indx]]
], 
Range[
                  r]
]; 
Select[Ladd, 
              Lstate[[#]] == 
                  0 &]
]; 

addcand[Lstate_, quo_, Ladd_] := 
          Module[
{ans = {}, Lstate2, quo2, len, 
              i}, 
len = 
              Length[Ladd]; 
For[i = 1, 
              i ≤ len, i++, 
Lstate2 = 
                ReplacePart[Lstate, 1, 
                  Ladd[[i]] ]; 
quo2 = 
                quo/(Lp[[Ladd[[i]] ]] - 
                      1); 
(Lstate2[[Ladd[[#]] ]] = -1) & /@ 
Range[
                  i - 1]; 
If[
                    Lstate2[[#]] == 
                        0 && 
Mod[quo2, 
                          Lp[[#]] - 1] ≠ 
                        0, 
Lstate2[[#]] = -1] & /@ 
Range[r]; 
AppendTo[
                ans, {Lstate2, 
                  quo2}]
]; 
ans
]; 

genans[L_] := 
          Module[
{ans = {}, Lstate, quo, res, add2, i, 
              j}, 
For[i = 1, 
              i ≤ 
                Length[L], i++, 
{Lstate, quo} = 
                L[[i]]; 
For[add2 = 0, 
                add2 ≤ 1, add2++, 
If[
                  add2 == 1, 
Lstate[[1]] = 
                    1
]; 
For[j = 1, 
                  j ≤ 
                    s, j++, 
If[((j ≤ 
                                r0 && 
Lstate[[j]] ≠ 1) || 
j > r0) && 
Mod[
                          quo, Lq[[j]] ] == 
                        0, 
Break[]
]
]; 
If[
                  j ≠ 
                    s + 1, 
Continue[]
]; 
res = 
                  Cases[
Transpose[{Lp, 
                        Lstate}], 
{x_, 1} -> 
                      x]; 
res = 
                  m Times @@ res/
                      Times @@ (res - 1); 
ans = {ans, 
                    res}
]
]; 
ans
]; 

Switch[m, 
0, 
          Return[{0}], 
1, 
          Return[{1, 
              2}], 
_?(OddQ[#] || Negative[#] &), 
          Return[{}]
]; 
init[]; 
main[]
]

Here finding all 847 solutions takes only a few seconds.

Timing[Length[invphi[(2^5) (3^5) (5^5)] ] ]

In this case, factoring each of the resulting integers would take about 150 seconds, so checking the answer takes many times longer than finding it.

Timing[
      invphi[(2^3) (31^3) (313^3) (619^3) (1511^3) 
          (1733^3)] ]

Copyright © 1999, Maxim Rytin