Mija miesiąc od ostatniego posta, czas na kolejny algorytm. Ale nie ma żadnego nowego mimo wytężonej pracy.
Przez ten czas przyglądałem się zastosowaniem wzorów skróconego mnożenia do przekształceń liczb. Ponownie odkryłem metodę Fermata a^2-b^2 w trzech odsłonach:
- klasyczną (np. u Knutha Sztuka programowania 4.5.4.C), kiedy odjemna jest zawsze kwadratem, zaś zwiększam b czasem zmniejszając a;
- odwróconą, kiedy odjemnik jest zawsze kwadratem, wtedy zmniejszamy b od czasu do czasu zwiększając a;
- mieszaną, kiedy naprzemiennie przenoszę jakąś wartość między odjemną a odjemnikiem pilnując, by zawsze co najmniej jedno z nich było kwadratem.
Algorytmy te zapisane za pomocą systemów niedziesiątkowych w obszarze z dużymi a oraz b sprowadzają się do prostych konwersji typu
1'2'a_p = 1'0'(a-1)_{p+1}
co oznacza, że dla bardzo dużych liczb działamy na małych wartościach. Nie zmienia to jednak złożoności.
Sprawdziłem także zachowanie się sumy sześcianów. Co prawda, nie każda liczba wyrazi się jako suma dwu sześcianów, ale przekształcenia wzięte z różnicy kwadratów działają i dla sześcianów.
Kiedy zbliżamy się do dzielników, jakiekolwiek przekształcenia wskazujące: 'jesteśmy coraz bliżej' zaczynają szwankować. Dlatego lepszą grupą są algorytmy szukajace dzielników mimochodem. Skaczące między różnymi podstawami.
W szczególności zastosowanie cechy podzielności przez (znany) dzielnik d w systemie o podstawie d+1 (cecha: różnica naprzemienna cyfr w zapisie liczby n jest 0 lub jest podzielna przez n) liczby trójcyfrowej n = a'b'c_p prowadzi do równania kwadratowego
ar^2 - br + c = 0
którego rozwiązania r czasem (nie zawsze) mogą przybliżać okolice dzielnika. W kolejnej iteracji przedstawiamy n w systemie o podstawie p+r albo p-r.
Kryterium: suma naprzemienna cyfr jest dodatnia, dzielnik może być mniejszy niż aktualna podstawa. Suma naprzemienna cyfr jest ujemna, dzielnik jest przy podstawie większej. Suma naprzemienna równa 0, dzielnik jest dokładnie o 1 mniejszy niż aktualna podstawa liczby.
Lecz zmian znaku sumy naprzemiennej cyfr liczby n zapisanej w różnych systemach pozycyjnych jako trójcyfrowa jest wiele. Wskazówką powodzenia może być fakt, że bardzo blisko dzielnika te zmiany są rzadkie. Tuż przy samym dzielniku równanie sugeruje daleki skok. Najdalsze propozycje są najbliżej dzielnika. Wtedy lepiej zastosować metodę siłową.
20 grudnia 2013
28 listopada 2013
Pierwiastek kwadratowy, przyspieszona wersja
Algorytm we wcześniejszym poście wymagał znalezienia pierwiastka kwadratowego z liczby, liczonego dwukrotnie podczas inicjacji.
Publikowałem już sposób na liczenie pierwiastka rok temu, w listopadzie 2012 roku, lecz teraz znalazłem nieco szybszą technikę.
Ustalmy liczbę n oraz drugą liczbę p położoną między pierwiastkiem sześciennym oraz kwadratowym z n. Nic więcej o liczbie p nie potrzeba wiedzieć, jest ona pierwszym przybliżeniem pierwiastka. Resztę załatwią konwersje, a być może także wzór skróconego mnożenia.
Dzielimy n z resztą przez p (dwukrotnie), by uzyskać postać
n = (a*p+b)*p+c (1)
Jest to zapis liczby trójcyfrowej.
Możemy sprawdzić, czy a, b, c stanowią współczynniki trójmianu kwadratowego d^2+2de+e^2+f, jeśli tak, pierwiastek można przybliżyć przez
(dp+e)^2+f.
Jeśli jednak b>2de, też możemy posłużyć się wzorem, do którego należy dodać poprawkę (b-2de)*p. Sposób też prowadzi do celu, lecz jest znacznie dłuższy. Zbieżność do pierwiastka jest wolniejsza.
Zatem pierwsza faza, sprowadzamy a do postaci podzielnej przez 4, b do liczby parzystej za pomocą formuł:
a-1 AND b+=p; b-1 AND c+=p. (2)
Np. mając (a,b,c) = (5, 4, 34) dla p=125 uzyskujemy z pomocą formuł (2)
(a,b,c) = (5-1, 4+125, 34) = (4, 129, 34) = (4, 129-1, 34+125) = (4, 128, 159)
Korzystamy z konwersji podwajania systemu
(a,b,c)_p = (a/4, b/2, c)_(2p) (3)
czyli w przykładzie mamy (a,b,c) = (1, 64, 159)_250 .
Kiedy już a<4, oraz nie dopasujemy wzoru skróconego mnożenia, dążymy do zmniejszenia b. W odróżnieniu od zeszłorocznego algorytmu najlepiej to zrobić za pomocą konwersji o floor(a*b/3). Konwersja ta zmniejsza b o co najmniej połowę, może także zmniejszyć a. Branie połowy (a*b/2) może doprowadzić do zmniejszenie b do wartości ujemnej, zaś chcemy zachować a=1 jako niezmiennik, oraz wartość liczby trójcyfrowej powinna pozostać większa niż pierwiastek kwadratowy.
Przykład dla n=8934053
Przyjmijmy p=1000, wtedy (a,b,c) = (8, 934, 53)_1000
Pierwsza konwersja (3) sprowadza do (a,b,c) = (2, 467, 53)_2000
konwersja o floor(a*b/3) = 311
2 467 53
2 467-2*311
2-1 2156 53 // dodaję do b (2000+311) = 2311
1 2156-311 53-2156*311
1 1845-291 -670463 // dodaję do c 291*p
1 1554 2038
oraz (a,b,c) = (1, 1554, 2038)_2311
konwersja o floor(a*b/3) = 518, nowa podstawa p = 2311+518 = 2829
1 1554 2038
1 1554-518
1 1036 2038
1 518 2038-1036*518 // dodajemy do c 189*p
1 518-189 -534610
1 329 71
konwersja o floor(a*b/3) = 109, nowa podstawa p = 2829+109 = 2938
1 329 71
1 329-109
1 220 71
1 220-109 71-220*109
1 111-9 -23909 // dodajamy do c 9*p
1 102 2533
Sprawdzamy wzór skróconego mnożenia (1+51)^2 = 1+102+2601 oraz porównujemy. Niezgodność na pozycji c, jest za dużo o 68. Ale jest dobrze.
Zatem (a,b,c)_p = (1*p+51)^2-68 = (2938+51)^2-68 = 2989^2-68.
Zmniejszamy kwadrat
2989^2-68 = 2988^2+5976+1-68 = 2988^2+5909.
Znaleziony został pierwiastek całkowitoliczbowy 2988 oraz reszty: 5909 z niedomiarem, 68 z nadmiarem.
Dokładniejszą wartość pierwiastka mozna policzyć za pomocą interpolacji 2988 + 5909/(5909+68) = 2988,988623; dokładniejsza wartość 2988,9886249365352845723814
Błąd dopiero na szóstej pozycji po przecinku.
Publikowałem już sposób na liczenie pierwiastka rok temu, w listopadzie 2012 roku, lecz teraz znalazłem nieco szybszą technikę.
Ustalmy liczbę n oraz drugą liczbę p położoną między pierwiastkiem sześciennym oraz kwadratowym z n. Nic więcej o liczbie p nie potrzeba wiedzieć, jest ona pierwszym przybliżeniem pierwiastka. Resztę załatwią konwersje, a być może także wzór skróconego mnożenia.
Dzielimy n z resztą przez p (dwukrotnie), by uzyskać postać
n = (a*p+b)*p+c (1)
Jest to zapis liczby trójcyfrowej.
Możemy sprawdzić, czy a, b, c stanowią współczynniki trójmianu kwadratowego d^2+2de+e^2+f, jeśli tak, pierwiastek można przybliżyć przez
(dp+e)^2+f.
Jeśli jednak b>2de, też możemy posłużyć się wzorem, do którego należy dodać poprawkę (b-2de)*p. Sposób też prowadzi do celu, lecz jest znacznie dłuższy. Zbieżność do pierwiastka jest wolniejsza.
Zatem pierwsza faza, sprowadzamy a do postaci podzielnej przez 4, b do liczby parzystej za pomocą formuł:
a-1 AND b+=p; b-1 AND c+=p. (2)
Np. mając (a,b,c) = (5, 4, 34) dla p=125 uzyskujemy z pomocą formuł (2)
(a,b,c) = (5-1, 4+125, 34) = (4, 129, 34) = (4, 129-1, 34+125) = (4, 128, 159)
Korzystamy z konwersji podwajania systemu
(a,b,c)_p = (a/4, b/2, c)_(2p) (3)
czyli w przykładzie mamy (a,b,c) = (1, 64, 159)_250 .
Kiedy już a<4, oraz nie dopasujemy wzoru skróconego mnożenia, dążymy do zmniejszenia b. W odróżnieniu od zeszłorocznego algorytmu najlepiej to zrobić za pomocą konwersji o floor(a*b/3). Konwersja ta zmniejsza b o co najmniej połowę, może także zmniejszyć a. Branie połowy (a*b/2) może doprowadzić do zmniejszenie b do wartości ujemnej, zaś chcemy zachować a=1 jako niezmiennik, oraz wartość liczby trójcyfrowej powinna pozostać większa niż pierwiastek kwadratowy.
Przykład dla n=8934053
Przyjmijmy p=1000, wtedy (a,b,c) = (8, 934, 53)_1000
Pierwsza konwersja (3) sprowadza do (a,b,c) = (2, 467, 53)_2000
konwersja o floor(a*b/3) = 311
2 467 53
2 467-2*311
2-1 2156 53 // dodaję do b (2000+311) = 2311
1 2156-311 53-2156*311
1 1845-291 -670463 // dodaję do c 291*p
1 1554 2038
oraz (a,b,c) = (1, 1554, 2038)_2311
konwersja o floor(a*b/3) = 518, nowa podstawa p = 2311+518 = 2829
1 1554 2038
1 1554-518
1 1036 2038
1 518 2038-1036*518 // dodajemy do c 189*p
1 518-189 -534610
1 329 71
konwersja o floor(a*b/3) = 109, nowa podstawa p = 2829+109 = 2938
1 329 71
1 329-109
1 220 71
1 220-109 71-220*109
1 111-9 -23909 // dodajamy do c 9*p
1 102 2533
Sprawdzamy wzór skróconego mnożenia (1+51)^2 = 1+102+2601 oraz porównujemy. Niezgodność na pozycji c, jest za dużo o 68. Ale jest dobrze.
Zatem (a,b,c)_p = (1*p+51)^2-68 = (2938+51)^2-68 = 2989^2-68.
Zmniejszamy kwadrat
2989^2-68 = 2988^2+5976+1-68 = 2988^2+5909.
Znaleziony został pierwiastek całkowitoliczbowy 2988 oraz reszty: 5909 z niedomiarem, 68 z nadmiarem.
Dokładniejszą wartość pierwiastka mozna policzyć za pomocą interpolacji 2988 + 5909/(5909+68) = 2988,988623; dokładniejsza wartość 2988,9886249365352845723814
Błąd dopiero na szóstej pozycji po przecinku.
25 listopada 2013
Różnica kwadratów w faktoryzacji
Przedstawienie liczby rozkładanej n w postaci dwu kwadratów leży u podstaw jednych z najszybszych metod faktoryzacji. Są to np. metoda sita kwadratowego, sita ciał liczbowych, a nawet ułamków łańcuchowych.
Wykorzystują one kongruencję Legendre'a
x*x = y*y (mod n) ,
gdzie wartości są dodatnie, x<y<n oraz x+y nie sumuje się do n.
Sposób stosowania tej kongruencji znajduje się w "Teorii liczb w informatyce" Song Y. Yana.
Teraz jednak chcę przedstawić inny algorytm, w którym wykorzystuję jedną z najbardziej prostych konwersji systemów niedziesiątkowych, a zarazem wzór skróconego mnożenia:
p^2 + 2*p + 1 = (p+1)^2 + 0*p + 0 = (p+1)^2 (1) .
Kiedy przedstawię liczbę n w postaci:
n = a*a - b*b - c , (2)
gdzie a jest najmniejsze możliwe a^2 < n < (a+1)^2, b jest największe możliwe. Wartość b^2+c można uzyskać wieloma sposobami, ale nie udało mi się jej zmniejszyć.
Przyrost między kolejnymi kwadratami p^2 oraz (p+1)^2 jest równy 2*p+1.
Dodając to do wyrażeń a^2, b^2+c ich różnica pozostanie równa n.
Po zwiększeniu c można dalej zwiększać b, co prowadzi do następującego algorytmu, który zatrzymuje się gdy c=0. Wtedy dzielnikami są:
n = a^2 - b^2 = (a+b)*(a-b) .
Drugim wyjściem jest zwiększenie b powyżej pierwiastka z n, co świadczy o braku dzielników.
inicjacja: a=deck(sqrt(n)); c=a^2-n; b=floor(sqrt(c)); c=c-b^2;
k=a; // kopia dla ograniczenia zakresu
while ( b<k ) {
a++;
c+=2*a+1; // dodanie niedomiaru od najbliższego większego kwadratu
while( 1 ) {
d = 2*b+1; // by nie liczyć dwukrotnie niedomiaru do kwadratu b
if( c<d ) break;
else { c-=d; b++; } // b wzrasta do kolejnego kwadratu kosztem c
if( 0==c ) return text("dzielniki: ", a-b, a+b);
}
}
return text("n jest pierwsza");
Funkcja text() wypisuje swoje argumenty na wyjściu. Pierwiastek z n powoli wzrasta po kolejnych kwadratach. Zaś b^2+c przyrastajace o tę samą wartość modyfikuje się do największego możliwego dla tej wartości kwadratu.
Ze wzrostem b spada liczność obliczeń w pętli wewnętrznej.
Złożoność.
Jeśli przyjmiemy, że pętla wewnętrzna wykona się dokładnie raz, pogorszymy złożoność. Jest to spowodowane tym, że przy każdym, nawet fikcyjnym zwiększeniu a bedziemy zwiększać b, zaś istotniejsze są zwiększenia b. Wyrażenia d=2*b+1 wykonują się w czasie liniowym, pozostałe przekształcenia to czas praktycznie stały. Pętla zewnętrzna wykona się nie więcej niż pierwiastek z n razy, zwiększając b o 1 w każdym kroku. Praktycznie mniej, gdyż b może nie być bardzo małe. Dodatkowo, w pierszych iteracjach zwiększanie b odbywa się bardzo szybko, zwalniając w czasie przebiegu do 1-2 zwiększeń b na jedno zwiększenie a. Podsumowując, złożoność tego algorytmu nie przekracza
O( n sqrt(n) ) dla pamięci,
max{ O( sqrt(n) ), O( sqrt ) } dla operacji.
Dla 18703 mamy wartości początkowe a=137, b^2+c = 66, zatem b=8, c=2.
Jest to złośliwy przypadek, gdyż dzielniki są dosyć odległe, zostaną znalezione dopiero gdy b=129, czyli w 51 na 56 iteracji. Podczas każdej iteracji pętli zewnętrznej wartość b zwiększała się od 1 do 8, najczęściej 2 lub 3.
W następnym kroku mamy:
a=138, jest dodane 2*137+1 = 275;
od c = 2+275 = 277 odejmujemy 2*8+1= 17 (wtedy b=9)
2*9+1 = 19 (oraz b=10)
2*10+1 = 21 (oraz b=11)
jeszcze parę razy i otrzymamy c=17 przy b=18,
Mamy do czynienia z wyrażeniem 138^2 - (18^2+17) = 19044 - 341 = 18703.
W kolejnej iteracji dodamy 2*138+1 = 277. Będziemy zabierać 37, 39, itd. dopóki b=24.
Wartości dodawane i zabierane tworzą proste ciągi arytmetyczne.
Algorytm dla liczb postaci n = p*q, gdzie p i q są pierwsze ma dodatkową własność. Wartość 2*b wskazuje minimalną odległość między tymi dzielnikami.
Wykorzystują one kongruencję Legendre'a
x*x = y*y (mod n) ,
gdzie wartości są dodatnie, x<y<n oraz x+y nie sumuje się do n.
Sposób stosowania tej kongruencji znajduje się w "Teorii liczb w informatyce" Song Y. Yana.
Teraz jednak chcę przedstawić inny algorytm, w którym wykorzystuję jedną z najbardziej prostych konwersji systemów niedziesiątkowych, a zarazem wzór skróconego mnożenia:
p^2 + 2*p + 1 = (p+1)^2 + 0*p + 0 = (p+1)^2 (1) .
Kiedy przedstawię liczbę n w postaci:
n = a*a - b*b - c , (2)
gdzie a jest najmniejsze możliwe a^2 < n < (a+1)^2, b jest największe możliwe. Wartość b^2+c można uzyskać wieloma sposobami, ale nie udało mi się jej zmniejszyć.
Przyrost między kolejnymi kwadratami p^2 oraz (p+1)^2 jest równy 2*p+1.
Dodając to do wyrażeń a^2, b^2+c ich różnica pozostanie równa n.
Po zwiększeniu c można dalej zwiększać b, co prowadzi do następującego algorytmu, który zatrzymuje się gdy c=0. Wtedy dzielnikami są:
n = a^2 - b^2 = (a+b)*(a-b) .
Drugim wyjściem jest zwiększenie b powyżej pierwiastka z n, co świadczy o braku dzielników.
inicjacja: a=deck(sqrt(n)); c=a^2-n; b=floor(sqrt(c)); c=c-b^2;
k=a; // kopia dla ograniczenia zakresu
while ( b<k ) {
a++;
c+=2*a+1; // dodanie niedomiaru od najbliższego większego kwadratu
while( 1 ) {
d = 2*b+1; // by nie liczyć dwukrotnie niedomiaru do kwadratu b
if( c<d ) break;
else { c-=d; b++; } // b wzrasta do kolejnego kwadratu kosztem c
if( 0==c ) return text("dzielniki: ", a-b, a+b);
}
}
return text("n jest pierwsza");
Funkcja text() wypisuje swoje argumenty na wyjściu. Pierwiastek z n powoli wzrasta po kolejnych kwadratach. Zaś b^2+c przyrastajace o tę samą wartość modyfikuje się do największego możliwego dla tej wartości kwadratu.
Ze wzrostem b spada liczność obliczeń w pętli wewnętrznej.
Złożoność.
Jeśli przyjmiemy, że pętla wewnętrzna wykona się dokładnie raz, pogorszymy złożoność. Jest to spowodowane tym, że przy każdym, nawet fikcyjnym zwiększeniu a bedziemy zwiększać b, zaś istotniejsze są zwiększenia b. Wyrażenia d=2*b+1 wykonują się w czasie liniowym, pozostałe przekształcenia to czas praktycznie stały. Pętla zewnętrzna wykona się nie więcej niż pierwiastek z n razy, zwiększając b o 1 w każdym kroku. Praktycznie mniej, gdyż b może nie być bardzo małe. Dodatkowo, w pierszych iteracjach zwiększanie b odbywa się bardzo szybko, zwalniając w czasie przebiegu do 1-2 zwiększeń b na jedno zwiększenie a. Podsumowując, złożoność tego algorytmu nie przekracza
O( n sqrt(n) ) dla pamięci,
max{ O( sqrt(n) ), O( sqrt ) } dla operacji.
Dla 18703 mamy wartości początkowe a=137, b^2+c = 66, zatem b=8, c=2.
Jest to złośliwy przypadek, gdyż dzielniki są dosyć odległe, zostaną znalezione dopiero gdy b=129, czyli w 51 na 56 iteracji. Podczas każdej iteracji pętli zewnętrznej wartość b zwiększała się od 1 do 8, najczęściej 2 lub 3.
W następnym kroku mamy:
a=138, jest dodane 2*137+1 = 275;
od c = 2+275 = 277 odejmujemy 2*8+1= 17 (wtedy b=9)
2*9+1 = 19 (oraz b=10)
2*10+1 = 21 (oraz b=11)
jeszcze parę razy i otrzymamy c=17 przy b=18,
Mamy do czynienia z wyrażeniem 138^2 - (18^2+17) = 19044 - 341 = 18703.
W kolejnej iteracji dodamy 2*138+1 = 277. Będziemy zabierać 37, 39, itd. dopóki b=24.
Wartości dodawane i zabierane tworzą proste ciągi arytmetyczne.
Algorytm dla liczb postaci n = p*q, gdzie p i q są pierwsze ma dodatkową własność. Wartość 2*b wskazuje minimalną odległość między tymi dzielnikami.
21 listopada 2013
Wariacja faktoryzacja przez proste dzielenie
Rozkład liczby n na czynniki za pomocą prostego dzielenia polega na kolejnym wstawianiu do ilorazu n/b kolejnych, najlepiej pierwszych, ewentualnie nieparzystych, liczb b. Reszta równa 0 oznacza, że b jest dzielnikiem.
Dla dużych liczb n b też jest duże, chociaż mniejsze niż pierwiastek z n. Wtedy dzielenie jest uciążliwe.
Istnieje sposób, by zmniejszyć argumenty dzielenia, by nie dzielić n/b, lecz jakieś c/b, gdzie c<n. A nawet dwa sposoby.
Pierwszy korzysta z przedstawienia liczby n jako 'liczby dwycyfrowej' jakiegoś systemu liczenia o dużej podstawie, oraz konwersje pozwalają na znajdowanie reszt - niestety, metoda wymaga liczenia ilorazów dla kolejnych liczb naturalnych przy zmniejszaniu dzielnej.
Metoda, którą opiszę dalej, pozwala przeskakiwać liczby parzyste. Zatem tworzymy ilorazy c/b, gdzie c jest nieparzyste, mniejsze od n.
Zapiszmy liczbę n w postaci wyrażenia:
n = a*b + c (1)
będziemy przekształcać wartości a, b, oraz c w taki sposób, by podczas przekształceń a tworzyło ciąg nierosnący, b było ciągiem rosnącym po wartościach nieparzystych, zaś c tworzyło ciąg przedziałami monotoniczny.
Algorytm jest rozgałęziony, a jego ogólny schemat wygląda następująco:
inicjacja: a = floor(n/3), b=3, c=n%3.
pętla dopóki a>b
czy b dzieli c (warunek 0 = c%b)? jeśli tak, b jest dzielnikiem, wyjście;
b = b+2;
rozgałęzienie, jeśli c<2a
c = c + (a%b)*(b-2);
a = b*floor(a/b) - 2*floor(a/b);
w przeciwnym razie (c>2a)
c = c-2a;
koniec pętli
Pierwsza część rozgałęzienia jest całkowito-liczbową operacją a*(b-2)/b, która zwiększa c oraz zmniejsza a. Mamy tu dodatkowo jedno dzielenie oraz mnożenie, co komplikuje algorytm. Można je przekształcić, by dzielić przez b oraz odjąć podwojony iloraz od a.
Z kolei druga część rozgałęzienia jest zwykłym odejmowaniem.
W początkowej fazie algorytmu przeważa pierwsza część, z dodatkowymi działaniami. Wykonują się jednak one na małych wartościach, dużo mniejszych niż pierwiastek z n. Pod koniec algorytmu najgorszym kawałkiem jest mnożenie (a%b)*(b-2), którego wartość może być stosunkowo blisko n. Zwiększa ona c do bardzo dużych wartości, umożliwiając stosowanie drugiej, prostszej odnogi.
Już po sprawdzeniu około 20% przypadków na możliwe wartości b do głosu dochodzi druga część rozgałęzienia, zaś przy 30% praktycznie dominuje.
Cały czas należy sprawdzać podzielność c przez b. Chociaż c (zwłaszcza przy dominacji drugiej odnogi) szybko maleje nawet do wartości bliskich b.
Fragment przykładu numerycznego, liczba 8 934 053 = 1087 * 8219.
inicjacja: a*b+c = 2 978 017 * 3 + 2
b = 3+2 = 5;
pierwsza odnoga: 2 978 017 - 2 = 2 978 015, bo 2 978 017 % 5 = 2,
c = 2 + 2*3 = 8
zmniejszanie a: 2 978 015 * 3 / 5 = 2 978 015 - 2*(595 603) = 1 786 809
nieco dalej:
a = 122 049; b = 73; c = 24 476;
sprawdzamy c%b = 21, nie jest zerem; pierwsza odnoga
b=75; reszta a%b = 24
c = 24 476 + 24*(75-2) = 26 228;
a = (122 049 - 24) - 2*(122 049 - 24)/75 = 118 771
jeszcze dalej:
a = 69 125; b = 127; c = 155 178;
sprawdzamy c%b = 111, druga odnoga
b = 129;
c = 155 178 - 2 * 69 125 = 16 928;
możemy znów sprawdzać c%b = 16 928 % 129
dla b>400 co parę przekształceń mamy kilka iteracji odnogą drugą, zaś blisko końca algorytmu mamy
a = 8685; b = 1027; c = 14 558;
sprawdzamy c%b = 180; pierwsza odnoga
b = 1029;
c = c + 453*1027 = 479 789;
a = (a-453) - 2*(a-453)/1029 = 8232-16 = 8216;
Teraz przez 29 iteracji powtarza się odejmowanie w drugiej odnodze, zanim znów zastosujemy pierwszą.
Zakończenie algorytmu
a*b+c = 8216 * 1087 + 3261, oraz 3261 = 3*1087.
Algorytm nie jest zbyt dobry. Może być nawet gorszy niż zwykłe dzielenia. Może być traktowany jako ciekawostka.
Dla dużych liczb n b też jest duże, chociaż mniejsze niż pierwiastek z n. Wtedy dzielenie jest uciążliwe.
Istnieje sposób, by zmniejszyć argumenty dzielenia, by nie dzielić n/b, lecz jakieś c/b, gdzie c<n. A nawet dwa sposoby.
Pierwszy korzysta z przedstawienia liczby n jako 'liczby dwycyfrowej' jakiegoś systemu liczenia o dużej podstawie, oraz konwersje pozwalają na znajdowanie reszt - niestety, metoda wymaga liczenia ilorazów dla kolejnych liczb naturalnych przy zmniejszaniu dzielnej.
Metoda, którą opiszę dalej, pozwala przeskakiwać liczby parzyste. Zatem tworzymy ilorazy c/b, gdzie c jest nieparzyste, mniejsze od n.
Zapiszmy liczbę n w postaci wyrażenia:
n = a*b + c (1)
będziemy przekształcać wartości a, b, oraz c w taki sposób, by podczas przekształceń a tworzyło ciąg nierosnący, b było ciągiem rosnącym po wartościach nieparzystych, zaś c tworzyło ciąg przedziałami monotoniczny.
Algorytm jest rozgałęziony, a jego ogólny schemat wygląda następująco:
inicjacja: a = floor(n/3), b=3, c=n%3.
pętla dopóki a>b
czy b dzieli c (warunek 0 = c%b)? jeśli tak, b jest dzielnikiem, wyjście;
b = b+2;
rozgałęzienie, jeśli c<2a
c = c + (a%b)*(b-2);
a = b*floor(a/b) - 2*floor(a/b);
w przeciwnym razie (c>2a)
c = c-2a;
koniec pętli
Pierwsza część rozgałęzienia jest całkowito-liczbową operacją a*(b-2)/b, która zwiększa c oraz zmniejsza a. Mamy tu dodatkowo jedno dzielenie oraz mnożenie, co komplikuje algorytm. Można je przekształcić, by dzielić przez b oraz odjąć podwojony iloraz od a.
Z kolei druga część rozgałęzienia jest zwykłym odejmowaniem.
W początkowej fazie algorytmu przeważa pierwsza część, z dodatkowymi działaniami. Wykonują się jednak one na małych wartościach, dużo mniejszych niż pierwiastek z n. Pod koniec algorytmu najgorszym kawałkiem jest mnożenie (a%b)*(b-2), którego wartość może być stosunkowo blisko n. Zwiększa ona c do bardzo dużych wartości, umożliwiając stosowanie drugiej, prostszej odnogi.
Już po sprawdzeniu około 20% przypadków na możliwe wartości b do głosu dochodzi druga część rozgałęzienia, zaś przy 30% praktycznie dominuje.
Cały czas należy sprawdzać podzielność c przez b. Chociaż c (zwłaszcza przy dominacji drugiej odnogi) szybko maleje nawet do wartości bliskich b.
Fragment przykładu numerycznego, liczba 8 934 053 = 1087 * 8219.
inicjacja: a*b+c = 2 978 017 * 3 + 2
b = 3+2 = 5;
pierwsza odnoga: 2 978 017 - 2 = 2 978 015, bo 2 978 017 % 5 = 2,
c = 2 + 2*3 = 8
zmniejszanie a: 2 978 015 * 3 / 5 = 2 978 015 - 2*(595 603) = 1 786 809
nieco dalej:
a = 122 049; b = 73; c = 24 476;
sprawdzamy c%b = 21, nie jest zerem; pierwsza odnoga
b=75; reszta a%b = 24
c = 24 476 + 24*(75-2) = 26 228;
a = (122 049 - 24) - 2*(122 049 - 24)/75 = 118 771
jeszcze dalej:
a = 69 125; b = 127; c = 155 178;
sprawdzamy c%b = 111, druga odnoga
b = 129;
c = 155 178 - 2 * 69 125 = 16 928;
możemy znów sprawdzać c%b = 16 928 % 129
dla b>400 co parę przekształceń mamy kilka iteracji odnogą drugą, zaś blisko końca algorytmu mamy
a = 8685; b = 1027; c = 14 558;
sprawdzamy c%b = 180; pierwsza odnoga
b = 1029;
c = c + 453*1027 = 479 789;
a = (a-453) - 2*(a-453)/1029 = 8232-16 = 8216;
Teraz przez 29 iteracji powtarza się odejmowanie w drugiej odnodze, zanim znów zastosujemy pierwszą.
Zakończenie algorytmu
a*b+c = 8216 * 1087 + 3261, oraz 3261 = 3*1087.
Algorytm nie jest zbyt dobry. Może być nawet gorszy niż zwykłe dzielenia. Może być traktowany jako ciekawostka.
08 listopada 2013
Modyfikacja heurystyki rho Pollarda
Kilkanaście lat temu zetknąłem się z heurystyką rho Pollarda. Polega ona na obliczaniu kolejnej wartości wielomianu W[x] modulo n startujacego z ziarna a. Metoda ładnie i dostępnie opisana w wikipedii, książkach do kryptografii.
Opiszę tu swoje (alternatywne) podejście, bo jak większość moich pomysłów, ciężko jest mi je publikować.
We wszystkich tych źródłach pojawia się informacja, że wielomianem nie może być wielomian x*x-2 ani x*x.
A właśnie drugi z tych wielomianów użyłem, kiedy zacząłem testować heurystykę. Występowała niejednoznaczność, którą można było bardzo łatwo usunąć biorąc kolejną wartość wielomianu jako
f(x) = min( x^2, n-x^2 ) mod n .
Zastosowanie tego wielomianu, nie dość, że zmniejsza o połowę krotność występujacych współczyników, to jeszcze ma kilka dodatkowych własności bardzo przyspieszajacych obliczenia. Można używać indeksów dla szacowania, kiedy i gdzie może pojawić się dzielnik.
Zmienia się także warunek stopu. Nie jest już jeden:
nwd( f(a)^m, n ),
gdzie f(a)^m oznacza m-krotne złożenie f liczone w a,
czy w modyfikacji p-1: nwd( f(a)^m -1, n ) .
Odpada także przygotowanie liczby do obliczeń przez zignorowanie kilku początkowch iteracji, jakie spotkałem w jednej z książek dotyczących algorytmów.
Wybieramy ziarno 1<a<(n-1)/2. Zapamiętujemy jego wartość (dla indeksu).
Obliczamy f(a), kiedy po raz pierwszy uzyskamy wartość większą niż n, wtedy zapamiętujemy dodatkowo b jako mniejszą z reszt z dzielenia przez n: b = min( b%n, (n-b)%n ).
Powtarzamy obliczenia zliczając złożenia, aż uzyskamy wartość a, b, 0 albo 1.
W przypadku uzyskania 0 dzielnikiem jest ziarno a.
Gdy uzyskamy ciąg okresowy z wyrazem a lub b, liczymy dla tego wyrazu:
np. d=nwd( a-1,n).
Gdy 0<d<n/2, d jest dzielnikiem. Przy wartości d=0 jeszcze nic straconego. Liczymy dla pewnego elementu x =f(a)^m:
e = nwd( x, f(x) ) ,
który jest kolejnym kandydatem na dzielnik.
Najciekawszy przypadek jest wtedy, gdy cyklem jest 1. W standardowym podejściu ten przypadek powoduje błąd heurystyki. Ale nie tutaj.
Zaczynamy algorytm z ziarnem 1<a+1<(n-1)/2. I teraz tylko liczby pierwsze mogą uzyskać ciąg złożony z 1. W sprawdzanych przeze mnie przypadkach pojawia się inny ciąg stały złożony z dzielnika n. Cykle jedynek i dzielnika przeplatają się wzajemnie dla kolejnych ziaren.
Gdy nie znajdziemy dzielnika, zwiększamy ziarno o 1, oraz puszczamy algorytm jeszcze raz, pamiętając dodatkowo poprzednie wartości a, b. Teraz istotne jest ich znalezienie w nowo powstającym ciągu.
Cykl dla sąsiedniego ziarna nie musi zawierać zapamiętanych wartości. Jest wtedy rozłączny, co pozwala oszacować, czy wogóle jest miejsce na dzielniki. Czy będą związane z jakąś liczbą pierwszą, czyli gdzie ich szukać.
Kilkanaście lat temu myślałem o pamiętaniu ich w tablicy, by startować z ziarna jeszcze nie występującego.
Ale po co zapamietywać liczność cyklu? Jest on potrzebny dla oczacowania, czy mamy do czynienia z liczbą pierwszą czy złożoną. Cykl długości (n-3)/4 powtarzający się dla dwu sąsiednich (względnie pierwszych) ziaren wskazuje liczbę pierwszą.
Uwaga, ziaren równych 0, 1, n-1 nie należy brać pod uwagę, uzyska się stały ciąg trywialnej wartości.
Przykłady.
Liczba 51.
Ziarno 2, ciąg wpada w cykl [1].
Ziarno 3, ciąg wpada w cykl [18], nwd(18-1, 51) = 17, oraz 51=3*17.
Przeplatane cykle [1] oraz [18].
Liczba 35.
Ziarno 2, ciąg [2, 4, 11, 16,... ] ogon [2,4], cykl (druga z zapamiętanych wartości to 11) [11, 16], liczymy -1+11 = 10 (mod 35), d = nwd( 10, 35 ) = 5 jest dzielnikiem.
Ziarno 5, ciąg [5, 10, 5...] z cyklem [5, 10]. Teraz d=nwd(4,35)=1. Drugi warunek stopu e=nwd(5,10)=5 wskazuje dzielnik 35.
Liczba pierwsza 47.
Dla wiekszości ziaren mamy cykle długości 11 elementów, ziarno 2 ma zapamiętaną wartość 2, ziarno 3 ma cykl rozłączny (nie zawierający 2), oba cykle wybierają 22 elementy z (47-3)/2 = 22 dostępnych. Jest to liczba pierwsza.
Opiszę tu swoje (alternatywne) podejście, bo jak większość moich pomysłów, ciężko jest mi je publikować.
We wszystkich tych źródłach pojawia się informacja, że wielomianem nie może być wielomian x*x-2 ani x*x.
A właśnie drugi z tych wielomianów użyłem, kiedy zacząłem testować heurystykę. Występowała niejednoznaczność, którą można było bardzo łatwo usunąć biorąc kolejną wartość wielomianu jako
f(x) = min( x^2, n-x^2 ) mod n .
Zastosowanie tego wielomianu, nie dość, że zmniejsza o połowę krotność występujacych współczyników, to jeszcze ma kilka dodatkowych własności bardzo przyspieszajacych obliczenia. Można używać indeksów dla szacowania, kiedy i gdzie może pojawić się dzielnik.
Zmienia się także warunek stopu. Nie jest już jeden:
nwd( f(a)^m, n ),
gdzie f(a)^m oznacza m-krotne złożenie f liczone w a,
czy w modyfikacji p-1: nwd( f(a)^m -1, n ) .
Odpada także przygotowanie liczby do obliczeń przez zignorowanie kilku początkowch iteracji, jakie spotkałem w jednej z książek dotyczących algorytmów.
Wybieramy ziarno 1<a<(n-1)/2. Zapamiętujemy jego wartość (dla indeksu).
Obliczamy f(a), kiedy po raz pierwszy uzyskamy wartość większą niż n, wtedy zapamiętujemy dodatkowo b jako mniejszą z reszt z dzielenia przez n: b = min( b%n, (n-b)%n ).
Powtarzamy obliczenia zliczając złożenia, aż uzyskamy wartość a, b, 0 albo 1.
W przypadku uzyskania 0 dzielnikiem jest ziarno a.
Gdy uzyskamy ciąg okresowy z wyrazem a lub b, liczymy dla tego wyrazu:
np. d=nwd( a-1,n).
Gdy 0<d<n/2, d jest dzielnikiem. Przy wartości d=0 jeszcze nic straconego. Liczymy dla pewnego elementu x =f(a)^m:
e = nwd( x, f(x) ) ,
który jest kolejnym kandydatem na dzielnik.
Najciekawszy przypadek jest wtedy, gdy cyklem jest 1. W standardowym podejściu ten przypadek powoduje błąd heurystyki. Ale nie tutaj.
Zaczynamy algorytm z ziarnem 1<a+1<(n-1)/2. I teraz tylko liczby pierwsze mogą uzyskać ciąg złożony z 1. W sprawdzanych przeze mnie przypadkach pojawia się inny ciąg stały złożony z dzielnika n. Cykle jedynek i dzielnika przeplatają się wzajemnie dla kolejnych ziaren.
Gdy nie znajdziemy dzielnika, zwiększamy ziarno o 1, oraz puszczamy algorytm jeszcze raz, pamiętając dodatkowo poprzednie wartości a, b. Teraz istotne jest ich znalezienie w nowo powstającym ciągu.
Cykl dla sąsiedniego ziarna nie musi zawierać zapamiętanych wartości. Jest wtedy rozłączny, co pozwala oszacować, czy wogóle jest miejsce na dzielniki. Czy będą związane z jakąś liczbą pierwszą, czyli gdzie ich szukać.
Kilkanaście lat temu myślałem o pamiętaniu ich w tablicy, by startować z ziarna jeszcze nie występującego.
Ale po co zapamietywać liczność cyklu? Jest on potrzebny dla oczacowania, czy mamy do czynienia z liczbą pierwszą czy złożoną. Cykl długości (n-3)/4 powtarzający się dla dwu sąsiednich (względnie pierwszych) ziaren wskazuje liczbę pierwszą.
Uwaga, ziaren równych 0, 1, n-1 nie należy brać pod uwagę, uzyska się stały ciąg trywialnej wartości.
Przykłady.
Liczba 51.
Ziarno 2, ciąg wpada w cykl [1].
Ziarno 3, ciąg wpada w cykl [18], nwd(18-1, 51) = 17, oraz 51=3*17.
Przeplatane cykle [1] oraz [18].
Liczba 35.
Ziarno 2, ciąg [2, 4, 11, 16,... ] ogon [2,4], cykl (druga z zapamiętanych wartości to 11) [11, 16], liczymy -1+11 = 10 (mod 35), d = nwd( 10, 35 ) = 5 jest dzielnikiem.
Ziarno 5, ciąg [5, 10, 5...] z cyklem [5, 10]. Teraz d=nwd(4,35)=1. Drugi warunek stopu e=nwd(5,10)=5 wskazuje dzielnik 35.
Liczba pierwsza 47.
Dla wiekszości ziaren mamy cykle długości 11 elementów, ziarno 2 ma zapamiętaną wartość 2, ziarno 3 ma cykl rozłączny (nie zawierający 2), oba cykle wybierają 22 elementy z (47-3)/2 = 22 dostępnych. Jest to liczba pierwsza.
04 listopada 2013
Faktoryzacja a system Fibonacciego
Liczby Fibonacciego tworzone są wzorem rekurencyjnym:
F[0]=0, F[1], F[n+2] = F[n+1]+F[n]
W wikipedii opisano system Fibonacciego, w którym w definicji systemu o cyfrach a[i] oraz podstawach p[i]
suma_{-1<i} a[i]*p[i]
przyjęto
p[i] = F[i-2], 1<i, a[i] in {0,1}.
Dodawanie i odejmowanie w tym systemie jest proste, mnożenie i dzielenie koszmarne - lepiej się trzymać z daleka. Dla usunięcia niejednoznaczności wprowadzono formułę łączącą w okienku kolejne trzy cyfry systemu Fibonacciego:
[0,1,1]_F = [1,0,0]_F (1) .
Nie musimy mieć cyfr równych tylko {0,1}. Możemy czasowo pozwolić sobie na dowolne 'cyfry' całkowite, a wtedy mnożenie i dzielenie nieco się ułatwi.
Przykład. Liczba [3,0,0,0,0,0]_F = 3*F[7] = 3*13 = 39.
Własność (1) można wykorzystać do faktoryzacji. Sprawdziłem trzy metody przekształcania liczby n na postać p*F[m].
Idea jest następująca, stosując prawa działające na liczbach Fibonacciego tworzę ciąg 'cyfr' ze zbioru {0,p}, do pozycji F[0] dodaję ewentualną resztę całkowitą r.
1) Liczba nieparzysta p zapisana liczbą Fibonacciego F jako p*F (użyta gramatyka (p|0)*r ), gdzie p*p<n<p*p+2.
zmniejszamy p = p-2.
Przechodzimy rekursywnie modyfikując 3-cyfrowe okienko [p+a,b,c] = [p,a+b,a+c] lub [a,b,c] = [0,a+b,a+c], a<p.
Okienko przesuwamy w kierunku 'cyfr' mniej znaczących. Kiedy się nie da, dla okienka [a,b,c] stosujemy przekształcenie d=2*b+c+r oraz wypełniamy: [a, d%(2*p), e*p+r'], gdzie r' jest nową resztą d%p, e in {0,1} jest dopasowane, by nie zmienić wartości liczby.
Przykład p=21, w ostatnim okienku mamy [21, 34, 4] = [21,0,72] = [21,21,21+9].
Teraz należy 'naprawić' liczbę Fibonacciego, wracamy z okienkiem posługując się uogólnieniem okienka (1) postaci:
[0,p,p] = [p,0,0] (2) .
Wielokrotnie to stosując (czasem z nawrotami dla ostatnich cyfr) uzyskamy nową wartość
(p+2)*F +r = p*F'+r'.
Kontynuując ze zmniejszajacym się nieparzystym p znajdziemy dzielnik lub dojdziemy do standardowej liczby Fibonacciego dla n. Wtedy kończymy.
W każdej iteracji wartość p się zmniejsza, F zwiększa. Za wyjątkiem obróbki końcowej cyfry mamy do czynienia tylko z dodawaniem, odejmowaniem i porównywaniem liczb nie większych niż pierwiastek z n.
2) Jest to algorytm dla liczb F mniejszych niż pierwiastek, wykorzystujący dzielenie. W pierwszym kroku dzielimy p=n/3, uzyskując postać [p,0]+r = p*10_F+r. Sprawdzamy, czy r=p albo r=0, wtedy kończymy mając dzielnik p.
W następnym kroku sprawdzamy, o ile możemy zmniejszyć 'cyfry' p za pomocą dzielenia: p/(F+1) <= x, dla nieparzystego p przyjmujemy za x najbliższą parzystą większą niż iloraz. W następnym kroku odejmujemy p=p-x, co zwiększy nam resztę r o x*F. Podobnie jak w 1), 'naprawiamy' liczbę za pomocą (2).
Porównujemy p z r oraz zmniejszamy p = p-2. Kończymy, gdy p stanie się mniejsze niż liczba Fibonacciego F lub znajdziemy dzielnik p.
Liczba Fibonacciego się zwiększa do długości osiągniętej przy pierwiastku z n, wartości p początkowo bardzo szybko maleją, póżniej słabiej, bo ilorazy dążą do 1, iloraz (całkowity) równy 0 oznacza, że p<F.
3) Jest odwróceniem algorytmu 1). Startujemy z liczby Fibonacciego dla wartości n oraz usiłujemy doprowadzić pierwszą 'cyfrę' okienka do wartości p albo do 0. Należy odpowiednio reagować, gdyż wartości w okienku mają silną tendencję do tworzenia ciągów rozbieżnych. Teraz wygodniej jest czasami nadawać 'cyfrom' wartości ujemne i czasami zciągać cyfrę z już ustalonej pozycji.
Zwiększamy p o 2. Są kłopoty z ustaleniem kryterium stopu, gdyż sprawdzanie, czy p<F za każdym zwiększeniem p jest bezcelowe. Z kolei konieczność 'naprawy' liczby po dotarciu okienka do cyfry najmniej znaczącej zachodzi stosunkowo rzadko.
Przykładowy fragment działania na okienku dla p=67: [65,0,65]=[67,-2,63]; przesunięcie [-2,63,0]=[0,61,-2]; przesunięcie [61,-2,0]=[0,59,61] itd.
Następuje zciąganie liczby Fibonacciego do krótszego ciągu. Operujemy na wzrastajacych wartościach dodawaniem, odejmowaniem i porównywaniem. Osiągniemy pierwiastek z n tylko dla liczb pierwszych.
Liczba w systemie Fibonacciego jest długa, dochodzi do 5 log n, ale przekształcenia i pętle (może z wyjątkiem napraw) są liniowe. Wszystkie przedstawione algorytmy mają porównywalną krotność przypadków do sprawdzenia. Bo albo p się zmienia, albo F.
Są to kolejni kandydaci na wielomianowe algorytmy faktoryzacji, gdyż znalezienie wartości cyfry liczby Fibonacciego można wykonać w czasie logarytmicznym.
Algorytm 1) sprawdza najpierw duże dzielniki, pozostałe najpierw małe.
F[0]=0, F[1], F[n+2] = F[n+1]+F[n]
W wikipedii opisano system Fibonacciego, w którym w definicji systemu o cyfrach a[i] oraz podstawach p[i]
suma_{-1<i} a[i]*p[i]
przyjęto
p[i] = F[i-2], 1<i, a[i] in {0,1}.
Dodawanie i odejmowanie w tym systemie jest proste, mnożenie i dzielenie koszmarne - lepiej się trzymać z daleka. Dla usunięcia niejednoznaczności wprowadzono formułę łączącą w okienku kolejne trzy cyfry systemu Fibonacciego:
[0,1,1]_F = [1,0,0]_F (1) .
Nie musimy mieć cyfr równych tylko {0,1}. Możemy czasowo pozwolić sobie na dowolne 'cyfry' całkowite, a wtedy mnożenie i dzielenie nieco się ułatwi.
Przykład. Liczba [3,0,0,0,0,0]_F = 3*F[7] = 3*13 = 39.
Własność (1) można wykorzystać do faktoryzacji. Sprawdziłem trzy metody przekształcania liczby n na postać p*F[m].
Idea jest następująca, stosując prawa działające na liczbach Fibonacciego tworzę ciąg 'cyfr' ze zbioru {0,p}, do pozycji F[0] dodaję ewentualną resztę całkowitą r.
1) Liczba nieparzysta p zapisana liczbą Fibonacciego F jako p*F (użyta gramatyka (p|0)*r ), gdzie p*p<n<p*p+2.
zmniejszamy p = p-2.
Przechodzimy rekursywnie modyfikując 3-cyfrowe okienko [p+a,b,c] = [p,a+b,a+c] lub [a,b,c] = [0,a+b,a+c], a<p.
Okienko przesuwamy w kierunku 'cyfr' mniej znaczących. Kiedy się nie da, dla okienka [a,b,c] stosujemy przekształcenie d=2*b+c+r oraz wypełniamy: [a, d%(2*p), e*p+r'], gdzie r' jest nową resztą d%p, e in {0,1} jest dopasowane, by nie zmienić wartości liczby.
Przykład p=21, w ostatnim okienku mamy [21, 34, 4] = [21,0,72] = [21,21,21+9].
Teraz należy 'naprawić' liczbę Fibonacciego, wracamy z okienkiem posługując się uogólnieniem okienka (1) postaci:
[0,p,p] = [p,0,0] (2) .
Wielokrotnie to stosując (czasem z nawrotami dla ostatnich cyfr) uzyskamy nową wartość
(p+2)*F +r = p*F'+r'.
Kontynuując ze zmniejszajacym się nieparzystym p znajdziemy dzielnik lub dojdziemy do standardowej liczby Fibonacciego dla n. Wtedy kończymy.
W każdej iteracji wartość p się zmniejsza, F zwiększa. Za wyjątkiem obróbki końcowej cyfry mamy do czynienia tylko z dodawaniem, odejmowaniem i porównywaniem liczb nie większych niż pierwiastek z n.
2) Jest to algorytm dla liczb F mniejszych niż pierwiastek, wykorzystujący dzielenie. W pierwszym kroku dzielimy p=n/3, uzyskując postać [p,0]+r = p*10_F+r. Sprawdzamy, czy r=p albo r=0, wtedy kończymy mając dzielnik p.
W następnym kroku sprawdzamy, o ile możemy zmniejszyć 'cyfry' p za pomocą dzielenia: p/(F+1) <= x, dla nieparzystego p przyjmujemy za x najbliższą parzystą większą niż iloraz. W następnym kroku odejmujemy p=p-x, co zwiększy nam resztę r o x*F. Podobnie jak w 1), 'naprawiamy' liczbę za pomocą (2).
Porównujemy p z r oraz zmniejszamy p = p-2. Kończymy, gdy p stanie się mniejsze niż liczba Fibonacciego F lub znajdziemy dzielnik p.
Liczba Fibonacciego się zwiększa do długości osiągniętej przy pierwiastku z n, wartości p początkowo bardzo szybko maleją, póżniej słabiej, bo ilorazy dążą do 1, iloraz (całkowity) równy 0 oznacza, że p<F.
3) Jest odwróceniem algorytmu 1). Startujemy z liczby Fibonacciego dla wartości n oraz usiłujemy doprowadzić pierwszą 'cyfrę' okienka do wartości p albo do 0. Należy odpowiednio reagować, gdyż wartości w okienku mają silną tendencję do tworzenia ciągów rozbieżnych. Teraz wygodniej jest czasami nadawać 'cyfrom' wartości ujemne i czasami zciągać cyfrę z już ustalonej pozycji.
Zwiększamy p o 2. Są kłopoty z ustaleniem kryterium stopu, gdyż sprawdzanie, czy p<F za każdym zwiększeniem p jest bezcelowe. Z kolei konieczność 'naprawy' liczby po dotarciu okienka do cyfry najmniej znaczącej zachodzi stosunkowo rzadko.
Przykładowy fragment działania na okienku dla p=67: [65,0,65]=[67,-2,63]; przesunięcie [-2,63,0]=[0,61,-2]; przesunięcie [61,-2,0]=[0,59,61] itd.
Następuje zciąganie liczby Fibonacciego do krótszego ciągu. Operujemy na wzrastajacych wartościach dodawaniem, odejmowaniem i porównywaniem. Osiągniemy pierwiastek z n tylko dla liczb pierwszych.
Liczba w systemie Fibonacciego jest długa, dochodzi do 5 log n, ale przekształcenia i pętle (może z wyjątkiem napraw) są liniowe. Wszystkie przedstawione algorytmy mają porównywalną krotność przypadków do sprawdzenia. Bo albo p się zmienia, albo F.
Są to kolejni kandydaci na wielomianowe algorytmy faktoryzacji, gdyż znalezienie wartości cyfry liczby Fibonacciego można wykonać w czasie logarytmicznym.
Algorytm 1) sprawdza najpierw duże dzielniki, pozostałe najpierw małe.
14 października 2013
Faktoryzacja przez iloczyn dzielników
Ten 'algorytm', a właściwie przepis powstał, kiedy badałem, czy można zmniejszyć krotność rozpatrywanych przypadków.
Pierwsze algorytmy mają po sqrt(n)/2 przypadków, publikowane w tym roku (2013) zmniejszyły krotność do 2^(i+1), gdzie 2^(2i)<n<2^(2i+2).
Opisywany dalej przepis może mieć mniej przypadków, ale nie można podać ile. Bardzo silnie zależy ona od podstawy systemu.
Zauważmy, że mając liczby a i b, zwiększając o 1 cyfrę a na pozycji i-tej, iloczyn zmienia się o wartość b*10^i. Podobnie zmniejszając o 1 cyfrę. Pozwala to wytworzyć ciąg zbieżny do wartości złożonej n. Kiedy iloczyn jest za duży, zmniejszamy jeden z dzielników, kiedy za mały, zwiększamy.
Ale należy to robić sprytnie. W systemie dziesiątkowym liczba złożona n mająca duże dzielniki pierwsze ma cyfrę jedności 1, 3, 7, 9. Wykluczone zostały dzielniki 2 oraz 5. Zatem mozemy dobrać a=11, b = floor(n/11), co spowoduje, że iloczyn a*b>n.
Dobierzmy i=1. Zmniejszamy b tak długo, dopóki wartość ilorazu nie będzie mniejsza niż n. Ale zatrzymajmy się, gdy cyfry najmniej znaczące iloczynu oraz n modulo 10^i będą równe. Wtedy zwiększamy i=i+1, oraz kontynujemy, biorąc pod uwagę dwie cyfry n. Tym razem zmniejszamy cyfrę dziesiątek liczby b. Możemy kontynuować z setkami itd.
Kiedy iloczyn stanie się odpowiednio mały: a*b<n, zwiększamy a do najbliższej wartości z cyfrą 1, 3, 7, 9 na pozycji jednostek.Powoduje to zwiększenie iloczynu.
Uzyskujemy ciąg wartości oscylujacy wokół wartości n. Kiedy wreszcie iloczyn stanie się równy n, nasze wartości a i b są dzielnikami.
Testy przekonały mnie, że sposób jest niezależny od systemu, w którym liczymy. Mamy do czynienia tylko z odejmowaniem i dodawaniem liczb w innym systemie, ewentualnie mnożonych przez małe stałe (nie większe niż podstawa systemu).
Zaś w różnych systemach liczność podejrzanych wartości a jest inna. Podejrzewam, że najgorzej jest w systemach o podstawach będących liczbami pierwszymi (każda wartość prócz 0 może wystąpić), najłatwiej w systemach będących liczbami złożonymi z różnych liczb pierwszych.
Przykładowo, system o podstawie 210 = 2*3*5*7 zawiera w cyfrze najmniej znaczącej informację o podzielności przez dowolną liczbę pierwszą pierwszej dziesiątki. Zaś poszukiwaną postacią liczby a jest liczba pierwsza (dziesiątkowa) mniejsza niż 210, albo jedna z pięciu złożonych wartości: 121, 143, 169, 187, 209.
Rozłożymy liczbę 6767 = 67*101. W systemie o podstawie 210 ma ona postać 32'47_{210}, przyjmijmy (dla wygody) a=10, b=3'46_{210}.
Iloczyn a*b = 10*3'46_{210} = 32'40_{210} = 6760<n, gdyż 10*46 = 460 = 2*210+40
Zwiększając a do 11, otrzymujemy:
a*b = 11*3'46_{210} = 32'40_{210}+ 3'46_{210} = 35'86_{210} = 7436.
Teraz zmniejszamy b, odejmując 11 z 'cyfry' dziesiątek, dopóki 86-k*11<47, albo 'cyfra jedności' b będzie równa 47. Pierwszy warunek będzie spełniony szybciej, po zmianie b o 61_{210}, zatem liczymy
a*b = 11*3'42_{210} = 35'86_{210} - 11*61_{210} = 35'86_{210} - 3'41_{210} = 32'45_{210} = 6765
Znów zwiększamy a do najbliższej liczby pierwszej 13, dodając do iloczynu 2*3'42_{210}, itd.
Ostatecznie uzyskamy a=67, b=101, skacząc po samych liczbach pierwszych dla a.
Pierwsze algorytmy mają po sqrt(n)/2 przypadków, publikowane w tym roku (2013) zmniejszyły krotność do 2^(i+1), gdzie 2^(2i)<n<2^(2i+2).
Opisywany dalej przepis może mieć mniej przypadków, ale nie można podać ile. Bardzo silnie zależy ona od podstawy systemu.
Zauważmy, że mając liczby a i b, zwiększając o 1 cyfrę a na pozycji i-tej, iloczyn zmienia się o wartość b*10^i. Podobnie zmniejszając o 1 cyfrę. Pozwala to wytworzyć ciąg zbieżny do wartości złożonej n. Kiedy iloczyn jest za duży, zmniejszamy jeden z dzielników, kiedy za mały, zwiększamy.
Ale należy to robić sprytnie. W systemie dziesiątkowym liczba złożona n mająca duże dzielniki pierwsze ma cyfrę jedności 1, 3, 7, 9. Wykluczone zostały dzielniki 2 oraz 5. Zatem mozemy dobrać a=11, b = floor(n/11), co spowoduje, że iloczyn a*b>n.
Dobierzmy i=1. Zmniejszamy b tak długo, dopóki wartość ilorazu nie będzie mniejsza niż n. Ale zatrzymajmy się, gdy cyfry najmniej znaczące iloczynu oraz n modulo 10^i będą równe. Wtedy zwiększamy i=i+1, oraz kontynujemy, biorąc pod uwagę dwie cyfry n. Tym razem zmniejszamy cyfrę dziesiątek liczby b. Możemy kontynuować z setkami itd.
Kiedy iloczyn stanie się odpowiednio mały: a*b<n, zwiększamy a do najbliższej wartości z cyfrą 1, 3, 7, 9 na pozycji jednostek.Powoduje to zwiększenie iloczynu.
Uzyskujemy ciąg wartości oscylujacy wokół wartości n. Kiedy wreszcie iloczyn stanie się równy n, nasze wartości a i b są dzielnikami.
Testy przekonały mnie, że sposób jest niezależny od systemu, w którym liczymy. Mamy do czynienia tylko z odejmowaniem i dodawaniem liczb w innym systemie, ewentualnie mnożonych przez małe stałe (nie większe niż podstawa systemu).
Zaś w różnych systemach liczność podejrzanych wartości a jest inna. Podejrzewam, że najgorzej jest w systemach o podstawach będących liczbami pierwszymi (każda wartość prócz 0 może wystąpić), najłatwiej w systemach będących liczbami złożonymi z różnych liczb pierwszych.
Przykładowo, system o podstawie 210 = 2*3*5*7 zawiera w cyfrze najmniej znaczącej informację o podzielności przez dowolną liczbę pierwszą pierwszej dziesiątki. Zaś poszukiwaną postacią liczby a jest liczba pierwsza (dziesiątkowa) mniejsza niż 210, albo jedna z pięciu złożonych wartości: 121, 143, 169, 187, 209.
Rozłożymy liczbę 6767 = 67*101. W systemie o podstawie 210 ma ona postać 32'47_{210}, przyjmijmy (dla wygody) a=10, b=3'46_{210}.
Iloczyn a*b = 10*3'46_{210} = 32'40_{210} = 6760<n, gdyż 10*46 = 460 = 2*210+40
Zwiększając a do 11, otrzymujemy:
a*b = 11*3'46_{210} = 32'40_{210}+ 3'46_{210} = 35'86_{210} = 7436.
Teraz zmniejszamy b, odejmując 11 z 'cyfry' dziesiątek, dopóki 86-k*11<47, albo 'cyfra jedności' b będzie równa 47. Pierwszy warunek będzie spełniony szybciej, po zmianie b o 61_{210}, zatem liczymy
a*b = 11*3'42_{210} = 35'86_{210} - 11*61_{210} = 35'86_{210} - 3'41_{210} = 32'45_{210} = 6765
Znów zwiększamy a do najbliższej liczby pierwszej 13, dodając do iloczynu 2*3'42_{210}, itd.
Ostatecznie uzyskamy a=67, b=101, skacząc po samych liczbach pierwszych dla a.
21 września 2013
Wzory na budowę dzielników, uzupełnienie w dziesiątkowym
Rozważany w ostatnim poście algorytm okazał się niezależny od podstawy systemu. Sprawdził się w binarnym, sprawdził się też w dziesiątkowym.
Funkcja jest postaci
(1) k(x,n,p,a,b) = -(p+b)/p + n/(p*(px+p+a)
dla potęgi podstawy systemu p, liczby rozkładanej n, odpowiednio przygotowanych wartości a, b spełniajacych warunek
(2) a*b%p = n%p.
Dzielniki zależne od punktów kratowych (x,k) spełniają zależności
p*(x+1)+a, p*(k+1)+b.
Aby szybciej wyznaczyć dzielniki, z użyciem mniejszych liczb, trzeba odpowiednio przygotować a oraz b.
Pokażę to na przykładzie n = 8934053 = 1087*8219. Wartość p = 10.
Najmniej znacząca cyfra n to 3, szukam w tabliczce mnożenia iloczynów, które kończą się na 3. W systemie dziesiątkowym znajduję 4 takie. Ze względu na symetrię dwa można odrzucić. Są to pary (1,3) oraz (7,9).
Jeśli podstawa systemu jest liczbą pierwszą, takich par jest tyle, ile niezerowych cyfr - połowa się powtarza ze względu na symetrię. Ale tylko w tym kroku. W dalszych symetria najczęściej jest zaburzona oraz trzeba sprawdzać wszystkie przypadki (10 w dziesiątkowym).
Mamy już pary: 1*3 = 3 oraz 7*9=63. W podanym we wcześniejszym poście przebiegu powinniśmy odjąć je od n. Lecz teraz odcinamy najmniej znaczącą cyfrę zarówno od n, jak i od tego iloczynu. Uzyskamy odpowiednio e=0 oraz e=6, oraz n=893405.
Modyfikujemy sposób wyznaczania a, b oraz e. Zmieni to nieco sposób wyznaczania punktów kratowych.
Załóżmy, że wykorzystujemy parę (a,b) = (7,9), oraz przechodzimy do następnej cyfry: p = p*10 = 100, n=89340. Dołączymy cyfry bardziej znaczące do a, b. Zanim to jednak zrobimy, zapamiętamy wartość dziesiątek e=6 z iloczynu 63.
Dokładamy do a jako cyfrę najbardziej znaczącą kolejno c=0, 1, ..., 9, dopasowywując d w taki sposób, aby powstałe liczby
a*b = (c##7) * (d##9) %100 = 53.
Interesuje nas tylko cyfra dziesiątek 5. Zatem obliczenie sprowadza się do prostej kongruencji liniowej: 9c+7d+e=5, czyli 9c+7d+1=0.
Na przykład dla c=0 powstaje a=07, dalej 0*9+7*d+1=0 skąd d=7 oznacza, że b=79. Rzeczywiście 7*79=553.
W ogólności mnożymy dokładaną cyfrę do a przez cyfrę najmniej znaczącą b, dodajemy doń iloczyn cyfry dokładanej do b z najmniej znaczącą cyfrą a, fragment iloczynu wcześniejszych cyfr trzymaną w e oraz iloczyn dokładanych cyfr, przesunięty o p. W tym przypadku mamy
0*9+7*7+6+0*7*100 = 56.
Odcinamy cyfrę najmniej znaczącą 6 przyjmując nową wartość e=5. Jest to całość z dzielenia a*c przez p.
W kolejnym przypadku dokładamy do a=7 cyfrę 1, uzyskujemy 1*9+7*d+1=0 skąd d=0 oraz wartości a=17, b=09, oraz e = (1*9+0*1+6)/10 = 1.
Zamiast wyznaczać iloczyn a*b, dodamy do siebie aktualne wartości
g = a+b+e+1*p
U nas to będzie g = 07+79+5+1*100 = 191.
Sprawdzanie przypadków zmieniło się, należy sprawdzić podzielność (n-g)%(1##a), (n-g)%(1##b). W naszym przypadku (89340-191)%109 oraz (89340-191)%179. Punkty całkowite są punktami kratowymi krzywej (1), które wkładamy do równań na dzielniki.
Do następnej iteracji z p=1000 wartość n przejdzie jako n=8934.
Wtedy z przybliżenia dzielników powstałych z a=87, b=19 dla c=0 uzyskamy
0*9+2*d+6=0, skąd d=2, b=219, e=(0*19+2*87+16+0*2*100)/10 = 19;
g = 087+219+19+1*1000 = 1325;
8934-1325 = 7*1087 wskazuje dzielniki, co kończy obliczenia.
Oczywiście, wartości a=07, b=79 nie przybliżają dzielników, robią to wartości a=87, b=19. Ale i tak po przerobieniu co najwyżej 220 przypadków (440 dzieleń) będziemy znali dzielniki.
Dzielenie przez wartości nieparzyste to 542 przypadki.
Funkcja jest postaci
(1) k(x,n,p,a,b) = -(p+b)/p + n/(p*(px+p+a)
dla potęgi podstawy systemu p, liczby rozkładanej n, odpowiednio przygotowanych wartości a, b spełniajacych warunek
(2) a*b%p = n%p.
Dzielniki zależne od punktów kratowych (x,k) spełniają zależności
p*(x+1)+a, p*(k+1)+b.
Aby szybciej wyznaczyć dzielniki, z użyciem mniejszych liczb, trzeba odpowiednio przygotować a oraz b.
Pokażę to na przykładzie n = 8934053 = 1087*8219. Wartość p = 10.
Najmniej znacząca cyfra n to 3, szukam w tabliczce mnożenia iloczynów, które kończą się na 3. W systemie dziesiątkowym znajduję 4 takie. Ze względu na symetrię dwa można odrzucić. Są to pary (1,3) oraz (7,9).
Jeśli podstawa systemu jest liczbą pierwszą, takich par jest tyle, ile niezerowych cyfr - połowa się powtarza ze względu na symetrię. Ale tylko w tym kroku. W dalszych symetria najczęściej jest zaburzona oraz trzeba sprawdzać wszystkie przypadki (10 w dziesiątkowym).
Mamy już pary: 1*3 = 3 oraz 7*9=63. W podanym we wcześniejszym poście przebiegu powinniśmy odjąć je od n. Lecz teraz odcinamy najmniej znaczącą cyfrę zarówno od n, jak i od tego iloczynu. Uzyskamy odpowiednio e=0 oraz e=6, oraz n=893405.
Modyfikujemy sposób wyznaczania a, b oraz e. Zmieni to nieco sposób wyznaczania punktów kratowych.
Załóżmy, że wykorzystujemy parę (a,b) = (7,9), oraz przechodzimy do następnej cyfry: p = p*10 = 100, n=89340. Dołączymy cyfry bardziej znaczące do a, b. Zanim to jednak zrobimy, zapamiętamy wartość dziesiątek e=6 z iloczynu 63.
Dokładamy do a jako cyfrę najbardziej znaczącą kolejno c=0, 1, ..., 9, dopasowywując d w taki sposób, aby powstałe liczby
a*b = (c##7) * (d##9) %100 = 53.
Interesuje nas tylko cyfra dziesiątek 5. Zatem obliczenie sprowadza się do prostej kongruencji liniowej: 9c+7d+e=5, czyli 9c+7d+1=0.
Na przykład dla c=0 powstaje a=07, dalej 0*9+7*d+1=0 skąd d=7 oznacza, że b=79. Rzeczywiście 7*79=553.
W ogólności mnożymy dokładaną cyfrę do a przez cyfrę najmniej znaczącą b, dodajemy doń iloczyn cyfry dokładanej do b z najmniej znaczącą cyfrą a, fragment iloczynu wcześniejszych cyfr trzymaną w e oraz iloczyn dokładanych cyfr, przesunięty o p. W tym przypadku mamy
0*9+7*7+6+0*7*100 = 56.
Odcinamy cyfrę najmniej znaczącą 6 przyjmując nową wartość e=5. Jest to całość z dzielenia a*c przez p.
W kolejnym przypadku dokładamy do a=7 cyfrę 1, uzyskujemy 1*9+7*d+1=0 skąd d=0 oraz wartości a=17, b=09, oraz e = (1*9+0*1+6)/10 = 1.
Zamiast wyznaczać iloczyn a*b, dodamy do siebie aktualne wartości
g = a+b+e+1*p
U nas to będzie g = 07+79+5+1*100 = 191.
Sprawdzanie przypadków zmieniło się, należy sprawdzić podzielność (n-g)%(1##a), (n-g)%(1##b). W naszym przypadku (89340-191)%109 oraz (89340-191)%179. Punkty całkowite są punktami kratowymi krzywej (1), które wkładamy do równań na dzielniki.
Do następnej iteracji z p=1000 wartość n przejdzie jako n=8934.
Wtedy z przybliżenia dzielników powstałych z a=87, b=19 dla c=0 uzyskamy
0*9+2*d+6=0, skąd d=2, b=219, e=(0*19+2*87+16+0*2*100)/10 = 19;
g = 087+219+19+1*1000 = 1325;
8934-1325 = 7*1087 wskazuje dzielniki, co kończy obliczenia.
Oczywiście, wartości a=07, b=79 nie przybliżają dzielników, robią to wartości a=87, b=19. Ale i tak po przerobieniu co najwyżej 220 przypadków (440 dzieleń) będziemy znali dzielniki.
Dzielenie przez wartości nieparzyste to 542 przypadki.
09 września 2013
Wzory na budowę dzielników
Kolejny algorytm faktoryzacji opisywany w tym blogu powstał dzięki odwróceniu kolejności w geometrycznym mnożeniu w systemie binarnym.
Mnożenie geometryczne polega na zapisaniu cyfr dzielników przy brzegach prostokąta, wewnątrz zapisujemy iloczyny tych cyfr, przy czym cyfra dziesiątek jest 'na innej przekątnej' niż cyfra jedności. Następnie sumujemy po przekątnych, np. wartość na cyfrze dziesziątek powstaje z sumy:
a_1*b_0 + a_0*b_1 + (a_0*b_0)/10.
Po uwzględnieniu przeniesień mamy wynik.
Sposób ten zastosowany w systemie binarnym wskazuje prostą zależność między bitami dzielników, co pozwala wyznaczyć te dzielniki, niestety, tylko na 50%. Ale wystarczy.
Połączyłem te informacje z moim pierwszym algorytmem, w którym zakładałem, że dzielniki sa odpowiedniej postaci (2^i*c+a) i sprawdzałem, czy jest tak rzeczywiście.
Niektóre informacje okazały sie nadmiarowe, i tak oto powstał wzór na funkcję dzielników k(x) liczby n, przy potędze dwójki i:
(1) k = -(i+b)/i + n/(i (ix+i+a) )
W tym wzorze n jest rozkładaną liczbą, i jest potęgą 2, szukamy dzielników większych niż i, wartości a oraz b są resztami z dzielenia n przez i, spełniającymi zależność bitową:
(2) (n - a*b) AND (i-1) = 0
Oznacza to, że a*b oraz n mają dokładnie takie same bity najmniej znaczące, czyli stanowią swoiste 'bitowe przybliżenie' dzielników.
Po przekształceniu wzoru (1) do wspólnego mianownika uzyskujemy funkcję holomorficzną postaci
(3) k = (-bx+q)/(ix+a),
w której i jest potęgą 2, q powstaje z ilorazu różnicy n oraz a*b, co jest liczbą naturalną bliską n*i^(-2). Dzielniki znajdujemy, gdy funkcja ta ma punkt kratowy (x,k), tzn. x, k są całkowite (interesują nas tylko nieujemne).
Wspomniane wzory na dzielniki:
i*(k+1)+b
oraz
i+x*i+a
Przykład: n=8934053, i=1024, znaleziona para (a,b) = (27,63) spełnia zależność (2), ich iloczyn w szesnastkowym to 0x2A5, oraz n jest postaci n = 0x400*d + 0x2A5.
funkcja (1) dla tego przypadku:
k = -(1024+63)/1024 + 8934053/(1024(1024x+1024+27))
przekształca się do
k = -1087(x-7)/(1024x+1061)
Punkt kratowy (x,k) = (7,0) wstawiamy do wzorów na dzielniki:
1024*(0+1)+63 = 1087
1024+7*1024+27 = 8219
Sprawdzamy, że rzeczywiście 8934053 = 1087*8219.
Mając postać (3), wystarczy przeglądać x naturalne w poszukiwaniu punktów kratowych, licznik maleje, mianownik rośnie, wartości są nie większe niż n*i^(-2). Im większe i, tym mniej wartości mamy do sprawdzenia, gdyż dla x rosnącego iloraz zmniejsza się wykładniczo. Kiedy zejdzie poniżej 0, nie znajdziemy dzielników. Ale trzeba spradzić parę początkowych wartości, żeby dzielnik nie umknął.
Największym kłopotem może być dopasowanie pary (a,b). Ale zaczynając od i=4, gdy bit na pozycji 'dziesiątek' liczby n jest ustawiony (n&2=1), a=1, b=3. Gdy bit ten jest zgaszony (n&2=0), mamy dwie możliwości a=b=1 albo a=b=3.
Teraz zwiększajac i = 2*i, sprawdzamy kombinacje par (a,b), (a+i,b), (a, b+i), (a+i,b+1). Dwie z nich mają takie same bity najmniej znaczace jak n - te zostawiamy, pozostałe dwie mają znak różny i je usuwamy. W ten sposób tworzy nam się drzewo binarne przypadków, w którym co najmniej jedna gałąź wskazuje rozkład. Dla liczb pierwszych jest to (a,b) = (1,n).
Mnożenie geometryczne polega na zapisaniu cyfr dzielników przy brzegach prostokąta, wewnątrz zapisujemy iloczyny tych cyfr, przy czym cyfra dziesiątek jest 'na innej przekątnej' niż cyfra jedności. Następnie sumujemy po przekątnych, np. wartość na cyfrze dziesziątek powstaje z sumy:
a_1*b_0 + a_0*b_1 + (a_0*b_0)/10.
Po uwzględnieniu przeniesień mamy wynik.
Sposób ten zastosowany w systemie binarnym wskazuje prostą zależność między bitami dzielników, co pozwala wyznaczyć te dzielniki, niestety, tylko na 50%. Ale wystarczy.
Połączyłem te informacje z moim pierwszym algorytmem, w którym zakładałem, że dzielniki sa odpowiedniej postaci (2^i*c+a) i sprawdzałem, czy jest tak rzeczywiście.
Niektóre informacje okazały sie nadmiarowe, i tak oto powstał wzór na funkcję dzielników k(x) liczby n, przy potędze dwójki i:
(1) k = -(i+b)/i + n/(i (ix+i+a) )
W tym wzorze n jest rozkładaną liczbą, i jest potęgą 2, szukamy dzielników większych niż i, wartości a oraz b są resztami z dzielenia n przez i, spełniającymi zależność bitową:
(2) (n - a*b) AND (i-1) = 0
Oznacza to, że a*b oraz n mają dokładnie takie same bity najmniej znaczące, czyli stanowią swoiste 'bitowe przybliżenie' dzielników.
Po przekształceniu wzoru (1) do wspólnego mianownika uzyskujemy funkcję holomorficzną postaci
(3) k = (-bx+q)/(ix+a),
w której i jest potęgą 2, q powstaje z ilorazu różnicy n oraz a*b, co jest liczbą naturalną bliską n*i^(-2). Dzielniki znajdujemy, gdy funkcja ta ma punkt kratowy (x,k), tzn. x, k są całkowite (interesują nas tylko nieujemne).
Wspomniane wzory na dzielniki:
i*(k+1)+b
oraz
i+x*i+a
Przykład: n=8934053, i=1024, znaleziona para (a,b) = (27,63) spełnia zależność (2), ich iloczyn w szesnastkowym to 0x2A5, oraz n jest postaci n = 0x400*d + 0x2A5.
funkcja (1) dla tego przypadku:
k = -(1024+63)/1024 + 8934053/(1024(1024x+1024+27))
przekształca się do
k = -1087(x-7)/(1024x+1061)
Punkt kratowy (x,k) = (7,0) wstawiamy do wzorów na dzielniki:
1024*(0+1)+63 = 1087
1024+7*1024+27 = 8219
Sprawdzamy, że rzeczywiście 8934053 = 1087*8219.
Mając postać (3), wystarczy przeglądać x naturalne w poszukiwaniu punktów kratowych, licznik maleje, mianownik rośnie, wartości są nie większe niż n*i^(-2). Im większe i, tym mniej wartości mamy do sprawdzenia, gdyż dla x rosnącego iloraz zmniejsza się wykładniczo. Kiedy zejdzie poniżej 0, nie znajdziemy dzielników. Ale trzeba spradzić parę początkowych wartości, żeby dzielnik nie umknął.
Największym kłopotem może być dopasowanie pary (a,b). Ale zaczynając od i=4, gdy bit na pozycji 'dziesiątek' liczby n jest ustawiony (n&2=1), a=1, b=3. Gdy bit ten jest zgaszony (n&2=0), mamy dwie możliwości a=b=1 albo a=b=3.
Teraz zwiększajac i = 2*i, sprawdzamy kombinacje par (a,b), (a+i,b), (a, b+i), (a+i,b+1). Dwie z nich mają takie same bity najmniej znaczace jak n - te zostawiamy, pozostałe dwie mają znak różny i je usuwamy. W ten sposób tworzy nam się drzewo binarne przypadków, w którym co najmniej jedna gałąź wskazuje rozkład. Dla liczb pierwszych jest to (a,b) = (1,n).
27 sierpnia 2013
Porównywanie ułamków zwykłych za pomocą ciągu Fareya
Porównywanie ułamków zwykłych jest zagadnieniem z klasy szóstej podstawówki.
Sprowadza się ułamki do wspólnego mianownika i porównuje liczniki. Komputerowo spradza się znak 'iloczynu na krzyż' (licznik 1 * mianownik 2 - licznik 2 * mianownik 1).
Można to jednak zrobić bez pomocy mnożenia, dzielenia.
Co to jest ciąg Fareya. Mając dwa ułamki a/b oraz c/d, tworzymy nowy ułamek, którego licznik powstaje przez dodanie liczników, zaś mianownik jako suma mianowników. Nowo powstały ułamek (a+c)/(b+d) rozdziela ułamki, z których powstał w porządku liczb wymiernych.
Ciąg ten może służyć jako uzyskanie ciągu wszystkich liczb wymiernych przedziału [0/1; 1/1] przez wstawianie powstałych ułamków między już dane:
0/1 < 1/2 < 1/1;
0/1 < 1/3 < 1/2 < 2/3 < 1/1; itd.
Zastosowanie ciagu Fareya do porównywania ułamków a/c, b/d można zastosować wtedy, gdy licznik i mianownik jednego z ułamków są większe niż drugiego a<c, b<d. Wtedy traktujemy je jako pierwsze elementy ciągu Fareya, oraz wyliczamy trzeci następująco (c-a)/(b-d). Ten nowy ułamek ma wartości stosunkowo małe, i może być porównany z a/b, czasem nawet z 1/1. Uzyskany porządek przenosi się jako porównanie początkowych ułamków.
Kiedy nie można zastosować podanego kryterium, stosujemy szacowania:
- ten ułamek o równym mianowniku jest większy, gdy jego licznik jest większy;
- ten ułamek o równym liczniku jest większy, gdy jego mianownik jest mniejszy;
W ten sposób znajdujemy porządek ułamków np. 7/18 oraz 8/13, 7<8 oraz 18>13, zatem 7/18 < 8/13.
Porównamy 371/1243 z 721/1571. Skorzystamy z ciągu Fareya obliczając następny element
(721-371) / (1571-1243) = 350/328 > 1, bo 350>328
Zatem porządek jest następujący:
371/1243 < 721/1571 < 350/328
Przekształcenie Fareya zmniejsza wartości porównywanych ułamków. Można zapisać to procedurą rekursywną, u1 oznacza ułamek właściwy u1(l1,m1), l1<m1:
int cmpfr( u1, u2 ) {
if( u1!= u2 && l1<=l2 && m1>=m2 ) return u1<u2;
if( l1<l2 && m1<m2 ) {
if( u3(l2-l1, m2-m1) ) return cmpfr( u1, u3 ); // ulamek niewlasciwy
else return u1<u2;
}
return u1>=u2;
};
Sprowadza się ułamki do wspólnego mianownika i porównuje liczniki. Komputerowo spradza się znak 'iloczynu na krzyż' (licznik 1 * mianownik 2 - licznik 2 * mianownik 1).
Można to jednak zrobić bez pomocy mnożenia, dzielenia.
Co to jest ciąg Fareya. Mając dwa ułamki a/b oraz c/d, tworzymy nowy ułamek, którego licznik powstaje przez dodanie liczników, zaś mianownik jako suma mianowników. Nowo powstały ułamek (a+c)/(b+d) rozdziela ułamki, z których powstał w porządku liczb wymiernych.
Ciąg ten może służyć jako uzyskanie ciągu wszystkich liczb wymiernych przedziału [0/1; 1/1] przez wstawianie powstałych ułamków między już dane:
0/1 < 1/2 < 1/1;
0/1 < 1/3 < 1/2 < 2/3 < 1/1; itd.
Zastosowanie ciagu Fareya do porównywania ułamków a/c, b/d można zastosować wtedy, gdy licznik i mianownik jednego z ułamków są większe niż drugiego a<c, b<d. Wtedy traktujemy je jako pierwsze elementy ciągu Fareya, oraz wyliczamy trzeci następująco (c-a)/(b-d). Ten nowy ułamek ma wartości stosunkowo małe, i może być porównany z a/b, czasem nawet z 1/1. Uzyskany porządek przenosi się jako porównanie początkowych ułamków.
Kiedy nie można zastosować podanego kryterium, stosujemy szacowania:
- ten ułamek o równym mianowniku jest większy, gdy jego licznik jest większy;
- ten ułamek o równym liczniku jest większy, gdy jego mianownik jest mniejszy;
W ten sposób znajdujemy porządek ułamków np. 7/18 oraz 8/13, 7<8 oraz 18>13, zatem 7/18 < 8/13.
Porównamy 371/1243 z 721/1571. Skorzystamy z ciągu Fareya obliczając następny element
(721-371) / (1571-1243) = 350/328 > 1, bo 350>328
Zatem porządek jest następujący:
371/1243 < 721/1571 < 350/328
Przekształcenie Fareya zmniejsza wartości porównywanych ułamków. Można zapisać to procedurą rekursywną, u1 oznacza ułamek właściwy u1(l1,m1), l1<m1:
int cmpfr( u1, u2 ) {
if( u1!= u2 && l1<=l2 && m1>=m2 ) return u1<u2;
if( l1<l2 && m1<m2 ) {
if( u3(l2-l1, m2-m1) ) return cmpfr( u1, u3 ); // ulamek niewlasciwy
else return u1<u2;
}
return u1>=u2;
};
23 sierpnia 2013
Przekształcenia przed dzieleniami - odczyt bitowy dzielników
Jeden z algorytmów mnożenia zastosowany w systemie binarnym wskazuje, że znajomość najmniej znaczących bitów dzielnika pozwala z odpowiedniego bitu iloczynu odczytać kolejny bit obu z dzielników. Są dwie możliwości, między innymi ze wględu na prawo przemienności.
Zatem podejście jest takie. Liczbę rozkładaną n zapisuję w postaci czwórki liczb (a, b, r, c), gdzie a oraz b są kandydatami na dzielniki, do których dołączam przewidywany najbardziej znaczący bit. Liczba r to rezerwa wskazująca, jak wartość jest jeszcze dostępna, gdy wykorzystam informacje o bieżących ustawieniach liczb a oraz b. Wartość c jest potęgą 2, wskazuje, który bit jest aktualnie dopasowywany.
Przebieg. Najpierw przekształcam wartości czwórki, następnie sprawdzam pary wartości oraz wywołuję rekursywnie obie możliwości.
Inicjacja czwórki jest następująca:
(a=1, b=1, r=n-1, c=1)
oraz oznacza, że n = r*c + a*b = (n-1)*1 + 1*1 = n-1 +1.
Teraz w zależności od najmniej znaczącego bitu r, mam dwie możliwości:
bit r&1 wyzerowany:
a'=a, b'=b, r'=r;
albo
a''=a+c, b''=b+c, r''=r-a-b-c;
bit r&1 ustawiony:
a'=a, b'=b+c, r'=r-a;
albo
a''=a+c, b''=b, r''=r-b;
Po tych przekształceniach przyglądam się własnościom par liczb (a',r') oraz (b',r'). Interesuje mnie zwłaszcza podzielność, np. a'|r'.
Wyodrębniłem następujące warunki, stosując parę (d,r):
- d=1, wskazuje rozkład trywialny n = 1*n;
- d|r, sprawdzam teraz, czy d|n. Jeśli tak, d jest dzielnikiem liczby n, jeśli nie, liczba n jest pierwsza. W obu przypadkach mogę zakończyć rozkład;
podobnie należy sprawdzić, kiedy r staje się zerem;
- d>r, rezerwy są zbyt małe, aby móc zwiększać kandydata na dzielnika, ta gałąż rekurencji jest bezużyteczna, można uciąć też drugą gałąź z tego węzła.
Wywołanie rekursywne dzieli rezerwę na 2 oraz przesuwa c na kolejny bit bardziej znaczący, co w praktyce oznacza wywołanie
(a', b', r'/2, 2*c) oraz (a'', b'', r''/2, 2*c).
Zachowanie się wartości.
Liczby a oraz b w kolejnych wywołaniach stanowią ciągi niemalejące, przyrost jest potęgą 2.
Liczba r stanowi ciąg malejący, ograniczony z góry przez n*2^(-i), i jest licznikiem wywołania. Dzielna jest też ograniczona przez tę wartość, zaś dzielnik dosyć szybko rośnie. Przypadek podzielności występuje parzystą krotność razy, dla każdego z dzielników, także dla liczb pierwszych i nie zawsze wskazuje kandydata na dzielnik trywialny.
Kiedy a stanie się większe niż r, rekursja się kończy, co powoduje, że kroność wywołania stanowi połowę logarytmu binarnego z n. Sprawdzamy 2^(1+1/2 lg n) możliwych ustawień bitów w dzielnikach.
Przykład: 143
czwórka inicjująca (1, 1, 142, 1),
pierwsze wywołanie rekursji zwraca (1, 1, 71, 2)
Ponieważ ~2|71, mamy przypadki czwórek (1, 3, 71-1=70, 2) oraz symetryczną (3, 1, 70, 2), 1|70 i jest przypadkiem trywialnym;
Kolejne wywołania rekursywne zwracają (1, 3, 35, 4) oraz (3, 1, 35, 4).
W obu przypadkach ~2|35, zajmijmy się pierwszym:
mamy (1, 7, 34, 4) oraz (5, 3, 32, 4);
wywołania rekursywne (1, 7, 17, 8) oraz (5, 3, 16, 8);
pierwsze z nich daje (1, 15, 16, 8) oraz (9, 7, 10, 8);
po kolejnym wywołaniu rekursywnym obydwu przypadków uzyskujemy r mniejsze niż któryś z kandydatów na dzielniki, należy zająć się (5, 3, 16, 8).
Mamy 2|16, czyli czwórkę (5, 3, 16, 8), albo (13, 11, 0, 8),
Wartość r=0, sprawdzamy, czy 13|143, i rzeczywiście, mamy rozkład 143=13*11.
Zatem podejście jest takie. Liczbę rozkładaną n zapisuję w postaci czwórki liczb (a, b, r, c), gdzie a oraz b są kandydatami na dzielniki, do których dołączam przewidywany najbardziej znaczący bit. Liczba r to rezerwa wskazująca, jak wartość jest jeszcze dostępna, gdy wykorzystam informacje o bieżących ustawieniach liczb a oraz b. Wartość c jest potęgą 2, wskazuje, który bit jest aktualnie dopasowywany.
Przebieg. Najpierw przekształcam wartości czwórki, następnie sprawdzam pary wartości oraz wywołuję rekursywnie obie możliwości.
Inicjacja czwórki jest następująca:
(a=1, b=1, r=n-1, c=1)
oraz oznacza, że n = r*c + a*b = (n-1)*1 + 1*1 = n-1 +1.
Teraz w zależności od najmniej znaczącego bitu r, mam dwie możliwości:
bit r&1 wyzerowany:
a'=a, b'=b, r'=r;
albo
a''=a+c, b''=b+c, r''=r-a-b-c;
bit r&1 ustawiony:
a'=a, b'=b+c, r'=r-a;
albo
a''=a+c, b''=b, r''=r-b;
Po tych przekształceniach przyglądam się własnościom par liczb (a',r') oraz (b',r'). Interesuje mnie zwłaszcza podzielność, np. a'|r'.
Wyodrębniłem następujące warunki, stosując parę (d,r):
- d=1, wskazuje rozkład trywialny n = 1*n;
- d|r, sprawdzam teraz, czy d|n. Jeśli tak, d jest dzielnikiem liczby n, jeśli nie, liczba n jest pierwsza. W obu przypadkach mogę zakończyć rozkład;
podobnie należy sprawdzić, kiedy r staje się zerem;
- d>r, rezerwy są zbyt małe, aby móc zwiększać kandydata na dzielnika, ta gałąż rekurencji jest bezużyteczna, można uciąć też drugą gałąź z tego węzła.
Wywołanie rekursywne dzieli rezerwę na 2 oraz przesuwa c na kolejny bit bardziej znaczący, co w praktyce oznacza wywołanie
(a', b', r'/2, 2*c) oraz (a'', b'', r''/2, 2*c).
Zachowanie się wartości.
Liczby a oraz b w kolejnych wywołaniach stanowią ciągi niemalejące, przyrost jest potęgą 2.
Liczba r stanowi ciąg malejący, ograniczony z góry przez n*2^(-i), i jest licznikiem wywołania. Dzielna jest też ograniczona przez tę wartość, zaś dzielnik dosyć szybko rośnie. Przypadek podzielności występuje parzystą krotność razy, dla każdego z dzielników, także dla liczb pierwszych i nie zawsze wskazuje kandydata na dzielnik trywialny.
Kiedy a stanie się większe niż r, rekursja się kończy, co powoduje, że kroność wywołania stanowi połowę logarytmu binarnego z n. Sprawdzamy 2^(1+1/2 lg n) możliwych ustawień bitów w dzielnikach.
Przykład: 143
czwórka inicjująca (1, 1, 142, 1),
pierwsze wywołanie rekursji zwraca (1, 1, 71, 2)
Ponieważ ~2|71, mamy przypadki czwórek (1, 3, 71-1=70, 2) oraz symetryczną (3, 1, 70, 2), 1|70 i jest przypadkiem trywialnym;
Kolejne wywołania rekursywne zwracają (1, 3, 35, 4) oraz (3, 1, 35, 4).
W obu przypadkach ~2|35, zajmijmy się pierwszym:
mamy (1, 7, 34, 4) oraz (5, 3, 32, 4);
wywołania rekursywne (1, 7, 17, 8) oraz (5, 3, 16, 8);
pierwsze z nich daje (1, 15, 16, 8) oraz (9, 7, 10, 8);
po kolejnym wywołaniu rekursywnym obydwu przypadków uzyskujemy r mniejsze niż któryś z kandydatów na dzielniki, należy zająć się (5, 3, 16, 8).
Mamy 2|16, czyli czwórkę (5, 3, 16, 8), albo (13, 11, 0, 8),
Wartość r=0, sprawdzamy, czy 13|143, i rzeczywiście, mamy rozkład 143=13*11.
07 sierpnia 2013
Faktoryzacja, przepis z prostymi przekształceniami
Skrzyżowałem faktoryzację sprawdzającą równocześnie małe i duże dzielniki z ostatnią faktoryzacją, w której sprawdzam zaledwie 2^i przypadków zamiast pierwiastka z liczby faktoryzowanej. Uzyskany algorytm cechują działania na zaledwie czterech liczbach osiągających wielkość pierwiastka z liczby rozkładanej.
Przepis jest następujący: najpierw dzielę liczbę n na dwa, w przybliżeniu równe kawałki, które dostarczą mi postać liczby trójcyfrowej liczby n. Mogę teraz zatrudnić trzy wątki. Jeden stosuje konwersję o 2, podczas której podstawa rośnie do pierwiastka z n, drugi stosuje konwersję o -2 sprawdzając dzielniki bliskie sobie, które są nieosiągalne przez pierwszy z wątków. Wreszcie trzeci przygląda się podstawie, oraz uwzględniając jej dzielnik trzy (powtarzający się co trzy iteracje w każdym z dwu wątków), eliminuje małe dzielniki przy współudziale wyrazu wolnego.
Ze względu na specyfikę obliczeń, program wymaga pamięci do korzystania z kopii niektórych wartości. Dostępne są następujące wartości:
a, b, c, p,
spełniające warunek n = (a*p+b)*p+c, indeksy podzielności p przez 3 (aby nie obliczać p%3 w każdej iteracji) oraz pewna liczba pomocnicza.
Kryterium zatrzymania stanowi wartość x, uzyskana pod koniec pracy pierwszego wątku. Zatrzymuje ona działanie drugiego wątku ograniczając liczność sprawdzanych przypadków. Trzeci wątek po znalezieniu dzielnika zatrzymuje cały program.
Przejdźmy do szczegółów, które przedstawię w taki sposób, jakby były omawiane po raz pierwszy.
NAPRAWA
Najpierw sposób naprawy postaci [a, b, c, p]. Spełniane są warunki: 0<a<4, -1<b<p, -1<c<p oraz n = (a*p+b)+c.
Jeśli b jest ujemne, dodajemy p do b oraz odejmujemy od a jedynkę:
if( 0>b ) { b+=p; a--; }
Jeśli c jest ujemne, dodajemy p do c oraz odejmujemy od b jedynkę:
if( 0>c ) { c+=p; b--; }
Jeśli b jest nie mniejsze niż p, zmniejszamy b o p oraz zwiększamy a:
if( b>p-1 ) { b-=p; a++; }
Podobnie c, jeśli jest zbyt duże, zmniejszamy je o p zwiększając b:
if( c>p-1 ) { c-=p; b--; }
Operacje te wystarczy powtórzyć do dwu razy w dowolnym z wątków przy uruchamianiu funkcji naprawa(), aby zapewnić wymagane warunki.
ALGORYTM
Przebieg algorytmu, liczbowe dane są z powietrza:
Podział liczby n spełniajacej warunek 2^(2i) < n < 2^(2n+2) na dwie liczby mniej więcej tej samej długości spełniające warunki:
p=2^i lub p=2^(i+1); b = n/p; c=n%p, b>p, n = b*p+c
np. dla uzyskanego p=1024, b nie przekracza 4096, zaś c jest ograniczone w przedziale [0, 1023).
Z liczby b wydobywamy a=b/p; b=b%p. Uzyskujemy wtedy wstępną czwórkę liczb postaci:
[a, b, c, p] , n = (a*p+b)*p+c, a<4, b<p, c<p
Wartość podstawy p powinna być nieparzysta, i najlepiej podzielna przez 3.
Jeśli p jest parzyste, zwiększymy podstawę p o 1 przekształceniami:
{
p++;
b-=a;
naprawa();
c-=b;
b-=a;
naprawa();
}
WĄTEK 1
pętla zwiększajaca podstawę p, zmniejszająca skokowo a. Kończymy, gdy a osiagnie 0, gdyż wtedy wartość podstawy p będzie większa niż pierwiastek kwadratowy z n. Ciąg [a,b,c,p] w tym momencie może wyglądać następująco: [1, 7, 482 953, 3 593 385],
[1, 3, 482 943, 3 593 387],
[0, 3 593 388, 482 941, 3 593 389]
Wyjście z programu następuje wtedy, gdy c będzie zerem. Znamy wtedy dzielniki: n = p * (a*p+b).
{
p+=2;
b -= 2*a;
naprawa();
c -= 2*b;
b -= 2*a;
naprawa();
}
jeśli p jest liczbą podzielną przez 3 (co sprawdzamy odpowiednim licznikiem) przesyłamy parę [c,p] do wątku 3. Wartość taka powtarza się co trzy iteracje, np. dla p=17 278 137, później dla p+6: 17 278 143 itd.
WĄTEK 2
pętla zmniejszająca podstawę p, zwiększajaca skokowo a. Kończymy, gdy trzeci z wątków nakaże zakończyć.
Wyjście z programu następuje jak przy wątku 1, gdy c będzie zerem. Mamy wtedy dzielniki: n = p * (a*p+b). Pętla:
{
p-=2;
b += 2*a;
naprawa();
c += 2*b;
b += 2*a;
naprawa();
}
Jeśli p jest liczbą podzielną przez 3 (znowu licznik), przesyłamy parę [c,p] do wątku 3 podobnie jak przy wątku piewszym.
WĄTEK 3
Wątek ten sprawdza małe dzielniki, oraz zatrzymuje algorytm, kiedy zostanie sprawdzone 2^i<sqrt(n) przypadków. Bazuje na obserwacji:
liczba n ma dzielnik d, gdy podstawa systemu p oraz wyraz wolny c są wielokrotnościami d.
Zatem najpierw zmniejszamy kopię podstawy p dzieląc ją przez 3, a następnie liczymy nwd(p,c) by znaleźć wspólny dzielnik. Gdyż jeśli dzielnikiem jest wielokrotność 3, to dzielnikiem jest też 3.
{
while ( 0==(p%3) ) p/=3;
dzielnik d = nwd(p,c); // jeśli 1==p, stosujemy p=3
}
jeśli przesłane zostanie p z wątku pierwszego, zapamiętywana jest wartość x = p/3. Jeśli wątek drugi prześle p mniejsze niż x, oznacza to, że ten przypadek został już sprawdzony. Kiedy pierwszy z wątków zakończy swoje działanie, jego ostatnia wartość p/3 stanowi równocześnie liczbę, do której wystarczy sprawdzać w wątku 2, aby wyeliminować wszystkie możliwe dzielniki.
Podsumowanie
Wątki 1 oraz 2 sprawdzają bardzo duże dzielniki, nieco mniejsze niż pierwiastek z rozkładanej liczby. Wątek 3 sprawdza dzielniki mniejsze niż 2^i< sqrt(n), oraz zapobiega ponownemu przeliczaniu. Wątek 2 przy zmniejszaniu podstawy od pewnego miejsca zaczyna podawać dzielniki sprawdzone już przez wątek 3, czemu należy zapobiec.
Wątek 1 powinien być dominujący nad drugim, na ogół ma najwięcej przypadków do sprawdzenia. Wątek 3 jest najbardziej czasochłonny, wymaga najgorszych operacji, oraz może korzystać z większej krotności procesorów.
Wartości liczbowe w wątku pierwszym: a maleje, b i c skokowo maleją, p rośnie jednostajnie.
Wartości liczbowe w wątku drugim: a rośnie, b i c skokowo rosną, p maleje jednostajnie.
Wartości liczbowe w wątku trzecim: jedyne co można przewidzieć, że p nie przekroczy pierwiastka kwadratowego, zaś c jest ograniczone przez p. Wartość c może być stosunkowo mała. Przy kolejnych dostawach p przez wątki, krotność dzieleń przez 3 zachowuje się jak grzebień: 1, 2, 1, 3, 1, 2, 1, 4, 1, 2, 1, 3, 1, itd.
Przepis jest następujący: najpierw dzielę liczbę n na dwa, w przybliżeniu równe kawałki, które dostarczą mi postać liczby trójcyfrowej liczby n. Mogę teraz zatrudnić trzy wątki. Jeden stosuje konwersję o 2, podczas której podstawa rośnie do pierwiastka z n, drugi stosuje konwersję o -2 sprawdzając dzielniki bliskie sobie, które są nieosiągalne przez pierwszy z wątków. Wreszcie trzeci przygląda się podstawie, oraz uwzględniając jej dzielnik trzy (powtarzający się co trzy iteracje w każdym z dwu wątków), eliminuje małe dzielniki przy współudziale wyrazu wolnego.
Ze względu na specyfikę obliczeń, program wymaga pamięci do korzystania z kopii niektórych wartości. Dostępne są następujące wartości:
a, b, c, p,
spełniające warunek n = (a*p+b)*p+c, indeksy podzielności p przez 3 (aby nie obliczać p%3 w każdej iteracji) oraz pewna liczba pomocnicza.
Kryterium zatrzymania stanowi wartość x, uzyskana pod koniec pracy pierwszego wątku. Zatrzymuje ona działanie drugiego wątku ograniczając liczność sprawdzanych przypadków. Trzeci wątek po znalezieniu dzielnika zatrzymuje cały program.
Przejdźmy do szczegółów, które przedstawię w taki sposób, jakby były omawiane po raz pierwszy.
NAPRAWA
Najpierw sposób naprawy postaci [a, b, c, p]. Spełniane są warunki: 0<a<4, -1<b<p, -1<c<p oraz n = (a*p+b)+c.
Jeśli b jest ujemne, dodajemy p do b oraz odejmujemy od a jedynkę:
if( 0>b ) { b+=p; a--; }
Jeśli c jest ujemne, dodajemy p do c oraz odejmujemy od b jedynkę:
if( 0>c ) { c+=p; b--; }
Jeśli b jest nie mniejsze niż p, zmniejszamy b o p oraz zwiększamy a:
if( b>p-1 ) { b-=p; a++; }
Podobnie c, jeśli jest zbyt duże, zmniejszamy je o p zwiększając b:
if( c>p-1 ) { c-=p; b--; }
Operacje te wystarczy powtórzyć do dwu razy w dowolnym z wątków przy uruchamianiu funkcji naprawa(), aby zapewnić wymagane warunki.
ALGORYTM
Przebieg algorytmu, liczbowe dane są z powietrza:
Podział liczby n spełniajacej warunek 2^(2i) < n < 2^(2n+2) na dwie liczby mniej więcej tej samej długości spełniające warunki:
p=2^i lub p=2^(i+1); b = n/p; c=n%p, b>p, n = b*p+c
np. dla uzyskanego p=1024, b nie przekracza 4096, zaś c jest ograniczone w przedziale [0, 1023).
Z liczby b wydobywamy a=b/p; b=b%p. Uzyskujemy wtedy wstępną czwórkę liczb postaci:
[a, b, c, p] , n = (a*p+b)*p+c, a<4, b<p, c<p
Wartość podstawy p powinna być nieparzysta, i najlepiej podzielna przez 3.
Jeśli p jest parzyste, zwiększymy podstawę p o 1 przekształceniami:
{
p++;
b-=a;
naprawa();
c-=b;
b-=a;
naprawa();
}
WĄTEK 1
pętla zwiększajaca podstawę p, zmniejszająca skokowo a. Kończymy, gdy a osiagnie 0, gdyż wtedy wartość podstawy p będzie większa niż pierwiastek kwadratowy z n. Ciąg [a,b,c,p] w tym momencie może wyglądać następująco: [1, 7, 482 953, 3 593 385],
[1, 3, 482 943, 3 593 387],
[0, 3 593 388, 482 941, 3 593 389]
Wyjście z programu następuje wtedy, gdy c będzie zerem. Znamy wtedy dzielniki: n = p * (a*p+b).
{
p+=2;
b -= 2*a;
naprawa();
c -= 2*b;
b -= 2*a;
naprawa();
}
jeśli p jest liczbą podzielną przez 3 (co sprawdzamy odpowiednim licznikiem) przesyłamy parę [c,p] do wątku 3. Wartość taka powtarza się co trzy iteracje, np. dla p=17 278 137, później dla p+6: 17 278 143 itd.
WĄTEK 2
pętla zmniejszająca podstawę p, zwiększajaca skokowo a. Kończymy, gdy trzeci z wątków nakaże zakończyć.
Wyjście z programu następuje jak przy wątku 1, gdy c będzie zerem. Mamy wtedy dzielniki: n = p * (a*p+b). Pętla:
{
p-=2;
b += 2*a;
naprawa();
c += 2*b;
b += 2*a;
naprawa();
}
Jeśli p jest liczbą podzielną przez 3 (znowu licznik), przesyłamy parę [c,p] do wątku 3 podobnie jak przy wątku piewszym.
WĄTEK 3
Wątek ten sprawdza małe dzielniki, oraz zatrzymuje algorytm, kiedy zostanie sprawdzone 2^i<sqrt(n) przypadków. Bazuje na obserwacji:
liczba n ma dzielnik d, gdy podstawa systemu p oraz wyraz wolny c są wielokrotnościami d.
Zatem najpierw zmniejszamy kopię podstawy p dzieląc ją przez 3, a następnie liczymy nwd(p,c) by znaleźć wspólny dzielnik. Gdyż jeśli dzielnikiem jest wielokrotność 3, to dzielnikiem jest też 3.
{
while ( 0==(p%3) ) p/=3;
dzielnik d = nwd(p,c); // jeśli 1==p, stosujemy p=3
}
jeśli przesłane zostanie p z wątku pierwszego, zapamiętywana jest wartość x = p/3. Jeśli wątek drugi prześle p mniejsze niż x, oznacza to, że ten przypadek został już sprawdzony. Kiedy pierwszy z wątków zakończy swoje działanie, jego ostatnia wartość p/3 stanowi równocześnie liczbę, do której wystarczy sprawdzać w wątku 2, aby wyeliminować wszystkie możliwe dzielniki.
Podsumowanie
Wątki 1 oraz 2 sprawdzają bardzo duże dzielniki, nieco mniejsze niż pierwiastek z rozkładanej liczby. Wątek 3 sprawdza dzielniki mniejsze niż 2^i< sqrt(n), oraz zapobiega ponownemu przeliczaniu. Wątek 2 przy zmniejszaniu podstawy od pewnego miejsca zaczyna podawać dzielniki sprawdzone już przez wątek 3, czemu należy zapobiec.
Wątek 1 powinien być dominujący nad drugim, na ogół ma najwięcej przypadków do sprawdzenia. Wątek 3 jest najbardziej czasochłonny, wymaga najgorszych operacji, oraz może korzystać z większej krotności procesorów.
Wartości liczbowe w wątku pierwszym: a maleje, b i c skokowo maleją, p rośnie jednostajnie.
Wartości liczbowe w wątku drugim: a rośnie, b i c skokowo rosną, p maleje jednostajnie.
Wartości liczbowe w wątku trzecim: jedyne co można przewidzieć, że p nie przekroczy pierwiastka kwadratowego, zaś c jest ograniczone przez p. Wartość c może być stosunkowo mała. Przy kolejnych dostawach p przez wątki, krotność dzieleń przez 3 zachowuje się jak grzebień: 1, 2, 1, 3, 1, 2, 1, 4, 1, 2, 1, 3, 1, itd.
24 lipca 2013
Faktoryzacja metodą (nie kolejnych) dzieleń
Liczbę złożoną nieparzystą n = p*q można przedstawić za pomocą przedstawienia liczby q = (q_m,..., q_1, q_0) w systemie binarnym następującą tablicą wartości:
n = (p*q_m, ..., p*q_1, p*q_0).
Postać ta charakteryzuje się zerami na miejscach, w których q_i=0, oraz wartościami p na miejscach, w których w zapisie binarnym q_i=1, i=0..m. Traktując poszczególne pola tablicy jako wielomian W[2], uzyskamy wartość liczby n.
Np. n = 125 możemy przedstawić jako n = (5, 5, 0, 0, 5), gdyż 125 = 5* 25, oraz 25 = 11001b. Wartość jest równa
W[2] = ((((5*2+5)*2+0)*2+0)*2+5 = 125.
Czy podana postać może służyć do znajdywania rozkładu liczby na czynniki?
W jaki sposób?
Okazuje się, że postać ta wskazuje kolejne sposoby faktoryzacji. Dysponując tabelą zapisu binarnego liczby n, za pomocą operacji niezmienniczych
n_i = 2 * n_{i-1}, i=1,...,m
można doprowadzić do postaci binarnej dzielnika.
Postępowanie zaproponowane nieco dalej sprowadza się do dzieleń przez kolejne liczby nieparzyste, lecz w nieco innej kolejności niż normalnie, z mniejszymi wartościami dzielnych. Pierwotnie próbowałem szukać, czy znajdę sposób na domyślenie się postaci binarnej dzielnika.
Pierwsze próby wyrównywania wartości w zapisie binarnym liczby n czasem skutkowały bardzo szybkim znajdywaniem dzielnika, czasem zaś wykazywały, że ten nie istnieje. Kłopot polegał na tym, że żaden algorytm nie jest w stanie przewidzieć, które pola tabeli należy wyzerować, a które zostawić ustawione, czyli przewidzieć budowę czynnika q. Należy to robić losowo, algorytmem genetycznym lub sprawdzać wszystkie możliwe przypadki.
Algorytm genetyczny ma duże szanse zapętlenia. Dodatkowo wprowadza bardzo dużo operacji wyrównujących poszczególne wartości tabeli.
Lepiej jest zastosować kombinatorykę, wtedy żegnamy się z dopasowywaniem bitów stosując wykładniczą metodę siłową.
Naszym celem jest przygotować tabelę wypełnioną równomiernie daną wartością dodatnią. Następnie zerować poszczególne pola tablicy, starając się zachować taką samą wartość na wszystkich ustawionych pozycjach.
Jeśli się to uda, dodatnia wartość w tablicy jest wartością jednego z dzielników, zaś budowa tablicy potraktowana jako liczba binarna drugim.
Praktycznie wprowadzamy komórkę reszty, do której odkładamy wszelkie nadmiary, powstały z dzielenia pierwotnego, jak i z zerowania pól tablicy.
Tabelę wypełnimy równomiernie jedną wartością b = floor(n : (2^i-1)). Dzielimy tu przez np, 3, 7, 15, 31, 63, ... w zależności od ewentualnej długości dzielnika, którego szukamy.
Pierwsze i ostatnie pole są niepuste. Przygotujemy szablon będący liczbą binarną s = (2^i-1).
Wyzerowanie bitu szablonu na pozycji k odpowiada wyzerowaniu pola tablicy, równocześnie wartość b*2^k jest dodawana do reszty. Można ją zapisać instrukcją w C++ dla maski bitowej t:
r += (s & ~t)*b.
Jeżeli teraz reszta r jest wielokrotnością szablonu (s & t), wartość reszty można bez kłopotu rozdzielić między niepuste pola tablicy, o co nam chodziło. Praktycznie sprawdzamy, czy wartość liczbowa dzieli resztę (s & t) | r. Jeśli tak, mamy do czynienia z dzielnikami (s&t) * (b+r/(s&t)) = n.
Powtarzamy dla dowolnej wariacji ustawienia pola tablicy z wyjątkiem wartości skrajnych, lub do znalezienia dzielników.
Ostatecznie dzielimy najpierw liczbę n przez s=(2^i-1) dla pewnego i, a następnie zerując odpowiednie bity szablonu dzielimy resztę przez szablon. Reszta jest liczbą nie przekraczającą połowy n, zaś liczby zadane przez schemat mają ściśle określoną długość binarną. Można też wykorzystać reszty modulo schemat.
Można się jeszcze pokusić o to, by liczby zadane szablonem nie przekroczyły pierwiastka kwadratowego z n. Zmniejszy to krotność obliczeń.
Pierwsze z dzieleń jest tak specyficzne, że można wykorzystać przyspieszone dzielenie przez odpowiedniki dziewiątki. Polega ono na pocięciu liczby na paczki długości i cyfr, na wynik składają się zapisane kolejno wartości sumy narastającej paczek.
Przykład: 1441 = 11 * 131
Dzielimy 1441 : 3 = 480 z resztą 1, nic dalej nie można zrobić
Dzielimy 1441 : 7 = 205 z resztą 6
schemat 101b każe sprawdzić dzielenie
(2*205+6) : 101b = 416 : 5 \equiv 1 (mod 5)
Dzielimy 1441 : 15 = 96 z resztą 1
Schemat generuje 3 dzielenia przez odpowiednio 1001b, 1011b oraz 1101b.
W pierwszym przypadku mamy dzielenie
((2+4)*96+1) : 1001b = 577 : 9 \equiv 1 (mod 9),
w drugim
(4*96+1) : 1011b = 385 : 11 \equiv 0 (mod 11)
oraz dzielniki 11 oraz 96+385/11 = 96+35 = 131; w trzecim
(2*96+1) : 1101b = 193 : 13 = 11 (mod 13) .
Następuje zamiana dzieleń:
1441 : 3 1441:3
1441 : 5 416 : 5
1441 : 7 1441:7
1441 : 9 577 : 9
1441 : 11 385 : 11
193 : 13
n = (p*q_m, ..., p*q_1, p*q_0).
Postać ta charakteryzuje się zerami na miejscach, w których q_i=0, oraz wartościami p na miejscach, w których w zapisie binarnym q_i=1, i=0..m. Traktując poszczególne pola tablicy jako wielomian W[2], uzyskamy wartość liczby n.
Np. n = 125 możemy przedstawić jako n = (5, 5, 0, 0, 5), gdyż 125 = 5* 25, oraz 25 = 11001b. Wartość jest równa
W[2] = ((((5*2+5)*2+0)*2+0)*2+5 = 125.
Czy podana postać może służyć do znajdywania rozkładu liczby na czynniki?
W jaki sposób?
Okazuje się, że postać ta wskazuje kolejne sposoby faktoryzacji. Dysponując tabelą zapisu binarnego liczby n, za pomocą operacji niezmienniczych
n_i = 2 * n_{i-1}, i=1,...,m
można doprowadzić do postaci binarnej dzielnika.
Postępowanie zaproponowane nieco dalej sprowadza się do dzieleń przez kolejne liczby nieparzyste, lecz w nieco innej kolejności niż normalnie, z mniejszymi wartościami dzielnych. Pierwotnie próbowałem szukać, czy znajdę sposób na domyślenie się postaci binarnej dzielnika.
Pierwsze próby wyrównywania wartości w zapisie binarnym liczby n czasem skutkowały bardzo szybkim znajdywaniem dzielnika, czasem zaś wykazywały, że ten nie istnieje. Kłopot polegał na tym, że żaden algorytm nie jest w stanie przewidzieć, które pola tabeli należy wyzerować, a które zostawić ustawione, czyli przewidzieć budowę czynnika q. Należy to robić losowo, algorytmem genetycznym lub sprawdzać wszystkie możliwe przypadki.
Algorytm genetyczny ma duże szanse zapętlenia. Dodatkowo wprowadza bardzo dużo operacji wyrównujących poszczególne wartości tabeli.
Lepiej jest zastosować kombinatorykę, wtedy żegnamy się z dopasowywaniem bitów stosując wykładniczą metodę siłową.
Naszym celem jest przygotować tabelę wypełnioną równomiernie daną wartością dodatnią. Następnie zerować poszczególne pola tablicy, starając się zachować taką samą wartość na wszystkich ustawionych pozycjach.
Jeśli się to uda, dodatnia wartość w tablicy jest wartością jednego z dzielników, zaś budowa tablicy potraktowana jako liczba binarna drugim.
Praktycznie wprowadzamy komórkę reszty, do której odkładamy wszelkie nadmiary, powstały z dzielenia pierwotnego, jak i z zerowania pól tablicy.
Tabelę wypełnimy równomiernie jedną wartością b = floor(n : (2^i-1)). Dzielimy tu przez np, 3, 7, 15, 31, 63, ... w zależności od ewentualnej długości dzielnika, którego szukamy.
Pierwsze i ostatnie pole są niepuste. Przygotujemy szablon będący liczbą binarną s = (2^i-1).
Wyzerowanie bitu szablonu na pozycji k odpowiada wyzerowaniu pola tablicy, równocześnie wartość b*2^k jest dodawana do reszty. Można ją zapisać instrukcją w C++ dla maski bitowej t:
r += (s & ~t)*b.
Jeżeli teraz reszta r jest wielokrotnością szablonu (s & t), wartość reszty można bez kłopotu rozdzielić między niepuste pola tablicy, o co nam chodziło. Praktycznie sprawdzamy, czy wartość liczbowa dzieli resztę (s & t) | r. Jeśli tak, mamy do czynienia z dzielnikami (s&t) * (b+r/(s&t)) = n.
Powtarzamy dla dowolnej wariacji ustawienia pola tablicy z wyjątkiem wartości skrajnych, lub do znalezienia dzielników.
Ostatecznie dzielimy najpierw liczbę n przez s=(2^i-1) dla pewnego i, a następnie zerując odpowiednie bity szablonu dzielimy resztę przez szablon. Reszta jest liczbą nie przekraczającą połowy n, zaś liczby zadane przez schemat mają ściśle określoną długość binarną. Można też wykorzystać reszty modulo schemat.
Można się jeszcze pokusić o to, by liczby zadane szablonem nie przekroczyły pierwiastka kwadratowego z n. Zmniejszy to krotność obliczeń.
Pierwsze z dzieleń jest tak specyficzne, że można wykorzystać przyspieszone dzielenie przez odpowiedniki dziewiątki. Polega ono na pocięciu liczby na paczki długości i cyfr, na wynik składają się zapisane kolejno wartości sumy narastającej paczek.
Przykład: 1441 = 11 * 131
Dzielimy 1441 : 3 = 480 z resztą 1, nic dalej nie można zrobić
Dzielimy 1441 : 7 = 205 z resztą 6
schemat 101b każe sprawdzić dzielenie
(2*205+6) : 101b = 416 : 5 \equiv 1 (mod 5)
Dzielimy 1441 : 15 = 96 z resztą 1
Schemat generuje 3 dzielenia przez odpowiednio 1001b, 1011b oraz 1101b.
W pierwszym przypadku mamy dzielenie
((2+4)*96+1) : 1001b = 577 : 9 \equiv 1 (mod 9),
w drugim
(4*96+1) : 1011b = 385 : 11 \equiv 0 (mod 11)
oraz dzielniki 11 oraz 96+385/11 = 96+35 = 131; w trzecim
(2*96+1) : 1101b = 193 : 13 = 11 (mod 13) .
Następuje zamiana dzieleń:
1441 : 3 1441:3
1441 : 5 416 : 5
1441 : 7 1441:7
1441 : 9 577 : 9
1441 : 11 385 : 11
193 : 13
04 lipca 2013
Sprawdzanie między potęgami 2, czy można przyspieszyć
Ostatni algorytm, którego nie podałem w postaci informatycznej przeszukuje 'stożek' wartości mniejszych niż liczba faktoryzowana n, którego podstawą jest 2^i, oraz 2^(2i)<n<2^(2i+1).
Zastanawiało mnie, czy muszę przeszukiwać cały ten 'stożek'. Może istnieje 'ścieżka', która doprowadzi bezpośrednio do dzielnika.
Niech n = a*b+c, gdzie a<b jest potęgą 2, oraz c<a.
Stosowana konwersja podwaja a oraz odpowiednio zmniejsza b. Wartość c zachowuje się jak chce, nie przekraczając a. Dla wartości skrajnych tworzy ciąg niemalejący.
Do dalszej iteracji wybierałem ten z przedziałów [a, a+1], [a+1, a+2] w którym spełniony był warunek b1 < 2^k < b2, k<i.
W ten sposób dochodziłem do tego, że b po którejś iteracji stawało się potęgą 2.
Często znajdowałem pierwiastek, ale tylko wtedy, gdy był on bliski potęgi 2. Kiedy był dalej, należało przeliczyć reszty w krotności podwojonej ostatniej znalezionej reszty. Była to często wartość bliska podstawy 'stożka'.
W szczególności dla liczb pierwszych należało sprawdzać wszystkie wartości 'stożka'. Podobny obraz dawała także liczba złożona, w której dzielniki były najdalej 'odsunięte' od potęg 2.
Zatem ten sposób nie pozwala przyspieszać, należy przerachować wszystkie wartości parzyste 'grzebienia' nie większe niż 2^(i+1).
Dlatego 'grzebienia', gdyż sprawdzam wszystkie liczby postaci 2^k*d przedziału [2^i, 2^(i+1)), gdzie -1<k<i oraz d jest liczbą nieparzystą. Liczb takich jest tyle ile wynosi k+1, oraz 'tną' one pionowo 'stożek' na warstwy.
Zresztą, kto popatrzy na liczność przypadków przykładu poprzedniego posta, będzie wiedział, skąd nazwa 'grzebień'.
Zastanawiało mnie, czy muszę przeszukiwać cały ten 'stożek'. Może istnieje 'ścieżka', która doprowadzi bezpośrednio do dzielnika.
Niech n = a*b+c, gdzie a<b jest potęgą 2, oraz c<a.
Stosowana konwersja podwaja a oraz odpowiednio zmniejsza b. Wartość c zachowuje się jak chce, nie przekraczając a. Dla wartości skrajnych tworzy ciąg niemalejący.
Do dalszej iteracji wybierałem ten z przedziałów [a, a+1], [a+1, a+2] w którym spełniony był warunek b1 < 2^k < b2, k<i.
W ten sposób dochodziłem do tego, że b po którejś iteracji stawało się potęgą 2.
Często znajdowałem pierwiastek, ale tylko wtedy, gdy był on bliski potęgi 2. Kiedy był dalej, należało przeliczyć reszty w krotności podwojonej ostatniej znalezionej reszty. Była to często wartość bliska podstawy 'stożka'.
W szczególności dla liczb pierwszych należało sprawdzać wszystkie wartości 'stożka'. Podobny obraz dawała także liczba złożona, w której dzielniki były najdalej 'odsunięte' od potęg 2.
Zatem ten sposób nie pozwala przyspieszać, należy przerachować wszystkie wartości parzyste 'grzebienia' nie większe niż 2^(i+1).
Dlatego 'grzebienia', gdyż sprawdzam wszystkie liczby postaci 2^k*d przedziału [2^i, 2^(i+1)), gdzie -1<k<i oraz d jest liczbą nieparzystą. Liczb takich jest tyle ile wynosi k+1, oraz 'tną' one pionowo 'stożek' na warstwy.
Zresztą, kto popatrzy na liczność przypadków przykładu poprzedniego posta, będzie wiedział, skąd nazwa 'grzebień'.
Subskrybuj:
Posty (Atom)