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)


14 marca 2014

Spacerkiem koło dzielników

Rozwijałem pomysł odwrócenia iloczynu
N = q1*q2
gdzie wszystkie wartości są traktowane schematem Hornera. Oznacza to dla wartości zapisanych trzycyfrowo
q1=(a2*p1+a1)*p0+a0,
dla dowolnych liczb p1, p2 i ciągach współczynników (a), (b) zapis:
N = ((a2*p1+a1)*p0+a0) * ((b2*p1+b1)*p0+b0) =
(( a2*b2*p1*p0 + p0*(a2*b1+a1*b2) + (a2*b0+a0*b2) )*p1 +
(a1*b1*p0 + (a1*b0+a0*b1) )*p0 + a0*b0 =
n*p1*p0 + w[a1,b1]*p0 + a0*b0

Jeśliby udało się jednoznacznie dopasowywać kolejne współczynniki ciągów 'cyfr' a, b, mamy do czynienia z algorytmem działającym na coraz mniejszych wartościach, w którym wyszukiwanie zostałoby przerzucone na znajdywanie odpowiednich par (a,b) stanowiących kolejne wyrazy.

Owszem, wielomian w[a,b] pojawiający się we wzorze jest łatwo wyznaczany, współczynniki znajdujemy z równań kongruencyjnych, których liczność zależy od wartości współczynnika danego poziomu p. Jest jednak duży kłopot z jednoznacznością.
Można temu częściowo zaradzić, biorąc za p niewielką liczbę pierwszą. Wtedy mamy do rozwiązania p-1 kongruencji, z których jedna jest sprzeczna, jak to pisałem w poprzednim poście.
Nauczyłem się także wskazywać, które z par (a,b) nie warto badać przy szukaniu dzielników. Zatem ciągi (a), (b) powinny z pomocą (p) skonwertować się na dzielniki. A wtedy pojawia się pułapka.

Oto heurystyka, która pozwala przybliżać orbity dzielników zapisanych jako a*P+A,
gdzie P jest iloczynem wartosci p0*p1*p2*..., A jest liczbą uzyskaną z obciętego ciągu (a0, a1, ..., am) wzorem A=a0+p0*(a1+p1*(...(am)...)).
Inicjujemy n=(N-1):2, A=1, B=1, P=2, w[0,0]=1. 

Dane są: liczba n, wielomian w[a,b] = abP+aA+bB o współczynnikach P, A i B
Wybieramy liczbę pierwszą p, dbając, by p nie dzieliło n.
Dzielimy z resztą n = p*q+r, 0<r<p.
Szukamy rozwiązań równania kongruencyjnego  jako par liczb (a,b)
w[a,b] = r (mod p)
po prostu wstawiając kolejno a=0, ..., a=p-1. Uzyskamy wartość b.
Tworzymy równanie  zmiennych C=A, D=B
q = (A+aP)*C+(B+bP)*D
które traktujemy modulo min(A,B) uzyskując równanie rekurencyjne jednej zmiennej. Jeżeli rozwiązanie tego równania istnieje i jest (w kolejności) mniejsze niż p, równe 0, wielokrotnością p, w takim razie mamy dużą szansę, że para (a,b) stanowi lepsze przybliżenie dzielników.
Modyfikujemy w kolejności: A=A+aP, B=B+bP, P=p*P, n = (n-w[a,b]):p

Powtarzamy tak długo, dopóki n jest dodatnie.
Wartość n=0 oznacza, że udało nam się dokładnie wyznaczyć ciągi współczynników dzielników, czyli A, B są dzielnikami N.

A oto pułapka. Może się zdarzyć, że któryś ze współczynników A+aP, B=bP jest dzielnikiem oraz n>0, wtedy algorytm się nie zatrzymuje, lecz szuka dalej.

Przykład numeryczny, 8934053
n = (8934053-1):2 = 4467026
p = 3, w[a,b] = 2ab+a+b
rozwiązaniem równania w[a,b]=2 są pary (2,0), (0,2), zaś (1,1) daje równanie sprzeczne. Przyjmujemy a=(2,0) dzieki równaniu (n-2):3 = 5A+B
nowy wielomian w[a,b] = 6ab+a+5b
Dalej p=5, równanie dla a=4 jest 4b+4=3 (5)
Podaję w[a,b] oraz p
w[a,b] = 30ab+7a+29b, p=7, rozwiązanie (0,1) z D=4
w[a,b] = 210ab+37a+29b, p=3, rozwiązanie (0,2) z C=21=7*3
w[a,b] = 1470ab+457a+29b, p=7, rozwiązanie (0,6) z D=7
i tu następuje ominięcie, gdyż właściwym rozwiązaniem jest (6,1), wtedy współczynnik przy D jest dzielnikiem 1087 liczby 8934053 oraz następne n jest równe n=0.
Biorąc liczbę pierwszą większą od 11, też uzyskamy rozkład na dzielniki po tych sześciu iteracjach. A=8219, B=1087.

28 lutego 2014

Próba przyspieszenia znajdowania dzielników


Każdą liczbę naturalną można jednoznacznie przedstawić jako sumę uporządkowanego ciągu iloczynów współczynnika i potęgi podstawy systemu 
\sum_{k>0} a_k*p^k
(budowa systemu pozycyjnego, dla p=10 dziesiątkowego). 

Nic nie stoi na przeszkodzie, by potęgi podstawy systemu (p_k)_k zastąpić losowym, lecz uporządkowanym ciągiem liczb, np. kolejnych liczb pierwszych.
Wtedy liczbę naturalną N jednoznacznie przedstawiamy wyrażeniem, np.
127 = ((1*7+1)*5+2)*3+1
 za pomocą ciągów [1,1,2,1] oraz ciągu liczb pierwszych [7,5,3].
Rozpatrując to wyrażenie w ogólnej postaci, kiedy N jest liczbą złożoną, uzyskujemy dla liczb p i q pierwszych oraz szukanych współczynników a,b,c,d,e,f następującą postać:
N = ((ap+c)q+e) * ((bp+d)q+f) =
[[abpq + (ad+cb)q + (af+eb)]p + (cdq+(cf+ed))]q + ef 

Zauważmy, że w każdym bardziej zewnętrznym nawiasie prostokątnym potrzeba coraz mniej współczynnikow. Zatem znalezienie pary (e,f), z jej pomocą pary (c,d) oraz ostatecznie pary (a,b) wyznaczy nam ciągi zapisu dzielników N. 
Wyrażenie jest przedstawione jednoznacznie (własności działań po utożsamieniu wyrażeń przemiennych ab = ba).
Kiedy wylosujemy dowolną liczbę q, wyznaczamy parę (e,f). Wstawiamy 
n_1 = (N- ef)/q 
oraz kontynuujemy z mniejszą wartościa n_2 = n_1*p+(cdq+cf+ed) = n_1*p+w(), w której nieznane są tylko c i d. Z niej wyznaczamy parę (c,d) kongruencjami 
w() = (n mod p) modulo p. 
Wielomiany w() dla kolejnych poziomów modyfikują się w bardzo specyficzny sposób. Są to równania kongruencyjne nieliniowe dla największego indeksu k, kiedy to mamy iloczyn obu współczynników oraz wszystkich liczb pierwszych ciągu. Do tego składnika dodawane są liniowe fragmenty (a_k*b_i + a_i*b_k) mnożone przez iloczyn liczb pierwszych na pozycjach nie przekraczajacych i<k. Współczynniki a_i oraz b_i zostały już wyznaczone wcześniej. 

Jeśli liczba p jest pierwsza, znajdziemy p-1 par (c,d) takich, że c*d = n%p (mod p) dla kolejnych wartości c. Jedna z kongruencji jest w praktyce sprzeczna.
Przemienność może nam popsuć wybór, dlatego reszta z dzielenia przez pierwszą z ciągu liczb pierwszych nie powinna być kwadratem a^2. W przeciwnym przypadku mamy dwie możliwości: (a,a) oraz (-a,-a). Zapewni to jednoznaczność wyboru par (c,d) dla różnych c, z których w przypadku liczby złożonej dokładnie jedna staje się częścią ciągów dzielników N. 

W klasycznej metodzie dzielenia przez kolejne liczby mamy przeszukiwanie liczb pierwszych. Tutaj ten proces został przerzucony na znajdywanie rozwiązań równań kongruencyjnych, co odbywa się na stosunkowo małych liczbach - zależnych od wyboru q.
W ogólności tworzymy drzewo z krotnością gałęzi danego poziomu nie wiekszą niż wybrana liczba pierwsza. Drzewo to przyda się przycinać przez wybór odpowiedniej pary.
Aby zapobiec kłopotom z wyborem, warto uzywać liczb pierwszych jeszcze nie występujacych w ciągu. W przeciwnym przypadku kongruencje będą się upraszczać kosztem większych kłopotów z doborem pary. Tak działał jeden z wcześniejszych algorytmów faktoryzacji.

Zauważyłem pewne kryterium, działające dla początkowego wyboru par.
Wartość n modulo kwadrat wybranej liczby pierwszej jest najbliższe reszcie z dzielenia różnicy n-w() przez iloczyn początkowych liczb pierwszych ciagu  z niedomiarem, lecz to kryterium w pewnym momencie staje się fałszywe. Nie udaje mi się znaleźć lepszego, wszystie w pewnym momencie zaczynają wskazywać inne pary niż właściwa.

Przykład numeryczny N = 8934053.
Ciąg podstaw w kolejności od cyfr bardziej znaczących: [13, 11, 7, 5, 3].
Ponieważ N = 2 (mod 3), przyjmujemy A=[2], B=[1], czyli para (a,b) = (2,1).
Do dalszego przebiegu używamy n = (8934053-2*1)/3.
Pobieramy nową liczbę pierwszą q=5, modyfikujemy wielomian 
w() = 3ab+(a+2b) 
i szukamy kolejnych par w() = n%5 (mod 5). Wybieramy tę, która jest najlepsza, tu okazała się nią para (4,2). Wstawiamy wartości (a,b) = (4,2) wyznaczając wartość wielomianu w(), którą odejmujemy od n. Kontynuujemy.

Podsumowując cały rozkład.
Tablice dla dzielników są postaci: 
A = 8219 = [7, 1, 1, 4, 2]; 
B = 1087 = [0, 10, 2, 2, 1].
Ostatnia para (a,b) = (7,0) (zapisana jako pierwsze elementy w tablicach ciągów) jako jedyna wsród kandydatek byla równa n przez wstawienie do wielomianu tego poziomu
w() = 11*7*5*3ab + 105(10a+b) + 15(2a+b) + 3(2a+4b) + (a+2b)
Zakończyła ona proces znajdowania ciągów przez wyzerowanie n.

Liczby pierwsze nie wyzerują n, oraz dla kongruencji nie znajdziemy rozwiązań w małych liczbach naturalnych, ograniczonych liczbami pierwszymi.

Nie potrafię policzyć złożoności, czy proces jest szybszy niż standardowe dzielenia. Czy przycinanie drzew jest wystarczająco dobre.
Jest to po prostu kolejny sposób spojrzenia na faktoryzację.

30 stycznia 2014

Metoda 'kolejnych dzieleń' od końca

Przyglądanie się zmianom wyglądu liczb zapisanych jako iloczyny małych wartości oraz liczby w systemie Fibonacciego spowodowały powstanie kolejnego algorytmu na poszukiwania dzielników.
Złożoność jest podobna jak przy przeskakiwaniu na kolejny system minimalizujący resztę, kiedy liczba jest dwucyfrowa. W odróżnieniu jednak od tamtego algorytmu można przeskakiwać część wartości, oraz sprawdzanych jest około połowa dopuszczalnych wartości - liczb od 1 do pierwiastka z n.

Założenie. Niech liczba nieparzysta n będzie przedstawiona wzorem
n = p*q+r,
gdzie p jest liczbą nieparzystą ograniczającą pierwiastek kwadratowy n z góry, analogicznie q jest liczbą bliską pierwiastkowi z n.

Algorytm jest następujący:
while (1<p) {
  p-=2;
  r+=2*q;
  q+= r/p;
  r = r%p;
  if( 0=r ) return dzielniki(p,q);
}
return liczba pierwsza; 

Wartość p maleje liniowo do 1, q rośnie skokami w zależności od części całkowitej ilorazu r/p. Ponieważ startujemy z wartości bliskich sobie, iloraz r/p początkowo jest równy 2 i bardzo powoli rośnie (niemonotonicznie). Dopiero kiedy sprawdzimy wszystkie dzielniki tej samej krotności cyfr dziesiętnych, iloraz ten przekroczy 20 (dzielenia rzędu 19/1). Reszta równa 0 oznacza, że mamy do czynienia z dzielnikami n = p*q. Przechodzimy tylko po wartościach p nieparzystych, gdyż dla parzystych nie ma dzielników przy nieparzystym n. 

Dla bardzo dużych wartości można dostrzec dodatkową własność, która pozwoli pominąć część przypadków. Z reszty ubywa elementów według ciągu rozbieżnego, który lokalnie daje się wyznaczać. 

Zatem sposób ten nadaje się jako konkurent rozkładu metodą dzielenia z resztą przez kolejne liczby nieparzyste (ta sama krotność przypadków), przy czym startuje od bardzo dużych wartości. Dodatkowo te ogromne są nieco prostsze w liczeniu.

Przykład numeryczny, 8934053 = 2989*2988+2921
p = 2989, q = 2988, r = 2921,  (r+2q)/p = 2
na początku nowego przebiegu pętli mamy wartości
p = 2987, q = 2990, r = 2923,  (r+2q)/p = 2
blisko pierwiastków
p = 1091, q = 8188, r = 945,  (r+2q)/p = 15
p = 1089, q = 8203, r = 986,  (r+2q)/p = 16
p = 1087, q = 8219, r = 0
rozkład 8934053 = 1087*8219.
W tym przypadku trzeba sprawdzać 952 przypadki zamiast 543 jak przy metodzie kolejnych dzieleń ze względu na położenie dzielników.

Wyjaśnienie (próba): liczba zapisana systemem Fibonacciego np. q=100100101
jest przekształcana do iloczynu (hybryda)
p*q = p00p00p0p.
Odjęcie 2 na każdej niezerowej pozycji powoduje zmniejszenie iloczynu o 2q, który jest przedstawiany jako nowa hybryda liczby Fibonacciego z taką samą wartością (p-2) na wartościach niezerowych oraz resztą r trzymaną oddzielnie.
Przechodząc z systemu Fibonacciego na dziesiątkowy powstaje ww algorytm.

07 stycznia 2014

Systemy niedziesiątkowe, wygląd liczby złożonej wskazującej dzielniki

Końcówka roku upłynęła na oglądaniu liczby złożonej przedstawianej w różnych systemach liczbowych, celem poszukiwania możliwości zmniejszenia krotności rozpatrywanych przypadków.

Kilka z algorytmów publikowanych na tym blogu korzysta z przedstawienia liczby w systemie o podstawach nieparzystych. Wtedy liczba złożona n=p*q w systemie o podstawie q ma cyfrę jedności równą 0.
Należy sprawdzać do połowy z pierwiastka kwadratowego z n przypadków.

Kiedy podstawy są parzyste, można oprócz konwersji na systemy p-2, p+2 łatwo przechodzić na systemy 2p, p/2. Miałem nadzieję, że manipulując tymi czteroma przekształceniami uda się dotrzeć do dzielników.
Przekształcenia działają, pojawiają się dodatkowe postacie błyskawicznie zwracające dzielnik, np:
n = a*(p-1)+a, tzn. naprzemienna suma cyfr dzieli p;
n = a*(p+1)+(p-a), tzn. suma cyfr dzieli p;
Wzory działają także dla liczb większych niż dwucyfrowe. Trudno je jednak znaleźć bez przeszukiwania.

Jeżeli mamy szczęście, to dodatkowo suma (naprzemienna) cyfr spełnia ten warunek także dla podstaw będących wielokrotnościami o potęgi całkowite z 2.
I na tym się kończy. Wiadomo jeszcze, że w pierwszej z tych postaci a jest nieparzyste (bo p-1 jest parzyste oraz wartość sumy ma być nieparzysta).
Myślałem, że dla kolejnych podstaw (p-1)*2^k wystarczy brać najmniejsze z możliwych a (nieparzyste), ale przeskakiwana jest wtedy szukana postać.

Drugie przypuszczenie: sumy i sumy naprzemienne majace w dzielnikach określoną postać w sąsiedztwie będą zbiegać do właściwych. Lecz nie do końca. Pojawiają się inne małe wartości dla sum bliskich szukanych z daleka od dzielników. Także przy samych dzielnikach wartości nie są monotoniczne.

Liczność przypadków może być nieco mniejsza, ale i tak wymagane jest sprawdzanie wielu przypadków (podstawy parzyste). Czasem widać dzielnik po przeliczeniu zaledwie połowy potrzebnych przypadków.

Nasuwa mi się możliwość zastosowania w tym przypadku specyficznego dodawania + mającego wiele argumentów. Jeden z argumentów jest dany relacją jednoargumentową '< r()', która zachowuje sumę tak długo, dopóki jest mniejsza od r, jeśli jest większa, blokuje dodawanie zwracajac wartość pustą NIL. Np. +(1, 2, 3, 4, 5, <7) zwraca NIL, bo 1+2+3+4>7. Z kolei +(1,2,3, <16, 4, 5) = 1+2+3+4+5=15<16 zwraca 15.