Sito liczb złożonych

Czasem pojawia się pomysł, jak zmniejszyć krotność przypadków faktoryzacji zachłannej. Pomysł ten często jest chybiony. Podobnie jak poniższy.

Nieparzyste dzielniki liczby N są jednej z postaci: 4p+1 oraz 4p+3.
Liczba N także, lecz w jej przypadku łatwo ją znajdziemy.
Wtedy możemy dopasować a, b w taki sposób, by N = a+4b, stąd b=(N-a)/4.
Wykorzystując cechy podzielności, N jest podzielna przez a wtedy i tylko wtedy, gdy N = a*(4b)+a w przypadku dzielnika postaci 4p+1 (b=p); lub N=a*(4b)-a w przypadku dzielnika 4p+3 (b=p+1).
Sugeruje to, że możemy badać pary (a,b), w których mamy nie połowę, ale czwartą część wartości mniejszych niż pierwiastek z N. Niestety, znajdujemy tylko dzielniki określonej postaci. Dzielnik występuje gdy a dzieli b.

Ostatecznie uzyskujemy następujący sposób rozkładu:
załóżmy, że  N = 4b+1
tworzymy ciągi a=1+4k, b=(N-1)/4-k dla k dodatnich. Inicjacja a=5, b=(N-5)/4.
Dla każdej wartości b przypisujemy zbiór jego dzielników {d}, początkowo pusty.
Przekształcamy: a=a+4, b=b-1. Korzystając ze zbioru dzielników {d} znajdujemy częściowy rozkład b na czynniki. W rozkładzie tym z b są wyłączone dzielniki nie większe niż a, mnożone przez wartość c. Wartość c jest 1, liczbą pierwszą lub iloczynem liczb pierwszych większych niż a. Gdy c =1, nie mamy dzielników, do zbiorów przypisanych wartościom b-d dodajemy dzielniki d.
Gdy c jest liczbą większą niż a, obliczamy resztę s z dzielenia c przez a. Pozwala nam ona znaleźć resztę r z dzielenia b przez a. Teraz do zbioru przypisanego wartości b-r dokładamy jego dzielnik a. Nie musimy się przejmować pominiętymi liczbami pierwszymi. Mogą się jeszcze pojawić, albo są w danej gałęzi nieistotne.
Kontynuujemy z kolejnymi k, uzyskując zbiór liczb, które na pewno nie są dzielnikami N.
Występujące dzielenia operują na mniejszych wartościach jak zwykłe dzielenie N przez a, co najmniej 4 razy w przypadku, gdy zarówno a, jak i b są pierwsze.
Kończymy, gdy wartość a przekroczy pierwiastek z N.

Niestety, nie wyłapane są wszystkie dzielniki, do reszty należy zmienić znak, oraz zmodyfikować a=3, teraz b zwiększamy zamiast zmniejszać (b+d), oraz w przypadku c większych od a do zbioru dzielników wartości (b-r+a) dołączamy a.

Dla N = 4b+3 wartości inicjujące a są odwrócone. Mamy 3 zamiast 5 i na odwrót. Chodzi o to, żeby po odjęciu a wartość była podzielna przez 4.

Skomplikowane, przykład może nieco rozjaśni.
Liczba 8934053
a=5; b=(8934053-5)/4 = 2233512. Rozkład częściowy b = 2^3*279189. Liczba c=279189, jej reszta z dzielenia przez 5 jest 4. Różna od 0, zatem 5 nie jest dzielnikiem. Liczbą pierwszą 3 się nie przejmujemy.  
Znajdujemy resztę z dzielenia b przez a: r = 4*8 = 2 (5).
Uzupełniamy tabele dzielników dla wartości: [b=2233512-2, {2}] = [2233510, {2}]. c>a, zatem [b=2233512-2, {a}] = [2233510, {2,5}]

W następnym kroku mamy a=5+4=9, b = c = 2233511, reszta r=s=8, czyli [b=2233511-8, {9}]  nie jest dzielnikiem N. Najgorszy przypadek z największymi wartościami.

Kolejny krok a=9+4=13, b=2233510 mamy zbiór dzielników {2,5}, czyli b = 2*5*223351. Nie mamy dzielników, ale korzystając z reszty r=6 modyfikujemy zbiory dla wartości [2233508, {2}], [2233505, {5}] i [2233504, {13}]

W tej części nie znajdziemy dzielników. Zajrzyjmy do drugiej:
a=3, b=2233514
Tu pojawia się czynnik 3.
Przeskoczmy nieco przypadków:
a=411, b=2233616 = 2^4*7^3*11*37 jest pełnym rozkładem, czyli c=1, wprowadzając wartości do zbiorów [2233618, {2}], [2233623, {7}], [2233627, {11}, [2233653, {37}] nie modyfikujemy zbioru dla wartości b=2233616-242+411 = 2233785 wartością {a=411}. Nie musimy mając już  wszystkie dzielniki.

Dzielnik pojawia się dla a=1087, b= (8934053+1087)/4 = 2233785 = 5*411*1087. jest on równy 1087.

Dodawanie z pomocą logiki

Swego czasu w czerwcu 2015 roku opublikowałem posta o tytule: "Operatory arytmetyczne wyrażone operatorami logicznymi".

Niedawno spróbowałem dodawać następującą metodą rekursywną: 

przechodząc po cyfrach argumentów, przy generowaniu kolejnej sumy biorę mniejszą z nich dla pierwszego argumentu a, zaś moduł różnicy cyfr jako argument drugiej b. Suma docelowa jest równa 2a+b. 

Np. 4283 + 1728 = 1223 * 2 + 3565
bo np. min(4,1)=1, 4-1=3 na pozycji tysięcy itd. Na pozycji jedności jest min(3,8)=3, 8-3=5. 

Sposób ten czasem się zapętla, lecz jest proste rozwiązanie. Przekształcam b=2*c+e, gdzie e jest zerem lub jedynką w zależności od parzystości b, oraz do dalszych przekształceń stosuję a+c, wartość zmniejszona co najmniej dwukrotnie.
4283 + 2728 = (1223*2 + 3564)+1 = (1223+1782)*2+1

Tym razem zapętlenie nie występuje. 

Metoda działa też w innych systemach liczbowych. Zaś w binarnym sprowadza się do opublikowanego w czerwcu algorytmu. 

Jest coś więcej. Algorytm ten wyznacza w jednej iteracji bit najmniej znaczący oraz z dużą dokładnością bit najbardziej znaczący. Chodzi o ten na pozycji, do której występuje przeniesienie przy zmianie długości liczby. Nazwę maską wartość, w niej tylko skrajne bity mają znaczenie.

Kiedy następuje podsumowanie - dodawanie można zastąpić OR-owaniem. 

Przykład: 15+9 = 1111b + 1001b

pierwsza iteracja: a=1001b, b=0110b, do dalszej części przechodzą 

(1001b+11b)*2 + 0

które przekształcamy do (001b+11b)*2 + 10000b (b jest parzyste, bit najmniej znaczący jest równy 0)

druga iteracja: a = 1, b=10b, 

po przekształceniu: (1+1), maska 0000

trzecia iteracja: a=1, b=0, tu wyznaczenie 2a = 10 jest tak duże, że zahacza o wyznaczony już bit maski w poprzedniej iteracji.

Podsumowanie, czyli orowanie wartości 10000b | 0000<<1 | 10<<2 = 11000b  = 24

jest to wartość 15+9 = 24. 

Inny przykład: 15+7, już szybciej. 1111b + 111b

a = 111b, b = 1000b, maska 00000b

suma 111b + 100b 

a = 100b, b = 11b, maska 1001b, bit najbardziej znaczący pobrany z 2a = 1000b

suma 0+1 = 1, kończymy wartością 01, co dajemy jako maskę.

Podsumowanie: 00000 | 10010 | 0100 = 10110b

Są 3 iteracje w porównaniu z przechodzeniem 4 razy po bitach. Wszystkie przechodzenia można traktować równocześnie jedną operacją logiczną. 

Inne podejście do programowania

W Świecie Nauki 3/2015 jest artykuł "Wystarczy dodać pamięć" o komputerach pamięciowych, w których procesorem obliczeniowym jest sama pamięć. Potrzebne jest nowe oprogramowanie sterujące.
Jest ono bliskie mojemu pomysłowi, dążącemu do obliczeń rekursywnych.

Paradygmat o zachowywaniu stałego programu jest bardzo silny. A widzę potrzebę jego zniesienia.

Program może być tablicą wskaźników na dane umieszczone w pamięci. Następują dwa przebiegi.
Pierwszy z nich porównuje dane ze sobą, sortuje i przygotowywuje trójki wskaźników do obliczeń. Wykorzystuje formuły programowania relacyjnego, takiego jak w Prologu.
Dopiero drugi przebieg przelicza dane wskazane w pierwszym przebiegu.

Funkcje dodawania i mnożenia są wielowartościowe. Potrzebny jest zatem znacznik zakończenia, tu średnik. Składniki sortowane są najpierw funkcjami, potem malejąco. Czynniki sortowane są najpierw funkcjami, potem malejąco.
Funkcje odejmowania i dzielenia są jednoargumentowe, tworząc formuły zapisane w ONP:
a - b = + - b a
a : b = * : b a
Formuł jest dużo, oto niektóre z nich:
+ a 0 = a
+ - a a = 0
+ * a b; * a c;; = * + b c; a;
Ostatnia jest formułą, gdyż kolejność nie ma znaczenia. Sortowanie i tak może ją zmienić.

Przy odejmowaniu, dzieleniu następuje próba wyłączenia wartości tworząc funkcję, tak
+ - 5 12;
jest przekształcane formułą w drugiej fazie na (12 = 5+7)
+ - 5 7 5;
kolejny krok w pierwszej fazie wyznacza do obliczeń + - 5 5, czyli w drugiej fazie uzyskujemy
+ 7 0;
pierwsza faza znów modyfikuje do wyniku
7

Cechy podzielności w systemach niedziesiątkowych

Znane są cechy podzielności przez 2, 3, 4, 5, 9, 10, 11 w systemie dziesiątkowym.

Okazuje się, że w systemie o podstawie n cechy podzielności są równie łatwe do wyprowadzenia.
W szczególności, traktując kolejne cyfry liczby w systemie o podstawie n jako współczynniki wielomianu, znajdujemy resztę z dzielenia przez n-x obliczając wartość tego wielomianu dla x modulo (n-x).
Wynika to z tych sposobów mnożenia i dzielenia, w których dowodach korzystamy z konwersji na inny system.   
Sposoby te są przenośne na wszystkie systemy, w których rozpisaniu kolejne cyfry są przy potęgach liczby naturalnej większej od 1.

Przykładowo mając liczbę w systemie o podstawie 20, np 13 7 2 8 10, daje ona w wyniku dzielenia przez 19 resztę 2, bo cecha podzielności to wartość wielomianu w 1:
13+7+2+8+10 = 22_{19} (dziesiętnie 40) = 2 (mod 19)

Faktoryzacja dwoma kandydatami równocześnie, szczegóły

W poprzednim poście wspomniałem o pomyśle równoczesnego badania dwu dzielników s*t, s<t, w którym s i t tworzą ciągi arytmetyczne o przyroście |2|. Sposób nie unika mnożeń lub dzieleń, lecz obecne są przez stosunkowo małe wartości. Dodatkowo wartości te łatwo jest szacować.

Wygodnie jest nie pamiętać liczb n, lecz ich trójki [a,p,b]: n = a*p+b.
Jest to na dwu poziomach, jeden przy podstawie p, obliczanej przez dodanie zmniejszającej się o 8 wartości w. Drugi poziom dotyczy iloczynu a*w, rozpisanego trójką a*w = [c,t,d].
Wartość b zapisujemy analogiczną trójką b = [e,t,f]. Przy przekształcaniu odejmujemy b-a*w = [e-c,t,f-d] starając się, by wartość f-d była dodatnia. W razie potrzeby dodajemy s*t zmniejszając a. Wartość e-c dla sprawdzanego zakresu jest zawsze dodatnia. Podobnie można zrobić z wartością s.
Pamiętamy stare wartości (to te z primami). 
Zaczynamy od s=3, t będącą nieparzystą wartością nieco mniejszą niż pierwiastek kwadratowy  z rozkładanej liczby n, przyjmujemy w = 2*t-10;
Przy sprawdzaniu kolejnych dwu dzielników postępujemy następująco:

w = w'-8;
p = p'+w = (s'+2)*(t'-2) = s*t;
a*w = [ct, t, dt] = [cs, s, ds]
b' - a*w + s*t*x = [e, t', f] - [ct, t, dt] + s*t*x = [e-ct+s*x, t, 2*e+f-dt]
Rozwiązujemy nierówność liniową e-ct+s*x>0 zmiennej x szukając najmniejszego możliwego x. Od pewnego momentu jest x=1 lub x=0.
Jeśli 2*e+f-dt= 0, mamy dzielnik t.
Podobnie można zrobić z kandydatem na dzielnik s, uwzględniając znaleziony x. Tym razem:
b'-a*w + s*t*x = [e, s', f] - [cs, s, ds] + s*t*x = [e-cs+t*x, s, f-2*e-ds]
Jednak ostatnia wartość f-2*e-ds bywa na tyle mała, że trzeba dodać kilka s, aby uczynić ją nieujemną. Krotność dodanych s zmniejsza (e-cs+t*x).
Po tych przekształceniach zmniejszamy a = a'-x

Przykład numeryczny.
Wartość w = 2164, p' = s'*t' = 955*2035, 'obliczamy' p = p'+w = (s'+2)*(t'-2):
n = [a, p, b] =
4 * 1 943 425 + [570, 2035, 403] = 4 * 1 943 425 + [1215, 955, 28]
Sprawdzamy kandydatów t=2033, s=957:
w = 2164-8 = 2156,
a*w = [4, 2033, 492] = [9, 957, 11]
nowe b = [ct, t, dt] = 570*2035+403 - 4*2033-482 + 2033*957x =
(570-4+957x)*2033 + (1140+403-492) 
w tym przypadku mamy x=0, bo 570-4>0
b = 566*2033+1051 = [566, 2033, 1051], 1051<>0, t=2033 nie jest dzielnikiem
dla s:
b = 1215*955+28 - 9*957-11 + 957*2033*0 = 1206*957-2413 = 1206*957+458 = [1206, 957, 458]
s też nie jest dzielnikiem
wartości a nie zmniejszamy, bo x=0. 

Następna para s*t = 959*2033 ma a*w = 4*2156 = [4, 2031, 468] = [8, 959, 920]
Przekształcamy dopóki s stanie się większe niż t albo trzeci element trójki będzie zerem.

Wydaje mi się, że obliczanie s w analogiczny sposób jak t jest bezcelowe, lepiej podzielić znaną wartość b = ct*t+dt przez s resztami modulo.
Jak już wspomniałem poprzednio, dla małych s 'a*w' są gigantyczne, ale dla nich x jest duże. Zatem w maleje, a maleje, wkrótce zaczynamy coraz rzadziej dodawać s*t, tylko przy zmianie a.

Równoczesne sprawdzanie dwu kandydatów na dzielniki

Dokończyłem standardowy przykład, w którym wyrażałem n w systemach o podstawach s*t, s<t oraz s i t są liczbami nieparzystymi sumującymi się do stałej wartości parzystej równej sufitowi z pierwiastka z n.

Wykorzystałem
  n = a*(s*t)+b,
gdzie s=3+2k, t = 2(1/2 int sqrt n)-1-2k (nieparzysta wartość naturalna ograniczajaca pierwiastek z n z dołu).

Liczba n jest równaniem kwadratowym zmiennej k, zatem druga różnica jest stałą równą 8. Ewentualny dzielnik liczny n dzieli także resztę b, przy czym należy sprawdzić obydwa przypadki dla s oraz dla t.

Startujemy z s=3, t odpowiednio dopasowane do n. Przekształcamy n zwiększając równocześnie s o 2 oraz zmniejszając t o 2. Oznacza to przechodzenie na system o podstawie (s+2)(t-2), który różni się od danego o stałą (drugą różnicę, czyli 8a). Wartość a dla małych s bardzo szybko maleje. Przy małych a możemy zapisywać współczynnik b w systemach o podstawach s oraz t. Te przekształcenia są liniowe.
Mając iloczyn st, iloczyn (s+2)(t-2) uzyskuje się przez dodanie 8a do zapamiętanej wartości, ta zaś jest dodawana do podstawy. Ze wzrostem s (suma s+t jest stała) wartość ta maleje.

Jest to nadal trial division. Lecz wydaje mi się, że pomijając skrajne przypadki małych s uzyskamy algorytm mający mniej operacji aniżeli algorytm konwersji na kolejne systemy liczenia. Zwłaszcza jeśli nie zmieni się współczynnik a, gdyż po zmianie s, t zmniejszamy współczynnik b o pewną stale malejącą o 8 wartość (modulo st).
Niestety, występuje dzielenie b/s oraz b/t. Brak części ułamkowej tych ilorazów wskazuje dzielnik. Nie można przewidzieć jego wystąpienia.

Algorytm kończy się, gdy s stanie się większe od t. Praktycznie wyglądało to tak, że zapowiadał się koniec dla współczynnika a=4, dla którego jest kilkadziesiąt procent przypadków. Jednak w przykładzie pojawiły się wartości iloczynów st, dla których a było równe 3.
Obliczenia pomijając początkowe przypadki są bardzo proste. Jeszcze je badam.

Inne wyszukiwanie siłowe zamiast 'trial division'

Zapragnałem już pisać w C++ prototyp podejścia wielowątkowego relacyjno-funkcyjnego pod DOSa, kiedy pomyślałem sobie obejrzeć liczbę w systemach będących kwadratami kolejnych liczb naturalnych. Spodziewałem się, że pojawią małe reszty. Myliłem się.

Pojawiły się za to inne proste zależności. Mianowicie wiedziałem od początku, że nie będę musiał szukać kolejnego kwadratu siłowo. Różnica między p^2 oraz (p+2)^2 wynosi 4(p+1). Kiedy p jest nieparzyste, różnica ta jest kolejną wielokrotnością 8. Wynika to z p=2k+1, stąd 4((2k+1)+1) = 8k+8. Wartość k ma znaczenie w przebiegu i jest wykorzystywana.

Dla liczby n, począwszy od wartości p spełniajacej warunek p^4>n, dostajemy liczby 'dwucyfrowe'  [a, b]_{p^2}:
n = a*p^2 + b  = a*(2k+1)^2 + b.
Zatem liczb mających więcej niż 3 cyfry jest zaledwie pierwiastek czwartego stopnia.

Zwiękzając p, wartość a z początku bardzo szybko maleje, lecz gdy 8ak<p, pojawia się prostsza zależność.
Normalnie obliczałem odległość do kolejnego pierwiastka z p+1 do p+3, oraz zwiększałem zapamiętaną podstawę systemu o 8. W ten sposób znałem nową podstawę systemu bez podnoszenia do kwadratu, lecz przez dodanie sumy narastajacej o stałą. Teraz wystarczyło przekonwertować liczbę do odpowiedniego systemu.
Kiedy 'cyfra dziesiątek' a była odpowiednio mała, w szczególności spełniała warunek 8ak<p, okazało się, że konwersja sprowadza się do odjęcia stałej wartości 8a od b. Zatem można ją było zastąpić zwykłym odejmowaniem oraz naprawą w postaci, gdy b było ujemne operacjami
b=b+p^2; a=a-1
oraz oczywiście p=p+2. 
Przy stałym a reszty b (wzrost p) tworzą ciąg malejący. Dla małych a wykres jest piłokształtny, a ciąg lokalnych maksimum b tworzy ciąg rosnący. Przejścia są proste, nie wymagają mnożenia ani dzielenia, jedyne co trzeba pamiętać od pewnego momentu i modyfikować to wartości
{a, b, k, p^2}.
Wartość a maleje jak 1/x, b ma charakter piłokształtny, k rośnie liniowo, p^2 rośnie kwadratowo, lecz wyznacza się je liniowo. 

Najgorsze operacje to mnożenie 8ak po zmianie a, gdy a jest duże, oraz odejmowanie tej wartości od b. Oczywiście wtedy b staje się ujemne, oraz trzeba je skorygować dodając odpowiednią krotność aktualnego p^2.

Kolejnym niemiłym elementem jest znajdowanie dzielnika d, że n = d*q. Otóż wtedy d dzieli także b. Zatem mamy dzielenie (często dużych wartości) b przez d=(2k+1). Czyli to samo co w sposobie 'trial division'.

Wartości b były duże, gdyż p^2 dążyło do n. Tak więc metoda nie jest lepsza niż trywialne dzielenie przez kolejne wartości. Może tylko nieco dłużej utrzymuje mniejsze liczby dla dzielenia.
Metodę można bardzo łatwo odwrócić zaczynając od pierwiastka p^2<n<(p+3)^2, p nieparzyste. Wtedy zmniejszamy p, obliczenia przez pewien czas są bardzo proste, dopóki nie zaczynamy dodawać (konwersja) wartości większych niż podstawa p^2. Wtedy lepiej zmienić sposób, lecz kandydatów na dzielniki jest już stosunkowo mało (szacuję na parędziesiąt procent).

Metoda zasugerowała jednak kolejny sposób, w którym można badać równocześnie dwóch kandydatów na dzielniki. Przypadek ten jest obecnie sprawdzany.

Szybkie przechodzenie z dużych liczb binarnych na dziesiątkowe (lub odwrotnie)

Ponieważ 125 dziesiętne jest 1111101 binarnie, można bardzo szybko przeskakiwać między tymi systemami dla większych wartości.

Zacznijmy od dziesiątkowego. Dzielimy liczby po trzy (rozdzielamy tysiące), następnie każdą z wydzielonych trójek mnożymy przez taką potęgę 8, jak daleko jest ona od trójki zawierającej cyfrę jedności. Reszty przenosimy. Uzyskamy w ten sposób liczbę w systemie o podstawie 125.
Zwiększamy podstawę systemu o 3 algorytmem, który już wielokrotnie opisywałem (iteracyjnie po coraz większej krotności 'cyfr' odejmujemy daną od potrojonej wcześniejszej i korygujemy przeniesieniami, by była nieujemna).
Uzyskamy wartość w systemie o podstawie 128, i teraz każdą wartość zapisujemy siedmioma pozycjami binarnymi. Może poza najbardziej znaczącą.

Przykład: 33 567
W systemie o podstawie 125 uzyskujemy 264 567, po przenosinach nadmiarów jest to liczba 'trójcyfrowa' 2 18 67.
Konwersja na system o podstawie 128
2 |18 67
2 12 |67
2 6 31
Zatem wartość binarna to 2=0000010, 6=0000110, 31=0011111.
Razem 10.0000110.0011111, w szesnastkowym 0x831F.

W drugą stronę, z binarnego na dziesiętny.  
Najpierw tworzymy siódemki cyfr, następnie każdą z nich zapisujemy dziesiętnie w systemie o podstawie 128. Konwertujemy o -3 na system o podstawie 125, dzielimy przez odpowiednie potęgi 8 przenosząc reszty oraz uzyskujemy liczbę dziesiętną zapisaną tysiącami (wartości mniejsze uzupełniane zerami z wyjątkiem najbardziej znaczącej).

Przykład
1000011110100001 = 10 , 0001111 , 0100001 = 2 15 33 przy podstawie 128
2 |15 33
2 21 |33
2 27 96
Dzielenia: 2/64=0 r 2; 2*125+27 = 277, 277/8 = 34 r 5; 5*125+96 = 721.
Zera dodane 034, lecz jest to pozycja najbardziej znacząca.
Jest to liczba dziesiętna 34 721.

Położenie dzielników statystycznie

Pozwoliłem sobie na przegląd reszt z dzielenia przez kolejne liczby naturalne.
Wybrałem liczbę, dzieliłem ją przez liczby poczynając od pierwiastka do coraz mniejszych wartości.
Skupiłem się na poszukiwaniach lokalnie najmniejszych wartości, tzn. takich, że dzieląc kolejno przez q+1, q, q-1, wartość reszty n%q jest mniejsza od sąsiednich.

Jeśli liczba n = (ap+b)p+c (p jest podstawą systemu między pierwiastkiem sześciennym i kwadratowym z n), to wyznaczone przez q lokalne minima układają się dosyć specyficznie.
Weźmy przedział dla stałego a. Wtedy dla wartości bliskich iloczynom a oraz kwadratowi pewnej liczby oraz (a+1) oraz kwadratowi innej liczby lokalne minima są bardzo rzadkie.
Równie rzadko, chociaż częściej występują one dla liczb n o wartościach bliskich t = (ap + p/2)p+c, czyli dokładnie w połowie przedziału, w którym iloraz n/p^2 ma stałą część całkowitą.
Dalej, kiedy znajdę resztę lokalnie najmniejszą, oraz nie jestem zbyt blisko krawędzi (lub środka) przedziału o stałym a, łatwo wyznaczę konwersję, którą należy zastosować, by dostać się do najbliższej lokalnie najmniejszej reszty. Trafiam dokładnie w ten dzielnik lub w sąsiednią wartość.
Dla wartości p większych niż t przekształcenie jest postaci b/q, oraz jeśli od q odejmę b/q znów mam wartość reszty lokalnie najmniejszą.
Dla wartości q mniejszych niż t przekształcenie jest postaci (q+c)/(q-b). Wtedy b jest większe niż połowa q, dlatego dzielę przez coraz miejszą wartość, zmniejszając podstawę q o tę wartość trafię blisko kolejnej lokalnie  minimalnej wartości.

Przekształcenia te nie działają przy krawędziach lub na środku. Jest to spowodowane zbyt małymi różnicami między kolejnymi resztami. Różnica między resztami dla p, p-1 oraz między resztami p-1, p-2 różni się o dokładnie 2a. Przy krawędziach (b bliskie 0 lub b bliskie p) jest to tak mała wartość, że można zastosować wzór na wyraz szeregu arytmetycznego i uzyskamy dobre przybliżenie. Po trafieniu w lokalnie najmniejsze minimum ilorazy b/q wolno maleją do 3 pozwalajac pominąć mnóstwo wartości p, dla których dzielniki są większe. W okolicy połowy przedziału, czyli blisko t, mamy skoki, jest coś w rodzaju bifurkacji, oraz widać układające się monotonicznie ciągi, ale co druga reszta. Możemy je traktować oddzielnie, np. licząc tylko dla wartości p nieparzystych. Przy przekraczaniu t następuje skok. Od tego momentu ciąg co drugiej reszty powoli zaczyna się zachowywać coraz mniej przewidywalnie, aż wreszcie lokalne minimum pojawia się co trzecia reszta. Od tego momentu można przyjąć drugie przekształcenie (q+c)/(q-b), które powoli rośnie wskazując kolejne lokalnie minimalne reszty, wreszcie dojdziemy blisko krawędzi, gdzie iloraz jest duży, ale i tak bywa za mały dla sprawdzenia monotoniczności reszt.

Dla większych a obszar w środku staje się coraz bardziej przewidywalny, tak że można cały czas korzystać ze podanych przekształceń uważając tylko, kiedy b urośnie powyżej połowy podstawy p.

W ten sposób dla liczby 2^(101) pomijam nawet miliony przypadków dzielenia przez wartości, dla których reszta z dzielenia będzie dodatnia.
Zaś dzielniki, zwłaszcza te najbardziej poszukiwane (RSA) układają się najbardziej gęsto w okolicach niezbyt bliskich t (b bliskie połowie systemu), ale zarazem bardzo daleko od granic (b małe albo niewiele mniejsze niż podstawa systemu), przy którach a jest równe 1, 2, 3, i tak do 9. Zwłaszcza na obszarze wartości nieco mniejszych niż t (b większe niż połowa podstawy systemu).

Najmniejsze wartości lokalnie najmniejszych reszt występują, gdy b jest małe. Reszty te są mniejsze nie tylko od p, ale i od b. Drugim takim obszarem są bardzo duże b, bliskie podstawy systemu. Wtedy lokalnie najmniejsze reszty też są mniejsze niż (p-b). Oraz sporadycznie gdzie indziej.

Dodatkowo zauważyłem, że w badanym przypadku, którym była liczba złożona, minima lokalne sąsiadujące z dzielnikiem (w którym reszta jest oczywiście 0) miały wartość większą niż połowa podstawy systemu p.

Można to wykorzystać. Wiadomo, gdzie szanse na znalezienie dzielnika są znikome, zwłaszcza po przebadaniu dzieleń przez kilka podejrzanych wartości.

Operatory arytmetyczne wyrażone operatorami logicznymi

Przeszukując źródła nie zauważyłem, by gdziekolwiek były podane operatory arytmetyczne za pomocą operatorów logicznych. Zauważyłem tylko, np.  w Kalejdoskopie Techniki 7/1978, że w układach cyfrowych odejmowanie, mnożenie i dzielenie są funkcjami dodawania. Zaś wykłady z elektroniki obejmowały konstrukcje sumatorów oraz komparatorów, cyfra po cyfrze.
O ile komparator musi działać liniowo, o tyle dodawanie można zapisać funkcjami logicznymi AND &, OR |, NOT ~. Wymagana jest tylko rekursja oraz operatory przesunięcia SAL (w lewo, czyli mnożenie przez 2) oraz SAR (w prawo, czyli dzielenie przez 2).
Prawdę mówiąc, nie potrafię zapamiętać, który z SAL, SAR to mnożenie, a który dzielenie. Przed każdym użyciem muszę to sprawdzać.

Algorytmy są rekursywne, w najgorszym przypadku przewijają cały obrabiany bajt. Przy przesuwaniu wstawiane są zera.
I tak inkrementacja INC A wyraża się za pomocą ( ? : ) z C następująco:
INC A = ( A&1 ? SAL INC SAR A : A|1 );

Funkcja dekrementacji DEC A stosując maskę -2 = 0xF...FE
DEC A = ( A&1 ? A&(-2) : 1|(SAL DEC SAR A) );

Funkcja porównania = zwraca A (równe) albo 0 (różne)
(A=B) = ( (A|B)&~(A&B) ? 0 : A ) = ( A XOR B ? 0 : A);

Funkcja porównania mniejszości < jest bardziej skomplikowana. Najpierw usuwamy pola z jednocześnie ustawionymi bitami argumentów
C = A&~(A&B); D = B&~(A&B);
Teraz w pętli zmniejszamy, dopóki jedno z nich się nie wyzeruje
while (C&D) { SAR C; SAR D; }
Po tych zabiegach można zwrócić wartość porównania: większe B lub zero
(A < B) = ( C ? 0 : B);

Wreszcie funkcja dodawania, która się zakończy zanim rekursja przejdzie cały bajt:
0+B = B; 
(A+B) = (SAL A&B) + ((A|B)&~(A&B)) = (SAL A&B) + (A XOR B);
najgorszym przypadkiem jest (-1)+1, kiedy rekursja przechodzi kolejno po każdym bicie.

Algorytmy te nie mają praktycznego znaczenia, gdyż jednostka arytmetyczno-logiczna procesorów jest dedykowana także do obliczeń arytmetycznych. Wykonuje je równie szybko i sprawnie jak logiczne.

Teoretycznie moduł arytmetyczny procesora nie jest potrzebny. Wystarczy logiczny. 

Schemat mnożenia, czyli funkcja jednoargumentowa iloczynu

Przy pisemnym mnożeniu binarnym mamy ciągłe inną dodawanie liczby i jej przesuwanie. Jeżeli mnożymy binarnie przez jeden, tworzy się pewien schemat. Zapamiętany schemat można odtworzyć biorąc jako argument dowolną wartość. Uzyskamy w ten sposób iloczyn argumentu oraz liczby danej schematem.

Kiedy występuje kilka jedynek pod rząd, np. '0111', zamiast trzech dodawań wygodniej jest dodać liczbę na pozycji poprzedzajacej jedynkę, oraz odjąć na pozycji ostatniej jedynki.

Opracowuje sposób uzyskania schematu. Zliczam długość binarnych ciągów zer i jedynek, oraz wprowadzam dodatkowe elementy + oraz - oznaczające dodanie, odjęcie argumentu.

Algorytm jest następujacy: dany jest ciąg symboli 0,1, poprzedzony zerami. Domyślnie jest to zapis binarny liczby od cyfr najmniej znaczących.
Mam flagę F, początkowo wyzerowaną oraz licznik C 
1) powtarzam 2)- 4), dopóki liczba dodatnia i flaga wyzerowana
2) jeśli ostatni bit jest wyzerowany i flaga ustawiona zapamiętuję '+'C,  zeruję C i kasuję flagę;
3) jeśli ostatnie dwa bity to '11' oraz flaga wyzerowana, zapamiętuję wartość licznika '-'C, zeruję licznik i ustawiam flagę; jeśli flaga była już ustawiona zwiększam tylko licznik
4) jeśli ostatnie dwa bity to '01' oraz flaga wyzerowana, zapisuję licznik '+'C, zeruję licznik; Gdy flaga jest ustawiona, zwiększam tylko licznik

Zapamiętany zostaje ciąg, np. dla liczby binarnej 100111101 jest nim +1+3-1+0. Zer nie musimy zapisywać, lecz należy je uwzględniać.
Sposób użycia schematu.
Startujemy od liczby lub +. W obu przypadkach jest to wartość argumentu. Napotykając +, mnożymy przesuwaniem liczbę o 2 i dodajemy argument. Kiedy napotykamy -, mnożymy przez 2 przesuwaniem i odejmujemy argument. Kiedy występuje liczba, przesuwamy mnożąc wartość przez 2 do potęgi tej liczby.
Schemat sprawdzamy, wstawiając za argument 1. Kiedy uzyskamy wartość  liczby, schemat jest prawidłowy.
Sprawdźmy zatem schemat +1+3-1+:
+1, 1 mnożymy przez 2^1 uzyskując 2;
+3, przesunięcie i dodanie 1 da 2*2+1 = 5, przesunięcie 3 powoduje pomnożenie przez 8 do 40;
-1, przesunięcie i odjęcie 1 powoduje uzyskanie 79, które mnożymy jeszcze raz przez 2 do 158
+[0], przesunięcie i dodanie 1 powoduje uzyskanie 317, ponieważ ostatnie przesunięcie jest zerowe, mamy wartość ostateczną 317 rowną binarnie liczbie 1 0011 1101.

Wartości liczbowe są zmniejszonymi o 1 krotnościami kolejnych zer, jedynek. Każdy zapis '+/- a' odpowiada następujacej instrukcji niskiego poziomu dla argumentu Y oraz wartości X: (X<<1 +/- Y)<<a .
Uzyskane liczby zawsze są dodatnie. 

Największy wspólny dzielnik podejściem wielowartościowym

W poprzednim poście podałem, że dla nieparzystej wartości a największy wspólny dzielnik nwd (ang. greatest common divisor GCD) spełnia warunek:
nwd(a,b) = nwd(a,2b).
Pozwala to dodać kilka dodatkowych warunków, by szybko liczyć nwd metodami wielowartościowymi.

Załóżmy roboczo, że a<b.
Metody do liczenia nwd są następujące: różnica b-a, reszta z dzielenia b/a, dopasowanie potęgi 2 tak, by
2^i*a <= b < 2^(i+1)*a
oraz całkowite
2^(j-1)*b <= a < 2^j*b .
Warunek całkowitości nie jest spełniony dla b nieparzystych, wtedy ten warunek sprowadza się do a<b.
Obliczamy moduł różnic uzyskanych wartości oraz odpowiedniej a, b oraz dokładamy dodatnią resztę z dzielenia. Uzyskamy wartość dodatnią c jako minimum uzyskanych wartości. Powtarzamy dla pary a, c.

Wartości zmniejszają się w kierunku wartości największego wspólnego dzielnika. Reszta z dzielenia może spaść do 0, wtedy zostaje zignorowana.

Złożoność postępowania jest złożonością najgorszej operacji: tu reszty z dzielenia. Dopasowanie potęg polega na porównaniu, przesuwaniu binarnym i odejmowaniu. Sprawdzana jest też parzystość, np. instrukcją 'b&1'.Wartości maleją szybciej niż logarytmicznie.


Przykład liczbowy, jak u Knutha Sztuka Programowania 4.5.2:
nwd( 20451, 6035 )
przesuwania: 2*6035 = 12070 < 20451 < 2^2*6035 = 24140;
20451/2 blokowany ułamek, pozostaje 6035 < 20451;
reszta z dzielenia: (2*(10000-6035)+451)%6035 = 2346;
lista: [8381, 3689, 14416, 2346]

Do dalszej iteracji przechodzi para ( 6035, 2346 ).
przesuwania: 2*2346 = 4692 < 6035 < 2^2*2346 = 9384;
6035/2 blokada, degeneracja do 2346 < 6035;
reszta z dzielenia: (6035 - 2*2346) = 1343;
lista: [1343, 3349, 3689, 1343]

Do dalszej iteracji rzechodzi para ( 2346, 1343 ).
przesuwania: 1343 < 2346 < 2*1343 = 2686;
2346/2 = 1173 < 1343 < 2346;
reszta z dzielenia: (2346-1343) = 1003;
lista: [1003, 340, 170, 1003, 1003]

Do dalszej iteracji przechodzi para ( 1343, 170 ).
przesuwania: 2^2*170 = 680 < 1343 < 2^3*170 = 1360;
1343/2 blokada, degeneracja do 170 < 1343;
reszta z dzielenia: 153;
lista: [663, 17, 1173, 153]

Do dalszej iteracji przechodzi para ( 170, 17 ).
przesuwania: 2^3*17 = 136 < 170 < 2^4*17 = 272;
170/2 = 85 >17 blokada na 17 < 85;
reszta z dzielenia: 0 blokada;
lista: [34, 102, 68]

Do ostatniej iteracji przechodzi para ( 34, 17 ).
Teraz obliczenia blokują się bardzo często zerem, zostawiając [17, 51] jako  wartości listy.  Po wymianie mamy parę (17, 17).
Zatem 17 =  nwd(20451, 6035).

U Knutha było koło 8 iteracji, tutaj 6.

Czy można pominąć część dzielników podczas faktoryzacji

Przed burzami miewam czasem ciekawsze pomysły. Połączyłem metodę faktoryzacji, która sprawdza od razu duże i małe dzielniki opublikowaną 7 sierpnia 2013 roku
'Faktoryzacja, przepis z prostymi przekształceniami'
z jedną z pierwszych faktoryzacji przez konwersję na inny system.
Okazało się wtedy, że nie muszę sprawdzać nie tylko liczb parzystych jako kandydatów na dzielniki (po wyłączeniu 2 na pewno nie są dzielnikami), ale mogę nawet pomijać wszystkich kandydatów mniejszych niż około 1/3 pierwiastka kwadratowego kosztem około log_3 n obliczeń nwd.

Łącznie można pominąć n/6 wartości szukając dzielników liczby n. Dzielnik znajdziemy przeliczając nie więcej niż trzykrotność jego wartości.
W poniższym sposobie dzielnik d pojawia się co d wartości, na kolejny ślad obecności trafimy nie później niż zmiejszając podstawę systemu o 3d.
Dodatkowo, dla bardzo dużych liczb jesteśmy w stanie łatwo oszacować na podstawie reszty, że nie ma blisko dzielników, lecz nie można tego łatwo sformalizować. Praktycznie można przyspieszyć stosując w niektórych przypadkach konwersję nie o 2, lecz o 6. Tutaj ignorowane.

Lemat 1. Jeśli k jest względnie pierwsze z a, to
  nwd(a,b) = nwd(a, k*b).
Jeśli jest wspólny dzielnik d = nwd(a,b), to d dzieli a oraz b, zatem dzieli także k*b.
Jeśli d = nwd(a, k*b), to d dzieli równocześnie a oraz k*b. Skoro a oraz k są względnie pierwsze, zatem d dzieli b. Stąd nwd(a,b) = d.

Lemat 2. Jeśli n = ap^2+bp+c, to dzielnika n można szukać jako nwd(p,c) dla pewnego p.
Zapiszmy n nieco inaczej n = (ap+b)p+c. Liczba się dzieli przez d, gdy c się dzieli przez d oraz iloczyn (ap+b)p się dzieli przez d. Ponieważ p jest modyfikowane, można się ograniczyć do szukania tylko nwd(p,c).

Ciekawe rezultaty są wtedy, gdy d jest małą liczbą, np. 3. Wtedy co trzecia podstawa jest podzielna przez 3, co piąta jest podzielna przez 5. Po wyeliminowaniu 2 jako dzielnika podstawy systemów są nieparzyste.

Konwersja na system p, kiedy podstawa systemu jest równa p+2 dla liczb 'trójcyfrowych' n = ap^2+bp+c wygląda następująco:
a  b+2a
a  b+4a  c+(b+2a)
co można zapisać równaniami:
a'=a;
b'=b+4a;
c'=c+(b+2a);
Jeśli którakolwiek z b', c' jest większa niż p, wykonujemy operacje (c'-=p; b'++;) oraz (b'-=p; a'++;) odpowiednio. Nie więcej niż dwukrotnie każda.

W praktyce, zaczynając od podstawy będącej liczbą nieparzystą mniejszą od pierwiastka z n, mamy tak długo do czynienia z trzema 'cyframi', że trzecia część pierwiastka kwadratowego z n jest jeszcze liczbą 'trójcyfrową'.  

             Algorytm:
Inicjacja: niech p będzie największą liczbą nieparzystą mniejszą niż pierwiastek kwadratowy z liczby nieparzystej n. Przedstawiamy liczbę n = p^2+bp+c, gdzie b jest mniejsze niż 4 oraz -1<c<p. Zapamiętujemy podłogę z p/3 jako x.
Powtarzamy w pętli dopóki nie znajdziemy dzielnika:
Jeżeli c jest zerem, dzielnikiem jest p.
Jeśli p jest podzielne przez 3, szukamy liczby q takiej że p=q*3^i, gdzie i jest całkowite oraz 3 nie dzieli q. Obliczamy d = nwd(c,q). Jeśli d>1, d jest dzielnikiem.
Zmniejszamy p o 2 i porównujemy z x. Jeśli p>x wracamy do pętli, w przeciwnym przypadku n jest pierwsza.

Ponieważ co szóste p jest nieparzyste podzielne przez 3, w algorytmie przyjmujemy za p tylko wartości nieparzyste, zatem co trzecie p jest podzielne, a na dodatek niech t=p/3. Wtedy co trzecią iterację mamy wartość t-2, zapamiętując je nie musimy dzielić przez 3 więcej niż raz, a następnie zastosować funkcję rekursywną:  t[0]=p/3, t[j+1] = 1/3 t[j],
zmienna wewnętrzna i odmierza co 3 wejście, inicjowana jest jako reszta i=2 (t[j] = 3r+1), i=1 (t[j] = 3r+2), i=0 (t[j] = 3r), r naturalne;
new f( t[j] ): { t[j+1] = floor t[j]/3; i=3-(t[j]%3)); }, t[j]>1
tworzy t[j+1].ty element. Funkcja zwraca dodatnią wartość nwd albo ujemną wiadomość 'jest reszta z dzielenia przez 3':
f (t[j])(i) = {
  i--;
  if( 0<i) return -i; else {
    i = 3;
    t[j] -= 2;
    k = f(t[j+1]);
    if(0<k) return k; else return nwd(t[j],c);
  }
};

Złożoność.
Konwersje mają złożoność arytmetyczną stałą, pamięciowo liniową.
Obliczanie nwd ma najgorszą złożoność, i to ono decyduje o złożoności algorytmu. Przy okazji, wartości dla nwd są ograniczone, jedna przez p, druga przez p/3. Obydwie mniejsze niż pierwiastek z n.


Kontrprzykład numeryczny, liczba pierwsza n = 1447 = 38^2+3
a  b  c   p   t[j]  nwd(d,c)
1  0  3   38
1  2  4   37    // algorytm
1  6  12 35
1 10 28 33 t[0]=11 nwd(11,28)=1
1 15 21 31
1 20 26 29
1 26 16 27 t[0]=9 t[1]=3 t[2]=1 nwd(1,16)=1
2  7  22 25
2 16 11 23
3  5  19 21 t[0]=7  nwd(7,19)=1
4  0  3  19
5  0  2  17
6  6  7  15  t[0]=5  nwd(5,7)=1
8  7  4  13
Wszystkie możliwe liczby nieparzyste uwzględnione, od 13 do 37 przy p, mniejsze przy liczeniu nwd. Jest przeliczonych 13 spośród 19 wartości nieparzystych (2/3 oraz nadmiar). Pozostałe uwzględnione są w czterech nwd. 

Wskażnik na funkcję, ale nie na każdą

Spróbowałem, czy można wymieniać ciało metody klasy za pomocą wskaźnika.

Pierwsze próby robił Koza programowaniem genetycznym w LISPie. (Algorytmy genetyczne + struktury danych = programy ewolucyjne Z. Michalewicza)
Kompilator gcc nie przepuścił. Traktował jako konwersję 'int (*)(int)' na 'int (obiekt::*)(int)' i nawet reinterpret_cast nie pomógł. Nie udało się wewnątrz klasy odwoływać się do jej metod wskaźnikiem. 

Kilka dni i poddałem się. Nie zamierzam rezygnować z obiektów. Ale mam coraz większą ochotę blokować lub fałszować sprawdzanie typów.
Czyli coraz dalej od standardów.

Na razie sposobem na zmianę ciała funkcji w czasie działania programu w C++ wydaje mi się interpretacja ciągu tekstowego funktorów.