Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm new to Mathematica and it was suggested to me to go through the Project Euler problems in order to learn it. However, I can't quite figure out why my solution to #5 is so slow.

The problem:

What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

My Mathematica solution:

n = 1; 
max = Fold[Times, 1, Range[20]]; 
While[Fold[Plus, 0, Table[Mod[n, x], {x, 20} ]] != 0 And n < max, 
    Null; 
    n++
]; 
Print[n]

This works but takes a long time to finish (to give you an idea I had time to eat supper and it still wasn't done).

For comparison I wrote a close analogue of the above code in C++:

#include <algorithm>
#include <array>
#include <iostream>
#include <numeric>

int main(int argc, char* argv[])
{
    std::array<int, 20> divisors;
    std::iota(std::begin(divisors), std::end(divisors), 1);
    auto isMultipleOf = [&](unsigned long long x)
    {
        return std::all_of(std::begin(divisors), std::end(divisors), [=](int i) { return (x % i) == 0; });
    };

    unsigned long long n = 1;
    while (!isMultipleOf(n)) { n++; }

    std::cout << n << std::endl;

    std::cin.get();
    return 0;
}

It follows the same structure, a while loop that stops when Modulo returns 0 for each divisor from 1 to 20.

To give you an idea of how much faster this is, I've benchmarked it using nonius with the following setup:

#define NONIUS_RUNNER
#include "Nonius.h++"

#include <algorithm>
#include <array>
#include <iostream>
#include <numeric>

NONIUS_BENCHMARK("Project Euler #5", [](nonius::chronometer meter)
{
    std::vector<unsigned long long> results(meter.runs());
    std::array<int, 20> divisors;
    std::iota(std::begin(divisors), std::end(divisors), 1);

    auto isMultipleOf = [&](unsigned long long x)
    {
        return std::all_of(std::begin(divisors), std::end(divisors), [=](int i) { return (x % i) == 0; });
    };

    std::vector<int> storage(meter.runs());
    meter.measure([&](int i)
    {
        unsigned long long n = 1;
        while (!isMultipleOf(n)) { n++; }
        results[i] = n;
    });
})

And these are the results (and in HTML for here): enter image description here

As you can see, the average is about 5.2s seconds per run. In fact the whole benchmark took less time than one run in Mathematica and I'm really curious as to why.

share|improve this question
    
Correct me if I'm wrong. Your max is around 2*10^18, you start with n=1 and increment once per iteration until max is exceeded? If that's what's happening then yes, it will be slow. –  Daniel Lichtblau Mar 3 at 15:52
    
@DanielLichtblau Actually no I stop way before then. It stops either if max is exceeded or it finds the answer. And the answer is 232792560. The only reason that check is there is because it was taking so long I thought I might have created an infinite loop somehow. –  Borgleader Mar 3 at 16:11
    
Ah. Makes more sense now. The method by @Szabolcs looks like a good translation of your approach. Adding CompilationTarget -> "C" might boost speed somewhat. –  Daniel Lichtblau Mar 3 at 16:30

3 Answers 3

up vote 9 down vote accepted

You've tried to use Fold, which is good, but the spirit of the algorithm is still very "procedural", in that you're not utilizing the Listable properties of certain functions and you're brute-forcing your way through the set of integers. For instance, Mod can take a list as a second argument, Fold[Times, 1, Range@20] is better written as Times @@ Range@20 and so on. Procedural programming is one of the worst possible styles of programming in Mathematica.

The PE question boils down to "implement a least common multiple function" and Mathematica has one built in — LCM. To understand why your solution is slow, let's get a feel for how long it would take. The answer to the problem is:

LCM @@ Range@20
(* 232792560 *)

Using your code with TimeConstrained, n progresses by about 60000 in one second on my machine. Assuming it takes the same time for all integers, we're talking about roughly little more than 1 hour on my machine to reach the answer.

This is one of the things you'll realize when solving PE problems – brute forcing is almost never useful. As you progress through the problems, you'll encounter situations where the memory and time constraints are astronomical. You'll probably even come across numbers that have more zeros than there are atoms in the known universe!


For this example, you can convert the question requirements directly into functional code:

Fold[#/#2 /. 
    {
        Rational[n_, d_] :> # d, (* not evenly divisible, so make it even *)
        d_Integer :> # (* evenly divisible, so leave it alone *)
    } &, 
    1, 
    Range@20 (* numbers from 1-20 *)
]
(* 232792560 *)

This is almost instantaneous and doesn't use LCM. It uses functional and rule-based programming, rather than procedural. Once you're more familiar with the syntax, defaults and special cases, you'll be able to shorten the above code some more by condensing the pattern, but that's not necessary here.

share|improve this answer
    
Technically sweet. +1 –  rasher Mar 3 at 2:17
    
Thank you for the in depth explanation. –  Borgleader Mar 3 at 2:33

This is pretty much instantaneous for me:

 LCM @@ Range[1, 20]
share|improve this answer
    
Using LCM is probably against the spirit of using PE to learn a language... –  rm -rf Mar 3 at 1:49
3  
@rm-rf I would disagree - I think that's the very purpose of using PE to learn Mathematica. It would be against the spirit of Mathematica to force procedural programming when a built-in function will do. –  VF1 Mar 3 at 1:50
    
I didn't say it had to be procedural... it certainly is possible to do without (see my answer). :) –  rm -rf Mar 3 at 2:12
    
@rm-rf I see how your method is more educational. Still though, if I was ever given a choice between both pretty and efficient patterns and a built-in, the built-in would still be a clear winner. –  VF1 Mar 3 at 2:37

Since this is about learning Mathematica, here are a few suggestions:

Fold[Times, 1, list] is better expressed as Times @@ list. You'd use FoldList if you needed to retain all partial results. If instead of Times you had Plus, these two tasks would be accomplished directly by Total and Accumulate.

Some functions are Listable. Table[Mod[n, x], {x, 20}] could be written as Mod[n, Range[20]].

Combining the two improvements from above, Fold[Plus, 0, Table[Mod[n, x], {x, 20}]] can be written as Total@Mod[n, Range[20]].

And cannot be used as infix. && can. And[a,b] is equivalent to a && b and is what you wanted. a And b is equivalent to Times[a, And, b] and is not what you wanted.

At this point we've simplified your code to

n = 1;
max = Times @@ Range[20];
While[
 Total@Mod[n, Range[20]] != 0 && n < max,
 n++
]

This looks better and it's easier to read for people used to Mathematica. But it's not (appreciably) faster.

It is however some fine procedural code which lends itself naturally for compilation. We can try to Compile it. Not everything can be compiled in Mathematica. Those parts of the code that can't be compiled will run the normal way (through the main evaluator). This type of code however can usually be compiled fine without modifications.

cf = Compile[{},
  Module[{n, max},
   n = 1;
   max = Times @@ Range[20];
   While[
    Total@Mod[n, Range[20]] != 0 && n < max,
    n++
    ];
   n
   ]
  ]

We can run it as cf[]. It's still very very slow, but it does finish in a couple of minutes on my machine.

Generally, I don't think it's possible to implement this specific algorithm so it will run much faster than this (in Mathematica), so I stop here. I hope these comments have helped improve your Mathematica.


Of course the proper solution to the PE problem is to use a better algorithm than sequentially checking every number. See the other answers for that. What I showed you here was how to write your own brute force algorithm better in Mathematica.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.