Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
add comment

4 Answers 4

up vote 7 down vote accepted

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 # &)].

share|improve this answer
    
Good stuff. I wonder if there's any odd perfect numbers. –  RunnyKine Oct 20 '13 at 4:05
    
@RunnyKine Many have asked this question for centuries. Check out MathWorld or oddperfect.org –  KennyColnago Oct 20 '13 at 17:26
    
I just started learning this information, but what exactly does @ do? –  asik Oct 20 '13 at 18:27
    
@asik f @ x is equivalent to f[x], another possible form is x//f. Read e.g. Special Ways to Input Expressions tutorial in the documentation. –  Artes Oct 20 '13 at 18:29
    
Thank you! Kept seeing it, but I wanted to make sure. –  asik Oct 21 '13 at 3:21
add comment

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.

share|improve this answer
    
A little faster yet: Pick[#, Unitize[DivisorSigma[1, #] - 2 #], 0] &@Range[10^6] –  Michael E2 Oct 20 '13 at 1:29
    
@MichaelE2. Nice, but almost the same time as the MapThread method on my machine. –  RunnyKine Oct 20 '13 at 1:29
    
On mine it's 4.7 (MapThread) to 4.2. –  Michael E2 Oct 20 '13 at 1:37
    
@MichaelE2, it's definitely a little bit faster. –  RunnyKine Oct 20 '13 at 1:40
    
@MichaelE2, You actually don't need Unitize and that makes it a tiny bit faster. –  RunnyKine Oct 20 '13 at 1:48
show 2 more comments

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)]
share|improve this answer
add comment

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}} *)
share|improve this answer
    
Do you know why the alternate line Do[dsum[[Range[n, lim, n]]] += n, {n, 2, lim}] is 3 times faster than Do[dsum[[k]] += n, {n, 2, lim}, {k, n, lim, n}] in an uncompiled version of divsum, and yet 7 times slower in a compiled version of divsum? –  KennyColnago Oct 30 '13 at 21:31
    
@KennyColnago The one is faster in the uncompiled version because Part is "vectorized", but I'm not sure why it's so much slower in the compiled version. Its compiled code is 75% longer than the other's, which may contribute some. It could be that the vectorized read/write (with ..[[list]]) is slower than single element read/write when compiled. I have to go, but perhaps if I think of something, I'll let you know. –  Michael E2 Oct 30 '13 at 22:15
    
@KennyColnago First, Do[r = Range[n, lim, n]; dsum[[r]] += n, {n, 2, lim}] is a little more efficient, with compiled code similar to my code, but still much slower. Examining the C code (CCodeGenerate) of each shows mine has rather simple, direct C array addressing -- probably impossible to beat for speed. In the other loop, library functions are called to deal with the list of parts, one to get the parts and another to set the parts. This is the main difference in the code, and I suppose must account for the difference in speed. –  Michael E2 Oct 31 '13 at 3:01
    
Wow, thanks for the sleuthing! It helps a lot with my never ending search for efficiency. –  KennyColnago Oct 31 '13 at 21:34
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.