Loren on the Art of MATLAB

When is a Number Perfect? 4

Posted by Loren Shure,

Answer: When the sum of all divisors of a positive integer, except the number itself, equals the number. Here's a reference with more information from Wikipedia.

Contents

Finding All the Divisors

You may think that all we need to do is find all the prime factors and then we can test the sum of those versus the number itself. Searching MATLAB documentation for prime factors yields the function factor.

help factor
 FACTOR Prime factors.
    FACTOR(N) returns a vector containing the prime factors of N.
 
    This function uses the simple sieve approach. It may require large
    memory allocation if the number given is too big. Technically it is
    possible to improve this algorithm, allocating less memory for most
    cases and resulting in a faster execution time. However, it will still
    have problems in the worst case, so we choose to impose an upper bound 
    on the input number and error out for n > 2^32. 
  
    Class support for input N:
       float: double, single
 
    See also PRIMES, ISPRIME.

    Overloaded methods:
       sym/factor

    Reference page in Help browser
       doc factor

Let's see if the number 6 is a perfect number (having divisors 1, 2, and 3).

n = 6;
facs = factor(n);
perfectOrNot = sum(facs) == n
perfectOrNot =
     0

Really? But 6 is a perfect number! What's going wrong? Let's first look at facs.

facs
facs =
     2     3

The number 1 is not a prime factor, but is a divisor and needs to be included in the calculation. We can either add it into the divisors or subtract it from n.

divisors = [1 factor(n)]
perfectOrNot = sum(divisors) == n
divisors =
     1     2     3
perfectOrNot =
     1

Try Some More Numbers

Make a function to do our comparisons now.

perfectOrNot = @(n) sum([1 factor(n)]) == n;

Check out some other numbers. First my favorite, 17.

n = 17;
perfectOrNot(n)
ans =
     0

Next, another known perfect number, 28.

n = 28;
perfectOrNot(n)
ans =
     0

What's Wrong Now?

Let's trace through the steps again.

divisors = [1 factor(n)]
sum28 = sum(divisors)
divisors =
     1     2     2     7
sum28 =
    12

Actually, these are not all the divisors of 28. Plus the number 2 is duplicated. We need to add all combinations of products of the prime factors, deduplicate the products, and omit the product of all of them (which is the number itself). For 28, that means including 2*2 and 2*7.

moreDivs = [2*2 2*7]
moreDivs =
     4    14

Combine prime factors with combinations, introduce the 1, and deduplicate.

perfect28 = sum(unique([divisors moreDivs])) == n;

Algorithm for Getting All Divisors

We next need an algorithm to calculate all the unique divisors of n. If you don't want to make up the algorithm yourself, you can check the File Exchange for entries that list divisors and not polynomials. You will find at least one candidate function there.

Or you could turn to the Symbolic Math Toolbox which includes some library functions for number theory - a likely place for us to look for some help. On the MuPAD side of the toolbox, we can take advantage of the MuPAD function numlib::divisors(). Since there is currently not a MATLAB version of this function available, I'm using a function from the MuPAD engine. To do that, I supply the library and name of the function, and the inputs.

symDivs = double(feval(symengine,'numlib::divisors', n))
symDivs =
     1     2     4     7    14    28

Since the divisors function returns n itself, we need to omit it from the sum calculation.

sum(symDivs(1:end-1)) == n
ans =
     1

So, 28 is perfect.

References and Links

Here are some links from the MATLAB newsgroup and pointers to some relevant contributions on the File Exchange.

Have You Had Some Fun with Number Theory in MATLAB?

I'd love to hear about your investigations using MATLAB for number theory. Tell me about them here.


Get the MATLAB code

Published with MATLAB® R2012b

4 CommentsOldest to Newest

Thanks, Loren. Fun blog today! You missed a function on the File Exchange that’s useful here, but that doesn’t require anything other than core MATLAB. I posted FACTOR2 back in 2003 to provide easy computation of all factors of x, including x itself:

factor2(28)
ans =
     1     2     4     7    14    28


So:


trimFactors = @(x) x(1:end-1);
isperfect = @(x) sum(trimFactors(factor2(x)))==x;
isperfect(5)
isperfect(6)
isperfect(28)
ans =
     0
ans =
     1
ans =
     1

Cheers!
Brett

Hi Loren,

A nice trick for finding the list of all divisors of an integer that is more efficient than the symengine solution is to use kron.

n = 10581480;

n has a lot of divisors, including 17, a very interesting number for some of us, and 1729, a number that others consider even more interesting. Including n itself, there are 384 distinct integer divisors. It takes a bit to compute the list of all divisors.

tic,double(feval(symengine,’numlib::divisors’, n));toc
Elapsed time is 0.024904 seconds.

Here is a simple solution in MATLAB.

tic
divs = 1;
for f = factor(n)
divs = unique(kron(divs,[1,f]));
end
toc
Elapsed time is 0.004309 seconds.

Interestingly, avoiding the multiple calls to unique is not hard, but it does not seem to save any time.

tic
f = factor(n);
k = [0,find(diff(f)),numel(f)];
divs = 1;
for i = 2:numel(k)
divs = kron(divs,f(k(i)).^(0:(k(i) – k(i-1))));
end
toc
Elapsed time is 0.005534 seconds.

Either case is faster than the symbolic solution though.

John

Thanks, John.

Very nice solution. I wonder if the unique once vs. every loop isn’t a gain because divs generally doesn’t grow that large, even when there are duplicates.

–Loren

Of course any method that uses factorization to find perfect numbers is about as efficient as sorting an array by generating all permutations and stopping when one is found that is in order.

These postings are the author's and don't necessarily represent the opinions of MathWorks.