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.

Advertisements

Pages: 1 2

8 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

    @ max, maybe try Google

  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[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);
      }
    }
    
  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

    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"
    
    $ ./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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: