Perfect numbers

Question

The question given to me:

a. Find the perfect numbers between $1$ and $10^6$

b: Find the abundant numbers between $1$ and $1000$

---

For a, I wrote

Select[Range[1, 10^6], All]

I think I'm suppose to use divisor in there somewhere, but I'm not sure where to put that since I think me selecting all the numbers between $1$ and $10^6$ won't give me the list of all the perfect numbers...

For b, I think it's similar to part a, which I'm in a dead end with.

---

I'm trying to find a simple input since I'm still a beginner and want to stay with the basics.

Answer 1

Perfect numbers:

Select[ Range[10^6], Total[Divisors @ #] == 2 # &]

{6, 28, 496, 8128}

abundant numbers:

Select[ Range[10^3], Total[ Most @ Divisors @ #] > # &]//Short

{ 12, 18, 20, 24, 30, 36, 40, 42, 48, <<228>>, 
     968, 972, 978, 980, 984, 990, 992, 996, 1000} 

I used Short to to get only a few since there are:

Count[ Range[10^3], _?(Total[Most@Divisors@#] > # &)]

246

of them.

Edit

As RunnyKine pointed out that using DivisorSigma[1, #] & is more efficient than Total @ Divisors @ # &. Another improvement might be exploiting the fact that there are no known odd perfect numbers, it was verified that there is none below 10^1500. The largest known perfect number (48-th) has only 34850340 digits i.e. IntegerLength[2^(57885161 - 1) (2^57885161 - 1)]. Taking the above into account we get 2 times speed up with:

Pick[ #, DivisorSigma[1, #] - 2 #, 0]& @ Range[2, 10^6, 2]

However we can observe there are odd abundant numbers, but they are sparsely distributed among even ones. Below 1000 there is only one ( while there are 245 even ones):

Pick[ #, UnitStep[ DivisorSigma[1, #] - 2 # - 1], 1]& @ Range[1, 10^3, 2]

{945}

Below 10^6 there are

Length @ Pick[ #, UnitStep[ DivisorSigma[1, #] - 2 # - 1], 1] & /@ 
{Range[ 2, 10^6, 2], Range[ 1, 10^6, 2]}

{245549, 1996}

even and odd abundant numbers respectively.

It is remarkable that Length @ Pick[ Range[10^6], UnitStep[ DivisorSigma[1, #] - 2 # - 1], 1]] is faster than : Count[ Range[10^6], _?(DivisorSigma[1, #] > 2 # &)].

Answer 2

A faster approach to finding Perfect numbers using DivisorSigma

 Select[Range[10^6], DivisorSigma[1, #] == 2 # &]

{6, 28, 496, 8128}

Here's an even faster approach:

Pick[#, MapThread[Equal, {DivisorSigma[1, #], 2 #}], True] &[Range[10^6]]

and a little bit faster:

Pick[#, DivisorSigma[1, #] - 2 #, 0] &@Range[10^6]

For Abundant numbers do:

Select[Range[10^6], (DivisorSigma[1, #] - #) > # &]

and a faster approach as above:

Pick[#, MapThread[Greater, {DivisorSigma[1, #] - #, #}], True] &[Range[10^6]]

By giving Greater the Listable Attribute we can squeeze out some more performance:

SetAttributes[Greater, Listable]
Pick[#, DivisorSigma[1, #] > 2 #, True] &@Range[10^6]

Edit :

For this particular problem since we know that the largest Perfect number is less than 10000, we can begin to hit miliseconds regime using ParallelMap:

ParallelMap[Pick[#, DivisorSigma[1, #] - 2 #, 0] &, {Range[2, 5*^3, 2], 
    Range[5*^3 + 2, 1*^4, 2]}] // Flatten // AbsoluteTiming

{0.015627, {6, 28, 496, 8128}}

Of course this will also give us speed up if we scan the entire range.

Answer 3

Or how about the connection between even perfect numbers and Mersenne primes?

With[{p = Prime[Range[20]]},
     Pick[p, PrimeQ[2^p - 1]] /. q_ -> 2^(q - 1) (2^q - 1)]

Update for perfect numbers only.

As noted in previous answers, DivisorSigma[1,n] is faster than summing Divisors, and Pick is faster than Select. So even perfect numbers may be found by using

RepeatedTiming[Pick[#, DivisorSigma[1, #] - 2 #, 0] &@Range[2, 10^6, 2]]

which gives a timing of 1.34 s on my machine. These approaches are slow compared to the clever divisor sum by @MichaelE2. The repeated timing of his function is a tiny 0.056 s.

@QuantumDot points out the new v10.4 functions PerfectNumber[n] and PerfectNumberQ[n]. However, the following takes a glacial 4 seconds! Why?

RepeatedTiming[Pick[#, Map[PerfectNumberQ, #]] &@Range[2, 10^6, 2]]

The Help page for the new v10.4 function MersennePrimeExponent shows how to instantly calculate even perfect numbers.

RepeatedTiming[2^(# - 1)*(2^# - 1) &[MersennePrimeExponent[Range[4]]]]

Answer 4

DivisorSigma is no doubt fast on individual integers, but the problem of finding the sum of divisors over a large contiguous range of integers can be solved with a fast a C-style program, perfect for compiling if the sums stay within the size of a machine integer.

The idea is to turn the problem of summing the divisors around. An integer n is a divisor of each of its multiples. We keep a list of sums, and for every integer n, we add it to each sum whose position in the list is a multiple of n. We do this up to some pre-specified limit lim. At the end, the list contains the sums of the divisors of each integer up to lim.

divsum = Compile[{{lim, _Integer}},
   Module[{dsum},
    dsum = Table[1, {lim}];
    Do[
     dsum[[k]] += n,
     {n, 2, lim}, {k, n, lim, n}];
    dsum],
   RuntimeOptions -> "Speed", CompilationTarget -> "C"
   ];

divsum[10^6]; // AbsoluteTiming
(* {0.119278, Null} *)

DivisorSigma[1, Range[10^6]]; // AbsoluteTiming
(* {3.985526, Null} *)

We can define functions to find perfect and abundant numbers, using Pick to pick out the numbers according to the divisor-sum property.

perfect[lim_] := Pick[#, divsum[lim] - 2 #, 0] &@Range[lim];
abundant[lim_] := Pick[#, UnitStep[2 # - divsum[lim]], 0] &@Range[lim]

perfect[10^6] // AbsoluteTiming
(* {0.159563, {6, 28, 496, 8128}} *)

abundant[10^6] // Short // AbsoluteTiming
(* {0.168468, {12, 18, 20, 24, <<247537>>, 999992, 999996, 999999, 1000000}} *)

Answer 5

With the newly released version 10.4, you can use the built-in function PerfectNumber.

I don't have the latest version of 10.4; it would be nice to see how it compares with the answers given here...

---

but... based on my reading of the documentation page (under the section Details), I can only conclude one thing: PerfectNumber is bloatware. Arguments contrary to my conclusion are welcome in the comments section.