# 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 – Display Name 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... :/ – Display Name 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 – Display Name 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

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