## Perfect Totient Numbers

### January 8, 2019

We’ve seen several methods for calculating totients in previous exercises. Here we calculate the totient using Euler’s product formula, factoring n using a 2,3,5-wheel; our totient function takes time O(sqrt n, and we memoize the result because of all the iterated sub-totients:

```(define-memoized (totient n)
(let ((wheel (vector 1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((n n) (t n) (f 2) (w 0))
(if (< n (* f f))
(if (< 1 n) (- t (/ t n)) t)
(if (zero? (modulo n f))
(begin (while (zero? (modulo n f))
(set! n (/ n f)))
(loop n (- t (/ t f))
(+ f (vector-ref wheel w))
(if (= w 10) 3 (+ w 1))))
(loop n t (+ f (vector-ref wheel w))
(if (= w 10) 3 (+ w 1))))))))```

Now it is easy to compute the iterated totient sum:

```(define (iterated-totient-sum n)
(let loop ((n n) (s 0))
(if (= n 1) s
(let ((t (totient n)))
(loop t (+ s t))))))```

And here is the requested result:

```> (filter (lambda (n) (= (iterated-totient-sum n) n)) (range 1 10000))
(3 9 15 27 39 81 111 183 243 255 327 363 471 729 2187 2199
3063 4359 4375 5571 6561 8751)```

You can run the program at https://ideone.com/TzNR0v.  A reasonable alternative to the method shown above is to calculate the totients by sieving, as in a previous exercise.

Pages: 1 2

### 9 Responses to “Perfect Totient Numbers”

1. matthew said

Sieving is nice, and we can roll calculating the totient sums into the sieving loop:

```def perfecttotient(N):
a = list(range(N))
for p in range(2,N):
if a[p] == p:
for q in range(p,N,p):
a[q] = a[q]//p*(p-1);
a[p] += a[a[p]]
if a[p]-1 == p: yield p

for n in perfecttotient(10000000): print(n)
```

gets up to 9056583 in 20 seconds or so.

2. max said

what’s a totient?

3. Paul said

4. programmingpraxis said

@max: The totatives of a positive integer n are those positive integers less than n that are coprime to n; that is positive integers m for which gcd(m, n) = 1. For instance, the totatives of 30 are 1, 7, 11, 13, 17, 19, 23 and 29. The totient of a positive integer n is the number of totatives of n, so φ(30) = 8.

5. matthew said

And this similar C++ program gets up to 918330183 in a couple of minutes. After that memory gets a bit tight on my laptop.

```#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <vector>

typedef uint32_t T;

int main(int argc, char *argv[]) {
T N = 10000;
if (argc > 1) N = strtoul(argv,0,0);
std::vector<T> a(N);
for (T i = 0; i < N; i++) a[i] = i;
for (T p = 2; p < N; p++) {
if (a[p] == p) {
for (T q = p; q < N; q += p) {
a[q] = a[q]/p*(p-1);
}
}
a[p] += a[a[p]];
if (a[p]-1 == p) printf("%u\n", p);
}
}
```
6. matthew said

One useful place the totient function appears is in Euler’s Theorem: https://en.wikipedia.org/wiki/Euler%27s_theorem, which is fundamental to the correctness of RSA encryption.

7. matthew said

Or maybe not, if the comment in that Wikipedia page is to be believed.

8. Globules said

```import Data.Maybe (mapMaybe)
import Math.NumberTheory.ArithmeticFunctions (totient)
import System.Environment (getArgs, getProgName)

-- Iterated totient.
iterTot :: Int -> Int
iterTot = foldr (+) 1 . takeWhile (>1) . iterate totient . totient

-- A list of perfect iterated totients.
perfIterTots :: Int -> [Int]
perfIterTots = eqIdx . map iterTot . enumFromTo 1
where eqIdx = map snd . filter (uncurry (==)) . zip [1..]

main :: IO ()
main = do
prog <- getProgName
ns <- mapMaybe readMaybe <\$> getArgs
case ns of
[n] -> mapM_ print \$ perfIterTots n
_   -> putStrLn \$ "Usage: " ++ prog ++ " n"
```
```\$ ./perfecttotient 9999 | fmt
1 3 9 15 27 39 81 111 183 243 255 327 363 471 729 2187 2199 3063
4359 4375 5571 6561 8751
```
9. Daniel said

Here’s a solution in C.

```#include <stdio.h>
#include <stdlib.h>

typedef unsigned long ulong;

void perfect_totients(ulong* array, ulong n) {
array = 0;
array = 0;
for (ulong i = 2; i <= n; ++i) {
array[i] = i;
}
for (ulong i = 2; i <= n; ++i) {
if (array[i] == i) {
for (ulong j = 1;; ++j) {
if (j * i > n) break;
array[j * i] *= i - 1;
array[j * i] /= i;
}
}
array[i] += array[array[i]];
if (array[i] == i) printf("%lu\n", i);
}
}

int main(int argc, char* argv[]) {
if (argc != 2) return EXIT_FAILURE;
ulong n = strtoul(argv, NULL, 10);
ulong* array = malloc((n + 1) * sizeof(array));
perfect_totients(array, n);
free(array);
return EXIT_SUCCESS;
}
```

Example:

```\$ ./a.out 10000
3
9
15
27
39
81
111
183
243
255
327
363
471
729
2187
2199
3063
4359
4375
5571
6561
8751
```