Algorytm faktoryzacji jeszcze szybszy

Jeden z algorytmów, które niedawno opracowałem, podczas przygotowywania do implementacji okazał się ciekawszy niż go szacowałem.

Algorytm polega na zapisie liczby faktoryzowanej n jako iloczynu dwu liczb z jak najmniejszą resztą c:
n = a*b+c
gdzie a, b są liczbami naturalnymi dodatnimi, a<=b oraz c jest mniejsze od mniejszej z liczb a, b.
Inicjacja algorytmu: n = 1*n+0.
Liczba ma dzielniki, kiedy c=0.

Przebieg algorytmu to bezpośrednie przedstawienie
n = a*b0+c0 = (a+1)*b1+c1,
najpierw wyznaczamy pomocniczą zmienną
x = b0+a-c0,
którą dzielimy z resztą r przez (a+1) otrzymując (q,r). Nowe wartości b1 oraz c1 uzyskujemy przez zwykłe odejmowanie
b1 = b0 - q;
c1 = a - r;

Przykładowo mając liczbę 155 = 1*155+0, otrzymamy w pierwszej iteracji
(155+1-0) : (1+1) = 156 : 2 = (78,0)
oraz nowe przedstawienie 2*(155-78)+(1-0) = 2*77+1.
Kolejna iteracja ma dzielenie
(77+2-1) : (2+1) = 78 : 3 = (26,0)
oraz dostajemy: 3*(77-26)+(2-0) = 3*51+2.
Kontynuując przekształcenia mamy
(51+3-2) : (3+1) = 52 : 4 = (13,0)
oraz dostajemy: 4*(51-13)+(3-0) = 4*38+3.
W następnej iteracji z dzielenia
(38+4-3) : (4+1) = 39 : 5 = (7,4)
dostajemy 5*(38-7)+(4-4) = 5*31+0, czyli dzielniki 5, 31.

Możemy kontynuować dalej:
(31+5-0) : (5+1) = 36:6 = (6,0) oraz 6*(31-6)+(5-0) = 6*25+5;
(25+6-5) : (6+1) = 26:7 = (3,5) oraz 7*(25-3)+(6-5) = 7*22+1;
(22+7-1) : (7+1) = 28:8 = (3,4) oraz 8*(22-3)+(7-4) = 8*19+3;
(19+8-3) : (8+1) = 24:9 = (2,6) oraz 9*(19-2)+(8-6) = 9*17+2;
(17+9-2) : (9+1) = 24:10 = (2,4) oraz 10*(17-2)+(9-4) = 10*15+5;
(15+10-5) : (10+1) = 20:11 = (1,9) oraz 11*(15-1)+(10-9) = 11*14+1;
(14+11-1): (11+1) = 24:12 = (2,0) oraz 12*(14-2)+(11-0) = 12*12+11
i na tym zakończyć.
Dzielenia są coraz prostsze, chociaż charakteryzują się sporym szumem.
Dla większych liczb wywołujemy algorytm ponownie dla większego z dzielników (31) mając pewność, że nie mamy dzielników mniejszych niż (5). Zatem w pierwszych iteracjach można zastosować szybkie konwersje a*(2b)+c = (2a)*b+c, a*(b+1)+c = a*b+(a+c).

Algorytm rachuje na coraz mniejszych wartościach, dążących do pierwiastka kwadratowego z rozkładanej liczby, i dlatego nie uważam go za najlepszy. Jest za to najprostszy. 

Kiedy oglądałem go szykując do implementacji, okazało się, że sum b+a-c nie trzeba liczyć!
Wystarczy skorzystać z dzielenia nie przeprowadzonego do końca, o wejściu (a,b0,c0), oraz uzyskać (b1,c1) znacznie mniejszym kosztem.
Polega to na tym, że możemy korzystać z warunku a< d = (a+1) oraz stosować metodę połowienia. Wynik w jest wtedy wartością zewnętrzną dzielenia.

Oto ten sposób: na wejściu dzielenia mamy a, b, c, dzielnik d, iloraz w=0, na wyjściu uzyskamy b1=w, c1=a+c jako niewykorzystaną resztę.
Przebieg algorytmu:
jeżeli b0 jest nieparzysta; c = c+a modulo d, b0-- i czasem modyfikacja w++
teraz b0 jest parzysta, dzielimy ją przez 2, zaś a mnożymy przez 2;
teraz a moze być a>d, bierzemy je wtedy modulo d oraz zwiększamy w o b0;
kończymy, gdy b0=1 lub a=d. Podczas uaktualniania współczynników pobieramy przechowywaną kopię a zwiększoną o 1.

Zastosowanie na tym samym przykładzie n = 155 = 1*155+0
a=1, b=155, c=0, d=a+1=2, iloraz w=0
b=155 nieparzysta, c = c+a = 1, oraz b=154
połowienie: a = 1*2 = 2, b=154:2 = 77
a modulo d jest 0, zatem w = b0 = 77 oraz c = 0+1 = 1
uaktualniamy współczynniki a=2, b=77, c=1, mamy postać 2*77+1!


Kolejna iteracja:
a=2, b=77, c=1, d=3, w=0
b=77 nieparzysta, c = c+a = 3 = 0 (3), b=76, w=1
połowienie: a = 2*2 = 4 = 1 (3), b = 76:2 = 38, w = 1+38 = 39
b=38 parzysta, dalej
połowienie: a = 1*2 = 2, b = 38:2 = 19
b=19 nieparzysta, c = c+a = 2, b=18
połowienie: a = 2*2 = 4 = 1 (3), b = 18:2 = 9, w = 39+9 = 48
b=9 nieparzysta, c = c+a = 3 = 0 (3), b=8, w = 48+1 = 49
połowienie: a = 1*2 = 2, b = 8:2 = 4
b=4 parzysta
połowienie: a = 2*2 = 4 = 1 (3),  b t= 4:2 = 2, w = 49+2 = 51
b=2 parzysta
połowienie: a = 1*2 = 2, b = 2:2 = 1
b=1 kończy iterację, uaktualniamy współczynniki a=3, b=51, c=2+0=2,
postać 3*51+2.

W tym przypadku dzielenie jest proste, zazwyczaj mamy tyle przekształceń, ile wynosi logarytm binarny z b. Złożoność tego dzielenia jest logarytmiczna.
Zatem algorytm ma złożoność O(n log n) arytmetycznie oraz
O(n^2 log n) pamięciowo.

Dzielenie przez sumowanie cyfr

Na stronie spryciarze.pl opublikowany został sposób dzielenia przez 9 wykorzystujący dodawanie narastające cyfr od najbardziej znaczących do cyfry jedności.

Sprawdziłem ten sposób, polega on na tym, że 10 = 1*9+1, oraz odpowiednio grupując ze sobą wyrazy rozwinięcia liczby
an*10^n + ... + a1*10 + a0
uzyska się postać
9*(an*10^{n-1} + ... + (suma _{k>i} ak)*10^i +... + (an+...+a0)) + (an+...+a0)
Z tej postaci można już oszacować iloraz oraz uzyskać resztę (an+...+a0) mod 9.

Sposób można łatwo uogólnić na inne wartości dzielnika. Jeżeli p jest podstawą systemu, oraz dzielimy przez liczbę p-r, to iloraz powstaje przez
sumę narastającą
cn = an
c{n-1} = r*an+a{n-1} = r*cn+a{n-1}
c{n-2} = r*(r*an+a{n-1})+a{n-2} = r*c{n-1} + a{n-2}
...
c0 = r*c1+a0
Najprostsze obliczenia są dla małych r.

Przykład: 3726 : 8 
p=10, r=2, a = 3 7 2 6
c3 = 3
c2 = 2*3+7 = 13
c1 = 2*13+2 = 28
c0 = 2*28+6 = 2*3*8 + 2*4+6 = 7*8+6
Interpretując teraz te liczby jako współczynniki rozwinięcia przy podstawie p liczby (cn...c1), przenosimy nadmiary do liczb bardziej znaczących, zaś po dodaniu jeszcze cześci całkowitej z c0:8 uzyskamy
3*10^2 + 13*10 + 28 + 7 = 465
Ostatecznie 3726 : 8 = 465 reszta 6, co jest prawdziwe.

Zwróćmy uwagę, że obliczenia przy c0 są postaci (d*e+f):g, przy których nie musimy liczyć wartości d*e+f, lecz można wykorzystać dzielenie chłopów rosyjskich aby wyłączać g. W ten sposób szybko zmniejszymy wartości i będziemy liczyć na stosunkowo małych liczbach.

A jak podzielić przez większą liczbę, np. 87?
Dokładnie tak samo, lecz podstawą będzie większa wartość, tu 100.
Przykład: 82045 : 87
p=100, r=13, a = 8 20 45
c2 =8
c1 = 13*8+20 = 124
c0 = 13*124+45 =19*87+4
Iloraz (8+1)*100+(24+19) = 943,
82045 : 87 = 943 reszta 4

To może inna liczba, np. podzielimy przez 11.
Teraz lepiej jest zastosować inne przekształcenie:
ci = -r*c{i+1}+ai
Uzyskana wartość może być ujemna, ale przekształcenia nie zmieniają się.
Przykład: 1751 : 11
p=10, r=1 (dokładniej -1), a = 1 7 5 1
c3 = 1
c2 = -1*1+7 = 6
c1 = -1*6+5 = -1
c0 = -1*(-1)+1 = 2
Iloraz 1*10^2 + 6*10 - 1 + 0 (z reszty) = 159
1751 : 11 = 159 reszta 2

W tym przypadku pojawiają się mniejsze wartości, które mogą oscylować wokół zera jako ciąg rozbieżny. W przypadku reszty dodatniej za pomocą pożyczek doprowadzamy ją do liczby dodatniej - cyfry. Dla reszty ujemnej algorytm tylko szacuje wartość ilorazu - wymaga dopracowania.

Pierwsze z przekształceń lepiej pasuje do liczb większych niż p/2, drugie dla liczb mniejszych od p-p/2. 
Pozostają przypadki dzielenia przez 3, 4. Dla nich najlepiej jest najpierw pomnożyć przez 3, 2 odpowiednio, po czym podzielić przez 9, 8 odpowiednio. Uzyskane reszty są odpowiednimi wielokrotnościami, potrójną, poczwórną właściwych reszt.
Zaś dzielenie przez 5 w systemie dziesiątkowym można zastąpić mnożeniem przez 2 i odcięciem cyfry jedności.  

Sposób ten dla części wielkich liczb jest równoważny z moim dzieleniem przez zmiany systemów (dla ilorazu dwucyfrowego są nawet te same przekształcenia). Lecz w ogólności przy liczbach mniejszych od podstawy sprowadza się do kolejnego dzielenia w c0, choć czasem przez mniejszą wartość. Dużo zależy od odległości dzielnika od podstawy. 
Nie znając szybkich sposobów konwersji, godny polecenia dla dzielników bliskich podstawy.