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.