Przyglądanie się zmianom wyglądu liczb zapisanych jako iloczyny małych wartości oraz liczby w systemie Fibonacciego spowodowały powstanie kolejnego algorytmu na poszukiwania dzielników.
Złożoność jest podobna jak przy przeskakiwaniu na kolejny system minimalizujący resztę, kiedy liczba jest dwucyfrowa. W odróżnieniu jednak od tamtego algorytmu można przeskakiwać część wartości, oraz sprawdzanych jest około połowa dopuszczalnych wartości - liczb od 1 do pierwiastka z n.
Założenie. Niech liczba nieparzysta n będzie przedstawiona wzorem
n = p*q+r,
gdzie p jest liczbą nieparzystą ograniczającą pierwiastek kwadratowy n z góry, analogicznie q jest liczbą bliską pierwiastkowi z n.
Algorytm jest następujący:
while (1<p) {
p-=2;
r+=2*q;
q+= r/p;
r = r%p;
if( 0=r ) return dzielniki(p,q);
}
return liczba pierwsza;
Wartość p maleje liniowo do 1, q rośnie skokami w zależności od części całkowitej ilorazu r/p. Ponieważ startujemy z wartości bliskich sobie, iloraz r/p początkowo jest równy 2 i bardzo powoli rośnie (niemonotonicznie). Dopiero kiedy sprawdzimy wszystkie dzielniki tej samej krotności cyfr dziesiętnych, iloraz ten przekroczy 20 (dzielenia rzędu 19/1). Reszta równa 0 oznacza, że mamy do czynienia z dzielnikami n = p*q. Przechodzimy tylko po wartościach p nieparzystych, gdyż dla parzystych nie ma dzielników przy nieparzystym n.
Dla bardzo dużych wartości można dostrzec dodatkową własność, która pozwoli pominąć część przypadków. Z reszty ubywa elementów według ciągu rozbieżnego, który lokalnie daje się wyznaczać.
Zatem sposób ten nadaje się jako konkurent rozkładu metodą dzielenia z resztą przez kolejne liczby nieparzyste (ta sama krotność przypadków), przy czym startuje od bardzo dużych wartości. Dodatkowo te ogromne są nieco prostsze w liczeniu.
Przykład numeryczny, 8934053 = 2989*2988+2921
p = 2989, q = 2988, r = 2921, (r+2q)/p = 2
na początku nowego przebiegu pętli mamy wartości
p = 2987, q = 2990, r = 2923, (r+2q)/p = 2
blisko pierwiastków
p = 1091, q = 8188, r = 945, (r+2q)/p = 15
p = 1089, q = 8203, r = 986, (r+2q)/p = 16
p = 1087, q = 8219, r = 0
rozkład 8934053 = 1087*8219.
W tym przypadku trzeba sprawdzać 952 przypadki zamiast 543 jak przy metodzie kolejnych dzieleń ze względu na położenie dzielników.
Wyjaśnienie (próba): liczba zapisana systemem Fibonacciego np. q=100100101
jest przekształcana do iloczynu (hybryda)
p*q = p00p00p0p.
Odjęcie 2 na każdej niezerowej pozycji powoduje zmniejszenie iloczynu o 2q, który jest przedstawiany jako nowa hybryda liczby Fibonacciego z taką samą wartością (p-2) na wartościach niezerowych oraz resztą r trzymaną oddzielnie.
Przechodząc z systemu Fibonacciego na dziesiątkowy powstaje ww algorytm.
30 stycznia 2014
07 stycznia 2014
Systemy niedziesiątkowe, wygląd liczby złożonej wskazującej dzielniki
Końcówka roku upłynęła na oglądaniu liczby złożonej przedstawianej w różnych systemach liczbowych, celem poszukiwania możliwości zmniejszenia krotności rozpatrywanych przypadków.
Kilka z algorytmów publikowanych na tym blogu korzysta z przedstawienia liczby w systemie o podstawach nieparzystych. Wtedy liczba złożona n=p*q w systemie o podstawie q ma cyfrę jedności równą 0.
Należy sprawdzać do połowy z pierwiastka kwadratowego z n przypadków.
Kiedy podstawy są parzyste, można oprócz konwersji na systemy p-2, p+2 łatwo przechodzić na systemy 2p, p/2. Miałem nadzieję, że manipulując tymi czteroma przekształceniami uda się dotrzeć do dzielników.
Przekształcenia działają, pojawiają się dodatkowe postacie błyskawicznie zwracające dzielnik, np:
n = a*(p-1)+a, tzn. naprzemienna suma cyfr dzieli p;
n = a*(p+1)+(p-a), tzn. suma cyfr dzieli p;
Wzory działają także dla liczb większych niż dwucyfrowe. Trudno je jednak znaleźć bez przeszukiwania.
Jeżeli mamy szczęście, to dodatkowo suma (naprzemienna) cyfr spełnia ten warunek także dla podstaw będących wielokrotnościami o potęgi całkowite z 2.
I na tym się kończy. Wiadomo jeszcze, że w pierwszej z tych postaci a jest nieparzyste (bo p-1 jest parzyste oraz wartość sumy ma być nieparzysta).
Myślałem, że dla kolejnych podstaw (p-1)*2^k wystarczy brać najmniejsze z możliwych a (nieparzyste), ale przeskakiwana jest wtedy szukana postać.
Drugie przypuszczenie: sumy i sumy naprzemienne majace w dzielnikach określoną postać w sąsiedztwie będą zbiegać do właściwych. Lecz nie do końca. Pojawiają się inne małe wartości dla sum bliskich szukanych z daleka od dzielników. Także przy samych dzielnikach wartości nie są monotoniczne.
Liczność przypadków może być nieco mniejsza, ale i tak wymagane jest sprawdzanie wielu przypadków (podstawy parzyste). Czasem widać dzielnik po przeliczeniu zaledwie połowy potrzebnych przypadków.
Nasuwa mi się możliwość zastosowania w tym przypadku specyficznego dodawania + mającego wiele argumentów. Jeden z argumentów jest dany relacją jednoargumentową '< r()', która zachowuje sumę tak długo, dopóki jest mniejsza od r, jeśli jest większa, blokuje dodawanie zwracajac wartość pustą NIL. Np. +(1, 2, 3, 4, 5, <7) zwraca NIL, bo 1+2+3+4>7. Z kolei +(1,2,3, <16, 4, 5) = 1+2+3+4+5=15<16 zwraca 15.
Kilka z algorytmów publikowanych na tym blogu korzysta z przedstawienia liczby w systemie o podstawach nieparzystych. Wtedy liczba złożona n=p*q w systemie o podstawie q ma cyfrę jedności równą 0.
Należy sprawdzać do połowy z pierwiastka kwadratowego z n przypadków.
Kiedy podstawy są parzyste, można oprócz konwersji na systemy p-2, p+2 łatwo przechodzić na systemy 2p, p/2. Miałem nadzieję, że manipulując tymi czteroma przekształceniami uda się dotrzeć do dzielników.
Przekształcenia działają, pojawiają się dodatkowe postacie błyskawicznie zwracające dzielnik, np:
n = a*(p-1)+a, tzn. naprzemienna suma cyfr dzieli p;
n = a*(p+1)+(p-a), tzn. suma cyfr dzieli p;
Wzory działają także dla liczb większych niż dwucyfrowe. Trudno je jednak znaleźć bez przeszukiwania.
Jeżeli mamy szczęście, to dodatkowo suma (naprzemienna) cyfr spełnia ten warunek także dla podstaw będących wielokrotnościami o potęgi całkowite z 2.
I na tym się kończy. Wiadomo jeszcze, że w pierwszej z tych postaci a jest nieparzyste (bo p-1 jest parzyste oraz wartość sumy ma być nieparzysta).
Myślałem, że dla kolejnych podstaw (p-1)*2^k wystarczy brać najmniejsze z możliwych a (nieparzyste), ale przeskakiwana jest wtedy szukana postać.
Drugie przypuszczenie: sumy i sumy naprzemienne majace w dzielnikach określoną postać w sąsiedztwie będą zbiegać do właściwych. Lecz nie do końca. Pojawiają się inne małe wartości dla sum bliskich szukanych z daleka od dzielników. Także przy samych dzielnikach wartości nie są monotoniczne.
Liczność przypadków może być nieco mniejsza, ale i tak wymagane jest sprawdzanie wielu przypadków (podstawy parzyste). Czasem widać dzielnik po przeliczeniu zaledwie połowy potrzebnych przypadków.
Nasuwa mi się możliwość zastosowania w tym przypadku specyficznego dodawania + mającego wiele argumentów. Jeden z argumentów jest dany relacją jednoargumentową '< r()', która zachowuje sumę tak długo, dopóki jest mniejsza od r, jeśli jest większa, blokuje dodawanie zwracajac wartość pustą NIL. Np. +(1, 2, 3, 4, 5, <7) zwraca NIL, bo 1+2+3+4>7. Z kolei +(1,2,3, <16, 4, 5) = 1+2+3+4+5=15<16 zwraca 15.
20 grudnia 2013
Dalsze badania
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ą.
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ą.
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.
Subskrybuj:
Posty (Atom)