Write a script to find "Perfect Numbers"

Inspired by this tweet:

Write a script to find "Perfect Numbers." Here's the list of numbers up to 19 digits 😈:

6, 28, 496, 8138, 33550336, 8589869056, 137438691328, 2305843008139952128

Look forward to seeing everyone's solutions!

Comments

[Meghan (she/her)]

const perfect = [6, 28, 496, 8128, 33550336, 8589869056, 137438691328, 
2305843008139952128];

perfect.map(v => v.toString(2));
/*
[ "110",
  "11100",
  "111110000",
  "1111111000000",
  "1111111111111000000000000",
  "111111111111111110000000000000000",
  "1111111111111111111000000000000000000",
  "1111111111111111111111111111111000000000000000000000000000000" ] */

perfect.map(v => v.toString(3));
/*
[ "20",
  "1001",
  "200101",
  "102011001",
  "2100010112102001",
  "211011122100120010101",
  "111010202100010220210001",
  "120100210022221112202020222210100212001" ] */

perfect.map(v => v.toString(4));
/*
[ "12",
  "130",
  "13300",
  "1333000",
  "1333333000000",
  "13333333300000000",
  "1333333333000000000",
  "1333333333333333000000000000000" ] */

I find base 2 and 4 particularly interesting because for each perfect number, the amount of 1s/3s and 0s is the same.

edit: fixed my sequence because 8138 is not perfect, its 8128.

[Shi Ling]

Yeah, I had a suspicion that there's a geometric solution to finding perfect numbers and was trying to formulate it in my head. Putting it in base2 makes it quite clear. Hm... 🤔

[Kasey Speakman]

Another interesting observation. In the binary form of each perfect number, the count of leading 1s is a prime number. In fact, they are all Mersenne Exponents.

[Itachi Uchiha]

I wrote this script 3 years ago.

def perfect_number(num):
    su = 0
    for i in range(1, num):
        if num % i == 0:
            su += i
    return su == num

def perfect_number_main(n):
    if perfect_number(n) == True:
        print("%s is perfect" % n)
    else:
        pass

for i in range(1,10000):
    perfect_number_main(i)

[chenge]

Translate yours in Ruby.

def sum_divisions(n)
  (1...n).select{|i| n % i == 0}.sum
end

def perfect_number?(n)
  n == sum_divisions(n)
end

(1..10_000).each do |i|
  puts "#{i} is perfect" if perfect_number?(i)
end

[chenge]

In JS, not so good cause I'm new to JS. If any good writing comment please.

const sum = (accumulator, currentValue) => accumulator + currentValue;

function range(size, startAt = 0) {
  return [...Array(size).keys()].map(i => i + startAt);
}

function sumDivisions(n){
  arr = range(n-1, startAt=1).filter(x => n % x == 0)
  return arr.reduce(sum, 0);
}

function is_perfect_number(n){
  return sumDivisions(n) == n
}

range(10000, 1).forEach(function(x){
  if(is_perfect_number(x)){
    console.log(x);
  }
})

[chenge]

Do you have a better version for Go?

package main

import (
  "fmt"
)

func sum_divisions(n int) int {
  sum := 0
  for i := 1; i < n; i++ {
    if n % i == 0 {
      sum += i
    }
  }
  return sum
}

func is_perfect(n int) bool {
  return n == sum_divisions(n)
}

func main() {
  for i := 1; i <= 10000; i++ {
    if is_perfect(i) {
      fmt.Println(i)
    }
  }
}

[kaelscion]

How efficient do you remember this article being?

[Kasey Speakman]

F#

let getDivisors n =
    let limit = int64 (sqrt (float n))
    [
        for i = 1L to limit do
            if n % i = 0L then
                yield i
                let result = n / i
                if result <> i && result <> n then
                    yield result
    ]


let isPerfect n =
    n = (n |> getDivisors |> List.sum)

Usage

(*
   per @nektro's observation on perfect numbers in binary
   * perfect numbers have an odd number of bits
   * the first half (rounded up) are 1s
   * the second half (rounded down) are 0s
*)

// Creates a perfect number candidate of the given bit length.
// Not all candidates are perfect numbers.
// Example: createCandidate 3 returns 110b or 6
//          createCandidate 5 returns 11100b or 28
let createCandidate bitLength =
    let zeroStart = bitLength / 2 + 1
    let update num i =
        let shifted = num * 2L
        if i < zeroStart then shifted + 1L
        else                  shifted
    List.init bitLength id |> List.fold update 0L

let oddBitLengths = // 1, 3, 5, 7, 9, ... , 61
    Seq.init 31 (fun i -> i * 2 + 1)
    // should be 32, but last one overflows on sum step

oddBitLengths
|> Seq.map createCandidate
|> Seq.filter isPerfect
|> Seq.iter (fun i ->
    printfn "%20i %s" i (Convert.ToString(i, 2))
)
//                   1 1
//                   6 110
//                  28 11100
//                 496 111110000
//                8128 1111111000000
//            33550336 1111111111111000000000000
//          8589869056 111111111111111110000000000000000
//        137438691328 1111111111111111111000000000000000000
// 2305843008139952128 1111111111111111111111111111111000000000000000000000000000000

To filter the numbers down to just likely candidates, I am using @nektro 's observation. I first generate a limited list of possible candidates that fit in a 64 bit integer, then filter down to only perfect ones. It takes 40s or so to run, but all of that time is spent between the last and next-to-last result.

[Guillaume Martigny]

How much time did this script took to run ? You approach is really nice, but I'm worried by |> Seq.filter isPerfect perf.

[Kasey Speakman]

Runtime is about 40 seconds in Release mode. The largest cost is finding divisors. I used the least expensive algorithm I could find. But it is still expensive. So to get the runtime down that low, I optimize down the number of inputs using nektro’s observation.

[Guillaume Martigny]

You could gain a little by only generating prime numbers for your candidates. Eratosthenes' sieve is not that hard to implement.

[Kasey Speakman]

I will post another solution using primes. Thanks for the tip!

[Kasey Speakman]

I posted the prime-based solution here. The actual perfect number calculation time is now negligible, and the code can be a one-liner. But it requires being fed Mersenne Exponents. Fortunately, very few/small Mersenne Exponents quickly turn into perfect numbers which are too large for traditional memory structures.

[Guillaume Martigny]

First dumb "look at every integer" solution:

const range = n => (new Array(+n)).fill().map((_, i) => i + 1);
return range(limit).filter((n) => {
    return range(Math.ceil(n ** 0.5))
        .filter(d => n % d === 0)
        .reduce((acc, val) => acc + val + (n / val), 0) / 2 === n;
});

This algo has trouble going further than 1e6.

Then, I dig your hypothesis that all Perfect number are in the form (n+1)"1"..(n)"0" in base 2. I needed to find the logical sequence of valid n.

[6, 28, 496, 8128, 33550336, 8589869056, 137438691328].map(n => n.toString(2).match(/0/g).length);
// => [1, 2, 4, 6, 12, 16, 18]

So I looked up this sequence on EOIS and luckily it found a match: "Numbers n such that 2ⁿ⁺¹ - 1 is prime".

Which lead to this way more powerful script:

const range = n => (new Array(+n)).fill().map((_, i) => i + 1);
const isPrime = (n) => {
    if (n <= 1) return false;
    if (n <= 3) return true;
    if (n % 2 === 0 || n % 3 === 0) return false;
    for (let i = 5, l = n ** 0.5; i <= l; i += 6) {
        if (n % i === 0 || n % (i + 2) === 0) return false;
    }
    return true;
};
return range(limit).filter(n => isPrime(2 ** (n + 1) - 1)).map((n) => {
    return parseInt("1".repeat(n + 1) + "0".repeat(n), 2);
});

This time, it can't go any higher than 1e3 before crashing, but yield more results. (However, JS precision let me down on the last result)

[Shi Ling]

Came to a same thesis and worked out a similar solution. Also ran into the limits at the 9th perfect number. But JavaScript now supports native BigInt (only in Chrome and I believe newer versions of NodeJs though!). Now we just need a way to generate large prime numbers fast to generate perfect numbers. 😉

[Guillaume Martigny]

Thanks to your input on BigInt, I made a new solution.
Thanks a lot !

[Shi Ling]

👏👏👏👏👏👏

I feel a heart is not enough. I'm so happy to see the new version.

[Guillaume Martigny]

Yeah, I liked that we came with two solution from the same hypothesis.

I'm looking into BigInt (thanks a lot I didn't know !) and it looks promising.

[Modaf]

Inspired from programming-algorithms.net/article...
Takes less than 1sec to get to 137438691328 but can't find easily 2305843008139952128

from math import sqrt
def perfect(n) :
    s = 1
    rac = int(sqrt(n))
    for k in range(2, rac+2) :
        if n%k == 0 :
            s += k + n//k
    return s == n
def f(x) :
    return (2**(x-1))*((2**x)-1)
n = 137438691329
m = 1
res = []
while True :
    val = f(m)
    print(val)
    if perfect(val) :
        res.append(val)
    m+=1
    if val > n :
        break
print(res)

[chenge]

site is great for algo, thanks.

[George]

Here's my solution in C

#include <stdio.h>

void  main(){
  int n, i, sum;
  int mn, mx;

  printf("Starting range: ");
  scanf("%d",&mn);

  printf("Ending range: ");
  scanf("%d",&mx);

  printf("Output: ");
  for(n = mn; n <= mx ; n++){
    i = 1;
    sum = 0;
    while(i < n){
      if(n % i == 0)
           sum = sum + i;
          i++;
    }
    if(sum== n)
      printf("%d ",n);
  }
      printf("\n");
}

[Andrei Gatej]

I think you can reduce the number of iterations by using “ while (i <=n/2)”.

[Stefan Rendevski]

You can reduce it even further by using "while (i <= sqrt(n))" and noticing that divisors come in pairs, for example: if 16 / 2 = 8, both 2 and 8 are divisors of 16. There is no need to go up to 8 and check all integers up until that point, you only need to perform the division.

[Florimond Manca]

A (lazy but slow) Python solution which hasn't even yielded 33550336 in 5 minutes, despite a subtle optimization allowing to find divisors in O(√n)… Probably could use some optimization to prevent trying all possible integers. :-)

def find_divisors(n: int):
    divisors = []
    m = 1
    while m * m < n:
        if n % m == 0:
            yield m
            if m != 1:
                yield n / m
        m += 1


def is_perfect(n: int) -> bool:
    return sum(find_divisors(n)) == n


if __name__ == '__main__':
    from itertools import cycle
    all_integers = cycle()
    for n in filter(is_perfect, all_integers):
        print(n)

[Philip Hallstrom]

Working on getting better at Clojure so thought I'd give it a shot.

(ns philip.perfect-numbers)

(defn brute-force-perfect-number? [num]
  (== num (apply + (filter #(zero? (mod num %)) (range 1 (inc (/ num 2)))))))

(time (println (take 4 (filter brute-force-perfect-number? (iterate inc' 2)))))
; (6 28 496 8128)
; "Elapsed time: 2649.69715 msecs"

(defn perfect-number? [num]
  (loop [divisor 2
         subtotal 1]
    (let [
          q (quot num divisor)
          evenly-divisible? (zero? (mod num divisor))]
    (if (>= divisor q)
      (== num subtotal)
      (recur (inc divisor) (if evenly-divisible? (+ subtotal divisor q) subtotal))))))

(time (println (take 4 (filter perfect-number? (iterate inc' 2)))))
; (6 28 496 8128)
; "Elapsed time: 66.942315 msecs"

[Kasey Speakman]

F#, ridiculously short/fast solution

// perfect num calculation is 2^(p-1)*(2^p - 1)
// source: https://oeis.org/A000396
let getPerfect p =
    pown 2L (p - 1) * (pown 2L p - 1L)

Usage

// prime calculation is out of scope for this problem
// (and this few can be looked up easily)
// so I simply list them here.
let mersenneExponents =
    [1; 2; 3; 5; 7; 13; 17; 19; 31]

List.map getPerfect mersenneExponents
|> printfn "%A"
// [1L; 6L; 28L; 496L; 8128L; 33550336L; 8589869056L; 137438691328L; 2305843008139952128L]

This runs in 0.35ms or ~3500 ticks in Release mode (after a JIT prerun). It basically turns the problem around to be just a matter of inputting Mersenne Exponents. Not many are needed before overflowing traditional memory structures like int64.