2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder. What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

A naive implementation of this in Matlab (especially with for() loops) exceeded the 1 minute time limit. After playing around some with the LCM of 1-10 to arrive at the provided answer of 2520, I came to the following algorithm (which has been previously posted in the Project Euler solution forum in similar but less clear forms).

All numbers have a prime factorization. Create a set from the prime factorization, but allow repeats of a prime factor to be distinct from each other. The LCM of a set of numbers has a prime factorization that is equal to the union of the prime factorization sets of each number. Pseudo-code implementation is as follows:

- Find all primes 1 to N inclusive, call this P. -This can be done most efficiently for large N using a sieve.
- Form a subset of P of primes < sqrt(N), call this S.
- For all primes in S find one less than the largest power that still occurs in 1-N inclusive. -This can be done easily by calculating floor( log(N)/log(S) )-1
- The LCM is the product of all primes in P times the product of all primes in S raised to their max_exponent-1 found in step 3.

*N.b.* We need only need to find the max exponents of the primes in S, not P, because any prime greater than the sqrt(N) will not have at least a square <N so we can skip it.

The Matlab code is quite straightforward. I use the Matlab function primes() to generate all primes to N. For N=20 the LCM fits within a 64-bit double or int. For N>42 the LCM is too large to fit in an int, and the accuracy starts to degrade as a double. To avoid these issues, I represent the result symbolically using the Symbolic Math Toolbox.

The code is very fast and able to calculate N=10^7 under the 1 minute time limit, the result is 4342311 digits long!

Check my Github repo for the code.

need to find your matlab program that goes with the lecture

1D Linear Solver Demo M2.11 – Intro to DG

thx

scott lucas

this code fails for argument ‘4’. i wonder why

Good catch, it looks like Matlab’s symbolic math has weird floating point issues in the following sub-statement:

floor(log(N)./log(A))

This should be 2 in the case of N=4, but it instead returns 1. If you just compute:

floor(log(4)./log(2))

it properly returns 2.

A simple solution is to just perturb the numerator up by the local floating point resolution:

floor((log(N)+eps(N))./log(A))

This returns the expected result.