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;   
};

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.



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.