Dziś przedstawię nieco kodów. Python, liczba na początku oznacza poziom wcięcia.
Pierwszy kod to przechodzenie po kolejnych wielokrotnościach kolejnych liczb nieparzystych będących większymi niż liczba rozkładana. Można je podejrzeć wypisując zawartość linii assert. Kod zwraca tablicę wszystkich dzielników pierwszych rozkładu.
0def factorize(n):
1 divisors = []
1 while 0==(n%2):
2 divisors.append(2)
2 n//=2
1 while 0==(n%3):
2 divisors.append(3)
2 n//=3
1 p = 3
1 d = 3-(n%3)
1 r = (n+d)//p
1 while 1<n:
2 a = d+(r<<1)
2 p += 2
2 d = a%p
2 r = r-(a-d)//p
2 assert(d,". ",p,"*",r,"=",n+d)
2 if 0==d:
3 while 0==(n%p):
4 divisors.append(p)
4 n//=p
3 if p*p>n: # no prime divisors less than sqrt n
4 break
2 if p>=r: # finish due p*r = p*p > n
3 break
1 if 1<n: divisors.append(n)
1 return divisors
Obserwując działanie (assert = print), początkowo następuje szybkie zbliżanie się wartości p oraz r. Z czasem jednak następuje tłumienie, pierwsza różnica ukryta w wyrażeniu (a-d)/p maleje geometrycznie do 2. Dla dużych wartości odległość między p oraz r potrafiła być gigantyczna, zaś ich zbliżanie się już było rzędu dziesiątek, czyli sprawdzanie, sprawdzanie, sprawdzanie... Tylko, operacje asymptotycznie dążą do działań na liczbach rzędu pierwiastka kwadratowego z liczby rozkładanej n.
Algorytm cały czas sprawdza kolejne wartości nieparzyste, nawet gdy wiem z własności zachowania się (lokalna monotoniczność przedziałami), że dla tych wartości w otoczeniu nie ma dzielników.
Zatem lepiej jest podmienić ten algorytm dla r<p^2 na inny, w którym można wykorzystać informacje o braku dzielników w sąsiedztwie. Wybór padł na konwersje systemów niedziesiątkowych. Konwersja poniższa jest uogólnieniem konwersji z pracy "Conversion of Number System and Factorization" dla liczby trójcyfrowej [a,b,c]_p.
#conversion number [a,b,c]_{p} into [x,y,z]_{p+k}
0 def conversion(a,b,c,p,k):
1 b = b-k*a
1 c = c-k*b
1 b = b-k*a
1 p += k
1 d=0
1 if 0>c:
2 d = 1+(-c//p)
2 c += d*p
1 elif c>p-1:
2 d = -(c/p)
2 c = (c%p)
1 if c==p:
2 d=-1
2 c=0
1 b -= d
1 d=0
1 if 0>b:
2 d = 1+(-b//p)
2 b += d*p
1 a -= d
1 return [a,b,c,p]
Jest tu jedna instrukcja if, która wydaje się być nie na miejscu, gdyż warunek jest spełniany if wcześniej. W praktyce u mnie bez tego ostatniego if c==p, wartość c okazywała się być równa p wbrew wcześniejszemu c = (c%p).
I okazuje się, że dla bardzo, bardzo dużych wartości można zmniejszyć krotność rozpatrywanych iteracji, i to nawet z dwu stron. Trzeba robić to jednak bardzo delikatnie. W pierwszym przypadku reszty z dzielenia przez kolejne wartości tworzą ciąg lokalnie monotonicznie rosnący, w drugim lokalnie monotonicznie malejący. Zaś możliwość wystąpienia dzielnika jest tylko przy tych wartościach, w których następuje przeniesienie.
Sprowadza się to do uważnego śledzenia, kiedy warunek wyrażenia b-2*a*k oraz wartość bk po konwersji o k rozjeżdzają się. Dokładniej, jeśli e = (b-2*a*k)-bk przy spełnieniu pewnych warunków (nie zawsze można wziąć b z marszu), to mamy dla zapisu liczby n w przedziale podstaw systemu [p, p+k] wystąpienie dokładnie e przeniesień. Czyli dla końcówki algorytmu -- z tysięcy wartości wystarczy zlokalizować i sprawdzić kilka. I to w obszarze bardzo obciążającym obliczeniowo procesory działające innymi metodami.
