Dalsze badania

Mija miesiąc od ostatniego posta, czas na kolejny algorytm. Ale nie ma żadnego nowego mimo wytężonej pracy.
Przez ten czas przyglądałem się zastosowaniem wzorów skróconego mnożenia do przekształceń liczb. Ponownie odkryłem metodę Fermata a^2-b^2 w trzech odsłonach:
- klasyczną (np. u Knutha Sztuka programowania 4.5.4.C), kiedy odjemna jest zawsze kwadratem, zaś zwiększam b czasem zmniejszając a;
- odwróconą, kiedy odjemnik jest zawsze kwadratem, wtedy zmniejszamy b od czasu do czasu zwiększając a;
- mieszaną, kiedy naprzemiennie przenoszę jakąś wartość między odjemną a odjemnikiem pilnując, by zawsze co najmniej jedno z nich było kwadratem.

Algorytmy te zapisane za pomocą systemów niedziesiątkowych w obszarze z dużymi a oraz b sprowadzają się do prostych konwersji typu
1'2'a_p = 1'0'(a-1)_{p+1}
co oznacza, że dla bardzo dużych liczb działamy na małych wartościach. Nie zmienia to jednak złożoności.

Sprawdziłem także zachowanie się sumy sześcianów. Co prawda, nie każda liczba wyrazi się jako suma dwu sześcianów, ale przekształcenia wzięte z różnicy kwadratów działają i dla sześcianów.

Kiedy zbliżamy się do dzielników, jakiekolwiek przekształcenia wskazujące: 'jesteśmy coraz bliżej' zaczynają szwankować. Dlatego lepszą grupą są algorytmy szukajace dzielników mimochodem. Skaczące między różnymi podstawami.
W szczególności zastosowanie cechy podzielności przez (znany) dzielnik d w systemie o podstawie d+1 (cecha: różnica naprzemienna cyfr w zapisie liczby n jest 0 lub jest podzielna przez n) liczby trójcyfrowej n = a'b'c_p prowadzi do równania kwadratowego
ar^2 - br + c = 0
którego rozwiązania r czasem (nie zawsze) mogą przybliżać okolice dzielnika. W kolejnej iteracji przedstawiamy n w systemie o podstawie p+r albo p-r.
Kryterium: suma naprzemienna cyfr jest dodatnia, dzielnik może być mniejszy niż aktualna podstawa. Suma naprzemienna cyfr jest ujemna, dzielnik jest przy podstawie większej. Suma naprzemienna równa 0, dzielnik jest dokładnie o 1 mniejszy niż aktualna podstawa liczby.

Lecz zmian znaku sumy naprzemiennej cyfr liczby n zapisanej w różnych systemach pozycyjnych jako trójcyfrowa jest wiele. Wskazówką powodzenia może być fakt, że bardzo blisko dzielnika te zmiany są rzadkie. Tuż przy samym dzielniku równanie sugeruje daleki skok. Najdalsze propozycje są najbliżej dzielnika. Wtedy lepiej zastosować metodę siłową.

Pierwiastek kwadratowy, przyspieszona wersja

Algorytm we wcześniejszym poście wymagał znalezienia pierwiastka kwadratowego z liczby, liczonego dwukrotnie podczas inicjacji.
 Publikowałem już sposób na liczenie pierwiastka rok temu, w listopadzie 2012 roku, lecz teraz znalazłem nieco szybszą technikę.

Ustalmy liczbę n oraz drugą liczbę p położoną między pierwiastkiem sześciennym oraz kwadratowym z n. Nic więcej o liczbie p nie potrzeba wiedzieć, jest ona pierwszym przybliżeniem pierwiastka. Resztę załatwią konwersje, a być może także wzór skróconego mnożenia.

Dzielimy n z resztą przez p (dwukrotnie), by uzyskać postać
n = (a*p+b)*p+c                     (1)
Jest to zapis liczby trójcyfrowej.

Możemy sprawdzić, czy a, b, c stanowią współczynniki trójmianu kwadratowego d^2+2de+e^2+f, jeśli tak, pierwiastek można przybliżyć przez
(dp+e)^2+f.
Jeśli jednak b>2de, też możemy posłużyć się wzorem, do którego należy dodać poprawkę (b-2de)*p. Sposób też prowadzi do celu, lecz jest znacznie dłuższy. Zbieżność do pierwiastka jest wolniejsza.

Zatem pierwsza faza, sprowadzamy a do postaci podzielnej przez 4, b do liczby parzystej za pomocą formuł:
a-1 AND b+=p;   b-1 AND c+=p.               (2)
Np. mając  (a,b,c) = (5, 4, 34) dla p=125 uzyskujemy z pomocą formuł (2)
(a,b,c) = (5-1, 4+125, 34) = (4, 129, 34) = (4, 129-1, 34+125) = (4, 128, 159)
Korzystamy z konwersji podwajania systemu
(a,b,c)_p = (a/4, b/2, c)_(2p)                       (3)
czyli w przykładzie mamy (a,b,c) = (1, 64, 159)_250 .
Kiedy już a<4, oraz nie dopasujemy wzoru skróconego mnożenia, dążymy do zmniejszenia b. W odróżnieniu od zeszłorocznego algorytmu najlepiej to zrobić za pomocą konwersji o floor(a*b/3). Konwersja ta zmniejsza b o co najmniej połowę, może także zmniejszyć a. Branie połowy (a*b/2) może doprowadzić do zmniejszenie b do wartości ujemnej, zaś chcemy zachować a=1 jako niezmiennik, oraz wartość liczby trójcyfrowej powinna pozostać większa niż pierwiastek kwadratowy. 

Przykład dla n=8934053
Przyjmijmy p=1000, wtedy (a,b,c) = (8, 934, 53)_1000
Pierwsza konwersja (3) sprowadza do (a,b,c) = (2, 467, 53)_2000
konwersja o floor(a*b/3) = 311
2    467              53
2    467-2*311
2-1 2156            53              // dodaję do b (2000+311) = 2311
1    2156-311     53-2156*311
1    1845-291    -670463      // dodaję do c 291*p
1    1554           2038
oraz (a,b,c) = (1, 1554, 2038)_2311
konwersja o floor(a*b/3) = 518, nowa podstawa p = 2311+518 = 2829
1    1554         2038
1    1554-518
1    1036        2038
1    518          2038-1036*518  // dodajemy do c 189*p
1    518-189   -534610
1    329         71 
konwersja o floor(a*b/3) = 109, nowa podstawa p = 2829+109 = 2938
1    329        71
1   329-109 
1   220         71
1   220-109  71-220*109
1   111-9     -23909              // dodajamy do c 9*p
1   102        2533
Sprawdzamy wzór skróconego mnożenia (1+51)^2 = 1+102+2601 oraz porównujemy. Niezgodność na pozycji c, jest za dużo o 68. Ale jest dobrze.
Zatem (a,b,c)_p = (1*p+51)^2-68 = (2938+51)^2-68 = 2989^2-68.
Zmniejszamy kwadrat
2989^2-68 = 2988^2+5976+1-68 = 2988^2+5909.
Znaleziony został pierwiastek całkowitoliczbowy 2988 oraz reszty: 5909 z niedomiarem, 68 z nadmiarem.

Dokładniejszą wartość pierwiastka mozna policzyć za pomocą interpolacji 2988 + 5909/(5909+68) = 2988,988623; dokładniejsza wartość 2988,9886249365352845723814
Błąd dopiero na szóstej pozycji po przecinku.

Różnica kwadratów w faktoryzacji

Przedstawienie liczby rozkładanej n w postaci dwu kwadratów leży u podstaw jednych z najszybszych metod faktoryzacji. Są to np. metoda sita kwadratowego, sita ciał liczbowych, a nawet ułamków łańcuchowych.
Wykorzystują one kongruencję Legendre'a
x*x = y*y (mod n) ,
gdzie wartości są dodatnie,  x<y<n  oraz x+y nie sumuje się do n.
Sposób stosowania tej kongruencji znajduje się w "Teorii liczb w informatyce" Song Y. Yana.

Teraz jednak chcę przedstawić inny algorytm, w którym wykorzystuję jedną z najbardziej prostych konwersji systemów niedziesiątkowych, a zarazem wzór skróconego mnożenia:
p^2 + 2*p + 1 = (p+1)^2 + 0*p + 0 = (p+1)^2             (1) .

Kiedy przedstawię liczbę n w postaci:
n = a*a - b*b - c ,                                                         (2)
gdzie a jest najmniejsze możliwe a^2 < n < (a+1)^2, b jest największe możliwe. Wartość b^2+c można uzyskać wieloma sposobami, ale nie udało mi się jej zmniejszyć.
Przyrost między kolejnymi kwadratami p^2 oraz (p+1)^2 jest równy 2*p+1.
Dodając to do wyrażeń a^2, b^2+c ich różnica pozostanie równa n.
Po zwiększeniu c można dalej zwiększać b, co prowadzi do następującego algorytmu, który zatrzymuje się gdy c=0. Wtedy dzielnikami są:
n = a^2 - b^2 = (a+b)*(a-b) .
Drugim wyjściem jest zwiększenie b powyżej pierwiastka z n, co świadczy o braku dzielników. 

inicjacja: a=deck(sqrt(n)); c=a^2-n; b=floor(sqrt(c)); c=c-b^2; 
k=a; // kopia dla ograniczenia zakresu
while ( b<k ) {
  a++;
  c+=2*a+1; // dodanie niedomiaru od najbliższego większego kwadratu
  while( 1 ) {
    d = 2*b+1;  // by nie liczyć dwukrotnie niedomiaru do kwadratu b
    if( c<d ) break;
    else { c-=d; b++; } // b wzrasta do kolejnego kwadratu kosztem c
    if( 0==c ) return text("dzielniki: ", a-b, a+b);
  }
}
return text("n jest pierwsza");

Funkcja text() wypisuje swoje argumenty na wyjściu. Pierwiastek z n powoli wzrasta po kolejnych kwadratach. Zaś b^2+c przyrastajace o tę samą wartość modyfikuje się do największego możliwego dla tej wartości kwadratu.
Ze wzrostem b spada liczność obliczeń w pętli wewnętrznej.

Złożoność.
Jeśli przyjmiemy, że pętla wewnętrzna wykona się dokładnie raz, pogorszymy złożoność. Jest to spowodowane tym, że przy każdym, nawet fikcyjnym zwiększeniu a bedziemy zwiększać b, zaś istotniejsze są zwiększenia b.  Wyrażenia d=2*b+1 wykonują się w czasie liniowym, pozostałe przekształcenia to czas praktycznie stały. Pętla zewnętrzna wykona się nie więcej niż pierwiastek z n razy, zwiększając b o 1 w każdym kroku. Praktycznie mniej, gdyż b może nie być bardzo małe. Dodatkowo, w pierszych iteracjach zwiększanie b odbywa się bardzo szybko, zwalniając w czasie przebiegu do 1-2 zwiększeń b na jedno zwiększenie a. Podsumowując, złożoność tego algorytmu nie przekracza
O( n sqrt(n) ) dla pamięci,
max{ O( sqrt(n) ), O( sqrt ) } dla operacji.

Dla 18703 mamy wartości początkowe a=137, b^2+c = 66, zatem b=8, c=2.
Jest to złośliwy przypadek, gdyż dzielniki są dosyć odległe, zostaną znalezione dopiero gdy b=129, czyli w 51 na 56 iteracji. Podczas każdej iteracji pętli zewnętrznej wartość b zwiększała się od 1 do 8, najczęściej 2 lub 3.
W następnym kroku mamy:
a=138, jest dodane 2*137+1 = 275;
od c = 2+275 = 277 odejmujemy 2*8+1= 17 (wtedy b=9)
2*9+1 = 19 (oraz b=10)
2*10+1 = 21 (oraz b=11)
jeszcze parę razy i otrzymamy c=17 przy b=18,
Mamy do czynienia z wyrażeniem 138^2 - (18^2+17) = 19044 - 341 = 18703.
W kolejnej iteracji dodamy 2*138+1 = 277. Będziemy zabierać 37, 39, itd. dopóki b=24.
Wartości dodawane i zabierane tworzą proste ciągi arytmetyczne.


Algorytm dla liczb postaci n = p*q, gdzie p i q są pierwsze ma dodatkową własność. Wartość 2*b wskazuje minimalną odległość między tymi dzielnikami. 

Wariacja faktoryzacja przez proste dzielenie

Rozkład liczby n na czynniki za pomocą prostego dzielenia polega na kolejnym wstawianiu do ilorazu n/b kolejnych, najlepiej pierwszych, ewentualnie nieparzystych, liczb b. Reszta równa 0 oznacza, że b jest dzielnikiem.

Dla dużych liczb n b też jest duże, chociaż mniejsze niż pierwiastek z n. Wtedy dzielenie jest uciążliwe.
Istnieje sposób, by zmniejszyć argumenty dzielenia, by nie dzielić n/b, lecz jakieś c/b, gdzie c<n. A nawet dwa sposoby.
Pierwszy korzysta z przedstawienia liczby n jako 'liczby dwycyfrowej' jakiegoś systemu liczenia o dużej podstawie, oraz konwersje pozwalają na znajdowanie reszt - niestety, metoda wymaga liczenia ilorazów dla kolejnych liczb naturalnych przy zmniejszaniu dzielnej.

Metoda, którą opiszę dalej, pozwala przeskakiwać liczby parzyste. Zatem tworzymy ilorazy c/b, gdzie c jest nieparzyste, mniejsze od n.

Zapiszmy liczbę n w postaci wyrażenia:
n = a*b + c      (1)
będziemy przekształcać wartości a, b, oraz c w taki sposób, by podczas przekształceń a tworzyło ciąg nierosnący, b było ciągiem rosnącym po wartościach nieparzystych, zaś c tworzyło ciąg przedziałami monotoniczny.
Algorytm jest rozgałęziony, a jego ogólny schemat wygląda następująco:

inicjacja: a = floor(n/3), b=3, c=n%3.
pętla dopóki a>b
czy b dzieli c (warunek 0 = c%b)? jeśli tak, b jest dzielnikiem, wyjście;
b = b+2;
rozgałęzienie, jeśli c<2a
c = c + (a%b)*(b-2);  
a = b*floor(a/b) - 2*floor(a/b);
w przeciwnym razie (c>2a)
c = c-2a; 
koniec pętli

Pierwsza część rozgałęzienia jest całkowito-liczbową operacją a*(b-2)/b, która zwiększa c oraz zmniejsza a. Mamy tu dodatkowo jedno dzielenie oraz mnożenie, co komplikuje algorytm. Można je przekształcić, by dzielić przez b oraz odjąć podwojony iloraz od a.
Z kolei druga część rozgałęzienia jest zwykłym odejmowaniem.

W początkowej fazie algorytmu przeważa pierwsza część, z dodatkowymi działaniami. Wykonują się jednak one na małych wartościach, dużo mniejszych niż pierwiastek z n. Pod koniec algorytmu najgorszym kawałkiem jest mnożenie (a%b)*(b-2), którego wartość może być stosunkowo blisko n. Zwiększa ona c do bardzo dużych wartości, umożliwiając stosowanie drugiej, prostszej odnogi.
Już po sprawdzeniu około 20% przypadków na możliwe wartości b do głosu dochodzi druga część rozgałęzienia, zaś przy 30% praktycznie dominuje.
Cały czas należy sprawdzać podzielność c przez b. Chociaż c (zwłaszcza przy dominacji drugiej odnogi) szybko maleje nawet do wartości bliskich b.

Fragment przykładu numerycznego, liczba 8 934 053 = 1087 * 8219.
inicjacja: a*b+c = 2 978 017 * 3 + 2
b = 3+2 = 5; 
pierwsza odnoga: 2 978 017 - 2 = 2 978 015, bo 2 978 017 % 5 = 2,
c = 2 + 2*3 = 8
zmniejszanie a: 2 978 015 * 3 / 5 = 2 978 015 - 2*(595 603) = 1 786 809

nieco dalej:
a = 122 049; b = 73; c = 24 476;
sprawdzamy c%b = 21, nie jest zerem; pierwsza odnoga
b=75; reszta a%b = 24
c = 24 476 + 24*(75-2) = 26 228;
a = (122 049 - 24) - 2*(122 049 - 24)/75 = 118 771

jeszcze dalej:
a = 69 125; b = 127; c = 155 178;
sprawdzamy c%b = 111, druga odnoga
b = 129;
c = 155 178 - 2 * 69 125 = 16 928;
możemy znów sprawdzać c%b = 16 928 % 129

dla b>400 co parę przekształceń mamy kilka iteracji odnogą drugą, zaś blisko końca algorytmu mamy
a = 8685; b = 1027; c = 14 558;
sprawdzamy c%b = 180; pierwsza odnoga
b = 1029;
c = c + 453*1027 = 479 789;
a = (a-453) - 2*(a-453)/1029 = 8232-16 = 8216;
Teraz przez 29 iteracji powtarza się odejmowanie w drugiej odnodze, zanim znów zastosujemy pierwszą.

Zakończenie algorytmu
a*b+c = 8216 * 1087 + 3261, oraz 3261 = 3*1087.   

Algorytm nie jest zbyt dobry. Może być nawet gorszy niż zwykłe dzielenia. Może być traktowany jako ciekawostka.

Modyfikacja heurystyki rho Pollarda

Kilkanaście lat temu zetknąłem się z heurystyką rho Pollarda. Polega ona na obliczaniu kolejnej wartości wielomianu W[x] modulo n startujacego z ziarna a. Metoda ładnie i dostępnie opisana w wikipedii, książkach do kryptografii.
Opiszę tu swoje (alternatywne) podejście, bo jak większość moich pomysłów, ciężko jest mi je publikować.

We wszystkich tych źródłach pojawia się informacja, że wielomianem nie może być wielomian x*x-2 ani x*x.
A właśnie drugi z tych wielomianów użyłem, kiedy zacząłem testować heurystykę. Występowała niejednoznaczność, którą można było bardzo łatwo usunąć biorąc kolejną wartość wielomianu jako
f(x) = min( x^2, n-x^2 ) mod n .

Zastosowanie tego wielomianu, nie dość, że zmniejsza o połowę krotność występujacych współczyników, to jeszcze ma kilka dodatkowych własności bardzo przyspieszajacych obliczenia. Można używać indeksów dla szacowania, kiedy i gdzie może pojawić się dzielnik.
Zmienia się także warunek stopu. Nie jest już jeden:
nwd( f(a)^m, n ),
gdzie f(a)^m oznacza m-krotne złożenie f liczone w a,
czy w modyfikacji p-1: nwd( f(a)^m -1, n ) .
Odpada także przygotowanie liczby do obliczeń przez zignorowanie kilku początkowch iteracji, jakie spotkałem w jednej z książek dotyczących algorytmów.

Wybieramy ziarno 1<a<(n-1)/2. Zapamiętujemy jego wartość (dla indeksu).
Obliczamy f(a), kiedy po raz pierwszy uzyskamy wartość większą niż n, wtedy zapamiętujemy dodatkowo b jako mniejszą z reszt z dzielenia przez n: b = min( b%n, (n-b)%n ).
Powtarzamy obliczenia zliczając złożenia, aż uzyskamy wartość a, b, 0 albo 1.
W przypadku uzyskania 0 dzielnikiem jest ziarno a.
Gdy uzyskamy ciąg okresowy z wyrazem a lub b, liczymy dla tego wyrazu:
np. d=nwd( a-1,n).
Gdy 0<d<n/2, d jest dzielnikiem. Przy wartości d=0 jeszcze nic straconego. Liczymy dla pewnego elementu x =f(a)^m:
e = nwd( x, f(x) ) ,
który jest kolejnym kandydatem na dzielnik.
Najciekawszy przypadek jest wtedy, gdy cyklem jest 1. W standardowym podejściu ten przypadek powoduje błąd heurystyki. Ale nie tutaj.
Zaczynamy algorytm z ziarnem 1<a+1<(n-1)/2. I teraz tylko liczby pierwsze mogą uzyskać ciąg złożony z 1. W sprawdzanych przeze mnie przypadkach pojawia się inny ciąg stały złożony z dzielnika n. Cykle jedynek i dzielnika przeplatają się wzajemnie dla kolejnych ziaren.
Gdy nie znajdziemy dzielnika, zwiększamy ziarno o 1, oraz puszczamy algorytm jeszcze raz, pamiętając dodatkowo poprzednie wartości a, b. Teraz istotne jest ich znalezienie w nowo powstającym ciągu. 
Cykl dla sąsiedniego ziarna nie musi zawierać zapamiętanych wartości. Jest wtedy rozłączny, co pozwala oszacować, czy wogóle jest miejsce na dzielniki. Czy będą związane z jakąś liczbą pierwszą, czyli gdzie ich szukać.
Kilkanaście lat temu myślałem o pamiętaniu ich w tablicy, by startować z ziarna jeszcze nie występującego.

Ale po co zapamietywać liczność cyklu? Jest on potrzebny dla oczacowania, czy mamy do czynienia z liczbą pierwszą czy złożoną. Cykl długości (n-3)/4 powtarzający się dla dwu sąsiednich (względnie pierwszych) ziaren wskazuje liczbę pierwszą.

Uwaga, ziaren równych 0, 1, n-1 nie należy brać pod uwagę, uzyska się stały ciąg trywialnej wartości.

Przykłady.
Liczba 51.
Ziarno 2, ciąg wpada w cykl [1].
Ziarno 3, ciąg wpada w cykl [18], nwd(18-1, 51) = 17, oraz 51=3*17.
Przeplatane cykle [1] oraz [18]. 

Liczba 35.
Ziarno 2, ciąg [2, 4, 11, 16,... ] ogon [2,4], cykl (druga z zapamiętanych wartości to 11) [11, 16], liczymy -1+11 = 10 (mod 35), d = nwd( 10, 35 ) = 5 jest dzielnikiem.

Ziarno 5, ciąg [5, 10, 5...] z cyklem [5, 10]. Teraz d=nwd(4,35)=1. Drugi warunek stopu e=nwd(5,10)=5 wskazuje dzielnik 35.

Liczba pierwsza 47.
Dla wiekszości ziaren mamy cykle długości 11 elementów, ziarno 2 ma zapamiętaną wartość 2, ziarno 3 ma cykl rozłączny (nie zawierający 2), oba cykle wybierają 22 elementy z (47-3)/2 = 22 dostępnych. Jest to liczba pierwsza.

Faktoryzacja a system Fibonacciego

Liczby Fibonacciego tworzone są wzorem rekurencyjnym:
F[0]=0, F[1], F[n+2] = F[n+1]+F[n]
W wikipedii opisano system Fibonacciego, w którym w definicji systemu o cyfrach a[i] oraz podstawach p[i]
suma_{-1<i} a[i]*p[i]
przyjęto
p[i] = F[i-2], 1<i, a[i] in {0,1}.

Dodawanie i odejmowanie w tym systemie jest proste, mnożenie i dzielenie koszmarne - lepiej się trzymać z daleka. Dla usunięcia niejednoznaczności wprowadzono formułę łączącą w okienku kolejne trzy cyfry systemu Fibonacciego:
[0,1,1]_F = [1,0,0]_F      (1) .

Nie musimy mieć cyfr równych tylko {0,1}. Możemy czasowo pozwolić sobie na dowolne 'cyfry' całkowite, a wtedy mnożenie i dzielenie nieco się ułatwi.
Przykład. Liczba [3,0,0,0,0,0]_F = 3*F[7] = 3*13 = 39.

Własność (1) można wykorzystać do faktoryzacji. Sprawdziłem trzy metody przekształcania liczby n na postać p*F[m].
Idea jest następująca, stosując prawa działające na liczbach Fibonacciego tworzę ciąg 'cyfr' ze zbioru {0,p}, do pozycji F[0] dodaję ewentualną resztę całkowitą r.

1) Liczba nieparzysta p zapisana liczbą Fibonacciego F jako p*F (użyta gramatyka (p|0)*r ),  gdzie p*p<n<p*p+2.
zmniejszamy p = p-2.
Przechodzimy rekursywnie modyfikując 3-cyfrowe okienko [p+a,b,c] = [p,a+b,a+c] lub [a,b,c] = [0,a+b,a+c], a<p.
Okienko przesuwamy w kierunku 'cyfr' mniej znaczących. Kiedy się nie da, dla okienka [a,b,c] stosujemy przekształcenie d=2*b+c+r oraz wypełniamy: [a, d%(2*p), e*p+r'], gdzie r' jest nową resztą d%p, e in {0,1} jest dopasowane, by nie zmienić wartości liczby.
Przykład p=21, w ostatnim okienku mamy [21, 34, 4] = [21,0,72] = [21,21,21+9].
Teraz należy 'naprawić' liczbę Fibonacciego, wracamy z okienkiem posługując się uogólnieniem okienka (1) postaci:
[0,p,p] = [p,0,0]   (2) .
Wielokrotnie to stosując (czasem z nawrotami dla ostatnich cyfr) uzyskamy nową wartość
(p+2)*F +r = p*F'+r'.
Kontynuując ze zmniejszajacym się nieparzystym p znajdziemy dzielnik lub dojdziemy do standardowej liczby Fibonacciego dla n. Wtedy kończymy.

W każdej iteracji wartość p się zmniejsza, F zwiększa. Za wyjątkiem obróbki końcowej cyfry mamy do czynienia tylko z dodawaniem, odejmowaniem i porównywaniem liczb nie większych niż pierwiastek z n. 

2) Jest to algorytm dla liczb F mniejszych niż pierwiastek, wykorzystujący dzielenie. W pierwszym kroku dzielimy p=n/3, uzyskując postać [p,0]+r = p*10_F+r. Sprawdzamy, czy r=p albo r=0, wtedy kończymy mając dzielnik p.
W następnym kroku sprawdzamy, o ile możemy zmniejszyć 'cyfry' p za pomocą dzielenia: p/(F+1) <= x, dla nieparzystego p przyjmujemy za x najbliższą parzystą większą niż iloraz. W następnym kroku odejmujemy p=p-x, co zwiększy nam resztę r o x*F. Podobnie jak w 1), 'naprawiamy' liczbę za pomocą (2).
Porównujemy p z r oraz zmniejszamy p = p-2. Kończymy, gdy p stanie się mniejsze niż liczba Fibonacciego F lub znajdziemy dzielnik p.

Liczba Fibonacciego się zwiększa do długości osiągniętej przy pierwiastku z n, wartości p początkowo bardzo szybko maleją, póżniej słabiej, bo ilorazy dążą do 1, iloraz (całkowity) równy 0 oznacza, że p<F.

3) Jest odwróceniem algorytmu 1). Startujemy z liczby Fibonacciego dla wartości n oraz usiłujemy doprowadzić pierwszą 'cyfrę' okienka do wartości p albo do 0. Należy odpowiednio reagować, gdyż wartości w okienku mają silną tendencję do tworzenia ciągów rozbieżnych. Teraz wygodniej jest czasami nadawać 'cyfrom' wartości ujemne i czasami zciągać cyfrę z już ustalonej pozycji.
Zwiększamy p o 2. Są kłopoty z ustaleniem kryterium stopu, gdyż sprawdzanie, czy p<F za każdym zwiększeniem p jest bezcelowe. Z kolei konieczność 'naprawy' liczby po dotarciu okienka do cyfry najmniej znaczącej zachodzi stosunkowo rzadko.
Przykładowy fragment działania na okienku dla p=67: [65,0,65]=[67,-2,63]; przesunięcie [-2,63,0]=[0,61,-2]; przesunięcie [61,-2,0]=[0,59,61] itd.

Następuje zciąganie liczby Fibonacciego do krótszego ciągu. Operujemy na wzrastajacych wartościach dodawaniem, odejmowaniem i porównywaniem. Osiągniemy pierwiastek z n tylko dla liczb pierwszych.

Liczba w systemie Fibonacciego jest długa, dochodzi do 5 log n, ale przekształcenia i pętle (może z wyjątkiem napraw) są liniowe. Wszystkie przedstawione algorytmy mają porównywalną krotność przypadków do sprawdzenia. Bo albo p się zmienia, albo F.

Są to kolejni kandydaci na wielomianowe algorytmy faktoryzacji, gdyż znalezienie wartości cyfry liczby Fibonacciego można wykonać w czasie logarytmicznym.
Algorytm 1) sprawdza najpierw duże dzielniki, pozostałe najpierw małe.

Faktoryzacja przez iloczyn dzielników

Ten 'algorytm', a właściwie przepis powstał, kiedy badałem, czy można zmniejszyć krotność rozpatrywanych przypadków.
Pierwsze algorytmy mają po sqrt(n)/2 przypadków, publikowane w tym roku (2013) zmniejszyły krotność do 2^(i+1), gdzie 2^(2i)<n<2^(2i+2).
Opisywany dalej przepis może mieć mniej przypadków, ale nie można podać ile. Bardzo silnie zależy ona od podstawy systemu.

Zauważmy, że mając liczby a i b, zwiększając o 1 cyfrę a na pozycji i-tej, iloczyn zmienia się o wartość b*10^i. Podobnie zmniejszając o 1 cyfrę. Pozwala to wytworzyć ciąg zbieżny do wartości złożonej n. Kiedy iloczyn jest za duży, zmniejszamy jeden z dzielników, kiedy za mały, zwiększamy.

Ale należy to robić sprytnie. W systemie dziesiątkowym liczba złożona n mająca duże dzielniki pierwsze ma cyfrę jedności 1, 3, 7, 9. Wykluczone zostały dzielniki 2 oraz 5. Zatem mozemy dobrać a=11, b = floor(n/11), co spowoduje, że iloczyn a*b>n.
Dobierzmy i=1. Zmniejszamy b tak długo, dopóki wartość ilorazu nie będzie mniejsza niż n. Ale zatrzymajmy się, gdy cyfry najmniej znaczące iloczynu oraz n modulo 10^i będą równe. Wtedy zwiększamy i=i+1, oraz kontynujemy, biorąc pod uwagę dwie cyfry n. Tym razem zmniejszamy cyfrę dziesiątek liczby b. Możemy kontynuować z setkami itd.
Kiedy iloczyn stanie się odpowiednio mały: a*b<n, zwiększamy a do najbliższej wartości z cyfrą 1, 3, 7, 9  na pozycji jednostek.Powoduje to zwiększenie iloczynu.
Uzyskujemy ciąg wartości oscylujacy wokół wartości n. Kiedy wreszcie iloczyn stanie się równy n, nasze wartości a i b są dzielnikami.

Testy przekonały mnie, że sposób jest niezależny od systemu, w którym liczymy. Mamy do czynienia tylko z odejmowaniem i dodawaniem liczb w innym systemie, ewentualnie mnożonych przez małe stałe (nie większe niż podstawa systemu).
Zaś w różnych systemach liczność podejrzanych wartości a jest inna.  Podejrzewam, że najgorzej jest w systemach o podstawach będących liczbami pierwszymi (każda wartość prócz 0 może wystąpić), najłatwiej w systemach będących liczbami złożonymi z różnych liczb pierwszych.

Przykładowo, system o podstawie 210 = 2*3*5*7 zawiera w cyfrze najmniej znaczącej informację o podzielności przez dowolną liczbę pierwszą pierwszej dziesiątki. Zaś poszukiwaną postacią liczby a jest liczba pierwsza (dziesiątkowa) mniejsza niż 210, albo jedna z pięciu złożonych wartości: 121, 143, 169, 187, 209.

Rozłożymy liczbę 6767 = 67*101. W systemie o podstawie 210 ma ona postać 32'47_{210}, przyjmijmy (dla wygody) a=10, b=3'46_{210}.
Iloczyn a*b = 10*3'46_{210} = 32'40_{210} = 6760<n, gdyż 10*46 = 460 = 2*210+40
Zwiększając a do 11, otrzymujemy:
a*b = 11*3'46_{210} = 32'40_{210}+ 3'46_{210} = 35'86_{210} = 7436.
Teraz zmniejszamy b, odejmując 11 z 'cyfry' dziesiątek, dopóki 86-k*11<47, albo 'cyfra jedności' b będzie równa 47. Pierwszy warunek będzie spełniony szybciej, po zmianie b o 61_{210}, zatem liczymy
a*b = 11*3'42_{210} = 35'86_{210} - 11*61_{210} = 35'86_{210} - 3'41_{210} = 32'45_{210} = 6765
Znów zwiększamy a do najbliższej liczby pierwszej 13, dodając do iloczynu 2*3'42_{210}, itd. 
Ostatecznie uzyskamy a=67, b=101, skacząc po samych liczbach pierwszych dla a.

Wzory na budowę dzielników, uzupełnienie w dziesiątkowym

Rozważany w ostatnim poście algorytm okazał się niezależny od podstawy systemu. Sprawdził się w binarnym, sprawdził się też w dziesiątkowym.

Funkcja jest postaci
(1)  k(x,n,p,a,b) = -(p+b)/p + n/(p*(px+p+a)
dla potęgi podstawy systemu p, liczby rozkładanej n, odpowiednio przygotowanych wartości a, b spełniajacych warunek
(2)  a*b%p = n%p.
Dzielniki zależne od punktów kratowych (x,k) spełniają zależności
p*(x+1)+a, p*(k+1)+b.

Aby szybciej wyznaczyć dzielniki, z użyciem mniejszych liczb, trzeba odpowiednio przygotować a oraz b.
Pokażę to na przykładzie n = 8934053 = 1087*8219. Wartość p = 10.
Najmniej znacząca cyfra n to 3, szukam w tabliczce mnożenia iloczynów, które kończą się na 3. W systemie dziesiątkowym znajduję 4 takie. Ze względu na symetrię dwa można odrzucić. Są to pary (1,3) oraz (7,9).
Jeśli podstawa systemu jest liczbą pierwszą, takich par jest tyle, ile niezerowych cyfr - połowa się powtarza ze względu na symetrię. Ale tylko w tym kroku. W dalszych symetria najczęściej jest zaburzona oraz trzeba sprawdzać wszystkie przypadki (10 w dziesiątkowym).

Mamy już pary: 1*3 = 3 oraz 7*9=63. W podanym we wcześniejszym poście przebiegu powinniśmy odjąć je od n. Lecz teraz odcinamy najmniej znaczącą cyfrę zarówno od n, jak i od tego iloczynu. Uzyskamy odpowiednio e=0 oraz e=6, oraz n=893405.
Modyfikujemy sposób wyznaczania a, b oraz e. Zmieni to nieco sposób wyznaczania punktów kratowych.

Załóżmy, że wykorzystujemy parę (a,b) = (7,9), oraz przechodzimy do następnej cyfry: p = p*10 = 100, n=89340. Dołączymy cyfry bardziej znaczące do a, b. Zanim to jednak zrobimy, zapamiętamy wartość dziesiątek e=6 z iloczynu 63.
Dokładamy do a jako cyfrę najbardziej znaczącą kolejno c=0, 1, ..., 9, dopasowywując d w taki sposób, aby powstałe liczby
a*b = (c##7) * (d##9) %100 = 53.
Interesuje nas tylko cyfra dziesiątek 5. Zatem obliczenie sprowadza się do prostej kongruencji liniowej: 9c+7d+e=5, czyli 9c+7d+1=0.
Na przykład dla c=0 powstaje a=07, dalej 0*9+7*d+1=0 skąd d=7 oznacza, że b=79. Rzeczywiście 7*79=553.
W ogólności mnożymy dokładaną cyfrę do a przez cyfrę najmniej znaczącą b, dodajemy doń iloczyn cyfry dokładanej do b z najmniej znaczącą cyfrą a, fragment iloczynu wcześniejszych cyfr trzymaną w e oraz iloczyn dokładanych cyfr, przesunięty o p. W tym przypadku mamy
0*9+7*7+6+0*7*100 = 56.
Odcinamy cyfrę najmniej znaczącą 6 przyjmując nową wartość e=5. Jest to całość z dzielenia a*c przez p.

W kolejnym przypadku dokładamy do a=7 cyfrę 1, uzyskujemy 1*9+7*d+1=0 skąd d=0 oraz wartości a=17, b=09, oraz e = (1*9+0*1+6)/10 = 1.

Zamiast wyznaczać iloczyn a*b, dodamy do siebie aktualne wartości
g = a+b+e+1*p
U nas to będzie g = 07+79+5+1*100 = 191.
Sprawdzanie przypadków zmieniło się, należy sprawdzić podzielność (n-g)%(1##a), (n-g)%(1##b). W naszym przypadku (89340-191)%109 oraz  (89340-191)%179. Punkty całkowite są punktami kratowymi krzywej (1), które wkładamy do równań na dzielniki.
Do następnej iteracji z p=1000 wartość n przejdzie jako n=8934.
Wtedy z przybliżenia dzielników powstałych z a=87, b=19 dla c=0 uzyskamy
0*9+2*d+6=0, skąd d=2, b=219, e=(0*19+2*87+16+0*2*100)/10 = 19;
g = 087+219+19+1*1000 = 1325;
8934-1325 = 7*1087 wskazuje dzielniki, co kończy obliczenia.

Oczywiście, wartości a=07, b=79 nie przybliżają dzielników, robią to wartości a=87, b=19. Ale i tak po przerobieniu co najwyżej 220 przypadków (440 dzieleń) będziemy znali dzielniki.
Dzielenie przez wartości nieparzyste to 542 przypadki.

Wzory na budowę dzielników

Kolejny algorytm faktoryzacji opisywany w tym blogu powstał dzięki odwróceniu kolejności w geometrycznym mnożeniu w systemie binarnym.
Mnożenie geometryczne polega na zapisaniu cyfr dzielników przy brzegach prostokąta, wewnątrz zapisujemy iloczyny tych cyfr, przy czym cyfra dziesiątek jest 'na innej przekątnej' niż cyfra jedności. Następnie sumujemy po przekątnych, np. wartość na cyfrze dziesziątek powstaje z sumy:
a_1*b_0 + a_0*b_1 + (a_0*b_0)/10.
Po uwzględnieniu przeniesień mamy wynik.
Sposób ten zastosowany w systemie binarnym wskazuje prostą zależność między bitami dzielników, co pozwala wyznaczyć te dzielniki, niestety, tylko na 50%. Ale wystarczy.
Połączyłem te informacje z moim pierwszym algorytmem, w którym zakładałem, że dzielniki sa odpowiedniej postaci (2^i*c+a) i sprawdzałem, czy jest tak rzeczywiście.

Niektóre informacje okazały sie nadmiarowe, i tak oto powstał wzór na funkcję dzielników k(x) liczby n, przy potędze dwójki i:
(1)    k = -(i+b)/i + n/(i (ix+i+a) )
W tym wzorze n jest rozkładaną liczbą, i jest potęgą 2, szukamy dzielników większych niż i, wartości a oraz b są resztami z dzielenia n przez i, spełniającymi zależność bitową:
(2)   (n - a*b) AND (i-1) = 0
Oznacza to, że a*b oraz n mają dokładnie takie same bity najmniej znaczące, czyli stanowią swoiste 'bitowe przybliżenie' dzielników.

Po przekształceniu wzoru (1) do wspólnego mianownika uzyskujemy funkcję holomorficzną postaci 
(3)   k = (-bx+q)/(ix+a),
w której i jest potęgą 2, q powstaje z ilorazu różnicy n oraz a*b, co jest liczbą naturalną bliską n*i^(-2). Dzielniki znajdujemy, gdy funkcja ta ma punkt kratowy (x,k), tzn. x, k są całkowite (interesują nas tylko nieujemne).

Wspomniane wzory na dzielniki:
    i*(k+1)+b
oraz
    i+x*i+a

Przykład: n=8934053, i=1024, znaleziona para (a,b) = (27,63) spełnia zależność (2), ich iloczyn w szesnastkowym to 0x2A5, oraz n jest postaci n = 0x400*d + 0x2A5.
funkcja (1) dla tego przypadku:
k = -(1024+63)/1024 + 8934053/(1024(1024x+1024+27))
przekształca się do
k = -1087(x-7)/(1024x+1061)
Punkt kratowy (x,k) = (7,0) wstawiamy do wzorów na dzielniki:
1024*(0+1)+63 = 1087
1024+7*1024+27 = 8219
Sprawdzamy, że rzeczywiście 8934053 = 1087*8219.

Mając postać (3), wystarczy przeglądać x naturalne w poszukiwaniu punktów kratowych, licznik maleje, mianownik rośnie, wartości są nie większe niż n*i^(-2). Im większe i, tym mniej wartości mamy do sprawdzenia, gdyż dla x rosnącego iloraz zmniejsza się wykładniczo. Kiedy zejdzie poniżej 0, nie znajdziemy dzielników. Ale trzeba spradzić parę początkowych wartości, żeby dzielnik nie umknął.

Największym kłopotem może być dopasowanie pary (a,b). Ale zaczynając od i=4, gdy bit na pozycji 'dziesiątek' liczby n jest ustawiony (n&2=1), a=1, b=3. Gdy bit ten jest zgaszony (n&2=0), mamy dwie możliwości a=b=1 albo a=b=3.
Teraz zwiększajac i = 2*i, sprawdzamy kombinacje par (a,b), (a+i,b), (a, b+i), (a+i,b+1). Dwie z nich mają takie same bity najmniej znaczace jak n - te zostawiamy, pozostałe dwie mają znak różny i je usuwamy. W ten sposób tworzy nam się drzewo binarne przypadków, w którym co najmniej jedna gałąź wskazuje rozkład. Dla liczb pierwszych jest to (a,b) = (1,n).

Porównywanie ułamków zwykłych za pomocą ciągu Fareya

Porównywanie ułamków zwykłych jest zagadnieniem z klasy szóstej podstawówki.
Sprowadza się ułamki do wspólnego mianownika i porównuje liczniki. Komputerowo spradza się znak 'iloczynu na krzyż' (licznik 1 * mianownik 2  - licznik 2 * mianownik 1).

Można to jednak zrobić bez pomocy mnożenia, dzielenia.

Co to jest ciąg Fareya. Mając dwa ułamki a/b oraz c/d, tworzymy nowy ułamek, którego licznik powstaje przez dodanie liczników, zaś mianownik jako suma mianowników. Nowo powstały ułamek (a+c)/(b+d) rozdziela ułamki, z których powstał w porządku liczb wymiernych.
Ciąg ten może służyć jako uzyskanie ciągu wszystkich liczb wymiernych przedziału [0/1; 1/1] przez wstawianie powstałych ułamków między już dane:
0/1 < 1/2 < 1/1;
0/1 < 1/3 < 1/2 < 2/3 < 1/1; itd.

Zastosowanie ciagu Fareya do porównywania ułamków a/c, b/d można zastosować wtedy, gdy licznik i mianownik jednego z ułamków są większe niż drugiego a<c, b<d. Wtedy traktujemy je jako pierwsze elementy ciągu Fareya, oraz wyliczamy trzeci następująco (c-a)/(b-d). Ten nowy ułamek ma wartości stosunkowo małe, i może być porównany z a/b, czasem nawet z 1/1. Uzyskany porządek przenosi się jako porównanie początkowych ułamków.

Kiedy nie można zastosować podanego kryterium, stosujemy szacowania:
- ten ułamek o równym mianowniku jest większy, gdy jego licznik jest większy;
- ten ułamek o równym liczniku jest większy, gdy jego mianownik jest mniejszy;
W ten sposób znajdujemy porządek ułamków np. 7/18 oraz 8/13, 7<8 oraz 18>13, zatem 7/18 < 8/13.

Porównamy 371/1243 z 721/1571. Skorzystamy z ciągu Fareya obliczając następny element
(721-371) / (1571-1243) = 350/328 > 1, bo 350>328
Zatem porządek jest następujący:
371/1243 < 721/1571 < 350/328

Przekształcenie Fareya zmniejsza wartości porównywanych ułamków. Można zapisać to procedurą rekursywną, u1 oznacza ułamek właściwy u1(l1,m1), l1<m1:

int cmpfr( u1, u2 ) {
  if( u1!= u2 && l1<=l2 && m1>=m2 ) return u1<u2;
  if( l1<l2 && m1<m2 ) {
    if( u3(l2-l1, m2-m1) ) return cmpfr( u1, u3 ); // ulamek niewlasciwy
    else return u1<u2;
  }
  return u1>=u2;   
};

Przekształcenia przed dzieleniami - odczyt bitowy dzielników

Jeden z algorytmów mnożenia zastosowany w systemie binarnym wskazuje, że znajomość najmniej znaczących bitów dzielnika pozwala z odpowiedniego bitu iloczynu odczytać kolejny bit obu z dzielników. Są dwie możliwości, między innymi ze wględu na prawo przemienności.

Zatem podejście jest takie. Liczbę rozkładaną n zapisuję w postaci czwórki liczb (a, b, r, c), gdzie a oraz b są kandydatami na dzielniki, do których dołączam przewidywany najbardziej znaczący bit. Liczba r to rezerwa wskazująca, jak wartość jest jeszcze dostępna, gdy wykorzystam informacje o bieżących ustawieniach liczb a oraz b. Wartość c jest potęgą 2, wskazuje, który bit jest aktualnie dopasowywany.

Przebieg. Najpierw przekształcam wartości czwórki, następnie sprawdzam pary wartości oraz wywołuję rekursywnie obie możliwości.

Inicjacja czwórki jest następująca:
(a=1, b=1, r=n-1, c=1)
oraz oznacza, że n = r*c + a*b = (n-1)*1 + 1*1 = n-1 +1.
Teraz w zależności od najmniej znaczącego bitu r, mam dwie możliwości:
bit r&1 wyzerowany:
a'=a, b'=b, r'=r;
albo
a''=a+c, b''=b+c, r''=r-a-b-c;
bit r&1 ustawiony:
a'=a, b'=b+c, r'=r-a;
albo
a''=a+c, b''=b, r''=r-b;

Po tych przekształceniach przyglądam się własnościom par liczb (a',r') oraz (b',r'). Interesuje mnie zwłaszcza podzielność, np. a'|r'.
Wyodrębniłem następujące warunki, stosując parę (d,r):
- d=1, wskazuje rozkład trywialny n = 1*n;
- d|r, sprawdzam teraz, czy d|n. Jeśli tak, d jest dzielnikiem liczby n, jeśli nie, liczba n jest pierwsza.  W obu przypadkach mogę zakończyć rozkład;
podobnie należy sprawdzić, kiedy r staje się zerem;
- d>r, rezerwy są zbyt małe, aby móc zwiększać kandydata na dzielnika, ta gałąż rekurencji jest bezużyteczna, można uciąć też drugą gałąź z tego węzła.

Wywołanie rekursywne dzieli rezerwę na 2 oraz przesuwa c na kolejny bit bardziej znaczący, co w praktyce oznacza wywołanie
(a', b', r'/2, 2*c) oraz (a'', b'', r''/2, 2*c).

Zachowanie się wartości.
Liczby a oraz b w kolejnych wywołaniach stanowią ciągi niemalejące, przyrost jest potęgą 2.
Liczba r stanowi ciąg malejący, ograniczony z góry przez n*2^(-i), i jest licznikiem wywołania. Dzielna jest też ograniczona przez tę wartość, zaś dzielnik dosyć szybko rośnie. Przypadek podzielności występuje parzystą krotność razy, dla każdego z dzielników, także dla liczb pierwszych i nie zawsze wskazuje kandydata na dzielnik trywialny.
Kiedy a stanie się większe niż r, rekursja się kończy, co powoduje, że kroność wywołania stanowi połowę logarytmu binarnego z n. Sprawdzamy 2^(1+1/2 lg n) możliwych ustawień bitów w dzielnikach.

Przykład: 143
czwórka inicjująca (1, 1, 142, 1),
pierwsze wywołanie rekursji zwraca (1, 1, 71, 2)
Ponieważ ~2|71, mamy przypadki czwórek (1, 3, 71-1=70, 2) oraz symetryczną (3, 1, 70, 2), 1|70 i jest przypadkiem trywialnym;
Kolejne wywołania rekursywne zwracają (1, 3, 35, 4) oraz (3, 1, 35, 4).
W obu przypadkach ~2|35, zajmijmy się pierwszym:
mamy (1, 7, 34, 4) oraz (5, 3, 32, 4);
wywołania rekursywne (1, 7, 17, 8) oraz (5, 3, 16, 8);
pierwsze z nich daje (1, 15, 16, 8) oraz (9, 7, 10, 8);
po kolejnym wywołaniu rekursywnym obydwu przypadków uzyskujemy r mniejsze niż któryś z kandydatów na dzielniki, należy zająć się (5, 3, 16, 8).
Mamy 2|16, czyli czwórkę (5, 3, 16, 8), albo (13, 11, 0, 8),
Wartość r=0, sprawdzamy, czy 13|143, i rzeczywiście, mamy rozkład 143=13*11.

Faktoryzacja, przepis z prostymi przekształceniami

Skrzyżowałem faktoryzację sprawdzającą równocześnie małe i duże dzielniki z ostatnią faktoryzacją, w której sprawdzam zaledwie 2^i przypadków zamiast pierwiastka z liczby faktoryzowanej. Uzyskany algorytm cechują działania na zaledwie czterech liczbach osiągających wielkość pierwiastka z liczby rozkładanej.

Przepis jest następujący: najpierw dzielę liczbę n na dwa, w przybliżeniu równe kawałki, które dostarczą mi postać liczby trójcyfrowej liczby n. Mogę teraz zatrudnić trzy wątki. Jeden stosuje konwersję o 2, podczas której podstawa rośnie do pierwiastka z n, drugi stosuje konwersję o -2 sprawdzając dzielniki bliskie sobie, które są nieosiągalne przez pierwszy z wątków. Wreszcie trzeci przygląda się podstawie, oraz uwzględniając jej dzielnik trzy (powtarzający się co trzy iteracje w każdym z dwu wątków), eliminuje małe dzielniki przy współudziale wyrazu wolnego.

Ze względu na specyfikę obliczeń, program wymaga pamięci do korzystania z kopii niektórych wartości. Dostępne są następujące wartości:
a, b, c, p,
spełniające warunek n = (a*p+b)*p+c, indeksy podzielności p przez 3 (aby nie obliczać p%3 w każdej iteracji) oraz pewna liczba pomocnicza. 
Kryterium zatrzymania stanowi wartość x, uzyskana pod koniec pracy pierwszego wątku. Zatrzymuje ona działanie drugiego wątku ograniczając liczność sprawdzanych przypadków. Trzeci wątek po znalezieniu dzielnika zatrzymuje cały program.

Przejdźmy do szczegółów, które przedstawię w taki sposób, jakby były omawiane po raz pierwszy.

  NAPRAWA
Najpierw sposób naprawy postaci [a, b, c, p]. Spełniane są warunki: 0<a<4, -1<b<p, -1<c<p oraz n = (a*p+b)+c.
Jeśli b jest ujemne, dodajemy p do b oraz odejmujemy od a jedynkę:
if( 0>b ) { b+=p; a--; }
Jeśli c jest ujemne, dodajemy p do c oraz odejmujemy od b jedynkę:
if( 0>c ) { c+=p; b--; }
Jeśli b jest nie mniejsze niż p, zmniejszamy b o p oraz zwiększamy a:
if( b>p-1 ) { b-=p; a++; }
Podobnie c, jeśli jest zbyt duże, zmniejszamy je o p zwiększając b:
if( c>p-1 ) { c-=p; b--; }
Operacje te wystarczy powtórzyć do dwu razy w dowolnym z wątków przy uruchamianiu funkcji naprawa(), aby zapewnić wymagane warunki.

  ALGORYTM
Przebieg algorytmu, liczbowe dane są z powietrza:
Podział liczby n spełniajacej warunek 2^(2i) < n < 2^(2n+2) na dwie liczby mniej więcej tej samej długości spełniające warunki:
p=2^i lub p=2^(i+1); b = n/p; c=n%p, b>p, n = b*p+c
np. dla uzyskanego p=1024, b nie przekracza 4096, zaś c jest ograniczone w przedziale [0, 1023).
Z liczby b wydobywamy a=b/p; b=b%p. Uzyskujemy wtedy wstępną czwórkę liczb postaci:
[a, b, c, p] , n = (a*p+b)*p+c, a<4, b<p, c<p

Wartość podstawy p powinna być nieparzysta, i najlepiej podzielna przez 3.
Jeśli p jest parzyste, zwiększymy podstawę p o 1 przekształceniami:
{
p++;
b-=a;
naprawa();
c-=b;
b-=a;
naprawa();
}

  WĄTEK 1
pętla zwiększajaca podstawę p, zmniejszająca skokowo a. Kończymy, gdy a osiagnie 0, gdyż wtedy wartość podstawy p będzie większa niż pierwiastek kwadratowy z n. Ciąg [a,b,c,p] w tym momencie może wyglądać następująco: [1, 7, 482 953, 3 593 385],
[1, 3, 482 943, 3 593 387],
[0, 3 593 388, 482 941, 3 593 389]
Wyjście z programu następuje wtedy, gdy c będzie zerem. Znamy wtedy dzielniki: n = p * (a*p+b).
{
p+=2;
b -= 2*a;
naprawa();
c -= 2*b;
b -= 2*a;
naprawa();
}
jeśli p jest liczbą podzielną przez 3 (co sprawdzamy odpowiednim licznikiem) przesyłamy parę [c,p] do wątku 3. Wartość taka powtarza się co trzy iteracje, np. dla p=17 278 137, później dla p+6: 17 278 143 itd.

  WĄTEK 2
pętla zmniejszająca podstawę p, zwiększajaca skokowo a. Kończymy, gdy trzeci z wątków nakaże zakończyć.
Wyjście z programu następuje jak przy wątku 1, gdy c będzie zerem. Mamy wtedy dzielniki: n = p * (a*p+b). Pętla:
{
p-=2;
b += 2*a;
naprawa();
c += 2*b;
b += 2*a;
naprawa();
}
Jeśli p jest liczbą podzielną przez 3 (znowu licznik), przesyłamy parę [c,p] do wątku 3 podobnie jak przy wątku piewszym.

  WĄTEK 3
Wątek ten sprawdza małe dzielniki, oraz zatrzymuje algorytm, kiedy zostanie sprawdzone 2^i<sqrt(n) przypadków. Bazuje na obserwacji:
liczba n ma dzielnik d, gdy podstawa systemu p oraz wyraz wolny c są wielokrotnościami d.
Zatem najpierw zmniejszamy kopię podstawy p dzieląc ją przez 3, a następnie liczymy nwd(p,c) by znaleźć wspólny dzielnik. Gdyż jeśli dzielnikiem jest wielokrotność 3, to dzielnikiem jest też 3.
{
while ( 0==(p%3) ) p/=3;
dzielnik d = nwd(p,c); // jeśli 1==p, stosujemy p=3
}
jeśli przesłane zostanie p z wątku pierwszego, zapamiętywana jest wartość x = p/3. Jeśli wątek drugi prześle p mniejsze niż x, oznacza to, że ten przypadek został już sprawdzony. Kiedy pierwszy z wątków zakończy swoje działanie, jego ostatnia wartość p/3 stanowi równocześnie liczbę, do której wystarczy sprawdzać w wątku 2, aby wyeliminować wszystkie możliwe dzielniki.

Podsumowanie
Wątki 1 oraz 2 sprawdzają bardzo duże dzielniki, nieco mniejsze niż pierwiastek z rozkładanej liczby. Wątek 3 sprawdza dzielniki mniejsze niż 2^i< sqrt(n), oraz zapobiega ponownemu przeliczaniu. Wątek 2 przy zmniejszaniu podstawy od pewnego miejsca zaczyna podawać dzielniki sprawdzone już przez wątek 3, czemu należy zapobiec.
Wątek 1 powinien być dominujący nad drugim, na ogół ma najwięcej przypadków do sprawdzenia. Wątek 3 jest najbardziej czasochłonny, wymaga najgorszych operacji, oraz może korzystać z większej krotności procesorów.

Wartości liczbowe w wątku pierwszym: a maleje, b i c skokowo maleją, p rośnie jednostajnie.
Wartości liczbowe w wątku drugim: a rośnie, b i c skokowo rosną, p maleje jednostajnie.
Wartości liczbowe w wątku trzecim: jedyne co można przewidzieć, że p nie przekroczy pierwiastka kwadratowego, zaś c jest ograniczone przez p. Wartość c może być stosunkowo mała. Przy kolejnych dostawach p przez wątki, krotność dzieleń przez 3 zachowuje się jak grzebień: 1, 2, 1, 3, 1, 2, 1, 4, 1, 2, 1, 3, 1, itd.

Faktoryzacja metodą (nie kolejnych) dzieleń

Liczbę złożoną nieparzystą n = p*q można przedstawić za pomocą przedstawienia liczby q = (q_m,..., q_1, q_0) w systemie binarnym następującą tablicą wartości:
n = (p*q_m, ..., p*q_1, p*q_0).
Postać ta charakteryzuje się zerami na miejscach, w których q_i=0, oraz wartościami p na miejscach, w których w zapisie binarnym q_i=1, i=0..m. Traktując poszczególne pola tablicy jako wielomian W[2], uzyskamy wartość liczby n.
Np. n = 125 możemy przedstawić jako n = (5, 5, 0, 0, 5), gdyż 125 = 5* 25, oraz 25 = 11001b. Wartość jest równa
W[2] = ((((5*2+5)*2+0)*2+0)*2+5 = 125.

Czy podana postać może służyć do znajdywania rozkładu liczby na czynniki?
W jaki sposób?
Okazuje się, że postać ta wskazuje kolejne sposoby faktoryzacji. Dysponując tabelą zapisu binarnego liczby n, za pomocą operacji niezmienniczych
n_i = 2 * n_{i-1}, i=1,...,m
można doprowadzić do postaci binarnej dzielnika.
Postępowanie zaproponowane nieco dalej sprowadza się do dzieleń przez kolejne liczby nieparzyste, lecz w nieco innej kolejności niż normalnie, z mniejszymi wartościami dzielnych. Pierwotnie próbowałem szukać, czy znajdę sposób na domyślenie się postaci binarnej dzielnika.

Pierwsze próby wyrównywania wartości w zapisie binarnym liczby n czasem skutkowały bardzo szybkim znajdywaniem dzielnika, czasem zaś wykazywały, że ten nie istnieje. Kłopot polegał na tym, że żaden algorytm nie jest w stanie przewidzieć, które pola tabeli należy wyzerować, a które zostawić ustawione, czyli przewidzieć budowę czynnika q. Należy to robić losowo, algorytmem genetycznym lub sprawdzać wszystkie możliwe przypadki.
Algorytm genetyczny ma duże szanse zapętlenia. Dodatkowo wprowadza bardzo dużo operacji wyrównujących poszczególne wartości tabeli.

Lepiej jest zastosować kombinatorykę, wtedy żegnamy się z dopasowywaniem bitów stosując wykładniczą metodę siłową.
Naszym celem jest przygotować tabelę wypełnioną równomiernie daną wartością dodatnią. Następnie zerować poszczególne pola tablicy, starając się zachować taką samą wartość na wszystkich ustawionych pozycjach.
Jeśli się to uda, dodatnia wartość w tablicy jest wartością jednego z dzielników, zaś budowa tablicy potraktowana jako liczba binarna drugim.

Praktycznie wprowadzamy komórkę reszty, do której odkładamy wszelkie nadmiary, powstały z dzielenia pierwotnego, jak i z zerowania pól tablicy.
Tabelę wypełnimy równomiernie jedną wartością b = floor(n : (2^i-1)). Dzielimy tu przez np, 3, 7, 15, 31, 63, ... w zależności od ewentualnej długości dzielnika, którego szukamy.
Pierwsze i ostatnie pole są niepuste. Przygotujemy szablon będący liczbą binarną s = (2^i-1).
Wyzerowanie bitu szablonu na pozycji k odpowiada wyzerowaniu pola tablicy, równocześnie wartość b*2^k jest dodawana do reszty. Można ją zapisać instrukcją w C++ dla maski bitowej t:
r += (s & ~t)*b.
Jeżeli teraz reszta r jest wielokrotnością szablonu (s & t), wartość reszty można bez kłopotu rozdzielić między niepuste pola tablicy, o co nam chodziło. Praktycznie sprawdzamy, czy wartość liczbowa dzieli resztę (s & t) | r. Jeśli tak, mamy do czynienia z dzielnikami (s&t) * (b+r/(s&t)) = n.

Powtarzamy dla dowolnej wariacji ustawienia pola tablicy z wyjątkiem wartości skrajnych, lub do znalezienia dzielników.


Ostatecznie dzielimy najpierw liczbę n przez s=(2^i-1) dla pewnego i, a następnie zerując odpowiednie bity szablonu dzielimy resztę przez szablon. Reszta jest liczbą nie przekraczającą połowy n, zaś liczby zadane przez schemat mają ściśle określoną długość binarną. Można też wykorzystać reszty modulo schemat.
Można się jeszcze pokusić o to, by liczby zadane szablonem nie przekroczyły pierwiastka kwadratowego z n. Zmniejszy to krotność obliczeń.

Pierwsze z dzieleń jest tak specyficzne, że można wykorzystać przyspieszone dzielenie przez odpowiedniki dziewiątki. Polega ono na pocięciu liczby na paczki długości i cyfr, na wynik składają się zapisane kolejno wartości sumy narastającej paczek.  

Przykład: 1441 = 11 * 131
Dzielimy 1441 : 3 = 480 z resztą 1, nic dalej nie można zrobić
Dzielimy 1441 : 7 = 205 z resztą 6
schemat 101b każe sprawdzić dzielenie
(2*205+6) : 101b = 416 : 5 \equiv 1 (mod 5)
Dzielimy 1441 : 15 = 96 z resztą 1
Schemat generuje 3 dzielenia przez odpowiednio 1001b, 1011b oraz 1101b.
W pierwszym przypadku mamy dzielenie
((2+4)*96+1) : 1001b = 577 : 9 \equiv 1 (mod 9),
w drugim
(4*96+1) : 1011b = 385 : 11 \equiv 0 (mod 11)
oraz dzielniki 11 oraz 96+385/11 = 96+35 = 131; w trzecim
(2*96+1) : 1101b = 193 : 13 = 11 (mod 13) .

Następuje zamiana dzieleń:
1441 : 3               1441:3
1441 : 5               416 : 5
1441 : 7               1441:7
1441 : 9               577 : 9
1441 : 11             385 : 11
                           193 : 13           

Sprawdzanie między potęgami 2, czy można przyspieszyć

Ostatni algorytm, którego nie podałem w postaci informatycznej przeszukuje 'stożek' wartości mniejszych niż liczba faktoryzowana n, którego podstawą jest 2^i, oraz 2^(2i)<n<2^(2i+1).

Zastanawiało mnie, czy muszę przeszukiwać cały ten 'stożek'. Może istnieje 'ścieżka', która doprowadzi bezpośrednio do dzielnika.
Niech n = a*b+c, gdzie a<b jest potęgą 2, oraz c<a.
Stosowana konwersja podwaja a oraz odpowiednio zmniejsza b. Wartość c zachowuje się jak chce, nie przekraczając a. Dla wartości skrajnych tworzy ciąg niemalejący.
Do dalszej iteracji wybierałem ten z przedziałów [a, a+1], [a+1, a+2] w którym spełniony był  warunek b1 < 2^k < b2, k<i.
W ten sposób dochodziłem do tego, że b po którejś iteracji stawało się potęgą 2. 

Często znajdowałem pierwiastek, ale tylko wtedy, gdy był on bliski potęgi 2. Kiedy był dalej, należało przeliczyć reszty w krotności podwojonej ostatniej znalezionej reszty. Była to często wartość bliska podstawy 'stożka'.

W szczególności dla liczb pierwszych należało sprawdzać wszystkie wartości 'stożka'. Podobny obraz dawała także liczba złożona, w której dzielniki były najdalej 'odsunięte' od potęg 2.

Zatem ten sposób nie pozwala przyspieszać, należy przerachować wszystkie wartości parzyste 'grzebienia' nie większe niż 2^(i+1). 
Dlatego 'grzebienia', gdyż sprawdzam wszystkie liczby postaci 2^k*d przedziału [2^i, 2^(i+1)), gdzie -1<k<i oraz d jest liczbą nieparzystą. Liczb takich jest tyle ile wynosi k+1, oraz 'tną' one pionowo 'stożek' na warstwy.
Zresztą, kto popatrzy na liczność przypadków przykładu poprzedniego posta, będzie wiedział, skąd nazwa 'grzebień'.

Ostatni algorytm faktoryzacji, jeszcze raz z przykładem i złożonością

-->

Opiszę jeszcze raz dokładnie ostatni algorytm faktoryzacji, z próbą policzenia złożoności. 

-->

Nieparzystą liczbę rozkładaną n spełniającą warunek 2²ⁱ < n <2²ⁱ⁺¹ dla pewnego i zapisujemy jako iloczyn dwu wartości n = a*b+c, gdzie a²=(2ⁱ)² < n, c<a.

Szukamy z założenia dzielników nieparzystych wykorzystując konwersję do podstawy parzystej o 2 większą, eliminując dzielnik 2 na początku algorytmu. Liczby złożone n mające dzielniki z przedziału (2ⁱ, 2ⁱ⁺¹) mają następującą postać:

a*b+a = a*(b+1)

dla podstaw p=b+1 parzystych. Mamy 2^i-1 przypadków złożoności konwersji zmiany systemu o 2.

Przygotowanie do faktoryzacji to przekształcenie liczby n w trójkę [a,b,c], w której a=2ⁱ jest największe możliwe, b jest częścią całkowitą n/a oraz c jest resztą z dzielenia n/a.

Jeden z procesorów konwertuje liczbę n na systemy o coraz większych podstawach, oraz przesyła postać [a, b, c] drugiemu procesorowi. Drugi dokonuje konwersji na systemy o coraz mniejszych podstawach.

Konwersja liczby n = [a,b,c], gdzie n = a*b+c na system o 2 większy. Praktycznie cyfra setek oraz cyfra dziesiątek są złączone w jedną liczbę dla zastosowania tożsamości

[a, b, c] = [a+2, b-e, c-2b+e(a+2)]

oraz a:=a+2; b:=b+e; c:=c-2b+e(a+2). Wartość e jest dobierana empirycznie w taki sposób, by c-2b+e(a+2) miała wartość dodatnią nie większą niż a+2. Złożoność konwersji zwiększania systemu o 2.

a = a+2, złożoność bitowa i = O(log n), gdyż a ma długość i oraz i=log n;

c = c+[(a+2)]^*-b-b, dwa odejmowania oraz do 6 dodawań, złożoność nie przekracza 6*log n, zwiększanie e jest w czasie stałym bitowo; dochodzą porównania możliwości odjęcia przy zachowaniu znaku, wykonywane w czasie liniowym bitowo, razem O(log n);

b = b-e, jedno odejmowanie małej liczby zawierającej co najwyżej 3 bity, prawie stałe;

a==c, porównania wykonywane w czasie liniowym O(log n).

Podsumowując, konwersja wykonywana jest O( log n * 2^{-1+log n}) = O(n log n) bitowo.

Konwersja liczby n = [a, b, c] gdzie n = a*b+c na system dwukrotnie mniejszy. Postać liczby a = 2^kd, gdzie 2 i d są względnie pierwsze.

[2*a, b, c] = [a, 2*b + f, c-f*a] .

Współczynnik sterujący f przyjmuje wartość 1, kiedy c>a, w przeciwnym przypadku 0. Procedurę kontynuujemy, dopóki a będzie parzyste. Krotność pętli jest równa k < i = log n.

Wyliczona wartość c jest porównywana z zerem, gdyż wtedy a oraz b są dzielnikami. Mamy do czynienia z dużą licznością przypadków, ale istnieje ograniczenie, gdyż szereg 1 + ½ + ¼ + ... +2^-i jest zbieżny do 2. Współczynniki szeregu wskazują liczność przypadków z k=1, k=2, ..., k=i odpowiednio dla wszystkich pobranych argumentów przedziału [2ⁱ, 2ⁱ⁺¹).

Złożoność konwersji podwajania podstawy.

a = a/2, przesunięcie bitowe liczby mniejszej niż pierwiastek z n, O(log n)

b = 2b + f, przesuniecie bitowe liczby zwiększanej do n, z ewentualnym dodaniem 1, O(n)

c = c-f*a, przy ustawionej fladze odejmowanie O(log n), w przeciwnym razie O(1).

Podsumowując, konwersja ta wykonywana jest k razy dla danych wejściowych, liczność z wykorzystaniem pierwszej z konwersji sumuje się do O( n 2^i-1 ) * O( 2*2^i-1*n ) = O( n² log n * 2^i-1*2*2^i-1 ) = O( n² log n 2^2i-1 ) = O(n² (log n)³ ) = O( 2^{2log n} (log n)³ ) .

Przykład n = 8934053, rozkład całkowity

Początek każdej linii: wyrażenie powstałe z konwersji z systemu o 2 większej, po ‘=’ konwersje z połowioną podstawą.

Inicjacja (znajdowanie najmniejszej potęgi 2ⁱ większej niż pierwiastek):

1*8934053+0 = 2*4467026+1 = 4*2233513+1 = 8*1116756+5 = 16*558378+5 = 32*279189+5 = 64*139594+37 = 128*69797+37 = 256*34898+165 = 512*17449+165 = 1024*8724+677 = 2048*4362+677

2048*4362+677 Start algorytmu

2050*4358+153 = 1025*8716+153

2052*4353+1697 = 1026*8707+671 = 513*17415+158

2054*4349+1207 = 1027*8699+180

2056*4345+733 = 1028*8690+733 = 514*17381+219 = 257*34762+219

2058*4341+275 = 1029*8682+275

2060*4336+1893 = 1030*8673+863 = 515*17347+348

2062*4332+1469 = 1031*8665+438

2064*4328+1061 = 1032*8657+29 = 516*17314+29 = 258*34628+29 = 129*69256+29

2066*4324+669 = 1033*8648+669

2068*4320+293 = 1034*8640+293 = 517*17280+293

2070*4315+2003 = 1035*8631+968

2072*4311+1661 = 1036*8623+625 = 518*17247+107 = 259*34494+107

2074*4307+1335 = 1037*8615+298

2076*4303+1025 = 1038*8606+1025 = 519*17213+506

2078*4299+731 = 1039*8598+731

2080*4295+453 = 1040*8590+453 = 520*17180+453 = 260*34361+193 = 130*68723+63 = 65*137446+63

2082*4291+191 = 1041*8582+191

2084*4286+2029 = 1042*8573+987 = 521*17157+466

2086*4282+1801 = 1043*8565+758

2088*4278+1589 = 1044*8557+545 = 522*17115+23 = 261*34230+23

2090*4274+1393 = 1045*8549+348

2092*4270+1213 = 1046*8541+167 = 523*17082+167

2094*4266+1049 = 1047*8533+2

2096*4262+901 = 1048*8524+901 = 524*17049+377 = 262*34099+115 = 131*68198+115

2098*4258+769 = 1049*8516+769

2100*4254+653 = 1050*8508+653 = 525*17017+128

2102*4250+553 = 1051*8500+553

2104*4246+469 = 1052*8492+469 = 526*16984+469 = 263*33969+206

2106*4242+401 = 1053*8484+401

2108*4238+349 = 1054*8476+349 = 527*16952+349

2110*4234+313 = 1055*8468+313

2112*4230+293 = 1056*8460+293 = 528*16920+293 = 264*33841+29 = 132*67682+29 = 66*135364+29 = 33*270728+29

2114*4226+289 = 1057*8452+289

2116*4222+301 = 1058*8444+301 = 529*16888+301

2118*4218+329 = 1059*8436+329

2120*4214+373 = 1060*8428+373 = 530*16856+373 = 265*33713+108

2122*4210+433 = 1061*8420+433

2124*4206+509 = 1062*8412+509 = 531*16824+509

2126*4202+601 = 1063*8404+601

2128*4198+709 = 1064*8396+709 = 532*16793+177 = 266*33586+177 = 133*67173+44

2130*4194+833 = 1065*8388+833

2132*4190+973 = 1066*8380+973 = 533*16761+440

2134*4186+1129 = 1067*8373+62

2136*4182+1301 = 1068*8365+233 = 534*16730+233 = 267*33460+233

2138*4178+1489 = 1069*8357+420

2140*4174+1693 = 1070*8349+623 = 535*16699+88

2142*4170+1913 = 1071*8341+842

2144*4167+5 = 1072*8334+5 = 536*16668+5

2146*4163+255 = 1073*8326+255

2148*4159+521 = 1074*8318+521 = 537*16636+521

2150*4155+803 = 1075*8310+803

2152*4151+1101 = 1076*8303+25 = 538*16606+25 = 269*33212+25

2154*4147+1415 = 1077*8294+338

2156*4143+1745 = 1078*8287+667 = 539*16575+128

2158*4139+2091 = 1079*8279+1012

2160*4136+293 = 1080*8272+293 = 540*16544+293 = 270*33089+23 = 135*66178+23

2162*4132+669 = 1081*8264+669

2164*4128+1061 = 1082*8256+1061 = 541*16513+520

2166*4124+1469 = 1083*8249+386

2168*4120+1893 = 1084*8241+809 = 542*16483+267 = 271*32966+267

2170*4117+163 = 1085*8234+163

2172*4113+617 = 1086*8226+617 = 543*16452+617

2174*4109+1087 = 1087*8219+0 

Dzielniki 8934053 = 1087*8219

Modyfikacje:

Najgorsze dla złożoności jest sprawdzanie przypadków dzielników mniejszych niż 2ⁱ. Dużo można zdziałać modyfikując ten fragment. Można to zrobić przez eliminację złożoności, kosztem utraty znajomości drugiego z dzielników, lub zmniejszając krotność przypadków do sprawdzenia.

1) mniejsza złożoność.

Nie musimy znać kolejnego dzielnika, gdyż znajdziemy go później dzieląc n przez znaleziony. Dlatego dla drugiego procesora nie przysyłamy trójki [a, b, c], ale parę [a, c]. Wewnątrz znika nam warunek mnożenia 2 przez b złożoności bitowej O(n) zastąpiony operacjami złożoności bitowej O( log n). Łącznie uzyskujemy złożoność bitową

O( n log n * (log n)³ ) = O( n (log n)⁴ ).

2) wielokrotność reszty daleko od dzielnika.

Aby uzyskać przypadek, że dla danego a uzyskamy c-a=0, wartości c zbyt małe lub zbyt duże będą skutkowały niepotrzebnym sprawdzaniem. Aby go chociaż częściowo eliminować, do procesora drugiego dołączamy informację o k, czyli czwórkę [a, b, c, k],

Kiedy mamy połowić podstawę systemu a=2^kd przy odpowiednio małej wartości c, sprawdzamy, czy a > 2^kc. Jeśli tak, pętlę można opuścić, c nie będzie na tyle duże, aby wskazać dzielnik.

Podobna sytuacja jest dla c=a-1. Tu w każdej iteracji c będzie malało, lecz za mało, by wskazać dzielnik.

3) szybkie wskazanie małych dzielników.

Podejście 2) potrafi zmniejszyć liczność rozpatrywanych przypadków. Lecz małe dzielniki można wykryć jeszcze szybciej, obliczając nwd(a,c) spośród wartości [a, b, c] przesyłanych do procesora drugiego. Algorytm ten mocno pogarsza złożoność, i działa dobrze zwłaszcza przy istnieniu małych dzielników. Czasem policzenie nwd(b,c) wskaże dzielnik jeszcze szybciej, lecz z uwagi na duże wartości nie jest zalecany.

4) wiele procesorów

Dotychczas rozważane było użycie dwu procesorów. Nic nie stoi na przeszkodzie, by przedział [2ⁱ, 2ⁱ⁺¹) podzielić na kilka przedziałów oraz rachować w każdym z nich niezależnie. Wystarczy odpowiednio wcześnie podczas przygotowania wyznaczyć podstawy dla granic przedziałów.

Algorytm faktoryzacji, jak on działa?

Kiedy nie wiadomo co robić, zaczynam szukać innych form, to było granie w gry na sucho, w edytorze tekstu. Gdyż myślałem, że z teorii liczb nic już nie wycisnę. Myliłem się.

Powstał nowy algorytm. Jest bardzo dobry na maszyny wieloprocesorowe, ma stosunkowo mało przypadków do sprawdzenia, i do końca nie wiem jak działa.

Ustalmy 3 procesory. Pierwszy delikatnie przekształca liczbę faktoryzowaną n zapisaną w postaci trójki
n = a*b+c ,
przy czym a rośnie od parzystej liczby ograniczajacej pierwiastek z n z góry do podwojonej tej wartości po liczbach parzystych. Ma zatem 2^i przypadków, gdzie 2^(2i) < n < 2^(2i+1).
Przekształcenie polega na konwersji na system o 2 wyższy, czyli dany tożsamościami:
a*b+c = (a+2)*d+e, gdzie d in {b-2, b-1, b} w zależności od znaku (c-2b).
Dokładniej:
a*b+c = (a+2)*b + (c-2b), -1 < c-2b;
a*b+c = (a+2)*(b-1) + (c-2b+a+2), -a-3 < c-2b <0;
a*b+c = (a+2)*(b-2) + (c-2b+2(a+2)), c-2b < -a-2.
Np. jeśli n jest postaci n = 12441^2+987, to a=12442, b=12440, c=987-12441+12442 = 988; następna wartość jest równa a=12444, zaś b i c są znajdowane według wzoru (tu jest b=12438).
Dla n postaci n=12442^2+c przyjmujemy a=12444, pozostałe są odpowiednio dopasowane, by b było bliskie a, zaś c jest resztą.
Wartości (a,b,c) są przekazywane drugiemu procesorowi.

Drugi z procesorów przyjmujący (a,b,c), gdzie a jest parzyste, usiłuje dokonać konwersji, w której a jest połowione, b jest podwajane, zaś c w zależności od znaku (c<a) jest albo pozostawiane, albo zmniejszane o a, z równoczesnym zwiększeniem b o 1. Wzorami
(2a)*b+c = a*(2b)+c, c<a;
(2a)*b+c = a*(2b+1)+(c-a), c>a-1.
Pętlę kończymy, gdy wartość zmniejszana a jest nieparzysta. Mamy do czynienia co najwyżej z logarytmem binarnym z a przypadków (b jest nie większe niż podwojony sufit pierwiastka kwadratowego z n). Ciąg wartości c jest ciągiem nierosnącym.
Watości (a,c) oraz (b,c) są przesyłane do trzeciego z procesorów. Jedna z tych wartości jest stosunkowo mała.

Trzeci procesor ma najwięcej roboty. Sprawdza największy wspólny dzielnik przesłanych wartości, jeśli jest on większy od 1, znaleziona wartość jest zarazem dzielnikiem liczby n.
Bardzo podobnie wygląda cecha podzielności przez 9 lub 11, która stała się prekursorem algorytmu. Jeśli cyfry liczby dwucyfrowej zapisanej w systemie niedziesiątkowej są równe, to podstawa systemu zwiększona o 1 jest dzielnikiem liczby.

W praktyce nie zawsze musimy przesyłać dane do trzeciego procesora. Ale musimy je przesłać, gdy następuje zmiana wartości c, które nie jest zbyt małe. Przy c=1 możemy zignorować ten krok, przy c pierwszych przyspieszyć sprawdzanie testując podzielność przez c.

Przykładowy fragment rozkładu znajdujący dzielnik: 
n = a*b+c = (2990 * 2987 + 2923) = 8934053
Procesor pierwszy: przelicza od a=2990 do a=5980:
(2990 , 2987 , 2923 ) = (2992 , 2985 , 2933) = ...
= (4348 , 2054 , 3261) = ... 
Procesor drugi (4348, 2054, 3261):
(4348 , 2054 , 3261) = (2174 , 4109 , 1087) = (1087 , 8219 , 0)
Procesor trzeci:
nwd(4348, 3261) = 1;
nwd(2054, 3261) = 1;
nwd(2174, 1087) = 1087, znaleziony dzielnik 1087
Zaś 8934053 = 1087 * 8219.

Pętla działania procesora drugiego tworzy 'grzebień'. Tzn. dla połowy przypadków jest jedna iteracja, każda dłuższa pętla sąsiaduje z dwoma najkrótszymi. Krotność iteracji tworzy podciągi (1, 2, 1, 4, 1, 2, 1, 8, 1, 2, 1, 4, ... ).

Algorytm powstał, kiedy połączyłem cechy podzielności przez liczby sąsiadujące z podstawą systemu, konwersje o potęgi dwójki oraz konwersje o 2.

Kryterium pierwszości liczby od strony logarytmów dyskretnych

Pierwszość liczby jest funkcją, która zwraca wartość: 'argument jest / nie jest liczbą pierwszą.  
Wydaje się idiotyczne, by sprawdzać pierwszość liczby za pomocą, wydawałoby się znacznie gorszego przypadku logarytmów dyskretnych. Przecież do liczenia logarytmów dyskretnych wymagana jest znajomość faktoryzacji liczby, aby obliczenia były prostsze. Nawet przy metodzie liczenia logarytmów dyskretnych publikowanych na tym blogu.

Lecz konstrukcja tworzenia podgrafu nie jest związana ze znajomością liczb pierwszych, i właśnie ta konstrukcja potrafi wskazać dzielniki. Sprawdziłem to między innymi dla liczb będących iloczynami małych liczb pierwszych typu 11*7.

Post może być nieco zagmatwany. Jeszcze nie znam sposobu przedstawienia tematu w sposób 'dla ludzi', zwłaszcza fragmentu eliminowania liczb pierwszych nie będących dzielnikami.

Mamy liczbę n, której dzielniki chcemy znaleźć.
Przygotujmy dwie tablice pomocnicze, w jednej trzymamy liczby pierwsze i niektóre złożone {2, 3, 5, 7, 11, 13, 17, 19, ... }.
Do drugiej będziemy odkładać liczby pierwsze, które na pewno nie są dzielnikami n.
Dobierzmy ziarno a, będzie to pierwsza wartość pierwszej z tablic, np. a=2.
Tworzymy ciąg  skończony (m_1, ..., m_i)
m_k  = a^k (mod n), k>0
w którym m_i jest jedną z trzech wartości: 0, 1 lub a, zaś i jest pierwszą napotkaną pozycją z jedną z tych wartości. Z małego twierdzenia Fermata i<n.
Licząc wartości m_k, od razu możemy modyfikować tablice pomocnicze: m_k jest liczbą złożoną, wtedy powinniśmy je zfaktoryzować. Ewentualnie sprawdzić podzielność przez a oraz liczby pierwsze występujące w drugiej z tablic. Jeśli wśród dzielników m_k nie ma a, przenosimy je do drugiej tablicy pomocniczej. Z pierwszej tablicy kasujemy a i wszystkie jego wielokrotności.
Dzieląc m_k przez liczby pierwsze drugiej tablicy możemy natrafić na kolejną liczbę pierwszą, ta na pewno nie jest dzielnikiem. Robimy z nią tak samo jak z a. Przenosimy ją do drugiej tablicy usuwając jej wielokrotnosci z pierwszej. Czasem się zdarza, że m_k jest złożona i nie widzimy jej dzielnika. Wtedy czeka na inne m_j takie, że nwd(m_k, m_j) jest pierwsza i robimy znów to samo co przy a.
Jest to najgorszy fragment, który należy dopracować. Przy czym przeliczenia pomocnicze wystarczy robić tylko przy spełnionym warunku m_{k-1}*a>m_k.

Odczyt wyniku, gdy już wyznaczymy długość ciągu i. Mamy przypadki:
m_i = 0 oznacza, że a było dzielnikiem, zaś n jest potęgą a (a jest pierwsza);
m_i = a oraz dla każdego k=1..i a| m_k. Wtedy a jest dzielnikiem n
wreszcie m_i = 1 oraz i= n-1 lub i = (n-1)/2 oznacza, że n jest liczbą pierwszą.
W pozostałych przypadkach reszta z dzielenia n przez i jest związana z dzielnikiem. Czasem jest to dzielnik, jak w przypadku 91, gdyż długość ciągu dla a=2 jest i=12 oraz 91 = 7*12+7. W tym przypadku mamy nawet gotowy rozkład.
Dla n=81 a=2 uzyskujemy ciąg zawierajacy 54 elementy, do drugiej tablicy zostały przeniesione liczby pierwsze oprócz 3, zatem 3 jest dzielnikiem, zaś 81-54 = 27.
Dla n=77 a=2 ciąg ma 31 elementów, w którym wśród m_k pojawiają się 2, 25, 9, 29 i inne liczby pierwsze inne niż 7 czy 11. Mamy 77 = 2*31+15, oraz suma długości pozostałych ciągów generowanych przez a=7 oraz a=11 sumuje się do 15. Przykładowo ciąg wyznaczony przez a=11 ma postać [11, 44, 22, 11].

Jeśli mamy dzielnik, to długość ciągu go wyznaczająca bywa drugim z dzielników. Tu dla n=77 długość ciągu wyznaczonego przez 7 jest 11.

Policzmy konkretny przykład n=91,
a=2
k  m_k   dzielniki do odstawienia
 1    2   
 2    4
 3    8
 4  16
 5  32
 6  64
 7  37   2, 37
 8  74
 9  57   złożona 57=3*19, nie widać dzielników, mogą się jeszcze pojawić
10 23  23
11 46 
12   1  wyznaczone i=k = 12
szacujemy n/i otrzymując 91 = 7*12 + 7, zauważamy rozkład 91 = 7*13.
Bez zauważenia, stosujemy kolejne ziarno a=3, domyślamy się, że dzielnik jest nie większy od uzyskanej reszty 7. 

Faktoryzacja przez dzielenie przez sume cyfr, implementacja

Faktoryzację polegającą na sprawdzaniu reszt z dzielenia przez kolejne liczby można zaimplementować mniej standardowo, sprawdzanych przypadków jest więcej, lecz złożoność arytmetyczna się polepsza.
Teoria jest taka, jak opisana w sąsiednim poście, kiedy szukałem ciągu konwersji wskazującego dzielnik.
Liczbę n przedstawiam w systemie o podstawie 2^i, oraz mając te m cyfr sprawdzam cechami podzielności przez liczby 2^i-k, k nieparzyste malejące od 2^i-1 do 1. Jeśli uzyskana suma cyfr dzieli / jest wielokrotnością podstawy, mamy dzielnik.

Algorytm wykorzystuje konwersję podwajającą podstawę. Oto jej implementacja.
Liczba dana jest tablicą unsigned int a[8], w której cyfra najmniej znacząca jest ostatnim elementem tablicy. Podstawą systemu jest unsigned int p.
[\code=cpp]
void divp( unsigned int * a, unsigned int p ) {
  unsigned int k=0x80; // 2^i
  unsigned int i=0; // iterator tablicy cyfr a
  unsigned int b=0; // zmienna pomocnicza
  while(i++<length) { // tu length = 8, tablica 8
    b = b*p + a[i];
    a[i] = b/k; // cyfra to czesc calkowita
    b %= k; // reszta przenoszona dalej
    k /= 2; // nastepny bit
  }
}
[\code]

Sam rozkład zaczyna się przedstawienia liczby w systemie szesnastkowej. Następuje sprawdzenie cech podzielności przez 2, 3, 5 (suma cyfr jest podzielna przez te wartości). Dalej podział na przedziały binarne, zaś w każdym z nich dla kolejnych wartości nieparzystych wykorzystywana cecha podzielności przez p.
[\code=cpp]
int rozklad() {
  unsigned int a[] = ; // inicjowanie
  unsigned int s = ; // suma dla inicjacji
  if(s) ; // sprawdzenie malych dzielnikow
  unsigned int p=9; // zaczynamy od systemu siodemkowego 16-9
  unsigned int k=7; // 16-9 = 7, do cechy podzielnosci
  while(1) {
    unsigned int i= a.length()-1; // na cyfrę najmniej znaczącą
    s = a[i]; // sumowanie, inicjacja
    unsigned int t=1; // zmienna wspomagająca znajdowanie wspolczynnikow
    while(--i) {
      t = (k*t)%p; // dla cyfry lub aktualnej lub wspolczynnika przy i-1
      s = s + (a[i]*t)%p; // aktualna cyfra
    }
    if( p==s ) return p; // dzielnikiem jest p
    if( 2>k ) {
      k=p; // przeliczenie postaci liczby a[]
      divp(a, k+1); // zmienia dlugosc a
    } else k -= 2; // do następnego dzielenia
    p += 2; // następna podstawa
    if() ;// wyjscie w przypadku liczby pierwszej
  }
}
[\code]

Złożoność modułu divp jest liniowa względem długości cyfr liczby n, czyli logarytmiczna. Pętla faktoryzacji rozdziela przypadki na przedziały [2^i, 2^(i+1)), których jest logarytm z n.
W każdym takim przedziale następuje sumowanie cyfr liczby modulo p, znowu mające krotność logarytmu z n. Pętla ma dokładnie 2^(i-1) przejść, zakończonych uruchomieniem funkcji divp() (liniowe, gdyż 2^(2i)<n<2^(2i+2)).
Najwięcej kosztują przypadki dla dużych k (ale wtedy a jest stosunkowo 'krótka'). Ogólnie złożoność pesymistyczna jest O(n ln n).
Z kolei iloczyny cząstkowe wykonywane modulo nie mają wartości większych niż potęga 2 przekraczająca n.
Jeden z czynników występujących iloczynów ma wykres piłokształtny (liniowo malejący o 2, gwałtowny wzrost do 2^i po osiągnięciu takiej podstawy).

Heurystyka pierwszości zamiast rozkładu

Testowanie rozkładu liczby na czynniki ma bardzo dużo przypadków do sprawdzenia. I cóż z tego, że procesory testują miliony przypadków na sekundę. Obliczenia trwają godziny, dni, miesiące.
Zatem należy poszukać takich zależności, które zmniejszą tę krotność.

Jednym z pomysłów było zastosowanie cech podzielności. Znamy podzielność przez 9 (suma cyfr cyfry dziesiętnej dzieli się przez 9), przez 11 (naprzemienna suma cyfr liczby dziesiętnej dzieli się przez 11).
Dla systemów o podstawie p parzystej, ich analogi brzmią następująco:
liczba n przedstawiona w systemie p dzieli się przez p-1, gdy suma jej cyfr dzieli się przez p-1;
liczba n przedstawiona w systemie p dzieli się przez p+1, gdy naprzemienna suma jej cyfr dzieli się przez p+1;

Przykłady liczbowe.
7 ' 595 ' 485 _ p=1088
w systemie o podstawie p ma sumę 'cyfr' 1087, naprzemienną sumę -103. Zatem istnieje dzielnik i jest równy 1088-1 = 1087.
1 ' 3904 ' 3903 _p=10280
ma sumę 'cyfr' 7808, naprzemienną sumę 0, zatem dzielnikiem jest 10281.

Aby znajdować dzielniki, wystarczy rozpatrywać podstawy parzyste(!) 4*m, co już zmniejszy krotność przypadków o połowę.

Uda się wyciągnąć coś jeszcze. Przechodząc po małych podstawach: 4, 8, 16, 20, 24, 30, 36 oznaczmy przez a sumę, zaś przez b sumę naprzemienną. Dla tak małych podstaw dla liczb złożonych okazuje się, że największy wspólny dzielnik d = nwd(p,a) jest równy mniejszej z p, a.
Dla sumy naprzemiennej wystarczy, że ta suma będzie 0, 1 lub -1, aby istniały dzielniki. Dla testowanej liczby pierwszej sumy od razu miały większe wartości.
Iloraz n/p dla takich p jest bardzo bliski liczbom całkowitym. Nie zawsze jest to dzielnik, nawet przy zerowaniu się sumy naprzemiennej.
Kolejna cecha namierzajaca to zmiana znaku sumy naprzemiennej. Towarzyszy rozmienieniu wartości na drobne, co występuje także przy dzielniku.
Przy czym nie możemy wchodzić w minima globalne, lecz lokalne sum naprzemiennych, które dla kolejnych podstaw mogą mieć wartości (4, -49, 3 -> dzielnik, -18, 1, 25, 0)

Czy można zatem zaprząc sumy do poszukiwania dzielnika?
Udaje się dla początkowych wartości namierzyć: w tym kierunku jest dzielnik. Lecz przybliżając się do niego za pomocą konwersji:
a[i]*(2^ip^i) + ... + a[1]*(2p) + a[0] = 2^ia[i]*p^i + ... + 2a[1]*p+a[0]
dla coraz większych p (za każdym razem p się podwaja) mamy coraz  większe trudności z utrzymaniem dzielnika 'na celowniku'.
Pojawiające się zerowania lub cechy wspomagające namierzanie dzielnika nie trzymają się jednej wartości, lecz przemieszczają się.

W jednym tescie został namierzony dzielnik dla p=136. Przechodząc do p=272 dzielnik 'zniknął'. Został znowu namierzony dla p=256, przesunięcie 20. Zaś praktycznie dzielnik w tym przypadku był równy 136*8.

Zatem heurystyka niskim kosztem wskazuje: mamy dzielnik. Powinien być bliski takiej to a takiej wartości pomnożonej przez potęgę 2.
Lecz ze znalezieniem należy przeszukiwać brutalnie - chociaż można po wartościach podzielnych przez 4.

Metoda jest hipotezą.

Iteracyjne przechodzenie między systemami

Mając faktoryzować liczbę postaci
n = (a*p+b)*p+c,
zastanawiałem się, czy można tak przekształcić tę liczbę, aby szybciej znaleźć jej ewentualne dzielniki. Bez użycia brute force dla poszczególnych liczb pierwszych lub liczb nieparzystych (jak nie znam wszystkich liczb pierwszych).

Pierwsze podejrzenie padło na postać:
n = e*p + (d*p+f)^2 = d*p^2 + (e+2*d*f)*p + f^2
Wzór sugeruje, że jeśli cyfra najmniej znacząca jest kwadratem, to wyłaczając odpowiedni kwadrat z n znajdę mniejszą wartość związaną z dzielnikiem. Pozostaje iloczyn 'e*p', który jednak podczas konwersji na odpowiednią podstawę modyfikuje wyraz wolny.

Kolejne pytanie, czy mogę tak przekształcać liczbę, aby mieć kontrolę nad zmianami wyrazu wolnego?
Wtedy odpowiedni ciąg konwersji zacznie przybliżać wyraz wolny do zera, dla którego istnieje dzielnik.

Podczas faktoryzacji przez zmianę systemów szybko uzyskuję liczby trzycyfrowe, zatem co robią konwersje z wyrazem wolnym liczb trzycyfrowych?
Okazało się, że przy konwersji z systemu o podstawie p na system o podstawie p+k na współczynniki liczb trzycyfrowych znalazłem wzory iteracyjne.
Załóżmy, że zwiększamy podstawę (k dodatnie):
Współczynnik a jest stały.
Współczynnik b zmienia się liniowo od a albo k: b = b-2*a*k. 
Współczynnik c zmienia się kwadratowo od k: c = a*k^2 - b*k + c.
Po wstawieniu nowych wartości
[a,b,c] = [a, b-2ak, akk-bk+c]
uzyskiwałem liczby, które nie były 'cyframi' przedziału [0,p), ale nie zawsze były miejsze, czasem także większe od p. Po 'naprawie' cyfr uzyskiwałem prawidłowe postaci cyfr w kolejnych systemach. 

Ale dobrze widać było tylko sąsiednie podstawy. Kiedy zwiększałem k, nawet tylko do 2, przy większych a bardzo często naprawa modyfikowała b, co powodowało szybko narastający błąd szacowań. Zachowania się wyrazu wolnego nie byłem w stanie przewidzieć.
Zatem ten sposób nie jest najlepszy. Może służyć jako zastąpienie wersji rekursywnej programu wersją iteracyjną (k=2), gdyż wtedy następuje do dwu pożyczek/przeniesień. Z kolei nie wiemy bez sprawdzenia, które z nich będzie.

Ślady faktoryzacji

Przez ostatnie tygodnie przygotowywałem ślad przebiegu dwu algorytmów faktoryzacji w języku angielskim. Ponieważ nie widzę, aby ten blog przyjmował załączniki tekstowe, wyciąłem nieco tekstu jako rysunki.

Faktoryzowana jest liczba

169,747,007 = 19*1087*8219 .

Narzędzia konwersji:

Użyte algorytmy:

Pierwszy z przykładowych algorytmów znajduje dzielnik następująco:

Po przekroczeniu pierwiastka sześciennego, coraz więcej przekształceń upraszcza się do postaci jak podana: 

W drugim z algorytmów startujemy od postaci kwadratu z liczby. Konwersje między systemami są tak proste, że zapisuję kolejne postaci jedną pod drugiej, prawdzając od czasu do czasu małe dzielniki. Początek algorytmu wygląda następujaco:

Po znalezieniu dzielnika19 następuje porządne zamieszanie, gdyż współczynniki a, b, c, p, s i t muszą być przeliczone, lecz nie zmienia to wartości q.
Zaś zakończenie algorytmu wygląda następująco:

W ostatnim kroku s=219, gdyż przy zwiększaniu wartości s reszta z dzielenia 1099 przez 5 jest coraz mniejsza. Miły szczegół upraszczający rachunki.

Liczba od Sierpińskiego

Sierpiński w swojej monografii "Co wiemy, a czego nie wiemy o liczbach pierwszych" wspomniał kilkakrotnie o wartości 2^{101}. Jest to liczba 101-bitowa, która oparła się próbom faktoryzacji przedwojennych matematyków. Dowiedziano się tylko, że jest złożona.
Heksadecymalnie ma ona wartość 0x1F,FFFF,FFFF,FFFF,FFFF,FFFF,FFFF.
Rozłożona została przez Cuningham Project, w którym podano rozkład, dwa dzielniki, oba mają kilkanaście cyfr, zaś sama liczba ma ich 31. W dodatku nie są bardzo bliskie siebie. Iloraz dzielników ma około pięciu cyfr dziesiętnych.

Zastanowiło mnie, czy moje algorytmy poradzą sobie z nią i przygotowałem liczbę do faktoryzacji pod najszybszy z nich. Zajęło mi to kilka dni.
Mam teraz do czynienia z liczbami mającymi piętnaście cyfr dziesiętnych (51 bitów). Lecz zastanawiam się ze startem.
Mianowicie czasem mylę się przy dodawaniu, a przy takich wartościach nie jestem w stanie przetestować poprawności wyniku.
Dalej, mam do sprawdzenia gigantyczną krotność postaci, są to setki biliardów, czyli 7e14.
Wyliczyłem, że w niektórych iteracjach mogę sobie darować nawet miliony przypadków, przeskakując podstawy, przy których najmniej znacząca cyfra tworzy ciąg, który mogę jawnie podać za pomocą równania kwadratowego. Lecz dla tego równania kwadratowego ze względu na wielkość współczynników nie mam gwarancji poprawnego rozwiązania.

Potrzebne jest wsparcie pewnego typu interaktywnego programu komputerowego, którego brak nakłonił mnie do szukania algorytmów publikowanych w tym blogu. Programu, który nie zawsze liczy, ale prezentuje wartości pośrednie, oczywiście z dopasowaną pode mnie grafiką.


Wracając jeszcze do obliczeń, do przygotowanie zastosowałem algorytm publikowany tutaj 5.11.2012 "Pierwiastek kwadratowy przy systemach niedziesiątkowych".
http://matformac.blogspot.com/search/label/pierwiastek%20kwadratowy

Kiedy miałem blisko połowy cyfr znaczących pierwiastka, każdy kolejny iloraz był coraz bliższy podstawie potęgi 2, którą w danym kroku używałem. Po przekroczeniu połowy cyfr znaczących, każdy kolejny iloraz był dokładnie potęgą 2. Zatem połowę dzieleń można było zastąpić mnożeniem przez odpowiednią potęgę 2. 
W dodatku, rachując w systemie szesnastkowym lub binarnym, podczas obliczeń część bitów z początku lub końca liczb nie była ruszana. Bity te były po prostu kopiowane. Pozwoliło to utworzyć 'ramkę' długości do 52 bitów, wewnątrz której były przeprowadzane obliczenia. Z jednym wyjątkiem - przeniesienie, gdy pierwsza cyfra ramki była maksymalną wartością.
Ostatecznie operacje mnożenia i dzielenia sprowadzały się do działań z argumentami będącymi liczbami nie większymi niż pierwiastek czwartego stopnia liczby rozkładanej lub potęgach 2.

Dzielenia podczas faktoryzacji

Testowałem kilka rodzajów dzieleń podczas faktoryzacji liczby dziewięciocyfrowej. Nie chodziło tu o sam rozkład, bo ten już znałem, lecz o zachowywanie się dzieleń.

Stosowałem algorytm, w którym liczba n jest traktowana jako iloczyn dwu liczb
n = a*b + c
gdzie a przyjmuje kolejne wartości naturalne, b jest największe możliwe oraz c<b jest resztą.

W algorytmie dzielimy (b+a-c):(a+1) = d:(a+1) = q z resztą r, oraz podstawiamy:
b = b-q;
r = a-r;
oraz zwiększamy a = a+1;

Sprawdzałem dzielenia w których z dzielnej d tworzyłem paczki długości krotności cyfr a+1, oraz dodawałem do siebie ich kombinację liniową. Np.
784532:87  miało paczki długości 2 oraz sumy
78
78*(100-87)+45 = 78*13+45 = 1059
1059*(100-87)+32 = 1059*13+32 = 13799
Rekursją do ostatniej liczby wyznaczana jest reszta z dzielenia 13799:87, pozostałe wartosci są dodawane.
q = 78*100+1059 + (13799/87)
Dla dzielenia przez liczby ponad paczkę stosowana modyfikacja:
784532:14 ma paczki długości 2 powstające następująco:
78
-78*(14-10)+45 = -267
--267*(14-10)+32 = 1100
i dalej tak samo. Metoda jest opłacalna do wartości 1,2 wielkości paczki.

Drugim dzieleniem jest bezpośrednie przechodzenie metodą chłopów rosyjskich
(a*b+c) = (a+1)*e+f
mnożąc a przez 2, dzieląc b przez 2, zas gdy b jest nieparzyste zwiększenie c = c+a, oraz przenoszenie nadmiarów do e. Operacji jest tyle, ile wynosi logarytm binarny z b dla każdego dzielenia.

Z kroku na krok zauważyłem, że pierwsza metoda dosyć szybko się stabilizuje, zaś kolejne ilorazy cząstkowe układają się blisko ciągu arytmetycznego. Lecz to jest tylko wrażenie. Dla dalszych iteracji pojawia się szum, który modyfikuje ilorazy, zmniejszając powoli przyrost ciągu arytmetycznego. Niemniej była do dosyć dobra aproksymacja kolejnego ilorazu. Wartości na ogół były mniejsze, choć zdarzał się wzrost większy od szacowanego.

W drugiej metodzie pierwsze współczynniki reszt w kolejnych iteracjach zwiększały się o 1, dopóki potęga 2 była mniejsza niż a. Później następowało rozbieganie się wartości. Odejmowane reszty od c były potęgami 2 modulo a+1. Ta regularność zachowywała się, kiedy przeskakiwałem do kolejnej liczby pierwszej (a*b+c) = (a+m)*e+f, gdzie a oraz a+m były liczbami pierwszymi. Była jednak znacznie mniej widoczna.
Występował też szum, był on jednak związany z zapisem binarnym a aniżeli operacjami.
Ten sposób pozwala zmniejszyć krotność przypadków do krotności liczb pierwszych nie większych niż pierwiastek faktoryzowanej liczby, ze złożonością wewnętrzną logarytmiczną. Dla człowieka jednak mało strawny. Szacuję go na O( log^2 n ) przy znajomości tablicy liczb pierwszych.

Podsumowując, ten sposób faktoryzacji dosyć szybko, przy pomocy programowania dynamicznego eliminuje większość wymaganych dzieleń, zostawiając dzielenia przez bardzo małe liczby. Z kolei dla sprawdzania szacowań pojawia się mnożenie przez czynniki rzędu czwartego stopnia rozkładanej liczby.
 

Czy istnieje wzór na dzielniki?

Przyglądajac się przekształcanej w różnych systemach pozycyjnych liczbie, zacząłem się zastanawiać, czy istnieje gotowy wzór, aby, mając postać liczby w jakimś systemie, dowiedzieć się, o ile należy dokonać konwersji, aby uzyskać dzielnik.
W teorii sprawa wydaje się prosta. Dla liczb co najmniej dwucyfrowych (czyli jednej z najczęściej spotykanych postaci) szukany wzór można podać jawnie (a0, a1 to 'cyfry' liczby w systemie o podstawie p):

przekształcić tak, aby cyfra najmniej znacząca a0 stała się zerem,
a0 - x * a1 = 0  (modulo p+x )
lub
a0 + x * a1 = 0 ( modulo p-x )

Równanie rekurencyjne liniowe, które za pomocą równania diofantycznego np.
x * a1 - py - xy = a0
w wielu przypadkach uda się rozwiązać mimo nieliniowości. Rozwiązanie zawsze istnieje. Istnieje dokładnie jedno dla liczb pierwszych.

Praktyka jest gorsza. Otóż wzór działa w bardzo wąskim zakresie, dopóki nie przeniesiemy nadmiaru lub pożyczki. Kiedy nastąpi takie zdarzenie, pojawia się modyfikacja.
Dla innego przypadku, takim niefortunnym zdarzeniem psujacym wynik jest brak przeniesienia nadmiaru lub pożyczki.
Zatem można go stosować w bardzo wąskich zakresach, nieekonomiczne.

Śledziłem zmiany cyfry najbardziej znaczącej a0 przy podstawach p dla liczb trójcyfrowych, kiedy zmniejszałem p o 2 sprawdzając nieparzystych kandydatów na dzielniki.
Okazało się, że szukany współczynnik x zachowuje się jak ciągła liczba rzeczywista, wzrastając od jakiejś wartości 2k do wartości 2k+2. Prawie płynnie. Wartość k zależała od cyfr bardziej znaczących, w szczególności cyfry setek.
Wartość cyfry najmniej znaczącej wzrastała w każdym kroku o 2k, 2k+2, itd.
Różnice między kolejnymi cyframi najmniej znaczącymi liczby postaci początkowej 1'0'2 przy malejącej podstawie były 1, 3, 5 itd.
Dlatego stosunkowo łatwo było wyznaczyć podstawę, w której następowało przeniesienie lub pożyczka, po czym, szczególnie przy większych k, należało modyfikować wzór.

Ślad cyfry najbardziej znaczącej brany modulo p ze zmniejszającym się bardzo dużym p przypomina ślad punktu na kole, które porusza się coraz szybciej. Najpierw widać przyrosty odległości, następnie ślad się rozmywa, później pojawiają się wolniej obracające się 'szprychy', zmienia się ich liczność. Przy dużych szybkościach zmienia się kierunek 'widocznego' obrotu i zdaje nam się, że koło się kręci 'do tyłu'. Tak jest między innymi, gdy cyfra dziesiątek będzie się zbliżała do podstawy p.

Mnożenie przez zmianę systemów, ulepszenie

Ćwicząc kolejne działania arytmetyczne wróciłem do algorytmu mnożenia przez zmianę systemów pozycyjnych. Opracowany algorytm wygląda w pseudokodzie następująco:
Dane: liczby a>=b.
Liczba a konwertowana do systemu o podstawie b, tu ma wygląd a'.
Mnożenie a' przez b wyglądające jak 10, iloczyn to liczba a'0.
Konwersja liczby a'0 na system początkowy, uzyskując a*b.
Wynik: iloczyn a*b. 

Tutaj cyfry są liczbami należącymi do przedziału [0,p), gdzie p jest podstawą systemu liczbowego.
Konwersja z systemu o podstawie r na system o podstawie b ma przebieg następujący, przy oznaczeniach:
a jest tablicą cyfr  [an ... a1 a0]
r=10^i jest potęga 10 spełniającą warunek r<= b < 10*r:

1. c =  najbardziej znacząca cyfra liczby a w systemie o podstawie r, c = [an]
2. Dopóki nie pobrane wszystkie cyfry liczby a powtarzaj 3.-5.
3. przenieś kolejną cyfrę z a do c: c = c + [ak] = [a'n ... ak]
4. dla wszystkich liczb tablicy c wykonaj  [cj] = [cj] - (b-r)*[c{j+1}], j>k-1
5. wszystkie liczby c sprowadź do cyfr za pomocą pożyczek lub przeniesień w systemie o podstawie b
6. a' = c
Liczba a' jest postacią liczby a w systemie o podstawie b

Niezmiennikiem konwersji jest fakt, że ciąg [an ... ak] przy podstawie pierwotnej p ma taką samą wartość jak ciąg [a'n ... a'k] przy podstawie p+(b-r).

Niezmiennik ten powoduje, że wszystkie operacje pętli 2.-5. z wyjątkiem ostatniej są operacjami do siebie odwrotnymi, i można je pominąć. Zostaje tylko jednorazowo przeprowadzona operacja 4. 5. dla liczby a'0.

Zatem mnożenie a*b w zależności od argumentów sprowadza się do kilku mnożeń na liczbach wielkości r<=b < a, których liczność jest wyznaczona przez a/b.

Pomnóżmy przykładowo 5.027.608.376 * 1008.
Wtedy b=1008, r=1000, i=3. Ciąg a to ciąg trójek cyfr dziesiętnych
[5, 27, 608, 376].
Z konwersji niesparowany jest przyrost o 8 ciągu z dopisanym zerem [5, 27, 608, 376, 0].
Krok 4. to dodawanie pomnożonych przez 8 wartości poprzedzających. Czyli mamy 4 mnożenia przez 8 i kilka dodawań. Uzyskujemy ciąg
[5, 67, 824, 5240, 3008].
Krok 5. to usunięcie nadmiarów, by wynik był w systemie o podstawie 1000:
[5, 67, 829, 243, 8].
Wreszcie ostatnia czynność, odczytanie wyniku iloczynu  w systemie dziesiątkowym:
5.067.829.243.008 


Mając liczby zawierające 480 oraz 500 bitów, możemy wziąć r = 2^i oraz system binarny. Ten sposób zamienia podany iloczyn na iloczyn liczb zawierających zaledwie około 480 bitów jedna (podstawa) oraz 500-480=20 bitów druga (cyfra 'dziesiątek' przy podstawie liczby 480-bitowej). Oznacza to, że zamieniamy  mnożenie dwu gigantycznych liczb na mnożenie przez małą wartość. Gdyż (b-r) jest liczbą mniejszą o co najmniej rząd wielkości od r.