StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<blockquote> <p>The next one should be 33550336.</p> </blockquote> <p>Your code (I fixed the indentation so that it does in principle what you want):</p> <pre><code>i = 2 a = open("perfect.txt", 'w') a.close() while True: sum = 0 for x in range(1, i+1): if i%x == 0: sum += x if sum / 2 == i: a = open("perfect.txt", 'a') a.write(str(i) + "\n") a.close() i += 1 </code></pre> <p>does <code>i</code> divisions to find the divisors of <code>i</code>.</p> <p>So to find the perfect numbers up to <code>n</code>, it does</p> <pre><code>2 + 3 + 4 + ... + (n-1) + n = n*(n+1)/2 - 1 </code></pre> <p>divisions in the <code>for</code> loop.</p> <p>Now, for <code>n = 33550336</code>, that would be</p> <pre><code>Prelude> 33550336 * (33550336 + 1) `quot` 2 - 1 562812539631615 </code></pre> <p>roughly 5.6 * 10<sup>14</sup> divisions.</p> <p>Assuming your CPU could do 10<sup>9</sup> divisions per second (it most likely can't, 10<sup>8</sup> is a better estimate in my experience, but even that is for machine <code>int</code>s in C), that would take about 560,000 seconds. One day has 86400 seconds, so that would be roughly six and a half days (more than two months with the 10<sup>8</sup> estimate).</p> <p>Your algorithm is just too slow to reach that in reasonable time.</p> <p>If you don't want to use number-theory (even perfect numbers have a very simple structure, and if there are any odd perfect numbers, those are necessarily huge), you can still do better by dividing only up to the square root to find the divisors,</p> <pre><code>i = 2 a = open("perfect.txt", 'w') a.close() while True: sum = 1 root = int(i**0.5) for x in range(2, root+1): if i%x == 0: sum += x + i/x if i == root*root: sum -= x # if i is a square, we have counted the square root twice if sum == i: a = open("perfect.txt", 'a') a.write(str(i) + "\n") a.close() i += 1 </code></pre> <p>that only needs about 1.3 * 10<sup>11</sup> divisions and should find the fifth perfect number in a couple of hours.</p> <p>Without resorting to the explicit formula for even perfect numbers (<code>2^(p-1) * (2^p - 1)</code> for primes <code>p</code> such that <code>2^p - 1</code> is prime), you can speed it up somewhat by finding the prime factorisation of <code>i</code> and computing the divisor sum from that. That will make the test faster for all composite numbers, and much faster for most,</p> <pre><code>def factorisation(n): facts = [] multiplicity = 0 while n%2 == 0: multiplicity += 1 n = n // 2 if multiplicity > 0: facts.append((2,multiplicity)) d = 3 while d*d <= n: if n % d == 0: multiplicity = 0 while n % d == 0: multiplicity += 1 n = n // d facts.append((d,multiplicity)) d += 2 if n > 1: facts.append((n,1)) return facts def divisorSum(n): f = factorisation(n) sum = 1 for (p,e) in f: sum *= (p**(e+1) - 1)/(p-1) return sum def isPerfect(n): return divisorSum(n) == 2*n i = 2 count = 0 out = 10000 while count < 5: if isPerfect(i): print i count += 1 if i == out: print "At",i out *= 5 i += 1 </code></pre> <p>would take an estimated 40 minutes on my machine.</p> <p>Not a bad estimate:</p> <pre><code>$ time python fastperf.py 6 28 496 8128 33550336 real 36m4.595s user 36m2.001s sys 0m0.453s </code></pre>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload