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.
Programming Praxis 
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.
what’s a totient?
@ max, maybe try Google
@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.
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[1],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); } }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.
Or maybe not, if the comment in that Wikipedia page is to be believed.
A brute force Haskell version.
import Data.Maybe (mapMaybe) import Math.NumberTheory.ArithmeticFunctions (totient) import System.Environment (getArgs, getProgName) import Text.Read (readMaybe) -- 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"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] = 0; array[1] = 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[1], NULL, 10); ulong* array = malloc((n + 1) * sizeof(array[0])); perfect_totients(array, n); free(array); return EXIT_SUCCESS; }Example: