Performance comparisons with Mathematica are likely to be very hard to reason about, since Mathematica isn’t, to my knowledge, much of an imperative language as much as a declarative language. For this example, I could imagine that Mathematica is building up an execution model that is very different from the one your Julia code implies. It would be pretty easy to write Julia code that would be much faster than the code you’ve written, but it would involve thinking explicitly about the execution model you’d like to employ.
— John On Mar 16, 2014, at 12:52 PM, Stefan Schwarz <[email protected]> wrote: > Hello, > > I started with Julia last week and to get into tracks I tried to implement > some of the Project Euler question, that I've already done in Mathematica. > > The specifics of this question is about problem #21: > > Let d(n) be defined as the sum of proper divisors of n (numbers less than n > which divide evenly into n). > If d(a) = b and d(b) = a, where a b, then a and b are an amicable pair and > each of a and b are called amicable numbers. > > For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, > 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, > 71 and 142; so d(284) = 220. > > Evaluate the sum of all the amicable numbers under 10000. > > > In Mathematica this was done very easily with: > > DivSum[n_] := DivisorSigma[1, n] - n > AmicableQ[n_Integer] := DivSum[n] != n && DivSum[DivSum[n]] == n > Plus @@ Select[Range[9999], AmicableQ] // Timing > > Yielding a tuple of time and the result: {0.178553, 31626} > > Which is correct. So I tried to achieve the same in Julia: > > function divisors(n::Int, prop::Bool=false) > if(!prop) > filter((x) -> rem(n, x) == 0, [1:n]) > else > filter((x) -> rem(n, x) == 0, [1:n-1]) > end > end > > function divsum(n::Int) > sum(divisors(n, true)) > end > > function amicable(n::Int) > divsum(n) != n && divsum(divsum(n)) == n > end > > This yielded the correct result, but awful timing: > > @time sum(filter(x -> amicable(x),[1:9999])) > > elapsed time: 24.025008574 seconds (4310013792 bytes allocated) > > So. I did something wrong, apparently. Rats. > > So I started to think about, where the hot spot was and > decided to pick the divisors function to be the suspect. > > So I rewrote it: > > function factors(n) > f = [one(n)] > for (p,e) in factor(n) > f = reduce(vcat, f, [f*p^j for j in 1:e]) > end > return f > end > > function sum_prop_factors(n) > sum(factors(n))-n > end > > I did even left out the sorting before return f back, since it isn't needed > at all. > > Executing this: > > @time sum(filter(x -> sum_prop_factors(x) != x && > sum_prop_factors(sum_prop_factors(x))==x, [2:9999])) > > yielded: > > elapsed time: 0.207550349 seconds (39478048 bytes allocated) >
