Lien vers le problème :
import IPython.display
IPython.display.IFrame(src='https://projecteuler.net/problem=500', width=800, height=300)
Explication pour 120 :
2**3 * 3 * 5
120
Le nombre de diviseurs différents est de $4 * 2 * 2$.
4 * 2 * 2
16
Si on avait une puissance de 2 avec autant de diviseurs ce serait :
2 ** 15
32768
Donc effectivement, c'est plus compact de rajouter des nombres premiers pour un total de diviseurs égal à 16. En fait, on peut toujours réduire le nombre de diviseurs en rajoutant une puissance d'un autre nombre premier plus grand.
2 ** 14 * 3
49152
En fait là on a maintenant 15 * 2 = 30 diviseurs. Du coup, si on en veut 16, il faut mettre 2 à la puissance 7.
2 ** 7 * 3
384
Effectivement, c'est déjà beaucoup plus petit. Et on vient d'apprendre la règle que quand on rajoute un nombre premier, on doit baisser le nombre de diviseurs de 2 d'un facteur 2 et rajouter un 1 sur le nombre suivant.
2 ** 3 * 3 * 5
120
2 ** 1 * 3 * 5 * 7
210
210 * 4 / 7
120.0
2 * 3 * 5 * 7 * 11
2310
Autre propriété : si les exposants sont tous 1 alors pour avoir 16 diviseurs, il y a maximum 4 nombres premiers dans l'écriture recherchée car $2 * 2 * 2 * 2 = 16$. Dans le cas de l'exposé, il y a au maximum $2^n = 2^{500500}$ c'est à dire $n = 500500$ nombres premiers dans l'écriture du nombre cherché. On pourrait donc déjà calculer une table de nombres premiers de cette taille là.
La stratégie pourrait donc être :
Formellement : le nombre cherché s'écrit à l'aide d'un certain nombre de nombres premiers $n$ (inconnu a priori)
$$ x = \prod_{i=1}^{n} p_i^{w_i} $$avec
$$ \prod_{i=1}^{n} (w_i + 1) = 2^{500500} $$On sait que $n$ vaut au maximum $500500$.
Les règles d'optimisation sont que si on baisse l'exposant du nombre le plus grand $n$, et que l'on incrémente l'exposant du nombre $k$ alors on obtient un nouveau nombre qui est multiplié par
$$\alpha = \frac{p_k^2}{p_n}$$Si $\alpha < 1$ alors l'opération est bénéfique.
Par contre, on veut conserver le nombre de diviseurs total, qui passe de $\prod_{i=1}^{n} (w_i + 1) = 2^{500500}$ à $\prod_{i=1 i \neq k}^{n-1} (w_i + 1) * (w_k + 3)$. Il est donc multiplié par
$$ \frac{w_k + 3}{(w_n + 1)(w_k + 1)} $$Si on veut que le changement du nombre de diviseurs soit égale à 1, c'est-à-dire qu'il y ait le même nombre de diviseurs au début comme à la fin, il faut que :
$$ w_n (w_k + 1) = 2 $$Donc en fait ce qu'il faudrait connaître c'est la décomposition en nombres premiers de 500500 pour savoir combien on peut incrémenter de nombre sur ce principe. Car on ne peut pas faire toutes les transformations que l'on veut !
Bonne nouvelle, Wolfram Alpha nous indique que :
$$500500 = 2^2 * 5^3 * 7 * 11 * 13$$Donc
$$2^{500500} = 2^{2^2 * 5^3 * 7 * 11 * 13}$$Si on redémarre dans le cas 120 avec
2 * 3 * 5 * 7
210
2**3 * 3**3
216
2**7 * 3
384
Le nombre de diviseurs est le bon : $2^4$.
L'idée est la suivante : on donne les poids et on essaye de modifier itérativement les poids au fur et à mesure.
Exemple, on boucle sur k en prenant l'indice le plus bas. Si la condition d'incrémentation est vérifiée, on calcule le facteur $\alpha$ obtenu. Si on gagne effectivement quelque chose, alors on le fait, sinon on passe.
L'algorithme sur les nombres premiers est tiré de http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n/3035188#3035188.
def rwh_primes1(n):
# http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a list of primes < n """
sieve = [True] * int(n/2)
for i in range(3,int(n**0.5)+1,2):
if sieve[int(i/2)]:
sieve[int(i*i/2)::i] = [False] * int((n-i*i-1)/(2*i)+1)
return [2] + [2*i+1 for i in range(1, int(n/2)) if sieve[i]]
print(rwh_primes1(100))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
from math import log
class euler500_solver():
def __init__(self, number_of_divisors, primes):
self.number_of_divisors = number_of_divisors
self.max_primes = int(log(number_of_divisors) / log(2))
self.primes = primes[:self.max_primes]
self.weights = [1] * self.max_primes
def step(self):
max_index = [(ind, v) for (ind, v) in enumerate(self.weights) if v!= 0][-1][0]
for i in range(max_index + 1):
if self.weights[max_index] * (self.weights[i] + 1) == 2:
if self.primes[i]**2 / self.primes[max_index] < 1:
self.weights[i] += 2
self.weights[max_index] -= 1
return True
return False
def compute_current_number(self, modulo=None):
factors = ([p**w for (p, w) in zip(self.primes, self.weights)])
if modulo == None:
product = 1
for f in factors:
product *= f
return product
else:
product = 1
for f in factors:
product = (product * f) % modulo
return product
e = euler500_solver(16, rwh_primes1(100))
e.primes
[2, 3, 5, 7]
e.weights
[1, 1, 1, 1]
e.compute_current_number()
210
e.step()
True
e.compute_current_number()
120
e.step()
False
e.compute_current_number()
120
On trouve la bonne solution. Comment étendre ça au cas qui nous intéresse ? Il faudrait idéalement un cas intermédiaire qui ne fasse pas exploser l'algorithme.
primes = rwh_primes1(10000000)
len(primes) / 500500
1.3278301698301698
Okay, on a assez de nombres premiers.
e = euler500_solver(2**500500, primes)
e.step()
True
e.step()
True
e.step()
True
last_step = e.step()
while last_step:
last_step = e.step()
e.compute_current_number(modulo=500500507)
493807392
%matplotlib inline
from pylab import plot, array, xlim
plot(e.weights)
xlim(0, 10000)
(0, 10000)
plot(array(e.weights).cumsum())
[<matplotlib.lines.Line2D at 0x10c9a3470>]
max(e.weights)
3
def compute_number_of_divisors(weights):
prod = 1
for weight in weights:
if weight != 0:
prod *= (weight + 1)
return prod
def compute_current_number(primes, weights):
product = 1
for p, w in zip(primes, weights):
if w != 0:
product *= p ** w
return product
On initialise :
p = primes[:5]
p
[2, 3, 5, 7, 11]
weights = [1, 1, 1, 1, 1]
Notre position de départ est :
compute_number_of_divisors(weights)
32
compute_current_number(p, weights)
2310
On veut améliorer ce résultat. On essaye donc de se débarasser du plus haut nombre premier :
weights[4] = 0
weights
[1, 1, 1, 1, 0]
Ceci réduit le nombre de diviseurs :
compute_number_of_divisors(weights)
16
compute_current_number(p, weights)
210
On peut augmenter n'importe lequel des facteurs premiers dans ce cas là :
weights[0] += 2
weights
[3, 1, 1, 1, 0]
compute_number_of_divisors(weights)
32
compute_current_number(p, weights)
840
On se retrouve au status quo. Et si on diminuait encore l'exposant le plus grand ?
weights[3] = 0
weights
[3, 1, 1, 0, 0]
compute_number_of_divisors(weights)
16
compute_current_number(p, weights)
120
On peut envisager deux possibilités :
Comparons les deux cas.
compute_number_of_divisors([7, 1, 1, 0, 0])
32
compute_current_number(p, [7, 1, 1, 0, 0])
1920
compute_number_of_divisors([3, 3, 1, 0, 0])
32
compute_current_number(p, [3, 3, 1, 0, 0])
1080
Donc là, on a plutôt intérêt à augmenter le facteur 3!
weights = [3, 3, 1, 0, 0]
On peut encore essayer d'aller plus loin.
compute_number_of_divisors([7, 3, 0, 0, 0])
32
compute_current_number(p, [7, 3, 0, 0, 0])
3456
Conclusions sur l'exemple avec 32 : j'ai appris que les poids ne pouvaient être que des puissances de 2 moins 1. Et que du coup pour choisir correctement, il fallait toujours vérifier le prochain voisin de poids différents des poids initiaux.
class euler500_solver():
def __init__(self, number_of_divisors, primes):
self.number_of_divisors = number_of_divisors
self.max_primes = int(log(number_of_divisors) / log(2))
self.primes = primes[:self.max_primes]
self.weights = [1] * self.max_primes
def step(self):
max_index = [(ind, v) for (ind, v) in enumerate(self.weights) if v!= 0][-1][0]
max_prime = self.primes[max_index]
for current_index in range(max_index + 1):
current_weight = self.weights[current_index]
current_prime = self.primes[current_index]
next_weight = 2 * current_weight + 1
delta_weight = next_weight - current_weight
if current_prime ** delta_weight / max_prime < 1:
self.weights[max_index] -= 1
self.weights[current_index] = next_weight
return True
break
return False
def compute_current_number(self, modulo=None):
factors = ([p**w for (p, w) in zip(self.primes, self.weights)])
if modulo == None:
product = 1
for f in factors:
product *= f
return product
else:
product = 1
for f in factors:
product = (product * f) % modulo
return product
On teste le nouvel algorithme sur le cas 16 :
e = euler500_solver(16, rwh_primes1(100))
e.step()
True
e.compute_current_number()
120
e.step()
False
Sur le cas 32 :
e = euler500_solver(32, rwh_primes1(100))
e.step()
True
e.compute_current_number()
840
e.weights
[3, 1, 1, 1, 0]
compute_number_of_divisors(e.weights)
32
Le résultat attendu semble marcher.
e = euler500_solver(128, rwh_primes1(100))
e.step()
True
e.step()
True
e.step()
False
e.weights
[3, 3, 1, 1, 1, 0, 0]
On peut maintenant tester le vrai cas :
e = euler500_solver(2**500500, primes)
%%time
last_step = e.step()
while last_step:
last_step = e.step()
CPU times: user 1min 23s, sys: 18.5 s, total: 1min 42s Wall time: 1min 50s
e.weights[:20]
[31, 15, 15, 15, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 3]
e.compute_current_number(modulo=500500507)
35407281
a = np.array([0, 1, 1, 0])
a.nonzero()[0][-1]
2
import numpy as np
from math import log
def rwh_primes1(n):
# http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a list of primes < n """
sieve = [True] * int(n/2)
for i in range(3,int(n**0.5)+1,2):
if sieve[int(i/2)]:
sieve[int(i*i/2)::i] = [False] * int((n-i*i-1)/(2*i)+1)
return [2] + [2*i+1 for i in range(1, int(n/2)) if sieve[i]]
class euler500_solver_numpy():
def __init__(self, number_of_divisors, primes):
self.number_of_divisors = number_of_divisors
self.max_primes = int(log(number_of_divisors) / log(2))
self.primes = np.array(primes[:self.max_primes])
self.weights = np.array([1] * self.max_primes, dtype=np.int16)
def step(self):
max_index = self.weights.nonzero()[0][-1]
max_prime = self.primes[max_index]
#current_weights = self.weights
#next_weights = 2 * self.weights + 1
#delta_weights = next_weights - current_weights
factors = np.power(self.primes, self.weights + 1) / max_prime
candidates = factors < 1
if factors[candidates].size > 0:
new_index = candidates.argmax()
self.weights[max_index] = 0
self.weights[new_index] = 2 * self.weights[new_index] + 1
return True
else:
return False
def iterate(self):
last_step = self.step()
while last_step:
last_step = self.step()
return last_step
def compute_current_number(self, modulo=None):
factors = ([p**w for (p, w) in zip(self.primes, self.weights)])
if modulo == None:
product = 1
for f in factors:
product *= f
return product
else:
product = 1
for f in factors:
product = (product * f) % modulo
return product
e = euler500_solver_numpy(32, rwh_primes1(100))
e.step()
True
e.compute_current_number()
840
primes = rwh_primes1(10000000)
e2 = euler500_solver_numpy(2**500500, primes)
%%time
e2.iterate()
CPU times: user 12 s, sys: 643 ms, total: 12.7 s Wall time: 12.7 s
False
result = 1
for p, w in zip(e2.primes, e2.weights):
result = (result * int(p) ** int(w)) % 500500507
result
35407281
e2.compute_current_number(modulo=500500507)
-c:55: RuntimeWarning: overflow encountered in long_scalars
19538647
e2.weights[:20]
array([31, 15, 15, 15, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 3], dtype=int16)