# Calculate practical numbers

### 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 xth 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
``````
-
(IMHO) this can be even move interesting if the fastest code (per language?) wins – Sarge Borsch Mar 15 '14 at 18:09
@SargeBorsch So you'll see tables of 250K entries all over the answers – Dr. belisarius Mar 15 '14 at 22:37
@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... :/ – Sarge Borsch Mar 15 '14 at 22:45
I have one algorithmic optimization in mind, but with current rules I'm too lazy to implement it :P – Sarge Borsch Mar 15 '14 at 22:46
@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? – Peter Taylor Mar 15 '14 at 23:19

# 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.

-

# 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.

-
I think you can get better performance at the expense of one char by bypassing odd integers – Dr. belisarius Mar 15 '14 at 16:05
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 Mar 15 '14 at 16:08
@belisarius Your suggestion cuts the time in half, as expected. – DavidC Mar 15 '14 at 20:49
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)&` – Dr. belisarius Mar 16 '14 at 6:19

# 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]
``````
-
When I try this in ghci it complains `parse error on input `='`. Do I need to use some flag or other? – Peter Taylor Mar 15 '14 at 12:41
@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` – mniip Mar 15 '14 at 12:51
@PeterTaylor `ghc --make -O3 [filename]`. You could also load it in ghci with `:l [filename]` but given the time constraints compiled is probably best. :) – Jonathan Van Matre Mar 15 '14 at 12:54
@JonathanVanMatre as seen in the above comment, `ghci` loads files specified in its arguments – mniip Mar 15 '14 at 12:55
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. – Peter Taylor Mar 15 '14 at 12:55

# 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/

-
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). – Peter Taylor Mar 16 '14 at 17:22
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. – user2846289 Mar 17 '14 at 22:43