19 marca 2015

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

26 lutego 2015

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)

03 grudnia 2014

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.

12 listopada 2014

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.

29 września 2014

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.

28 sierpnia 2014

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.

04 sierpnia 2014

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.

25 czerwca 2014

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. 

16 czerwca 2014

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. 

02 czerwca 2014

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.

30 maja 2014

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. 



21 maja 2014

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.

10 kwietnia 2014

Algorytm sortowania quicksort ze stanem

Algorytm sortowania szybkiego ma kilka kiepskich przypadków. Należy do nich już posortowana tablica (P. Wróblewski, Algorytmy, struktury danych i techniki programowania).
Wykorzystując tablicę stanów można tak zmodyfikować ten algorytm, że łączna krotność zamian elementów może sie zbliżyć do sortowania przez selekcję. Zaś sortowanie przez selekcję ma najmniejszą możliwą liczność zamian (R. Sedgewick, Algorytmy w C++). Ta wersja quicksortu nie zawsze osiągnie tę granicę, lecz liczność zmian mocno się zmniejsza.

Algorytm  quicksort dzieli tablicę na trzy części. Pobiera jeden z elementów oraz szuka dla niego docelowego miejsca w tablicy. Elementy poprzedzające są mniejsze niż rozpatrywany, następne są większe. Następnie quicksort wywołuje sam siebie ponownie dla tych części.

Moja modyfikacja sprawdza, czy tablica nie jest już posortowana. Dla wybranego elementu sprawdza, czy nie ma przypadkiem wartości równych danej.
Wprowadzam funkcję, która dla zadanej wartości zwraca funktor <, = albo >, związany ze stanem. Parametrem funktora jest element tablicy. Ponowne wywołanie tej funkcji usuwa funktor wydobywając zawartość - sortowany element.
Oczywiście dla tablicy długości poniżej 3 elementów użycie funktora jest niepotrzebne - sortowanie następuje od razu. 

Oto stany algorytmu quicksort oraz ich interpretacja:
1 - wszystkie elementy są równe,
2 - elementy tablicy tworzą ciąg niemalejący, porównanie jest z elementem poprzedzającym, nie z wybranym,
3 - elementy tablicy tworzą ciąg nierosnący, porównanie jest z elementem następnym, nie z wybranym,
4 - elementy tablicy nie są posortowane, oraz począwszy od pewnego miejsca nie ma elementów równych,  porównanie z elementem wybranym,
5 - elementy tablicy nie są posortowane, występują elementy równe elementowi aktualnie rozpatrywanemu, porównanie z elementem wybranym.

Dla pierwszych trzech stanów, wystarczy do jednego liniowego przejścia, by mieć posortowaną tablicę. Przerobienie tablicy w stanie 3 na listę posortowaną rosnącą wymaga użycia funkcji reverse, odwracajacej kolejność tablicy.
Stan czwarty odpowiada standardowemu podejściu do quicksort. Natomiast stan 5 wymaga dodatkowego przejścia zcalającego elementy równe sobie, też złożoności liniowej.   

Tablica stanów, kiedy porównanie dwu elementów a, b zwraca relację <, =, >. Funktor > jest zwracany, jeśli bieżący element jest większy niż porównywany lub równy w stanie 2, funktor < dla elementu mniejszego lub równego w stanie 3.
Stan  <    =    >, drugi z symboli to wartość (funktor)
1       3<  1=  2>
2       ?     2>  2>
3       3<  3<  ?
4       4<  5=  4>
5       5<  5=  5>
Obszary oznaczone pytajnikiem są oprócz zmiany stanu na 4 albo 5 są  powieleniem porównania, raz dla elementu poprzedzajacego, docelowo dla wybranego. Jeśli w ciagu słabo rosnącym jest element mniejszy od przedostatniego, lecz większy od wybranego, następuje zwrot funktora >, dla równego elementowi wybranemu zwrot funktora =. Jeśli porównywany element jest mniejszy zarówno od elementu poprzedzającego, jak i od wybranego, zwracany jest funktor <. Podobna sytuacja występuje dla ciągu słabo malejącego, wartość i stan jest nadawana w zależności od drugiego porównania.
Teraz w tabeli mamy ciąg funktorów. Jeśli mamy stan 5, blokujemy ponowne wywołanie funkcji do następnego przejścia. Teraz poruszając się od krańców tablicy zamieniamy ze sobą funktory > z początku tablicy oraz < z końca tablicy. W stanie 4 od razu od razu stosujemy ponownie funkcję, by odzyskać wartości elementów. Ostatnim krokiem jest zamiana elementu wybranego z ostatnim elementem o wartości <, wskazywanym indeksem, który służy do podziału tablicy do kolejnych iteracji.
W stanie 5 należy jeszcze raz przejść tablicę wymieniając elementy oznaczone relacją = z przyległymi do docelowego miejsca, modyfikując zakres tablicy. Elementy w środkowej części tablicy są wtedy na swoich miejscach i nie będą dalej ruszane.
Mamy tablicę, w której element wybrany jest na docelowej pozycji, mniejsze go poprzedzają, większe następują po rozpatrywanym elemencie.
Następuje wywołanie rekursywne dla obu części.

Przykład:
A S O R T I N G E X A M P L E
pierwszym, porównywanym elementem jest A, porównujemy w stanie 1:
S>A, stan 2, zwracamy >,
O<S, porównujemy O>A, stan 4, zwracamy >,
R>A, zwrot >
jeszcze parę wyrazów
A=A, stan 5, zwrot =
dokańczając porównania uzyskujemy tablicę w stanie 5
= > > > > > > > > > = > > > >
wyszukanie relacji < daje wartość pustą, zcalamy wartości z relacją = stosując funkcję. Uzyskujemy tablicę (jedna zamiana), z której można odciąć jeszcze jeden element A:
A | A O R T I N G E X S M P L E
wywołujemy quicksort ponownie. Po jednej stronie ciąg pusty czyli nie ma co robić, po drugiej
O R T I N G E X S M P L E
porównujemy w stanie 1:
R>O, stan 2, zwrot >
T>R, zwrot >
I<T, porównujemy I<O, stan 4, zwrot <
N<O, zwrot <
G<O, zwrot <
E<O, zwrot <
X>O, zwrot >
S>O, zwrot >
M<O, zwrot <
P>O, zwrot >
L<O, zwrot <
E<O, zwrot <
nowa tablica w stanie 4
= > > < < < < > > < > < <
zamieniamy ze sobą skrajne relacje (>,<) czyli: (R,E), (T,L), (X,M) uzyskując
O E L I N G E M | S X P T R
ostatnia zamiana elementów (O,M)
M E L I N G E | O | S X P T R
oraz ponowne wywołanie dla MELINGE oraz SXPTR itd.

26 marca 2014

Faktoryzacja metodą kolejnych dzieleń z cechami podzielności

Wymieniłem sposób dzielenia w faktoryzacji metodą kolejnych dzieleń.
Zastosowałem dzielenie, a które naprowadziły mnie systemy niedziesiątkowe.
Dla bardzo dużych liczb n, długości 2i lub 2i-1, dzielę n na dwie paczki o długościach nie przekraczających i:
n = a ## b = a*10^i+b
Cecha podzielności przez odpowiednio duże liczby jest postaci:

"liczba n jest podzielna przez p=10^(i+1)-r wtedy i tylko wtedy gdy wartość a*r+b jest podzielna przez p";

istnieje też podobna cecha, gdy p=10^i+r, cecha lepiej działa, gdy r nie jest większe niż 2*10^(i-1):

"liczba n jest podzielna przez p=10^i+r wtedy i tylko wtedy, gdy wyrażenie  -a*r+b jest podzielne przez p".

Co się wtedy dzieje dla bardzo dużych liczb, typu RSA?
Otóż w zależności od przedziału p \in (10^i/(k+1) , 10^i/k) cecha podzielności może zmniejszyć wartości biorące udział w obliczeniach. Mianowicie przyrost r można zastąpić wartością 10^i % r, co powoduje zmniejszenie czynnika, czasem bardzo duże.

Przyjmujemy p nieparzyste. Jeśli p>5*10^(i-1), to cecha pozwala liczyć wyrażenie ar+b, zaś dla liczby p-2 już tylko a(r+2)+b. Oznacza to, że w tym przedziale nie trzeba wykonywać iloczynu a*(r+2), tylko do poprzedniego dodać 2a.
Gdy p \in (0,(3)*10^i , 0,5*10^i ), dla porachowania wartości w (p-2) należy do wcześniejszego iloczynu dodać 4a. I tak dalej, każde kolejne zmniejszenie poniżej 10^i/k powoduje konieczość dokładania 2k*a.

Dodatkowo, możemy jeszcze liczyć kongruencje modulo p iloczynów.
Nie trzeba także sprawdzać wartości 10^i/k. Pojawia się ona automatycznie, kiedy wartości r przekroczą p.
Dla tych największych wartości nie warto stosować drugiej z cech (suma naprzemienna) ze względu na pogorszenie sytuacji obniżającej wykładnik potęgi. Wtedy warto używać działań modulo dla iloczynów, choć można wygenerować też równanie kwadratowe, którego wartość jest dodawana. 
Jednak, nie znalazłszy dzielnika dla tych największych wartości, można kontynuować, stosujac cechę podzielności w innym obszarze:
n = a ## b ## c = a*10^(2j)+b*10^j+c
"liczba n jest podzielna przez p = 10^(j+1)+r, r = (10^(j+1) % p) wtedy i tylko wtedy gdy (ar+b)r+c jest podzielne przez p"
"liczba n jest podzielna przez p = 10^j+r, r<10^j wtedy i tylko wtedy gdy -(-ar+b)r+c jest podzielne przez p"
To nie jest już tak wygodne jak wcześniejszy przypadek, ale nadal prosty.

Modyfikacja ta pozwala upraszczać od 1% do 90% przypadków (związane z wartością liczby n) w najgorszym obliczeniowo obszarze. Dodatkowe wykorzystanie podziału liczby n na trzy części pozwala sprawdzić około 90-99% przypadków.

Przykład numeryczny: n = 8934053 = 893*10^4 + 4053,
podzielność przez p=2501: r = 10^4 % 2501 = 2497
reszta n%p = (893*2497+4053)%p = (2233874)%p = (1430+1552)%p = 481
kolejne dodanie 6+2497 > 2500, następuje zmiana wielokrotności dodawanych a
podzielność przez p=2499: r = 10^4 % 2499 = 4
reszta n%p = (893*4+4053)%p = 7625%2499 = 128
podzielność przez p=2497: r = 10^4 % 2497 = 12 (1/4 k: dodać 2*4a=8a)
reszta n%p = (893*12+4053)%p = 2284  = (8*893+7625)%2497

podzielność przez 1001:  r = 10^4 % 1001 = 991
reszta (893*991+4053)%1001, ale także
-(-8+934)+53 = -873 = 128 (1001)