W tym artykule wracamy do rozwiązywania problemów algorytmicznych i podejmujemy problem #3 z zasobów ProjectEuler.
Problem 3: Largest prime factor
The prime factors of 13195 are 5, 7, 13 and 29.
What is the largest prime factor of the number 600851475143?
Zakres
W celu rozwiązania zadania musimy umieć określić czy dana liczba N jest liczbą pierwszą.
Wbrew pozorom nie jest to problem trywialny, a same liczby pierwsze to bardzo ciekawy i złożony temat.
Wykorzystamy ten problem do zgłębienia kilku podstawowych aspektów liczb pierwszych, choćby dlatego, że wiedza ta będzie przydatna do rozwiązania przyszłych problemów.
Przydatna praca dotycząca ewaluacji algorytmów sprawdzających czy dana liczba jest liczbą pierwszą:
deepai.org.
Liczby pierwsze
Wyznaczanie liczb pierwszych może być zrealizowane w różnoraki sposób, nie istnieje jednak żaden wzór. Obecne osiągnięcia matematyki są coraz, coraz bliżej, jednak wciąż za daleko.
Sito Eratostenesa
Sito Eratostenesa (~ -300 BCE): iteracyjny algorytm, które eliminuje liczby poprzez oznaczanie wielokrotności kolejnych liczb - począwszy od liczby 2 - jako liczby posiadające inne dzielniki. Po przeprowadzeniu tej eliminacji, pozostałe liczby będą liczbami pierwszymi.
Krótki kod algorytmu sita w Python:
def sieve(n: int):
"""Seive of Erathostenes"""
assert n > 1
limit = n ** (1/2) # sqrt
limit = math.ceil(limit)
A = [1] * n
for i in range(2, limit):
if A[i]:
for j in range(i * 2, n, i):
A[j] = 0
return [idx for idx, val in enumerate(A) if val and idx >= 2]
Liczby sprawdzamy w przedziale [2, sqrt(n)].
W liście A oznaczamy na False indeksy, które
otrzymujemy z wyliczenia j jako kolejno: 2i, 3i, 4i, 5i... aż do n. Skrypt tutaj.
Uruchamiamy:
$ time python eratho.py 50
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
real 0m0.023s
user 0m0.023s
sys 0m0.000s
Mierzymy czas wykonania dla liczby wynoszącej 1 million.
$ time python eratho.py 1000000
real 0m0.223s
user 0m0.215s
sys 0m0.008s
Dla liczby 10 millionów algorytm zabiera jednak sporo czasu.
$ time python eratho.py 10000000
real 0m2.238s
user 0m2.209s
sys 0m0.024s
Naiwnie proste i sprytne, do małych przedziałów jak znalazł.
Sito Sundarama
Sito Sundarama (1934): wariant sita Eratostenesa. Skrypt tutaj.
def sieve(n):
"""Sieve of Sundaram"""
k = (n - 2) // 2
A = [1] * (k + 1)
for i in range(1, k + 1):
j = i
while i + j + 2 * i * j <= k:
A[i + j + 2 * i * j] = 0
j += 1
primes = [2 * i + 1 for i in range(1, k + 1) if A[i]]
return ([2] if n > 2 else []) + primes
$ time python sundaram.py 10000000
real 0m7.850s
user 0m7.812s
sys 0m0.020s
Benchmark (Erathostenes vs Sundaram)
Wyniki uzyskujemy poprawne, ale sito Sundarama jest znacznie wolniejsze od oryginalnego wariantu.
Już przy zakresie 10 milionów wykonuje się 3 razy wolniej.
Krótki benchmark:
→ eratho
→ eratho-10000 1 0.00133794
→ eratho-100000 1 0.01448907
→ eratho-1000000 1 0.18316913
→ eratho-10000000 1 2.24281044
--------------------------------------------------------------------------------
→ sundaram
→ sundaram-10000 1 0.00336341
→ sundaram-100000 1 0.04481575
→ sundaram-1000000 1 0.5550501
→ sundaram-10000000 1 6.88032754
================================================================================
RESULTS
================================================================================
GROUP | CASE | N_TIMES | TOTAL TIME | N/SEC
--------------------------------------------------------------------------------
eratho | eratho-10000 | 1 | 0.00133794 | 747.42
eratho | eratho-100000 | 1 | 0.01448907 | 69.02
eratho | eratho-1000000 | 1 | 0.18316913 | 5.46
eratho | eratho-10000000 | 1 | 2.24281044 | 0.45
--------------------------------------------------------------------------------
sundaram | sundaram-10000 | 1 | 0.00336341 | 297.32
sundaram | sundaram-100000 | 1 | 0.04481575 | 22.31
sundaram | sundaram-1000000 | 1 | 0.5550501 | 1.8
sundaram | sundaram-10000000 | 1 | 6.88032754 | 0.15
--------------------------------------------------------------------------------
Wikipedia wspomina również o sicie Atkinsa, jednak nie będziemy go tutaj już implementować.
Metody
Wygenerowanie liczb pierwszych pozwoliłoby nam na szybsze wykonanie algorytmu. Co jednak, kiedy nie chcemy przechowywać wcześniej znanych liczb pierwszych?
Istnieją algorytmy sprawdzające czy dana liczba jest pierwszą, lecz każdy z nich ma właściwe sobie niuanse. Część jest niedeterministyczna (nie rozpoznaje wszystkich liczb pierwszych), część ma zastosowanie tylko dla liczb pierwszych charakteryzujących się szczególną właściwością, itd. Wrócimy do tego tematu w innej odsłonie tej serii.
W matematyce istnieje wciąż nierozwiązany problem Riemanna, który jest powiązany z liczbami pierwszymi: (Wikipedia).
W tym artykule zastosujemy po prostu podejście naiwne, czyli liczby sprawdzać będziemy po kolei, aż do zakresu pierwiastka sprawdzanej liczby.
Python
#!/usr/bin/env python
import sys
import math
def is_prime(n):
assert n >= 2
if n == 2:
return True
top = math.ceil(n ** (1/2))
for i in range(2, top + 1):
if n % i == 0:
return False
return True
def next_prime(n):
for i in range(n, 2 * n):
if is_prime(i):
return i
raise Exception("No prime.")
def max_prime_factor(n):
prime = 2
max_prime = 0
while n > 1:
if n % prime == 0:
n = n / prime
else:
prime = next_prime(prime + 1)
if prime > max_prime:
max_prime = prime
return max_prime
if __name__ == "__main__":
n = int(sys.argv[1]) if len(sys.argv) > 1 else 600851475143
print(max_prime_factor(n))
Podczas uruchomienia tego programu możemy podać inną liczbę, dla której chcemy znaleźć najwyższy czynnik pierwszy, a domyślnie program znajdzie czynnik dla liczby 600851475143 z przedmiotowego problemu.
$ time python prog.py
6857
real 0m0.029s
user 0m0.021s
sys 0m0.008s
Zapamiętujemy wynik, rozwiązania w innych językach mają dać wynik identyczny.
Sprawdźmy jeszcze wykonanie dla innej liczby (123123123):
$ time python code/prog.py 123123123
333667
real 0m0.989s
user 0m0.989s
sys 0m0.000s
C
Analogiczne rozwiązanie w języku C demonstruje poniższy listing.
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
int is_prime(int n)
{ /* naive */
if(n < 2) { return 0; }
else if(n == 2) { return 1; }
int top = ceil(sqrt(n)) + 1;
for(int i = 2; i < top; i++) {
if(n % i == 0)
return 0;
}
return 1;
}
int next_prime(int n)
{ /* naive prime gen. */
for(int i = n; i < (2 * n); i++) {
if(is_prime(i)) {
return i;
}
}
printf("NO PRIME! %d\n", n);
assert(0);
}
void test()
{
assert(next_prime(2) == 2);
assert(next_prime(3) == 3);
assert(next_prime(5) == 5);
assert(next_prime(7) == 7);
assert(next_prime(9) == 11);
assert(next_prime(17) == 17);
assert(next_prime(18) == 19);
}
int main(int argc, char *argv[])
{
test();
unsigned long number = (argc > 1 ? strtoul(argv[1], 0, 10) : 600851475143);
int prime = 2;
int max_prime = 0;
while(number > 1) {
if(number % prime == 0) {
number /= prime;
}
else {
prime = next_prime(prime + 1);
}
if(prime > max_prime) {
max_prime = prime;
}
}
printf("Max prime: %d\n", max_prime);
return 0;
}
Uruchamiamy go dla przedmiotowej liczby (bez argumentu):
$ time ./prog
Max prime: 6857
real 0m0.001s
user 0m0.000s
sys 0m0.001s
A sprawdźmy jeszcze dla liczby 123123123:
$ time ./prog 123123123
Max prime: 333667
real 0m0.050s
user 0m0.046s
sys 0m0.004s
Na koniec tej sekcji warto wspomnieć o tym ciekawym zasobie przybliżającym różne implementacje pierwiastkowania w C: codeproject.