Calculate practical numbers

Question

Definition

A positive integer n is a practical number (OEIS sequence A005153) iff all smaller positive integers can be represented as sums of distinct divisors of n.

For example, 18 is a practical number: its divisors are 1, 2, 3, 6, 9, and 18, and the other positive integers smaller than 18 can be formed as follows:

 4 =わ 1 +たす 3 5 =わ 2 +たす 3 7 =わ 1 +たす 6
 8 =わ 2 +たす 6 10 =わ 1 +たす 9 11 =わ 2 +たす 9
12 =わ 3 +たす 9 =わ 1 +たす 2 +たす 9 =わ 1 +たす 2 +たす 3 +たす 6
13 =わ 1 +たす 3 +たす 9 14 =わ 2 +たす 3 +たす 9 15 =わ 6 +たす 9
16 =わ 1 +たす 6 +たす 9 17 =わ 2 +たす 6 +たす 9

But 14 is not a practical number: its divisors are 1, 2, 7, and 14, and there's no subset of these which adds to 4, 5, 6, 11, 12, or 13.

Challenge

Write a program, function, or verb which takes as input a positive integer x and either returns or prints the x^th practical number, indexed from 1 for consistency with OEIS. Your code must be sufficiently efficient that it can handle inputs up to 250000 in less than two minutes on a reasonable desktop computer. (NB my reference implementation in Java manages 250000 in less than 0.5 seconds, and my reference implementation in Python manages it in 12 seconds).

Test cases

Input Expected output
1 1
8 18
1000 6500
250000 2764000
1000000 12214770
3000000 39258256

• Display Name: (IMHO) this can be even move interesting if the fastest code (per language?) wins
• Dr. belisarius: @SargeBorsch So you'll see tables of 250K entries all over the answers
• Display Name: @belisarius good point. but I think such cheating can be easily banned. Or the problem may require correct answers for any number, but then there would be difficulties when doing it in a language with no big integers in the standard library... :/
• Display Name: I have one algorithmic optimization in mind, but with current rules I'm too lazy to implement it :P
• Peter Taylor: @SargeBorsch, if you don't want to golf your code feel free to upload it to something like gist.github.com and drop a link in a comment here or in chat. FWIW I prefer code golf with generous performance constraints to fastest code for two reasons: firstly, the length of the code is more objectively measurable; secondly, it introduces an element of tradeoff: which speed optimisations can be left out in order to shorten the code without ruining the performance?

Answer 1

Mathematica, (削除) 126 (削除ここまで) 121 chars

Thanks to belisarius.

Using the formula on wikipedia.

f=(i=j=1;While[j<#,If[And@@Thread[#[[;;,1]]<2+Most@DivisorSum[FoldList[#Power@@#2&,1,#],#&]&@FactorInteger@++i],j++]];i)&

Examples:

f[1]

1

f[8]

18

f[250000]

2764000

It took 70s to compute f[250000] on my computer.

• Dr. belisarius: I think you can get better performance at the expense of one char by bypassing odd integers
• DavidC: In reducing the code from the OEIS submission, you slowed down the execution 10-fold. Just wondering, "why do you think your code runs so much slower than the OEIS example?"
• DavidC: @belisarius Your suggestion cuts the time in half, as expected.
• Dr. belisarius: The same in 119 chars: (i=j=1;While[j<#,If[And@@Thread[#[[;;,1]]<2+Most@DivisorSum[FoldList[#Power@@#2&,1,#],#&]&@FactorInteger@++i],j++]];i)&

Answer 2

J (99 chars)

f=:3 :0
'n c'=.0 1
while.c<y do.
'p e'=.__ q:n=.n+2
c=.c+*/(}.p)<:}:1+*/\(<:p^e+1)%<:p
end.
n+n=0
)

Since the problem statement asks for a "program, function or verb", someone had to make a J submission. J people will notice I didn't really golf (!) or optimize this. Like the other entries, I used Stewart's theorem, mentioned at the OEIS link, to test whether each even number is practical or not.

I don't have ready access to a "reasonable desktop computer" with J installed. On my six year old netbook f 250000 computes in 120.6 seconds, which is not quite under two minutes, but presumably on any slightly more reasonable computer this finishes in time.

Answer 3

Haskell - 329

s 1=[]
s n=p:(s$div n p)where d=dropWhile((/=0).mod n)[2..ceiling$sqrt$fromIntegral n];p=if null d then n else head d
u=foldr(\v l@((n,c):q)->if v==n then(n,c+1):q else(v,1):l)[(0,1)]
i z=(z<2)||(head w==2)&&(and$zipWith(\(n,_)p->n-1<=p)(tail n)$scanl1(*)$map(\(n,c)->(n*n^c-1)`div`(n-1))n)where w=s z;n=u w
f=((filter i[0..])!!)

Examples:

> f 1
1
> f 13
32
> f 1000
6500

Here's a small testing suite (prepend to the above):

import Data.Time.Clock
import System.IO
test x = do
 start <- getCurrentTime
 putStr $ (show x) ++ " -> " ++ (show $ f x)
 finish <- getCurrentTime
 putStrLn $ " [" ++ (show $ diffUTCTime finish start) ++ "]"
main = do
 hSetBuffering stdout NoBuffering
 mapM_ test [1, 8, 1000, 250000, 1000000, 3000000]

Test results after being compiled with ghc -O3:

1 -> 1 [0.000071s]
8 -> 18 [0.000047s]
1000 -> 6500 [0.010045s]
250000 -> 2764000 [29.084049s]
1000000 -> 12214770 [201.374324s]
3000000 -> 39258256 [986.885397s]

• Peter Taylor: When I try this in ghci it complains parse error on input `='. Do I need to use some flag or other?
• mniip: @PeterTaylor You cannot paste function definitions into ghci like that. Simplest you can do is save it to asdf.hs and run ghci asdf.hs, then from there you would be able to access f
• Jonathan Van Matre: @PeterTaylor ghc --make -O3 [filename]. You could also load it in ghci with :l [filename] but given the time constraints compiled is probably best. :)
• mniip: @JonathanVanMatre as seen in the above comment, ghci loads files specified in its arguments
• Peter Taylor: Ah, ok. In the meantime I've got it running with your test framework and ghc. Your computer's faster than mine, but it's still soundly inside the performance criterion on my computer at 98 seconds.

Answer 4

Javascript, (削除) 306 307 (削除ここまで) 282B

function y(r){for(n=r-1,k=1;n;k++)if(p=[],e=[],c=0,P=s=1,!((x=k)%2|1==x)){while(x>1){for(f=x,j=2;j<=Math.sqrt(f);j++)if(f%j==0){f=j;break}f!=p[c-1]?(p.push(f),e.push(2),c++):e[c-1]++,x/=f}for(i=0;c>i;i++){if(p[i]>P+1){s=0;break}P*=(Math.pow(p[i],e[i])-1)/(p[i]-1)}s&&n--}return k-1}

250k in approx. 6s on my laptop.

Commented un-golfed code: http://jsfiddle.net/82xb9/3/ now with better sigma-testing and a better if condition (thank you comments)

Pre-edit versions: http://jsfiddle.net/82xb9/ http://jsfiddle.net/82xb9/1/

• Peter Taylor: The question does ask for a function or program (JS doesn't have verbs), so rather than not counting the first line you should wrap the second line in a function declaration and replace the final k--; with return k-1. Although that increases your byte count slightly, you can save a few with things like replacing p[i]>=P+2 with p[i]>P+1 (and probably by removing the internal function call and using a break instead).
• user2846289: I think "testing sigma" part can be re-written for both size and speed: jsfiddle.net/3DTSa . Though this JS solution is very fast as it is.