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.

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.    

Projektowanie prototypu, wstęp

Aby sprawdzić siłę wymyślanych algorytmów, należy je przetestować praktycznie. 

Próbuję ożenić podejście funkcyjne z obiektowym. Muszą występować silne kompromisy.
Przykładowo obiekty odpowiadają relacjom bez argumentów, i są wskaźnikami do właściwej struktury. Jako takie mogą pojawiać się i znikać w dowolnym momencie.
Od strony podejścia funkcyjnego wprowadzane są liczne efekty uboczne oraz stan charakteryzowany licznością i typem argumentów. Eliminowane są też wycieki pamięci - właśnie użyty argument zostaje skasowany, co zmusza do specyficznej gospodarki pamięci 'z pamięcią'.

Kartoteka trzymająca dane została odpowiednio przygotowana, mam nadzieję. Wymaga jeszcze co najmniej dwu procedur: miotły do sprzątania pustych obszarów oraz podsumowania zawartości.
Po zbudowaniu prototypu kartoteki zacząłem testować moduł dodawania dwu liczb.

Moduł dodawania bezbłędnie przetrzymał dostawę 0, 1, 2, 3 argumentów  na wejściu przy wartościach sięgających 2^{65} z przekroczeniem zakresu slowa maszynowego (slowo maszynowe mam wielkości 2^{32}).
Nie tylko przetrzymał, ale i w odpowiednich przypadkach wyznaczył prawidłową wartość sumy.
Koszt: zostało użytych kilkanaście zmiennych, na jedną strukturę wymagane są co najmniej trzy plus kilka wspomagających.

Pierwiastek kwadratowy przy systemach niedziesiątkowych

W ostatnim poście pokazałem sposób faktoryzacji dużych liczb. Wejściową daną była poostać bliska pierwiastkowi kwadratowego z rozkładanej liczby n, czyli
    n = a b r _ p
gdzie a=1, b=0 lub b=1, r jest resztą z dzielenia n przez p, zaś p jest podstawą systemu. Parametry r,  p, b można liczyć klasycznie, ale istnieje też sposób mający złożoność arytmetyczną logarytmiczną.

Będziemy  wykorzystywać binarną postać liczby n.
Jako pierwszy krok potrzebujemy podział liczby n na dwie części, mniej więcej na połowy, które także oznaczam przez b, r. Dalsze modyfikacje sprowadzą te wartości do odpowiednich parametrów. Oto warunki podziału: podstawa jest potęgą 2 nieco większą niż połowa n:
  p = 2^k,
  n < p^2,
  b < p, r < p,
  b*p+r = n.
Warunki te zapewnia przyjęcie k jako połowy logarytmu binarnego z n,
  b = (n >> k) ;
  r = n % (2^k); 

Mając  już tę postać, posługujemy się konwersjami
    d e _ {p+f} = d {e+d*f} _ p
dla f będącymi potęgami 2. Oczywiście liczby e+d*f często są większe niż podstawa systemu p, zatem nadmiary są przenoszone do cyfry d.

Konwersje te stosujemy tylko wtedy, gdy
  b+2^{k+1} < p.
Chodzi o to, by przybliżać się do pierwiastka z nadmiarem - będziemy mieli wtedy ciągle do czynienia z liczbą dwucyfrową b r _ p. 
Algorytm to prosta pętla po zmniejszającym się k. Kiedy wspomniany warunek jest spełniony, stosowana jest konwersja. Krotność wykonania konwersji jest związana z krotnością zer w binarnym przedstawieniu pierwiastka. 

Zmniejszając iteracyjnie k = (k>>1), parametry b oraz p zbiegają do wartości pierwiastka kawadtowego, kiedy k=0, co najmniej jeden z nich osiąga wartość pierwiastka. Kiedy są różne, jest to mniejszy z nich.
Na zakończenie stosujemy konwersję na system z podstawą o 1 mniejszą, co doprowadza do potrzebnej trójcyfrowej postaci wejściowej.


Przykład liczbowy: szukanie pierwiastka 169 747 007 < 2^{28}.
Tutaj za k zamiast wykładnika 28/2 = 14 przyjęta jest cała potęga 2^{14} = 16384. W kolumnach podane są wartości b, r, p, k, spełnienie warunku b+2k < p oznaczane jest przez p-=k, co oznacza zastosowanie konwersji.

Przygotowanie: p = k = 16384, b = n*2^{-14} = 10360, r = n%k = 8767
b          r        p          k        konwersja o
10360  8767  16384  16384
10360  8767  16384  8192
10360  8767  16384  4096
10360  8767  16384  2048  p-=2048
11840  8767  14336  1024  p-=1024
12751  5695  13312  512
12751  5695  13312  256    p-=256
13001  5951  13056  128
13001  5951  13056  64
13001  5951  13056  32
13001  5951  13056  16      p-=16
13017  5327  13040  8        p-=8
13025  5207  13032  4
13025  5207  13032  2       p-=2
13027  5197  13030  1        p-=1
13028  5195  13029  0
Wiemy zatem, że pierwiastek całkowity to 13028, zakończenie to konwersja na system o tej podstawie (wystarczy wymienic b z p). Uzyskamy postać:
1  1  5195 _ {13028}
Wynika ona z przeniesienia nadmiaru: 13029 = 13028 +1 = p+1.

Faktoryzacja spradzająca duże i małe dzielniki

-->

Faktoryzacja liczb sprawdzająca prawie równocześnie największe i najmniejsze dzielniki jest możliwa do przedstawienia, chociaż algorytm jest bardziej skomplikowany. Przedstawię szkic za pomocą pseudokodu.

Ponieważ jednak podejście jest niestandardowe, najpierw funkcje biorące udział w algorytmach.

Liczba jest listą swoich cyfr, zaś cyfra to liczba z przedziału [0,p), gdzie p jest podstawą systemu liczbowego.

Last lista

zwraca ostatni element listy, w przypadku listy będącej liczbą, jej najmniej znaczącą cyfrę.

First lista

podobnie jak last, lecz zwraca pierwszy element listy, w przypadku liczb cyfrę najbardziej znaczącą.

Succ element

zwraca element poprzedzający dany iterator na liście, w braku takiego zwraca 0.

Move element lista

przenosi wskazany Element dołączając go do innej listy na koniec.

Rep [warunek] instrukcja

oznacza prostą pętlę powtarzającą się tak długo, dopóki warunek jest prawdziwy.

Fi instrukcje warunek

oznacza funkcję, która najpierw liczy instrukcje, a potem w zależności od warunku albo przywraca stare wartości, albo zostawia nowe, jest bardzo podobna do funkcji If, która robi to samo, ale w innej kolejności.

Konwersja [parametr]

tu oznacza konwersję zmniejszającą podstawę systemu o 2, chyba, że podany jest parametr. Wtedy zwiększa podstawę o ten parametr.

Corr [podstawa]

jest funkcją działającą na liczbach (listach) sprawdzającą, czy liczba jest poprawnie przedstawiona, w szczególności testuje warunki na cyfry, czy są liczbami całkowitymi w odpowiednim przedziale. Jeśli nie, nadmiary przenosi do cyfr bardziej znaczących, części ułamkowe jako liczby całkowite do cyfr mniej znaczących. Działa rekursywnie kończąc działanie, kiedy liczba jest we właściwej postaci.

Foreach lista : instrukcja

dla każdego elementu listy wykonuje podaną instrukcję.

Algorytm konwersji Konwersja [r=2]

Wejście: Liczba n pusta, Liczba m przy podstawie p+r, podstawa p

Rep [ m niepuste ] {

Move m.First n

Foreach n : element += r * Succ element

n.Corr [ p ]

}

Wyjście: Liczba n przy podstawie p, jako skutek uboczny może zmieniać p o r.

Algorytm faktoryzacji korzystający tylko z tej konwersji ma wadę: po znalezieniu dzielników nie wiemy, czy są one pierwsze czy złożone. Należy uruchomić algorytm ponownie dla każdego z dzielników, lub zapamiętać sprawdzonych kandydatów na dzielniki.

Jeśli dzielnik x dzieli podstawę systemu, aby dzielić liczbę, potrzebuje być także dzielnikiem wyrazu wolnego. W sprawdzaniu małych dzielników ograniczymy się tylko do iteracji, kiedy x jest dzielnikiem podstawy systemu.

Potrzebujemy trzech dodatkowych wartości: dzielnika minimalnego xmin, maksymalnego xmax oraz aktualnego x, inicjowanych x = xmin = 2, xmax=3.

Zaś dla większych liczb, często łatwiej i szybciej jest testować nieco większe wartości niż najbliższego kandydata na dzielnik. Mogą się one układać w pewien wzorek, w którym zaczynamy sprawdzać od większego kandydata, po czym w regularny sposób zmniejszamy do aktualnej najmniejszej podejrzanej wartości.

Przystąpmy do samego algorytmu faktoryzacji:

1. przygotowanie postaci liczby najmniejszej z trzycyfrowych n = [a, b, c], znalezienie podstawy p

2. Konwersja do najbliższej podstawy p nieparzystej i sprawdzenie dzielników xmin, x, xmax

3. dopasowanie licznika iteracji d

4. Rep [ istnieją kandydaci na dzielniki ] {

5. if 0= (d--) sprawdzenie i ponowne przeliczenie x, xmin, xmax, d

5a przy znalezionym dzielniku jego wyłączenie i restart algorytmu

6. Konwersja, modyfikacja podstawy p -=2

7. Corr [p]

8. if 0=n.Last czyli wyraz wolny równy 0, znaleziony dzielnik p; iteracyjny restart dla liczby z usuniętym Last, oraz z tą samą wartością xmin dla podstawy p

}

Ad 1. Chodzi o doprowadzenie do postaci b*p+c, w której b,c<p oraz p, b są największe z możliwych. W tym celu można przyjąć p jako sufit z pierwiastka kwadratowego liczby. Mając już tę postać stosujemy Konwercja [1], dzięki czemu uzyskujemy najmniejszą z możliwych liczbę trzycyfrową z a równym 1, b równym 0 albo 1 przy danej podstawie p;

Ad 2. W tym kroku sprawdzamy warunek p<xmin. Gdy 2=xmin, sprawdzamy podzielność przez 2 oraz stosujemy Konwersja [1], aby nie szukać niepotrzebnie dzielników po wartościach parzystych. Znalezione dzielniki usuwamy, i powtarzamy algorytm ze zmodyfikowaną liczbą. Gdy xmin jest większe, dobieramy x, xmax oraz d jak w 8;

Ad 3. Mała liczba naturalna d wskazująca, ile iteracji pętli będziemy czekać, aż x będzie dzielnikiem podstawy p. Jeśli nie jest wyznaczona w kroku 2, przyjmuje wartość 1 dla liczb postaci 3k+2, 2 dla liczb postaci 3k+4, wartość d=0 wykrywa podzielność przez 3;

Ad 4. Pętla 4-8, z której istnieją wyjścia na dwa sposoby, dla liczb pierwszych xmin>p, dla złożonych następuje ponowne wywołanie algorytmu dla poszczególnych dzielników z tymi samymi ustawieniami xmin;

Ad 5. W tym kroku d jako 0 wskazuje, że x dzieli p. Sprawdzamy, czy x dzieli też n.Last, czyli cyfrę najmniej znaczącą. Jeśli tak, mamy dzielnik x, zaś cała liczba ulega zniekształceniu, tak że należy zrestartować algorytm. Nie można podzielić podstawy oraz wyrazu wolnego przez x, gdyż wtedy mogą 'uciec' inne dzielniki. Jeżeli x>xmin, zmniejszamy x-=2 oraz liczymy d z pomocą y:

y = p modulo x; d = ( y%2 ? (y+x)/2 : y/2 )

Gdy x=xmin wstawiamy xmin = xmax oraz sprawdzamy wartości d liczone ww wzorem dla kolejnych elementów ciągu (p modulo (x+2*k) ). Interesuje nas podciąg malejący z wyjątkiem ostatniego elementu, który jest większy niż przedostatni. Taki podciąg cechuje się budową zbliżoną do (24, 22, 15, 8, 1, 17).

Przyjmujemy za xmax = (x+2*k) ostatniego elementu, oraz x = xmax-2.

Ad 8. Liczba jest iloczynem k*p, lecz nie wiemy, czy dzielniki są pierwsze czy złożone. Uruchamiamy algorytm dla tych wartości z aktualnym ustawieniem xmin;

Przykład liczbowy 169 747 007 = 19 * 1087 * 8219, szkic

Liczb niepodzielna przez 2. Do algorytmu wchodzi postać 

[a, b, c] = [1, 1, 5195] przy podstawie p=13028.

Następuje konwersja, by podstawa p była nieparzysta p=13027 = 3k+1, liczba [1, 3, 5197] , dwie iteracje doprowadzą do wartości p podzielnej przez 3. Obliczamy d dla p modulo 5 uzyskując 1, dla 7 jest większe. Mamy zatem d=1, xmin = 3, x = 5, xmax=7.

W pierwszej iteracji liczba przyjmuje postać 

[1, 7, 5207] przy p=13025. 

Podstawa p jest podzielna przez 5, lecz najmniej znacząca cyfra 5207 = 5*1041+2. Niezerowa reszta wskazuje, że 5 nie jest dzielnikiem. Zmniejszamy x do 3, d=1.

W drugiej iteracji liczba 

[1, 11, 5225] przy p=13023, 

p jest podzielne przez 3, sprawdzamy podzielność przez 3 cyfry najmniej znaczącej 5225 = 3*1741+2.

Teraz xmin=7, dane p jest podzielne przez 9, lecz x=9 nie jest dzielnikiem, jest 5 iteracji d=5 do sprawdzenia 7, reszta z dzielenia przez 11 jest większa, zatem xmax=11, x=7.

Po kolejnych kilku iteracjach okazuje się, że liczba postaci 

[1, 103, 7843] przy p=12977 

ma podstawę oraz najmniej znaczącą cyfrę podzielną przez x = xmin = 19.

Następuje restart algorytmu z xmin=19. Teraz należy powtórzyć krok 1. uzyskując postać liczby 169 747 007 : 19 = 8 934 053 jako 

[1, 3, 2923] przy p=2987. 

Teraz należy na nowo wyliczyć d=2, zaś do następnego dzielnika 21 jest 13 iteracji.

Kolejną postacią liczby jest [1, 7, 2933] przy p=2985.

I tak sprawdzamy kolejne duże dzielniki mniejsze od 2985, zarazem większe od 19.

Kiedy p opadnie do 1215, wiemy już z xmin=135, że dzielniki są w przedziale [xmax=157, 1215) lub sprawdzane właśnie 135, do sprawdzenia 173 dzielą nas dwie iteracje, zaś kolejne kandydatki: 171, 169, 167 i dalej dzieli taka sama krotność iteracji.

Osiągamy wreszcie postać [7, 580, 986] przy p=1089, dla której znowu następuje oczekiwanie na sprawdzenie dzielnika 217 (xmin=175) przy p=1085, lecz kolejna iteracja przekształca liczbę do postaci 

[7, 610, 0] przy p=1087. 

Mamy kolejne dzielniki: 1087 oraz 7*1087+610.

Obliczając jednak postać trzycyfrową tych wartości okazuje się, że są one mniejsze niż xmin=175, dzielniki są zatem liczbami pierwszymi.

Rozkład dokonany. 

W algorytmie tym rachujemy na coraz mniejszych wartościach, występujące dzielenie jest przez stosunkowo bardzo małe liczby. 

Logarytmy dyskretne - odczytywanie

W poprzednim poście z pierwszego października opisana została konstrukcja grafu, z której można odczytywać wartosci logarytmów dyskretnych. Czas na ćwiczenia praktyczne.

Porachujemy kilka logarytmów dyskretnych dla systemu Z_{17}. Został wybrany, gdyż dla tej liczby pierwszej pojawia się niejednoznaczność znajdowanych rozwiązań.
Dla przypomnienia podaję graf, którego lista jest wypisywana pionowo za numerem kolejnym wierzchołka. Wierzchołki listy są drzewami wysokości 1 o dwu liściach rozdzielonych '|'.
lista trywialna
0
lista długości 8:
1. 2 - 6 | 11
2. 4
3. 8 - 5 | 12
4. 16
5. 15 - 7 | 10
6. 13
7. 9 - 3 | 14
8. 1

Numer jest potrzebny, gdyż odczyt jest podobny jak w metodzie zliczania indeksów.
Rozwiązując równanie x^k = a (17) wyszukujemy pozycję a w grafie, tę pozycję listy dzielimy przez k. Jeśli nie jest to liczba całkowita dodajemy odpowiednią wielokrotność długości listy. Jeśli nie znajdujemy rozwiązania całkowitego, równanie nie ma rozwiązań. Wartość logarytmu odczytujemy z tego samego poziomu, co a.

Policzymy dla próby x^3 = 16 (17)
Szukamy 16 na listach. Znajdujemy na pozycji czwartej listy. Teraz znajdujemy położenie rozwiązania jako iloraz 4/3. Jest to liczba niewymierna, dodajemy długość listy (4+8)/3 = 12/3 = 4. Odczytujemy wartość: 16, oraz 16^3 = 16 (17).
Jest to proste, gdyż 16 = -1 (17) oraz (-1)^3 = (-1) = 16 (17).

Teraz coś trudniejszego, też z listy głównej:
x^5 = 9 (17)
Wartość 9 stoi na pozycji siódmej, zatem obliczamy pozycję wyniku (7+8)/5 = 3. Znajdujemy wartość 8, którą sprawdzamy: 8^5 = 9 (17). Zgadza się.

Kiedy wartość jest znaleziona w liściu, rozwiązanie też znajduje się na tym poziomie:
x^3 = 6 (17)
Wartość 6 jest na pozycji pierwszej listy, następuje tymczasowe 'podniesienie' do 2. Obliczamy pozycję wyniku (1+8)/3 = 3. Teraz 'schodzimy' poziom niżej znajdując dwie kandydatki: 5 oraz 12. Sprawdzamy je:
5^3 = 6 (17) oraz 12^3 = 11 (17).
Znaleźliśmy rozwiązanie 5, mając do sprawdzenia wybór jednej z dwu wartości.

Zwróćmy uwagę, że rachowane były wyższe potęgi niż druga. Dla tego sposobu wartość potęgi nie jest przeszkodą. W wikipedii większośc przykładów ograniczała się do potęgi drugiej, która przy tym podejściu jest często znajdowana 'z marszu'.


17 jest liczbą pierwszą. Sposób ten działa także dla liczb złożonych, lecz graf może być znacznie bardziej skomplikowany.
Przykładowo dla liczby 493 niektóre wartości mają do 8 pierwiastków, które trudno uporządkować liniowo. Dołączane do listy o długości 112 łańcuchy miały wysokość 1 lub 3. Dla 28 wartości przyłaczane łańcuchy długości 4 miały wspólny środkowy wierzchołek oraz maksimum lokowane na liście.
Wybieramy wtedy za jedną z list dowolny ze znalezionych łańcuchów długości 112, łańcuchy fragmentami są rozdzielne, lecz łączą się przy wartościach mających większą krotność pierwiastków. Przypomina mi to nieco budowę białek, kiedy poskręcane włókna wzajemnie się przeplatają. Pozostałe łańcuchy miały długości o 1 mniejsze niż dzielniki 493.
W przypadku doczepiania łańcuchów niektóre wartości musiały pojawiać się dwukrotnie, gdyż były sześcianami innych wartości, których łańcuchy nie obejmowały. Wartości te raz rozpoczynały łańcuch, za drugim razem sąsiadowały z korzeniem.