19/02/2009

Project Euler : meilleure solution au problème 12, en Python

J'avais posté lundi une solution très "premier jet" et absolument brutale pour le problème 12 de Project Euler. Evidemment, ceux qui auraient essayé d'exécuter ce programme auront remarqué à quel point le force brute mène à des temps d'exécution gigantesques. En tout cas d'un ordre de grandeur beaucoup trop loin de la "règle de la minute" pour être honnête.

Mais j'avais promis de faire plus propre plus tard, donc chose promise chose due.

L'idée ici se base sur une propriété amusante, c'est que si D(n) est le nombre total de diviseurs d'un naturel (1 et lui-même compris, bien sûr), alors D(n) = (a1 + 1) + (a2 + 2) + ... où les a1 a2... sont les exposants des facteurs premiers de n. Pratique.

Donc ce programme-ci, après avoir établi une liste de nombre premier assez grande compte tenu de la taille du problème (évaluation à la louche, j'avoue), va chercher pour chaque candidat le nombre de ses diviseurs en utilisant la propriété énoncée ci-dessus. Et ça va nettement plus vite, puisqu'on passe de plusieurs heures de calcul à... 35 secondes.

Encore une belle leçon sur l'utilité de trouver un algorithme plus performant que le premier jet.

#! /usr/bin/env python

# 2009/02/19 - euler012-2.py
# Meilleure solution au Probleme 12 de Project Euler
# http://projecteuler.net/index.php?section=problems&id=12
# Jean Karim Bockstael - jkb@jkbockstael.be

def trianglenumber(n):
    return (n * (n - 1) / 2)

def numberofdivisors(n, primes):
    # Recherche des facteurs premiers de n
    factors = []
    i = 0
    while i < len(primes) and primes[i] <= n / 2:
        if n % primes[i] == 0:
            factors.append(primes[i])
        i = i + 1
    # Pour chaque facteur premier, recherche de son exposant
    exponents =[]
    for i in range(0, len(factors)):
        exponent = 0;
        while n % factors[i] == 0:
            exponent = exponent + 1
            n = n / factors[i]
        exponents.append(exponent)
    # Calcul du nombre total de diviseurs
    cnt = 1
    for i in range(0, len(exponents)):
        cnt = cnt * (exponents[i] + 1)
    return cnt

def getfirstprimes(n):
    primes = [2]
    count = 1
    cand = 3
    while count < n:
        i = 0
        while i < len(primes):
            if cand % primes[i] == 0:
                break
            i = i + 1
        if i == len(primes):
            primes.append(cand)
            count = count + 1
        cand = cand + 2
    return primes
    
def euler12(n):
    primes = getfirstprimes(n * 10) # Largement assez
    rank = 1
    tri = trianglenumber(rank)
    while numberofdivisors(tri, primes) <= n:
        rank = rank + 1
        tri = trianglenumber(rank)
    return tri
    
print euler12(500)

Tags:

cc-by-nc | code (python) | project euler


Project Euler : solution au problème 13, en Python (20/02/2009)Project Euler : solution au problème 11, en Python (19/02/2009)