20 grudnia 2012

Algorytm faktoryzacji jeszcze szybszy

Jeden z algorytmów, które niedawno opracowałem, podczas przygotowywania do implementacji okazał się ciekawszy niż go szacowałem.

Algorytm polega na zapisie liczby faktoryzowanej n jako iloczynu dwu liczb z jak najmniejszą resztą c:
n = a*b+c
gdzie a, b są liczbami naturalnymi dodatnimi, a<=b oraz c jest mniejsze od mniejszej z liczb a, b.
Inicjacja algorytmu: n = 1*n+0.
Liczba ma dzielniki, kiedy c=0.

Przebieg algorytmu to bezpośrednie przedstawienie
n = a*b0+c0 = (a+1)*b1+c1,
najpierw wyznaczamy pomocniczą zmienną
x = b0+a-c0,
którą dzielimy z resztą r przez (a+1) otrzymując (q,r). Nowe wartości b1 oraz c1 uzyskujemy przez zwykłe odejmowanie
b1 = b0 - q;
c1 = a - r;

Przykładowo mając liczbę 155 = 1*155+0, otrzymamy w pierwszej iteracji
(155+1-0) : (1+1) = 156 : 2 = (78,0)
oraz nowe przedstawienie 2*(155-78)+(1-0) = 2*77+1.
Kolejna iteracja ma dzielenie
(77+2-1) : (2+1) = 78 : 3 = (26,0)
oraz dostajemy: 3*(77-26)+(2-0) = 3*51+2.
Kontynuując przekształcenia mamy
(51+3-2) : (3+1) = 52 : 4 = (13,0)
oraz dostajemy: 4*(51-13)+(3-0) = 4*38+3.
W następnej iteracji z dzielenia
(38+4-3) : (4+1) = 39 : 5 = (7,4)
dostajemy 5*(38-7)+(4-4) = 5*31+0, czyli dzielniki 5, 31.

Możemy kontynuować dalej:
(31+5-0) : (5+1) = 36:6 = (6,0) oraz 6*(31-6)+(5-0) = 6*25+5;
(25+6-5) : (6+1) = 26:7 = (3,5) oraz 7*(25-3)+(6-5) = 7*22+1;
(22+7-1) : (7+1) = 28:8 = (3,4) oraz 8*(22-3)+(7-4) = 8*19+3;
(19+8-3) : (8+1) = 24:9 = (2,6) oraz 9*(19-2)+(8-6) = 9*17+2;
(17+9-2) : (9+1) = 24:10 = (2,4) oraz 10*(17-2)+(9-4) = 10*15+5;
(15+10-5) : (10+1) = 20:11 = (1,9) oraz 11*(15-1)+(10-9) = 11*14+1;
(14+11-1): (11+1) = 24:12 = (2,0) oraz 12*(14-2)+(11-0) = 12*12+11
i na tym zakończyć.
Dzielenia są coraz prostsze, chociaż charakteryzują się sporym szumem.
Dla większych liczb wywołujemy algorytm ponownie dla większego z dzielników (31) mając pewność, że nie mamy dzielników mniejszych niż (5). Zatem w pierwszych iteracjach można zastosować szybkie konwersje a*(2b)+c = (2a)*b+c, a*(b+1)+c = a*b+(a+c).

Algorytm rachuje na coraz mniejszych wartościach, dążących do pierwiastka kwadratowego z rozkładanej liczby, i dlatego nie uważam go za najlepszy. Jest za to najprostszy. 

Kiedy oglądałem go szykując do implementacji, okazało się, że sum b+a-c nie trzeba liczyć!
Wystarczy skorzystać z dzielenia nie przeprowadzonego do końca, o wejściu (a,b0,c0), oraz uzyskać (b1,c1) znacznie mniejszym kosztem.
Polega to na tym, że możemy korzystać z warunku a< d = (a+1) oraz stosować metodę połowienia. Wynik w jest wtedy wartością zewnętrzną dzielenia.

Oto ten sposób: na wejściu dzielenia mamy a, b, c, dzielnik d, iloraz w=0, na wyjściu uzyskamy b1=w, c1=a+c jako niewykorzystaną resztę.
Przebieg algorytmu:
jeżeli b0 jest nieparzysta; c = c+a modulo d, b0-- i czasem modyfikacja w++
teraz b0 jest parzysta, dzielimy ją przez 2, zaś a mnożymy przez 2;
teraz a moze być a>d, bierzemy je wtedy modulo d oraz zwiększamy w o b0;
kończymy, gdy b0=1 lub a=d. Podczas uaktualniania współczynników pobieramy przechowywaną kopię a zwiększoną o 1.

Zastosowanie na tym samym przykładzie n = 155 = 1*155+0
a=1, b=155, c=0, d=a+1=2, iloraz w=0
b=155 nieparzysta, c = c+a = 1, oraz b=154
połowienie: a = 1*2 = 2, b=154:2 = 77
a modulo d jest 0, zatem w = b0 = 77 oraz c = 0+1 = 1
uaktualniamy współczynniki a=2, b=77, c=1, mamy postać 2*77+1!


Kolejna iteracja:
a=2, b=77, c=1, d=3, w=0
b=77 nieparzysta, c = c+a = 3 = 0 (3), b=76, w=1
połowienie: a = 2*2 = 4 = 1 (3), b = 76:2 = 38, w = 1+38 = 39
b=38 parzysta, dalej
połowienie: a = 1*2 = 2, b = 38:2 = 19
b=19 nieparzysta, c = c+a = 2, b=18
połowienie: a = 2*2 = 4 = 1 (3), b = 18:2 = 9, w = 39+9 = 48
b=9 nieparzysta, c = c+a = 3 = 0 (3), b=8, w = 48+1 = 49
połowienie: a = 1*2 = 2, b = 8:2 = 4
b=4 parzysta
połowienie: a = 2*2 = 4 = 1 (3),  b t= 4:2 = 2, w = 49+2 = 51
b=2 parzysta
połowienie: a = 1*2 = 2, b = 2:2 = 1
b=1 kończy iterację, uaktualniamy współczynniki a=3, b=51, c=2+0=2,
postać 3*51+2.

W tym przypadku dzielenie jest proste, zazwyczaj mamy tyle przekształceń, ile wynosi logarytm binarny z b. Złożoność tego dzielenia jest logarytmiczna.
Zatem algorytm ma złożoność O(n log n) arytmetycznie oraz
O(n^2 log n) pamięciowo.






08 grudnia 2012

Dzielenie przez sumowanie cyfr

Na stronie spryciarze.pl opublikowany został sposób dzielenia przez 9 wykorzystujący dodawanie narastające cyfr od najbardziej znaczących do cyfry jedności.

Sprawdziłem ten sposób, polega on na tym, że 10 = 1*9+1, oraz odpowiednio grupując ze sobą wyrazy rozwinięcia liczby
an*10^n + ... + a1*10 + a0
uzyska się postać
9*(an*10^{n-1} + ... + (suma _{k>i} ak)*10^i +... + (an+...+a0)) + (an+...+a0)
Z tej postaci można już oszacować iloraz oraz uzyskać resztę (an+...+a0) mod 9.

Sposób można łatwo uogólnić na inne wartości dzielnika. Jeżeli p jest podstawą systemu, oraz dzielimy przez liczbę p-r, to iloraz powstaje przez
sumę narastającą
cn = an
c{n-1} = r*an+a{n-1} = r*cn+a{n-1}
c{n-2} = r*(r*an+a{n-1})+a{n-2} = r*c{n-1} + a{n-2}
...
c0 = r*c1+a0
Najprostsze obliczenia są dla małych r.

Przykład: 3726 : 8 
p=10, r=2, a = 3 7 2 6
c3 = 3
c2 = 2*3+7 = 13
c1 = 2*13+2 = 28
c0 = 2*28+6 = 2*3*8 + 2*4+6 = 7*8+6
Interpretując teraz te liczby jako współczynniki rozwinięcia przy podstawie p liczby (cn...c1), przenosimy nadmiary do liczb bardziej znaczących, zaś po dodaniu jeszcze cześci całkowitej z c0:8 uzyskamy
3*10^2 + 13*10 + 28 + 7 = 465
Ostatecznie 3726 : 8 = 465 reszta 6, co jest prawdziwe.

Zwróćmy uwagę, że obliczenia przy c0 są postaci (d*e+f):g, przy których nie musimy liczyć wartości d*e+f, lecz można wykorzystać dzielenie chłopów rosyjskich aby wyłączać g. W ten sposób szybko zmniejszymy wartości i będziemy liczyć na stosunkowo małych liczbach.

A jak podzielić przez większą liczbę, np. 87?
Dokładnie tak samo, lecz podstawą będzie większa wartość, tu 100.
Przykład: 82045 : 87
p=100, r=13, a = 8 20 45
c2 =8
c1 = 13*8+20 = 124
c0 = 13*124+45 =19*87+4
Iloraz (8+1)*100+(24+19) = 943,
82045 : 87 = 943 reszta 4

To może inna liczba, np. podzielimy przez 11.
Teraz lepiej jest zastosować inne przekształcenie:
ci = -r*c{i+1}+ai
Uzyskana wartość może być ujemna, ale przekształcenia nie zmieniają się.
Przykład: 1751 : 11
p=10, r=1 (dokładniej -1), a = 1 7 5 1
c3 = 1
c2 = -1*1+7 = 6
c1 = -1*6+5 = -1
c0 = -1*(-1)+1 = 2
Iloraz 1*10^2 + 6*10 - 1 + 0 (z reszty) = 159
1751 : 11 = 159 reszta 2

W tym przypadku pojawiają się mniejsze wartości, które mogą oscylować wokół zera jako ciąg rozbieżny. W przypadku reszty dodatniej za pomocą pożyczek doprowadzamy ją do liczby dodatniej - cyfry. Dla reszty ujemnej algorytm tylko szacuje wartość ilorazu - wymaga dopracowania.

Pierwsze z przekształceń lepiej pasuje do liczb większych niż p/2, drugie dla liczb mniejszych od p-p/2. 
Pozostają przypadki dzielenia przez 3, 4. Dla nich najlepiej jest najpierw pomnożyć przez 3, 2 odpowiednio, po czym podzielić przez 9, 8 odpowiednio. Uzyskane reszty są odpowiednimi wielokrotnościami, potrójną, poczwórną właściwych reszt.
Zaś dzielenie przez 5 w systemie dziesiątkowym można zastąpić mnożeniem przez 2 i odcięciem cyfry jedności.  

Sposób ten dla części wielkich liczb jest równoważny z moim dzieleniem przez zmiany systemów (dla ilorazu dwucyfrowego są nawet te same przekształcenia). Lecz w ogólności przy liczbach mniejszych od podstawy sprowadza się do kolejnego dzielenia w c0, choć czasem przez mniejszą wartość. Dużo zależy od odległości dzielnika od podstawy. 
Nie znając szybkich sposobów konwersji, godny polecenia dla dzielników bliskich podstawy.    

27 listopada 2012

Projektowanie prototypu, wstęp

Aby sprawdzić siłę wymyślanych algorytmów, należy je przetestować praktycznie. 

Próbuję ożenić podejście funkcyjne z obiektowym. Muszą występować silne kompromisy.
Przykładowo obiekty odpowiadają relacjom bez argumentów, i są wskaźnikami do właściwej struktury. Jako takie mogą pojawiać się i znikać w dowolnym momencie.
Od strony podejścia funkcyjnego wprowadzane są liczne efekty uboczne oraz stan charakteryzowany licznością i typem argumentów. Eliminowane są też wycieki pamięci - właśnie użyty argument zostaje skasowany, co zmusza do specyficznej gospodarki pamięci 'z pamięcią'.

Kartoteka trzymająca dane została odpowiednio przygotowana, mam nadzieję. Wymaga jeszcze co najmniej dwu procedur: miotły do sprzątania pustych obszarów oraz podsumowania zawartości.
Po zbudowaniu prototypu kartoteki zacząłem testować moduł dodawania dwu liczb.

Moduł dodawania bezbłędnie przetrzymał dostawę 0, 1, 2, 3 argumentów  na wejściu przy wartościach sięgających 2^{65} z przekroczeniem zakresu slowa maszynowego (slowo maszynowe mam wielkości 2^{32}).
Nie tylko przetrzymał, ale i w odpowiednich przypadkach wyznaczył prawidłową wartość sumy.
Koszt: zostało użytych kilkanaście zmiennych, na jedną strukturę wymagane są co najmniej trzy plus kilka wspomagających.

05 listopada 2012

Pierwiastek kwadratowy przy systemach niedziesiątkowych

W ostatnim poście pokazałem sposób faktoryzacji dużych liczb. Wejściową daną była poostać bliska pierwiastkowi kwadratowego z rozkładanej liczby n, czyli
    n = a b r _ p
gdzie a=1, b=0 lub b=1, r jest resztą z dzielenia n przez p, zaś p jest podstawą systemu. Parametry r,  p, b można liczyć klasycznie, ale istnieje też sposób mający złożoność arytmetyczną logarytmiczną.

Będziemy  wykorzystywać binarną postać liczby n.
Jako pierwszy krok potrzebujemy podział liczby n na dwie części, mniej więcej na połowy, które także oznaczam przez b, r. Dalsze modyfikacje sprowadzą te wartości do odpowiednich parametrów. Oto warunki podziału: podstawa jest potęgą 2 nieco większą niż połowa n:
  p = 2^k,
  n < p^2,
  b < p, r < p,
  b*p+r = n.
Warunki te zapewnia przyjęcie k jako połowy logarytmu binarnego z n,
  b = (n >> k) ;
  r = n % (2^k); 

Mając  już tę postać, posługujemy się konwersjami
    d e _ {p+f} = d {e+d*f} _ p
dla f będącymi potęgami 2. Oczywiście liczby e+d*f często są większe niż podstawa systemu p, zatem nadmiary są przenoszone do cyfry d.

Konwersje te stosujemy tylko wtedy, gdy
  b+2^{k+1} < p.
Chodzi o to, by przybliżać się do pierwiastka z nadmiarem - będziemy mieli wtedy ciągle do czynienia z liczbą dwucyfrową b r _ p. 
Algorytm to prosta pętla po zmniejszającym się k. Kiedy wspomniany warunek jest spełniony, stosowana jest konwersja. Krotność wykonania konwersji jest związana z krotnością zer w binarnym przedstawieniu pierwiastka. 

Zmniejszając iteracyjnie k = (k>>1), parametry b oraz p zbiegają do wartości pierwiastka kawadtowego, kiedy k=0, co najmniej jeden z nich osiąga wartość pierwiastka. Kiedy są różne, jest to mniejszy z nich.
Na zakończenie stosujemy konwersję na system z podstawą o 1 mniejszą, co doprowadza do potrzebnej trójcyfrowej postaci wejściowej.


Przykład liczbowy: szukanie pierwiastka 169 747 007 < 2^{28}.
Tutaj za k zamiast wykładnika 28/2 = 14 przyjęta jest cała potęga 2^{14} = 16384. W kolumnach podane są wartości b, r, p, k, spełnienie warunku b+2k < p oznaczane jest przez p-=k, co oznacza zastosowanie konwersji.

Przygotowanie: p = k = 16384, b = n*2^{-14} = 10360, r = n%k = 8767
b          r        p          k        konwersja o
10360  8767  16384  16384
10360  8767  16384  8192
10360  8767  16384  4096
10360  8767  16384  2048  p-=2048
11840  8767  14336  1024  p-=1024
12751  5695  13312  512
12751  5695  13312  256    p-=256
13001  5951  13056  128
13001  5951  13056  64
13001  5951  13056  32
13001  5951  13056  16      p-=16
13017  5327  13040  8        p-=8
13025  5207  13032  4
13025  5207  13032  2       p-=2
13027  5197  13030  1        p-=1
13028  5195  13029  0
Wiemy zatem, że pierwiastek całkowity to 13028, zakończenie to konwersja na system o tej podstawie (wystarczy wymienic b z p). Uzyskamy postać:
1  1  5195 _ {13028}
Wynika ona z przeniesienia nadmiaru: 13029 = 13028 +1 = p+1.

26 października 2012

Faktoryzacja spradzająca duże i małe dzielniki

-->
Faktoryzacja liczb sprawdzająca prawie równocześnie największe i najmniejsze dzielniki jest możliwa do przedstawienia, chociaż algorytm jest bardziej skomplikowany. Przedstawię szkic za pomocą pseudokodu.
Ponieważ jednak podejście jest niestandardowe, najpierw funkcje biorące udział w algorytmach.

Liczba jest listą swoich cyfr, zaś cyfra to liczba z przedziału [0,p), gdzie p jest podstawą systemu liczbowego.

Last lista
zwraca ostatni element listy, w przypadku listy będącej liczbą, jej najmniej znaczącą cyfrę.
First lista
podobnie jak last, lecz zwraca pierwszy element listy, w przypadku liczb cyfrę najbardziej znaczącą.
Succ element
zwraca element poprzedzający dany iterator na liście, w braku takiego zwraca 0.
Move element lista
przenosi wskazany Element dołączając go do innej listy na koniec.
Rep [warunek] instrukcja
oznacza prostą pętlę powtarzającą się tak długo, dopóki warunek jest prawdziwy.
Fi instrukcje warunek
oznacza funkcję, która najpierw liczy instrukcje, a potem w zależności od warunku albo przywraca stare wartości, albo zostawia nowe, jest bardzo podobna do funkcji If, która robi to samo, ale w innej kolejności.
Konwersja [parametr]
tu oznacza konwersję zmniejszającą podstawę systemu o 2, chyba, że podany jest parametr. Wtedy zwiększa podstawę o ten parametr.
Corr [podstawa]
jest funkcją działającą na liczbach (listach) sprawdzającą, czy liczba jest poprawnie przedstawiona, w szczególności testuje warunki na cyfry, czy są liczbami całkowitymi w odpowiednim przedziale. Jeśli nie, nadmiary przenosi do cyfr bardziej znaczących, części ułamkowe jako liczby całkowite do cyfr mniej znaczących. Działa rekursywnie kończąc działanie, kiedy liczba jest we właściwej postaci.
Foreach lista : instrukcja
dla każdego elementu listy wykonuje podaną instrukcję.

Algorytm konwersji Konwersja [r=2]
Wejście: Liczba n pusta, Liczba m przy podstawie p+r, podstawa p
Rep [ m niepuste ] {
Move m.First n
Foreach n : element += r * Succ element
n.Corr [ p ]
}
Wyjście: Liczba n przy podstawie p, jako skutek uboczny może zmieniać p o r.

Algorytm faktoryzacji korzystający tylko z tej konwersji ma wadę: po znalezieniu dzielników nie wiemy, czy są one pierwsze czy złożone. Należy uruchomić algorytm ponownie dla każdego z dzielników, lub zapamiętać sprawdzonych kandydatów na dzielniki.
Jeśli dzielnik x dzieli podstawę systemu, aby dzielić liczbę, potrzebuje być także dzielnikiem wyrazu wolnego. W sprawdzaniu małych dzielników ograniczymy się tylko do iteracji, kiedy x jest dzielnikiem podstawy systemu.

Potrzebujemy trzech dodatkowych wartości: dzielnika minimalnego xmin, maksymalnego xmax oraz aktualnego x, inicjowanych x = xmin = 2, xmax=3.
Zaś dla większych liczb, często łatwiej i szybciej jest testować nieco większe wartości niż najbliższego kandydata na dzielnik. Mogą się one układać w pewien wzorek, w którym zaczynamy sprawdzać od większego kandydata, po czym w regularny sposób zmniejszamy do aktualnej najmniejszej podejrzanej wartości.

Przystąpmy do samego algorytmu faktoryzacji:
1. przygotowanie postaci liczby najmniejszej z trzycyfrowych n = [a, b, c], znalezienie podstawy p
2. Konwersja do najbliższej podstawy p nieparzystej i sprawdzenie dzielników xmin, x, xmax
3. dopasowanie licznika iteracji d
4. Rep [ istnieją kandydaci na dzielniki ] {
5. if 0= (d--) sprawdzenie i ponowne przeliczenie x, xmin, xmax, d
5a przy znalezionym dzielniku jego wyłączenie i restart algorytmu
6. Konwersja, modyfikacja podstawy p -=2
7. Corr [p]
8. if 0=n.Last czyli wyraz wolny równy 0, znaleziony dzielnik p; iteracyjny restart dla liczby z usuniętym Last, oraz z tą samą wartością xmin dla podstawy p
}

Ad 1. Chodzi o doprowadzenie do postaci b*p+c, w której b,c<p oraz p, b są największe z możliwych. W tym celu można przyjąć p jako sufit z pierwiastka kwadratowego liczby. Mając już tę postać stosujemy Konwercja [1], dzięki czemu uzyskujemy najmniejszą z możliwych liczbę trzycyfrową z a równym 1, b równym 0 albo 1 przy danej podstawie p;
Ad 2. W tym kroku sprawdzamy warunek p<xmin. Gdy 2=xmin, sprawdzamy podzielność przez 2 oraz stosujemy Konwersja [1], aby nie szukać niepotrzebnie dzielników po wartościach parzystych. Znalezione dzielniki usuwamy, i powtarzamy algorytm ze zmodyfikowaną liczbą. Gdy xmin jest większe, dobieramy x, xmax oraz d jak w 8;
Ad 3. Mała liczba naturalna d wskazująca, ile iteracji pętli będziemy czekać, aż x będzie dzielnikiem podstawy p. Jeśli nie jest wyznaczona w kroku 2, przyjmuje wartość 1 dla liczb postaci 3k+2, 2 dla liczb postaci 3k+4, wartość d=0 wykrywa podzielność przez 3;
Ad 4. Pętla 4-8, z której istnieją wyjścia na dwa sposoby, dla liczb pierwszych xmin>p, dla złożonych następuje ponowne wywołanie algorytmu dla poszczególnych dzielników z tymi samymi ustawieniami xmin;
Ad 5. W tym kroku d jako 0 wskazuje, że x dzieli p. Sprawdzamy, czy x dzieli też n.Last, czyli cyfrę najmniej znaczącą. Jeśli tak, mamy dzielnik x, zaś cała liczba ulega zniekształceniu, tak że należy zrestartować algorytm. Nie można podzielić podstawy oraz wyrazu wolnego przez x, gdyż wtedy mogą 'uciec' inne dzielniki. Jeżeli x>xmin, zmniejszamy x-=2 oraz liczymy d z pomocą y:
y = p modulo x; d = ( y%2 ? (y+x)/2 : y/2 )
Gdy x=xmin wstawiamy xmin = xmax oraz sprawdzamy wartości d liczone ww wzorem dla kolejnych elementów ciągu (p modulo (x+2*k) ). Interesuje nas podciąg malejący z wyjątkiem ostatniego elementu, który jest większy niż przedostatni. Taki podciąg cechuje się budową zbliżoną do (24, 22, 15, 8, 1, 17).
Przyjmujemy za xmax = (x+2*k) ostatniego elementu, oraz x = xmax-2.
Ad 8. Liczba jest iloczynem k*p, lecz nie wiemy, czy dzielniki są pierwsze czy złożone. Uruchamiamy algorytm dla tych wartości z aktualnym ustawieniem xmin;

Przykład liczbowy 169 747 007 = 19 * 1087 * 8219, szkic
Liczb niepodzielna przez 2. Do algorytmu wchodzi postać 
[a, b, c] = [1, 1, 5195] przy podstawie p=13028.
Następuje konwersja, by podstawa p była nieparzysta p=13027 = 3k+1, liczba [1, 3, 5197] , dwie iteracje doprowadzą do wartości p podzielnej przez 3. Obliczamy d dla p modulo 5 uzyskując 1, dla 7 jest większe. Mamy zatem d=1, xmin = 3, x = 5, xmax=7.
W pierwszej iteracji liczba przyjmuje postać 
[1, 7, 5207] przy p=13025. 
Podstawa p jest podzielna przez 5, lecz najmniej znacząca cyfra 5207 = 5*1041+2. Niezerowa reszta wskazuje, że 5 nie jest dzielnikiem. Zmniejszamy x do 3, d=1.
W drugiej iteracji liczba 
[1, 11, 5225] przy p=13023, 
p jest podzielne przez 3, sprawdzamy podzielność przez 3 cyfry najmniej znaczącej 5225 = 3*1741+2.
Teraz xmin=7, dane p jest podzielne przez 9, lecz x=9 nie jest dzielnikiem, jest 5 iteracji d=5 do sprawdzenia 7, reszta z dzielenia przez 11 jest większa, zatem xmax=11, x=7.
Po kolejnych kilku iteracjach okazuje się, że liczba postaci 
[1, 103, 7843] przy p=12977 
ma podstawę oraz najmniej znaczącą cyfrę podzielną przez x = xmin = 19.
Następuje restart algorytmu z xmin=19. Teraz należy powtórzyć krok 1. uzyskując postać liczby 169 747 007 : 19 = 8 934 053 jako 
[1, 3, 2923] przy p=2987. 
Teraz należy na nowo wyliczyć d=2, zaś do następnego dzielnika 21 jest 13 iteracji.
Kolejną postacią liczby jest [1, 7, 2933] przy p=2985.
I tak sprawdzamy kolejne duże dzielniki mniejsze od 2985, zarazem większe od 19.
Kiedy p opadnie do 1215, wiemy już z xmin=135, że dzielniki są w przedziale [xmax=157, 1215) lub sprawdzane właśnie 135, do sprawdzenia 173 dzielą nas dwie iteracje, zaś kolejne kandydatki: 171, 169, 167 i dalej dzieli taka sama krotność iteracji.
Osiągamy wreszcie postać [7, 580, 986] przy p=1089, dla której znowu następuje oczekiwanie na sprawdzenie dzielnika 217 (xmin=175) przy p=1085, lecz kolejna iteracja przekształca liczbę do postaci 
[7, 610, 0] przy p=1087. 
Mamy kolejne dzielniki: 1087 oraz 7*1087+610.
Obliczając jednak postać trzycyfrową tych wartości okazuje się, że są one mniejsze niż xmin=175, dzielniki są zatem liczbami pierwszymi.
Rozkład dokonany. 

W algorytmie tym rachujemy na coraz mniejszych wartościach, występujące dzielenie jest przez stosunkowo bardzo małe liczby. 


03 października 2012

Logarytmy dyskretne - odczytywanie

W poprzednim poście z pierwszego października opisana została konstrukcja grafu, z której można odczytywać wartosci logarytmów dyskretnych. Czas na ćwiczenia praktyczne.

Porachujemy kilka logarytmów dyskretnych dla systemu Z_{17}. Został wybrany, gdyż dla tej liczby pierwszej pojawia się niejednoznaczność znajdowanych rozwiązań.
Dla przypomnienia podaję graf, którego lista jest wypisywana pionowo za numerem kolejnym wierzchołka. Wierzchołki listy są drzewami wysokości 1 o dwu liściach rozdzielonych '|'.
lista trywialna
0
lista długości 8:
1. 2 - 6 | 11
2. 4
3. 8 - 5 | 12
4. 16
5. 15 - 7 | 10
6. 13
7. 9 - 3 | 14
8. 1

Numer jest potrzebny, gdyż odczyt jest podobny jak w metodzie zliczania indeksów.
Rozwiązując równanie x^k = a (17) wyszukujemy pozycję a w grafie, tę pozycję listy dzielimy przez k. Jeśli nie jest to liczba całkowita dodajemy odpowiednią wielokrotność długości listy. Jeśli nie znajdujemy rozwiązania całkowitego, równanie nie ma rozwiązań. Wartość logarytmu odczytujemy z tego samego poziomu, co a.

Policzymy dla próby x^3 = 16 (17)
Szukamy 16 na listach. Znajdujemy na pozycji czwartej listy. Teraz znajdujemy położenie rozwiązania jako iloraz 4/3. Jest to liczba niewymierna, dodajemy długość listy (4+8)/3 = 12/3 = 4. Odczytujemy wartość: 16, oraz 16^3 = 16 (17).
Jest to proste, gdyż 16 = -1 (17) oraz (-1)^3 = (-1) = 16 (17).

Teraz coś trudniejszego, też z listy głównej:
x^5 = 9 (17)
Wartość 9 stoi na pozycji siódmej, zatem obliczamy pozycję wyniku (7+8)/5 = 3. Znajdujemy wartość 8, którą sprawdzamy: 8^5 = 9 (17). Zgadza się.

Kiedy wartość jest znaleziona w liściu, rozwiązanie też znajduje się na tym poziomie:
x^3 = 6 (17)
Wartość 6 jest na pozycji pierwszej listy, następuje tymczasowe 'podniesienie' do 2. Obliczamy pozycję wyniku (1+8)/3 = 3. Teraz 'schodzimy' poziom niżej znajdując dwie kandydatki: 5 oraz 12. Sprawdzamy je:
5^3 = 6 (17) oraz 12^3 = 11 (17).
Znaleźliśmy rozwiązanie 5, mając do sprawdzenia wybór jednej z dwu wartości.

Zwróćmy uwagę, że rachowane były wyższe potęgi niż druga. Dla tego sposobu wartość potęgi nie jest przeszkodą. W wikipedii większośc przykładów ograniczała się do potęgi drugiej, która przy tym podejściu jest często znajdowana 'z marszu'.


17 jest liczbą pierwszą. Sposób ten działa także dla liczb złożonych, lecz graf może być znacznie bardziej skomplikowany.
Przykładowo dla liczby 493 niektóre wartości mają do 8 pierwiastków, które trudno uporządkować liniowo. Dołączane do listy o długości 112 łańcuchy miały wysokość 1 lub 3. Dla 28 wartości przyłaczane łańcuchy długości 4 miały wspólny środkowy wierzchołek oraz maksimum lokowane na liście.
Wybieramy wtedy za jedną z list dowolny ze znalezionych łańcuchów długości 112, łańcuchy fragmentami są rozdzielne, lecz łączą się przy wartościach mających większą krotność pierwiastków. Przypomina mi to nieco budowę białek, kiedy poskręcane włókna wzajemnie się przeplatają. Pozostałe łańcuchy miały długości o 1 mniejsze niż dzielniki 493.
W przypadku doczepiania łańcuchów niektóre wartości musiały pojawiać się dwukrotnie, gdyż były sześcianami innych wartości, których łańcuchy nie obejmowały. Wartości te raz rozpoczynały łańcuch, za drugim razem sąsiadowały z korzeniem.

01 października 2012

Logarytmy dyskretne - graf

W sierpniu podałem, że logarytmy dyskretne można znajdować za pomocą odpowiednio skonstruowanego grafu. Teraz podaję budowę tego grafu.

Najpierw nieco oznaczeń grafu zawierającego elementy zbioru Z_p:
ziarno - dowolna liczba 1< a < p; 
łańcuch - ciąg elementów dany resztami z a^k z dzielenia przez p, gdzie a jest ziarnem, p pochodzi z Z_p, k jest naturalne; z małego twierdzenia Fermata ciąg ten jest okresowy;
długość łańcucha - krotność elementów łańcucha liczonych od ziarna a do 0 lub 1 włącznie, albo do ponownego pojawienia się wartości a (łańcuch nie zawiera elementów 0, 1);  
lista - szczególny przypadek łańcucha, którego wierzchołki są grafami; listę tworzą najdłuższe łańcuchy (zazwyczaj);
struktura - zbiór częściowo uporządkowany podzbioru liczb Z_p relacją wziętą z łańcuchów; tworzy las czyli zbiór klasycznych drzew.

Grafy będące wierzchołkami listy to maksimum dla łańcuchów. Są one często
drzewami ukorzenionymi, czyli wierzchołek listy jest korzeniem drzewa.
Dla liczb złożonych mogą pojawić się także zbiory uporządkowane relacją z łańcuchów w sposób nieliniowy (istnieją rozgałęzienia łańcuchów). Taki przypadek spotkałem gdy p jest iloczynem nieco większych liczb, np p=493 = 17*29.

Tworzenie list.
Mamy dostępny zbiór wszystkich liczb Z_p. Będziemy je kolejno przenosić do grafu.
Wybieramy ziarno np. a=2. Tworzymy z niego łańcuch (a^k mod p), przemieszczając odpowiednie liczby. Przy okazji znajdujemy długość tego łańcucha n.
Jeśli 2n = p-1 lub n=p-1, p jest liczbą pierwszą (hipoteza sprawdzona dla p<50 i dla kilku większych p pierwszych).
Gdy n jest mniejsze, dzielimy p/n z resztą, szacując wartość reszty. Pozwala to przybliżyć wygląd grafów listy.
Za listy przyjmujemy łańcuchy. Jeśli z innego ziarna wygenerujemy łańcuch o wspólnym wierzchołku z pierwszym, zatrzymujemy generowanie, chyba że długość pierwszego uważamy za nieodpowiednią. Nowy łańcuch będzie miał więcej wierzchołków wspólnych, i czasem okazuje się lepszym kandydatem na listę.
Np. dla Z_{11} ziarno a=2 wygeneruje listę długości 6: a = (2^k) = (2 4 8 5 10 1), ale ziarno a=3 wygeneruje łańcuch odpowiedniej dla liczb pierwszych  długości 5: a = (3 9 5 4 1). 

Teraz mamy listę, oraz zbiór nieużytych wartości z Z_p. Wartości te usiłujemy przenieść doczepiając do listy jako grafy skierowane.
Wybieramy kolejne ziarno a nie należące do struktury, z którego generujemy łańcuch. Łańcuch tniemy na pierwszej wartości należącej do listy, oraz doczepiamy do odpowiedniego wierzchołka. Ostatecznie wszystkie wartości Z_p prócz ewentualnie 0 powinny zostać użyte.
Dla liczb złożonych mamy kilka niezależnych list. Mają one znacznie bardziej złożoną, i jak na razie niesklasyfikowaną budowę. Staram się tworzyć jak najbardziej symetryczne grafy.

Dla liczb pierwszych mamy 3 podstawowe listy (hipoteza prawdziwa dla liczb pierwszych mniejszych niż 75):
n=p-1, wszystkie elementy grypy multiplikatywnej Z_p^* są uporządkowane liniowo. Najłatwiejszy przypadek znajdowania, jak w metodzie zliczania indeksów (index calculi algorithm);
2n=p-1, wszystkie elementy listy mają doczepiony dokładnie jeden łańcuch, czyli dwuwierzchołkowe drzewo, którego korzeń jest na liście;
2n=p-1, co drugi element listy ma doczepione dwa łańcuchy, czyli pojawia się trójwierzchołkowe drzewo o korzeniu na liście, wysokości 1.

Liczby złożone mają bardziej skomplikowany wygląd, czasem nawet zamiast drzew pojawiają się grafy uporządkowane.

W przkładach poszukamy kwadratów, czyli takich wartości, które powstały z pomnożenia innych przez siebie.

Przykłady grafów dla rozwiązywania logarytmów dyskretnych:
Z_{19}, lista długości 18:
(2 4 8 16 13 7 14 9 18 17 15 11 3 6 12 5 10 1)
lista trywialna (0)
kwadraty: liczby przy wierzchołkach parzystych listy a[i] = {4, 16, 7, 9, 17, 11, 6, 5, 1} są kwadratami liczb na pozycjach a[i/2] oraz oraz a[(i+18)/2], np. 4 jest kwadratem 2 oraz a[(2+18)/2] = a[9] = 17.
Oraz oczywiście 0. 

Z_{14},
lista długości 3 zawierająca drzewa o 2 wierzchołkach, zapisywane (2 . 10):
( (2 . 10) (4 . 12) (8 . 6 ) );
lista liniowa długości 6: (3 9 13 11 5 1);
lista liniowa długości 1 zawierająca dzielnik: (7) i lista trywialna: (0)
kwadraty: kwadraty są na pozycjach parzystych list oraz w korzeniach doczepionych drzew, zbiór kwadratów: {2, 4, 8} + {9, 11, 1} + {7} + {0}

Z_{17}, oprócz listy trywialnej (0) lista mająca doczepione po dwa łańcuchy oddzielone przecinkami np. (2 . 6,11) oznacza złączenie łańcuchów (6-2-) oraz (11-2-) w jeden wierzchołek listy '2':
( (2 . 6,11)  4  (8 . 5,12)  16  (15 . 10,7)  13  (9 . 3,14)  1 )
kwadratami jak wcześniej są parzyste wierzchołki listy oraz te wierzchołki, które są korzeniami, tu wszystkie wierzchołki listy {2, 4, 8, 16, 15, 13, 9, 1}.

Sposób odczytywania logarytmów dyskretnych jest w poście z 3 października tego roku. 


31 sierpnia 2012

Mnożenie przez zmianę systemów - mnożenie chłopów rosyjskich


Sprawdziłem kilka algorytmów dotyczących mnożenia, korzystając z postaci uzyskanej przez zmianę systemów. Prawie wszystkie okazują się nieefektywne, mając duży narzut obliczeniowy. Dzielenie jest znacznie prostsze.

Powrotna wersja przekształceń liczby "a 0_p" do systemu dziesiątkowego wymaga podczas konwersji mnożenia przez stosunkowo duże liczby.

Wywoływanie iteracyjne konwersji zmieniających podstawę o 2^k dla kolejnych k jest szybkie komputerowo, lecz ma tyle wywołań, ile jedynek zawiera komputerowy zapis mniejszego z czynników. 

Pozostaje jeszcze jeden algorytm, bazujący na konwersjach
2a * p + b = a * (2p) + b
(2a+1) * p + b = a * (2p) + (p+b).
Liczbie "a 0_p" wartość a jest połowiona, po jakimś czasie ze względu na przenoszenie jedynek z nieparzystych wartości a do cyfry mniej znaczącej cała zawartość a przeniesie się do cyfry jednostek. Wtedy mamy wartość dziesiątkową iloczynu zapisaną w cyfrze jednostek.
Nie jest to nic innego jak tzw. mnożenie chłopów rosyjskich, i to jeszcze przy słabym uporządkowaniu czynników.

Ale mnożenie chłopów rosyjskich można zastosować od razu do czynników uporządkowanych malejąco, bez bawienia się konwersjami.
Mamy tylko tyle przekształceń, ile wynosi log_2 z mniejszego z czynników.

Oto przykładowe zastosowanie algorytmu mnożenia egipskiego. Pomnożymy 5*129.
Mniejsza jest 5, zatem będziemy dzielić 5 przez 2, zaś drugi czynnik 129 będziemy mnożyć przez dwa oraz dodawać do wyniku, kiedy reszta będzie niezerowa. 
pierwszy czynnik z dzieleniem, drugi czynnik p mnożony przez 2, wynik w
5/2=2 r 1  p=129   w = 129
2/2=1 r 0  p=258
1/2=0 r 1  p=516   w = 129+516 = 645
koniec, uzyskujemy wynik 5*129 = 645.
Bez bawienia się konwersjami, że 129 ma postać 1004_5 w piątkowym, zaś iloczyn 5*129 przedstawiamy jako 10040_5, jak przy mnożeniu przez zmianę podstaw.


22 sierpnia 2012

Publikacje

Międzynarodowa wymiana listów owocuje powstawaniem abstraktu opartego na moim pomyśle. Artykuł sprawdził i praktycznie napisał od nowa Djilali Behloul z Algerii.
Nie możemy zastosować standardowych modeli pamięciowych - prowadzą do  użytkowania stałych zasobów pamięci przez cały czas działania.
Z kolei model arytmetyczny, szacowany przeze mnie na
O ( n log^2 n )
jest jeszcze mniej złożony
O ( sqrt n log^2 n )
co dla bardzo dużych liczb jest poniżej złożoności liniowej.

Jest też konkurs prowadzony przez firmę Securitum na artykuł o szeroko pojętym bezpieczeństwie sieci. Wysłałem tam swoje najnowsze informacje dotyczące faktoryzacji oraz logarytmów dyskretnych. 
Jeżeli znajdę las dopasowany do danego logarytmu dyskretnego, znalezienie dowolnej z wartości staje się równoważne przeszukaniu tegoż, co jest wykonywane bardzo szybko.Zastanawiam się, czy można przewidzieć typ lasu.
Bliższe informację zamieszczę po zakończeniu konkursu.

18 lipca 2012

Dzielenie przez przedstawienie w innym systemie

Kończę dopieszczanie algorytmu dzielenia przez przedstawienie w innym systemie. Uzyskujemy resztę, a jak zastosujemy jeszcze raz przekształcenia od końca, także iloraz.

Zasada jest następująca:
mamy podzielić liczby a/b , a<b;
znajdujemy długość binarną liczby b:  2^k <= b < 2^{k+1};
dzielimy a na paczki mające dokładnie k bitów począwszy od cyfr najmniej znaczących, jest to przedstawienie a w systemie o podstawie 2^k;
stosujemy algorytm konwersji na system o podstawie b, przyrost  r = b-2^k ; algorytm ten jest podany w innym poście. Polega na tym, że dla kolejnych ciągów cyfr najbardziej znaczących stosujemy przekształcenie nie zmieniające wartości liczby a*(b-r) + c= a*p + (c+ar), oraz 'naprawiamy' za pomocą pożyczek lub nadmiarów tak, aby liczba na każdej pozycji stała się cyfrą;
cyfrą najmniej znaczącą jest reszta, odcinamy ją;
aby uzyskać iloraz, stosujemy algorytm konwersji na system o podstawie 2^k, przyrost r = -(b-2^k) ;
sklejamy paczki k-bitowe a w jedną wartość, będącą  ilorazem.

Mnożenie możemy przygotować w podobny sposób, bliższe informacje w następnym poście: "Mnożenie przez zmianę systemów"
Zatem ten sposób dzielenia jest łatwiejszy i mniej pracochłonny niż mnożenie.

Zachowanie algorytmu.
Gdy liczba a >> b, mamy dużo cyfr, oraz algorytm konwersji jest kwadratowy ze względu na krotność paczek mających k bitów. Dla dużych liczb stosunek ten się zmniejsza, co prowadzi do 'zmniejszania się' liczności cyfr , to z kolei bardzo upraszcza i przyspiesza konwersję. Jest to mile widziane dla bardzo dużych liczb.
Przyrost jest liczbą nie przekraczającą 2^k. Przypadkiem najgorszym są liczby nieco poniżej 2^k-1, dla nich przyrost jest największy.

Algorytm konwersji można zmodyfikować stosując iteracyjnie przyrost 2^m dla ustawionego bitu b[m], oraz malejących m= k ..  0.


Przykład: dzielimy liczby (zapis szesnastkowy dla zwiększenia czytelności)
-->
0x9927E4 : 0x7C  ( 10 037 137 : 124) 
liczba 0x7C = 111 1100b ma długość 7, podstawą jest 2^6 = 64 
liczbę 0x9927E4 = 1001 1001 0010 0111 1110 0100b dzielimy na paczki oddzielane kropkami
1001 10.01 0010 . 0111 11.10 0100 : 1.11 1100 = 38 18 31 36_{64} : 1 60_{64}
gdyż 100110b = 38, 010010 = 18 itd. 
liczbę tę konwertujemy na system większy o 60, czyli iteracyjnie będziemy 'naprawiać' w systemie  o podstawie 124 wyrażenia: 38*64 +18 = 38*124 + (18-60*38) =  38*124 + (-2262) = 19*124+94, i tak dalej ((19*124)+94)*64+31 = (19*124 + (94-60*19))*124 + (31-60*94) = ... 
liczba po konwersji : 5 32 97 40_{124} , reszta 40
mamy iloraz zapisany w systemie o podstawie 124: 5 32 97_{124}, należy go przerzucić do systemu o podstawie 64. Konwersja zaczyna się od 5*124 + 32 = 5*64 + (32+60*5) = ... 
liczba po konwersji : 19 48 49_{64} = 0x13C31 (80 945) 
tu znowu rozpisujemy 19 = 010011, 48 = 110000, 49 = 110001, sklejamy i odczytujemy binarnie.  
Podsumowanie dzielenia:  10 037 137 : 124 = 80 945 r. 40 

30 maja 2012

Faktoryzacja w praktyce

Wyskrobalem program do faktoryzacji w C++.
Niestety, na listach dynamicznych, przez co jest bardzo spawalniany. Ale i tak, kiedy już się rozpędzi i przekroczy pierwiastek sześcienny z liczby, w ciągu każdej minuty sprawdza dzielniki z przedziału zawierającego ok. 9 000 000 liczb, minus narzuty spowodowane zarządzaniem pamięcią. Wszystko to na procesorze o zegarze 1400 MHz.
Narzuty te zmniejszają  krotność sprawdzanych wartości  o parę procent. Można je łatwo zmniejszyć.

Zatem, swoim algorytmem do sprawdzenia, czy liczba n-cyfrowa jest pierwsza, mój algorytm będzie działał nieco dłużej niż n/9 10^{-6} minut, pomijając początkowy etap małych dzielników. Wtedy ze względu na krotność przekształceń czas ten jest kilka razy dłuższy. Trwa to jednak stosunkowo krótko.

15 maja 2012

Szukanie najwiekszego wspolnego dzielnika

Pomysły ze zmianami systemów liczbowych skutkują, że dzielenie staje się łatwiejsze i szybsze niż niektóre algorytmy mnożenia. Wystarczy, by w ilorazie
iloraz = dzielna / dzielnik
zapisać dzielną w systemie o podstawie dzielnik. Cyfra najmniej znacząca to reszta z dzielenia, pozostałe to wartość ilorazu. Wystarczy tylko jeszcze przekonwertować.

Technika ta wspomaga także szukanie  największego wspólnego dzielnika.  Przy czym do znajdowania złożoności lgorytmu wygodniejszy jest system binarny.

Podzielmy dzielną na paczki długości o jeden mniejszej niż długość dzielnika, jest ich a = sufit (dzielna/ dzielnik). Tak przygotowana dzielna jest liczbą a-cyfrową w systemie o podstawie b=2^i, gdzie
2^i <= dzielnik < 2^{i+1}
Przechodząc po cyfrach binarnych dzielnika, jeśli na danej pozycji 2^j występuje 1, konwertujemy do systemu o podstawie b = b+2^j.  Złożoność  tego algorytmu jest  kwadratowa  ze względu na  liczność  cyfr  równą  a. Korzysta z odejmowania złożoności liniowej względem cyfr. W najgorszym przypadku mamy 1 na wszystkich p pozycjach dzielnika. Powtarzamy p razy.
Uzyskamy część całkowitą i resztę.
Jeżeli reszta jest równa 1, nwd jest równe 1, kiedy zaś reszta jest równa 0, znaleźliśmy nwd równe aktualnej podstawie systemu.
W innych przypadkach powtarzamy postępowanie dla dzielnika i reszty. Uzyskujemy ciąg wartości zmniejszających się co najmniej dwukrotnie. Zatem jest ich co najwyżej tyle to cyfr dzielnika p.

Podsumowując, złożoność liczb n<m
nwd ( n, m ) = O( n(ap)^2 )  = O( n(m/n log n)^2 )
dla dużych n nie przekracza złożoności kwadratowej O(m^2), którą się przyjmuje klasycznie, gdyż najlepsze algorytmy działają nieco szybciej.

Przykład,
nwd( 1517, 1073 )
zapisujemy 1517 = 101 1110 1101_2 w systemie o podstawie 1073 = 100 0011 0001_2
1517 = 1 . (1 1110 1101)_{1024} = 1 493_{1024}
konwersja łącznie  o  przyrost 11 0001_2  = 49, generuje liczbę 1 . 444_{1073}
powtarzamy dla 1073 oraz 444 = 1 1011 1100_2
1073 = 4 49_{256} = 2 305_{384} = 2 185_{444}
powtarzamy dla 444 oraz 185 = 1011 1001_2
444 = 3 60_{128} = 2 74_{185}
powtarzamy dla 185 oraz 74 = 100 1010_2
185 = 2 57_{64} = 2 37_{74}
powtarzamy dla 74 oraz 37 = 10 0101_2
74 = 2 10_{32} = 2 0_{37}
ostatnia cyfra 0, zatem dzielnikiem jest aktualna podstawa 37
nwd( 1517, 1073 ) = 37

26 kwietnia 2012

Faktoryzacja o zlożonosci sześciennej

Znam dokładną złożoność algorytmu faktoryzacji publikowanego na tym blogu z początku roku. Algorytm posługuje się konwersjami  na  kolejny system.

Jego złożoność  tworzą dla liczby faktoryzowanej n:
Pętla sqrt n systemów;
konwersje wymagające (log n)^2 kroków wymagające operacji odejmowania i dodawania;
operacje działań na liczbach rzędu  sqrt n .
Łącznie algorytm ma  złożoność:
O( n log^2 n )

Dla porównania, najlepszy algorytm faktoryzacji podany u Knutha ma złożoność
O( e^{(1+epsilon) (log n)^{1/3} (log log n)^{2/3} )


30 marca 2012

Konwersja o połowę podstawy

Bardzo miły sposób mnożenia, dzielenia występuje wtedy, kiedy oba wyrażenia zapiszemy przy podstawie mniejszego z nich. Jest to wtedy odpowiednik mnożenia przez '10'.

Kłopotliwy jest powrót do systemu początkowego. Ale wykorzystując równe liczby zapisane inaczej:
a_n p^n + ... + a_1 p + a_0 przy podstawie 2n
2^n a_n p^n + ... + 2 a_1 p + a_0 przy podstawie n
uzyskujemy kolejną, stosunkowo szybką metodę konwersji.

Przykładowo, liczba trójcyfrowa w systemie o podstawie 288 mająca następującą postać:
1.54.53
po konwersji do systemu o podstawie połowę mniejszej, 144 przyjmuje postać:
4.108.53
Przekształcenia: [1*4, 54*2(<144), 53] = [4, 108, 53]. Nie trzeba używać przeniesień, które występują, gdy iloczyn na miejscu cyfry będzie większy niż podstawa systemu.

Podobnie przy zwiększaniu podstawy systemu należy dzielić przez potęgi 2 w takiej potędze, jak daleko są one od cyfry najmniej znaczącej. Oczywiście cyfry po każdej zmianie systemu warto naprawiać.

Sposób ten oraz konwersje na sąsiedni system pozycyjny p -> p+1 opisane kilka postów wcześniej pozwalają błyskawicznie pomnożyć, i wrócić do systemu początkowego. Może to być wygodniejsze niż mnożenie pisemne.

Do pełnej implementacji pozostaje jeszcze sprowadzenie liczby do systemu, w którym będziemy mnożyć, dzielić. Możemy to zrobić w czasie O( log^2 ).
Jeśli chcemy zapisać liczbę a w systemie o podstawie b, najpierw wpuszczamy a między dwie sąsiednie potęgi
b^n < a < b^(n+1).
Następnie połowimy b zacieśniając przedział występowania b jak w wyszukiwaniu binarnym. Kiedy z b dojdziemy do 0, znajdziemy najbardziej znaczącą cyfrę a w systemie o podstawie b.
Zagęszczeń danej cyfry mamy tyle, jak długa jest liczba b. Powtarza się tyle razy, ile cyfr ma docelowa postać liczby a.

18 lutego 2012

mnożenie wysokiej precyzji

Mam porównać swój algorytm faktoryzacji liczb ze standardowym. Liczby, które zadaję są za małe.
Zatem należy poznać bibliotekę rachunków na dużych liczbach - gmp.

Nie byłbym sobą, gdybym wpierw nie przejrzał dokumentacji i zastanowił się nad algorytmami.
Oczywiście dla dodawania i odejmowania algorytmy nie są wspomniane. Są liniowe, jak to opisuje Knuth w Sztuce programowania rozdział 4.3.1.
Większość algorytmów używanych przez gmp jest też podana przez Knutha. Najprostszy ma złożoność O(NM), gdzie N i M są licznościami cyfr czynników. Lepsze mają złożoność związaną z licznością cyfr N wzorem O(N 2^(sqrt(2 lg N)) log N) .

Lecz Knuth nie zauważył jeszcze jednego podejścia, związanego z podawaną wielokrotnie przez niego równoważną postacią liczby.
Liczba binarna a_n a_{n-1}... a_0 może być przedstawiona następująco:
A = (...(a_n*2+a_{n-1})*2+...)*2 + a_0
Oznaczmy A_i nawias zawierający jako najmniej znaczącą cyfrę a_i.
Obliczając dwie takie liczby A oraz B uzyskujemy iloczyn
A*B = (A_1*2 + a_0) * (B_1*2 + b_0 ) = ((A_1*B_1)*2 + A_1*b_0 + a_0*B_1)*2 + a_0*b_0
Ponieważ mnożenie przez 2 jest równoważne przesunięciu, mamy 3 mnożenia przez cyfry 0, 1; do wyniku dodajemy 3 sumy cząstkowe oraz mamy wywołanie rekurencyjne dla liczb mających mniej cyfr.

Złożoność tego algorytmu to 3 mnożenia, 3 dodawania oraz 3 przesunięcia na każdą iterację, za każdym razem wyznaczamy 2 cyfry. Do wyznaczonych cyfr najmniej znaczących już nie musimy wracać. Iteracji jest tyle, by któraś z liczb A_i, B_i stała się jednocyfrowa.
Ponieważ w przypadku binarnym mnożenia mozemy zastąpić dodawaniami liczb o długości co najwyżej max(N,M), cała złożoność nie przekracza złożoności sum coraz krótszych liczb. Jest zatem co najwyzej kwadratowa dla wiekszej z nich. Dokladniej porownujac, uzyskuje zlozonosc taka sama jak u Knutha O(NM), lecz na coraz mniejszych cyfrach.

28 stycznia 2012

Sposób konwersji

W ostatnim poście pokazałem, jak można stosować cechę podzielności, ale nie wskazałem szybkiego sposobu przechodzenia z systemu na system. Za to użyłem takiego przy konwertowaniu z systemu dziesiątkowego na ósemkowy.

Zatem dodatek: szybki sposób konwersji systemów pozycyjnych.

Zanim go jednak przedstawię, oto różnica między liczbą a cyfrą.
CYFRA danego systemu liczbowego to  liczba mieszcząca się między 0 włącznie oraz podstawą systemu liczbowego p, czyli w przedziale [0, p). Stanowi zarazem podstawowy składnik konstrukcji liczb. Dla małych podstaw cyfry mają własne symbole graficzne.
LICZBA [w systemie pozycyjnym] to [uporządkowany] ciąg cyfr a_n .. a_1 a_0.

W algorytmie poniższym cyfry są traktowane jako liczby, oraz w czasie początkowych przekształceń ignorujemy nałożone nań ograniczenia na wielkość.

Mając liczbę a_n...a_0 w systemie p, chcemy ją zapisać w systemie o podstawie p+r (p+r>1), stosujemy oparty na schemacie Hornera zapisu liczb
(...(a_n*p+a_{n-1})...)*p+a_0
następujący algorytm:
inicjacja: sprawdzamy przyrost r, załóżmy ze jest r>0 Bierzemy a_n, pozostałe cyfry czekają w kolejce [a_{n-1} , ... , a_0]
1) pobieramy cyfrę a_j z kolejki;
2) każdą cyfrę a_i (zaczynając od i=j, zwiększając do i=n ) modyfikujemy do liczby a_i := a_i - r*a_{i+1}, kiedy nie ma wartości i+1, traktujemy ją jako 0; 
3) powstałe liczby sprowadzamy do cyfr dodając t razy (p+r) - pożyczamy od cyfr bardziej znaczących, aby liczba stojąca na miejscu cyfry była w odpowiednim przedziale;
4) poprzedzającą cyfrę a_{i+1} zmniejszamy o znalezione w 3) t;
powtarzamy 3) 4) tak długo, aż liczba będzie się składała tylko z prawidłowych w danym systemie cyfr;
5) kiedy kolejka cyfr jest pusta, uzyskujemy postać liczby w systemie (p+r), w przeciwnym przypadku wracamy do 1)
Na zakończenie modyfikujemy aktualną krotność cyfr, gdyż wartość ta może się zmienić. 

Kiedy r jest dodatnie, w 2) dodajemy zamiast odejmować, zaś liczbę sprowadzamy do cyfry przez przekazywanie nadmiarów (p+r) cyfrom bardziej znaczącym.

Sposób nie jest tak szybki jak konwersje między systemami o podstawach będących potęgami jednej liczby a, w których liczbę zapisujemy przy podstawie a, a następnie cyfra powstaje z paczki odpowiedniej długości, np. konwersja z systemu o podstawie 27 na system o podstawie 9 przez system o podstawie 3:
22 2 4_{27} = 211 002 011_3 = 2 11 00 20 11_3 = 24064_9,
ale i tak znacznie szybszy niż stosowanie definicji systemów.

Przykład: przekonwertujemy szesnastkowe 0xBACA na system dziesiętny.
cyfry początkowe: [ 0xB=11, 0xA=10, 0xC=12, 0xA=10] 
Początek to lista 11, mamy przyrost 10-16=-6, czyli do poszczególnych cyfr będziemy dodawać iloczyny 6
zaczynamy pierwszą iterację:
11  10+6*11
11  76
na pozycji jedności mamy liczbę! 76, którą należy zmniejszyć, aby na powrót stała się cyfrą systemu, tym razem dziesiątkowego
(11+t) (76-t*10)    t=7
18  6
teraz cyfra dziesiątek jest za duża. Przenosimy nadmiar tworząc następną cyfrę, cyfrę setek.
1  8  6
Dla sprawdzenia, 11*16+10 = 186 = 1*10^2 + 8*10 + 6,
Kolejna iteracja, dołączamy cyfrę 0xC=12
1          8            6           12
1[+0]   8+6*1   6+6*8   12+6*6
1          14         54          48
przenosimy nadmiary (od razu wynik)  (11*16+10)*16+12 = 2988
2          9           8            8
Ostatnia cyfra, czyli ostatnia iteracja
2          9           8            8            10
2[+0]   9+6*2   8+6*9   8+6*8   10+6*8
2          21         62          56          58
Przenosimy nadmiary
4          7           8           1              8
Zatem szesnastkowe 0xBACA to dziesiętne 47818.


W drugą stronę podobnie, lecz mamy 5 iteracji, przyrost 6. Pierwsza jest postaci
4[-0]  7-6*4
4       -17
naprawa polega na dodawaniu docelowej podstawy 16, czyli
4-2   -17+2*16
2      15
Kolejna iteracja
2      15         8
2      15-6*2  8-6*15
2      3           -82
1      13         14
Kolejna iteracja
1      13         14          1
1      13-6*1  14-6*13  1-6*14
1      7           -64         -83
1      2           10          13
Ostatnia iteracja
1      2           10          13           8
1      2-6*1    10-6*2   13-6*10  8-6*13
1      -4          -2           -47         -70
0      11         10          12          10
Ostatecznie przekształcamy te szesnastkowe cyfry na liczbę 0xBACA.


23 stycznia 2012

Podzielnosc w systemach niedziesiatkowych

Ułamki egipskie oraz algorytm faktoryzacji zwrócily moja uwagę na systemy niedziesiątkowe. Na razie system binarny oraz systemy o podstawie będące potęgą 2.

Występują bardzo proste cechy podzielności dla liczb będących postaci 2^n +- 1.

Dzieląc liczbę w systemie o podstawie 2^n na paczki binarne długości n-1 znaków, liczba jest podzielna przez 2^n-1 gdy suma wartości poszczególnych paczek binarnych też jest podzielna. (por. cecha podzielności przez 9).
Dla liczb postaci 2^n+1 reszta jest równa reszcie sumy naprzemiennej tych paczek (por. cecha podzielności przez 11.

Do odpowiednich liczb należą 3, 5, 7, 9, 15, 17, 31, 33 i wiele innych.

Wystarczy przekształcić liczbę do postaci np. oktagonalnej by już mieć postać równoważną binarnej i zastosować cechę.

Przykład: sprawdzimy resztę z dzielenia 3523432 przez 17 równe 10001 binarnie.
Przekształcamy 3523432 do systemu ósemkowego, dodając do każdej cyfry z wyjątkiem pierwszej przed '|' podwojoną wcześniejszą i 'naprawiając cyfry' przez przenoszenie nadmiarów 8:
3 5 | 2 3 4 3 2
4 3 2 | 3 4 3 2
5 4 0 3 | 4 3 2
6 7 0 3 4 | 3 2
1 0 4 6 4 2 3 | 2
1 2 6 0 1 2 7 2 |
1 5 3 4 1 5 5 0
Zatem w ósemkowym liczba wyglada 15341550 = binarnie 001 101 011 100 001 101 101 000 = 0011 0101 1100 0011 011 01000 = heksadecymalnie 3 5 12 3 6 8.
Według cechy podzielności, naprzemienna suma -3+5-12+3-6+8 = -5 = 12 (17) ma taką samą resztę. Istotne są tu znaki, ale rzeczywiście, wartość ta w dziesiatkowym przedstawi się jako 207260*17+12.

Algorytm dotyczący sposobu przechodzenia między systemami jest opublikowany w następnym poście.