Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

Here's my attempt at Project Euler Problem #5, which is looking quite clumsy when seen the first time. Is there any better way to solve this? Or any built-in library that already does some part of the problem?

'''
    Problem 5: 
    2520 is the smallest number that can be divided by each of the numbers from 1 to 10

    What is the smallest number, that is evenly divisible by each of the numbers from
    1 to 20?
'''

from collections import defaultdict

def smallest_number_divisible(start, end):    
    '''
        Function that calculates LCM of all the numbers from start to end
        It breaks each number into it's prime factorization, 
        simultaneously keeping track of highest power of each prime number
    '''
    # Dictionary to store highest power of each prime number.
    prime_power = defaultdict(int)

    for num in xrange(start, end + 1):
        # Prime number generator to generate all primes till num
        prime_gen = (each_num for each_num in range(2, num + 1) if is_prime(each_num))

        # Iterate over all the prime numbers
        for prime in prime_gen:
            # initially quotient should be 0 for this prime numbers
            # Will be increased, if the num is divisible by the current prime
            quotient = 0

            # Iterate until num is still divisible by current prime
            while num % prime == 0:
                num = num / prime
                quotient += 1

            # If quotient of this priime in dictionary is less than new quotient,
            # update dictionary with new quotient
            if prime_power[prime] < quotient:
                prime_power[prime] = quotient

    # Time to get product of each prime raised to corresponding power  
    product = 1

    # Get each prime number with power
    for prime, power in prime_power.iteritems():
        product *= prime ** power

    return product 


def is_prime(num):
    '''
        Function that takes a `number` and checks whether it's prime or not
        Returns False if not prime
        Returns True if prime
    '''
    for i in xrange(2, int(num ** 0.5) + 1):
        if num % i == 0:
            return False
    return True


if __name__ == '__main__':
    print smallest_number_divisible(1, 20)

    import timeit
    t = timeit.timeit

    print t('smallest_number_divisible(1, 20)', 
             setup = 'from __main__ import smallest_number_divisible', 
             number = 100)

While I timed the code, and it came out with a somewhat ok result. The output came out to be:

0.0295362259729  # average 0.03

Any inputs?

share|improve this question
    
I think you can adapt Erathosthenes' prime sieve to compute is_prime plus the amount of factors in one go. –  flup Feb 6 '13 at 12:07
add comment

4 Answers 4

up vote 7 down vote accepted

You are recomputing the list of prime numbers for each iteration. Do it just once and reuse it. There are also better ways of computing them other than trial division, the sieve of Eratosthenes is very simple yet effective, and will get you a long way in Project Euler. Also, the factors of n are all smaller than n**0.5, so you can break out earlier from your checks.

So add this before the num for loop:

prime_numbers = list_of_primes(int(end**0.5))

And replace prime_gen with :

prime_gen =(each_prime for each_prime in prime_numbers if each_prime <= int(num**0.5))

The list_of_primes function could be like this using trial division :

def list_of_primes(n)
    """Returns a list of all the primes below n"""
    ret = []
    for j in xrange(2, n + 1) :
        for k in xrange(2, int(j**0.5) + 1) :
            if j % k == 0 :
                break
        else :
            ret.append(j)
    return ret

But you are better off with a very basic sieve of Erathostenes:

def list_of_primes(n) :
    sieve = [True for j in xrange(2, n + 1)]
    for j in xrange(2, int(sqrt(n)) + 1) :
        i = j - 2
        if sieve[j - 2]:
            for k in range(j * j, n + 1, j) :
                sieve[k - 2] = False
    return [j for j in xrange(2, n + 1) if sieve[j - 2]]

There is an alternative, better for most cases, definitely for Project Euler #5, way of going about calculating the least common multiple, using the greatest common divisor and Euclid's algorithm:

def gcd(a, b) :
    while b != 0 :
        a, b = b, a % b
    return a

def lcm(a, b) :
    return a // gcd(a, b) * b

reduce(lcm, xrange(start, end + 1))

On my netbook this gets Project Euler's correct result lightning fast:

In [2]: %timeit reduce(lcm, xrange(1, 21))
10000 loops, best of 3: 69.4 us per loop
share|improve this answer
    
I agree, using gcd is the easiest way here. I am quite sure it had to be written before or anywhere else anyway. (For very large numbers, this approach has some problems, but well, 20 is not a very large number) –  tb- Feb 6 '13 at 14:46
    
Hello @Jaime. Thanks for your response. I tried to use your list_of_prime method, but the problem is: - end ** 0.5 range generates just [2, 3] as prime numbers below 20. So, it's not considering 5, which is a valid prime diviser. Same in the case of 10. It's not considering 5. And hence I'm loosing some factors. Also, for numbers like 5, or 7, again it is not considering 5 and 7 respectively, which we should consider right? So, does that mean that range(2, int(end**0.5)) is not working well? –  Rohit Jain Feb 7 '13 at 3:27
    
However, if I replace the range with range(2, end), and in list comprehension also: - each_prime <= num, it's giving me correct result. –  Rohit Jain Feb 7 '13 at 3:30
add comment

You can use the Sieve of Erathosthenes just once AND count the factors while you filter the primes:

def computeMultiplicity(currentPrime, n):
    result = 0
    while n % currentPrime == 0:
        n = n / currentPrime
        result = result + 1
    return result

def euler5(n):
    result = 1
    isPrime = [True for _ in range(n + 1)]

    for currentPrime in range(2, n + 1):
        if isPrime[currentPrime]:
            multiplicity = 1      
            for multiple in range(2 * currentPrime, n + 1, currentPrime):
                isPrime[multiple] = False
                multiplicity = max(multiplicity, computeMultiplicity(currentPrime, multiple))
            result *= currentPrime ** multiplicity

    return result

if __name__ == '__main__':
    print(euler5(20))

    from timeit import timeit
    print(timeit('euler5(20)', setup='from __main__ import euler5', number=100))

Prints:

0.00373393183391
share|improve this answer
    
This code does not appear to be in Python. –  Gareth Rees Nov 2 '13 at 15:12
    
@GarethRees You are most observant ;) It is a Java implementation of a (way) more efficient algorithm. Shouldn't be too hard to port to python since it's all simple arithmetic computations. –  flup Nov 2 '13 at 15:46
    
@GarethRees Ported to Python. Better? :) –  flup Nov 2 '13 at 16:42
    
Yes, that's better. –  Gareth Rees Nov 2 '13 at 18:10
    
isPrime = [True for _ in range(n + 1)] can be written as isPrime = [True] * (n + 1). I would also add a comment that isPrime[0] and isPrime[1] are erroneous but inconsequential. –  200_success Jan 4 at 20:07
show 2 more comments

For a "small" problem like this, performance is a non-issue, and clarity should be the primary consideration. I would decompose the problem into well recognized and testable primitives.

Furthermore, you can eliminate the is_prime() test, as it is just unnecessary extra work. (Why?)

def prime_factors(n):
    '''
    Returns the prime factorization of n as a dict-like object whose keys
    are prime factors and whose values are their powers.

    >>> [(k, v) for k, v in prime_factors(360).iteritems()]
    [(2, 3), (3, 2), (5, 1)]
    '''
    # This is a simplistic, non-optimized implementation.  For better
    # performance with larger inputs, you could terminate the loop earlier, use
    # memoization, or implement the Sieve of Eratosthenes.
    factors_and_powers = defaultdict(int)
    for i in range(2, n+1):
        while n % i == 0:
            factors_and_powers[i] += 1
            n //= i
    return factors_and_powers

def product(factors_and_powers):
    '''
    Given a dict-like object whose keys are factors and whose values
    are their powers, returns the product.

    >>> product({})
    1
    >>> product({3: 1})
    3
    >>> product({2: 3, 3: 2, 5: 1})
    360
    '''
    return reduce(lambda mult, (factor, power): mult * (factor ** power),
                  factors_and_powers.iteritems(),
                  1)

def least_common_multiple(numbers):
    '''
    Returns the least common multiple of numbers, which is an iterable
    that yields integers.

    >>> least_common_multiple([1])
    1
    >>> least_common_multiple([2])
    2
    >>> least_common_multiple([6])
    6
    >>> least_common_multiple([4, 6])
    12
    '''
    lcm_factors = defaultdict(int)
    for n in numbers:
        for factor, power in prime_factors(n).iteritems():
            lcm_factors[factor] = max(lcm_factors[factor], power)
    return product(lcm_factors)

def smallest_number_divisible(start, end):
    '''
    Calculates the LCM of all the numbers from start to end, inclusive.  It
    breaks each number into its prime factorization, keeping track of highest
    power of each prime factor.

    >>> smallest_number_divisible(1, 10)
    2520
    '''
    return least_common_multiple(xrange(start, end + 1))

if __name__ == '__main__':
    import doctest
    doctest.testmod()
    print smallest_number_divisible(1, 20)
share|improve this answer
    
By the way, this version completes in about 1/3 of the time of the original code. –  200_success Nov 3 '13 at 19:04
add comment

There is no least common multiple function in the built-in library. However, there is a greatest common divisor function. You can take advantage of the fact that LCM(a, b) × GCD(a, b) = a b.

from fractions import gcd

def lcm(a, b):
    return (a * b) // gcd(a, b)

reduce(lcm, range(1, 20 + 1))

reduce() calls lcm() repeatedly, as in lcm(… lcm(lcm(lcm(1, 2), 3), 4), … 20).

share|improve this answer
add comment

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.