Wymieniłem sposób dzielenia w faktoryzacji metodą kolejnych dzieleń.
Zastosowałem dzielenie, a które naprowadziły mnie systemy niedziesiątkowe.
Dla bardzo dużych liczb n, długości 2i lub 2i-1, dzielę n na dwie paczki o długościach nie przekraczających i:
n = a ## b = a*10^i+b
Cecha podzielności przez odpowiednio duże liczby jest postaci:
"liczba n jest podzielna przez p=10^(i+1)-r wtedy i tylko wtedy gdy wartość a*r+b jest podzielna przez p";
istnieje też podobna cecha, gdy p=10^i+r, cecha lepiej działa, gdy r nie jest większe niż 2*10^(i-1):
"liczba n jest podzielna przez p=10^i+r wtedy i tylko wtedy, gdy wyrażenie -a*r+b jest podzielne przez p".
Co się wtedy dzieje dla bardzo dużych liczb, typu RSA?
Otóż w zależności od przedziału p \in (10^i/(k+1) , 10^i/k) cecha podzielności może zmniejszyć wartości biorące udział w obliczeniach. Mianowicie przyrost r można zastąpić wartością 10^i % r, co powoduje zmniejszenie czynnika, czasem bardzo duże.
Przyjmujemy p nieparzyste. Jeśli p>5*10^(i-1), to cecha pozwala liczyć wyrażenie ar+b, zaś dla liczby p-2 już tylko a(r+2)+b. Oznacza to, że w tym przedziale nie trzeba wykonywać iloczynu a*(r+2), tylko do poprzedniego dodać 2a.
Gdy p \in (0,(3)*10^i , 0,5*10^i ), dla porachowania wartości w (p-2) należy do wcześniejszego iloczynu dodać 4a. I tak dalej, każde kolejne zmniejszenie poniżej 10^i/k powoduje konieczość dokładania 2k*a.
Dodatkowo, możemy jeszcze liczyć kongruencje modulo p iloczynów.
Nie trzeba także sprawdzać wartości 10^i/k. Pojawia się ona automatycznie, kiedy wartości r przekroczą p.
Dla tych największych wartości nie warto stosować drugiej z cech (suma naprzemienna) ze względu na pogorszenie sytuacji obniżającej wykładnik potęgi. Wtedy warto używać działań modulo dla iloczynów, choć można wygenerować też równanie kwadratowe, którego wartość jest dodawana.
Jednak, nie znalazłszy dzielnika dla tych największych wartości, można kontynuować, stosujac cechę podzielności w innym obszarze:
n = a ## b ## c = a*10^(2j)+b*10^j+c
"liczba n jest podzielna przez p = 10^(j+1)+r, r = (10^(j+1) % p) wtedy i tylko wtedy gdy (ar+b)r+c jest podzielne przez p"
"liczba n jest podzielna przez p = 10^j+r, r<10^j wtedy i tylko wtedy gdy -(-ar+b)r+c jest podzielne przez p"
To nie jest już tak wygodne jak wcześniejszy przypadek, ale nadal prosty.
Modyfikacja ta pozwala upraszczać od 1% do 90% przypadków (związane z wartością liczby n) w najgorszym obliczeniowo obszarze. Dodatkowe wykorzystanie podziału liczby n na trzy części pozwala sprawdzić około 90-99% przypadków.
Przykład numeryczny: n = 8934053 = 893*10^4 + 4053,
podzielność przez p=2501: r = 10^4 % 2501 = 2497
reszta n%p = (893*2497+4053)%p = (2233874)%p = (1430+1552)%p = 481
kolejne dodanie 6+2497 > 2500, następuje zmiana wielokrotności dodawanych a
podzielność przez p=2499: r = 10^4 % 2499 = 4
reszta n%p = (893*4+4053)%p = 7625%2499 = 128
podzielność przez p=2497: r = 10^4 % 2497 = 12 (1/4 k: dodać 2*4a=8a)
reszta n%p = (893*12+4053)%p = 2284 = (8*893+7625)%2497
podzielność przez 1001: r = 10^4 % 1001 = 991
reszta (893*991+4053)%1001, ale także
-(-8+934)+53 = -873 = 128 (1001)
14 marca 2014
Spacerkiem koło dzielników
Rozwijałem pomysł odwrócenia iloczynu
N = q1*q2
gdzie wszystkie wartości są traktowane schematem Hornera. Oznacza to dla wartości zapisanych trzycyfrowo
q1=(a2*p1+a1)*p0+a0,
dla dowolnych liczb p1, p2 i ciągach współczynników (a), (b) zapis:
N = ((a2*p1+a1)*p0+a0) * ((b2*p1+b1)*p0+b0) =
(( a2*b2*p1*p0 + p0*(a2*b1+a1*b2) + (a2*b0+a0*b2) )*p1 +
(a1*b1*p0 + (a1*b0+a0*b1) )*p0 + a0*b0 =
n*p1*p0 + w[a1,b1]*p0 + a0*b0
Jeśliby udało się jednoznacznie dopasowywać kolejne współczynniki ciągów 'cyfr' a, b, mamy do czynienia z algorytmem działającym na coraz mniejszych wartościach, w którym wyszukiwanie zostałoby przerzucone na znajdywanie odpowiednich par (a,b) stanowiących kolejne wyrazy.
Owszem, wielomian w[a,b] pojawiający się we wzorze jest łatwo wyznaczany, współczynniki znajdujemy z równań kongruencyjnych, których liczność zależy od wartości współczynnika danego poziomu p. Jest jednak duży kłopot z jednoznacznością.
Można temu częściowo zaradzić, biorąc za p niewielką liczbę pierwszą. Wtedy mamy do rozwiązania p-1 kongruencji, z których jedna jest sprzeczna, jak to pisałem w poprzednim poście.
Nauczyłem się także wskazywać, które z par (a,b) nie warto badać przy szukaniu dzielników. Zatem ciągi (a), (b) powinny z pomocą (p) skonwertować się na dzielniki. A wtedy pojawia się pułapka.
Oto heurystyka, która pozwala przybliżać orbity dzielników zapisanych jako a*P+A,
gdzie P jest iloczynem wartosci p0*p1*p2*..., A jest liczbą uzyskaną z obciętego ciągu (a0, a1, ..., am) wzorem A=a0+p0*(a1+p1*(...(am)...)).
Inicjujemy n=(N-1):2, A=1, B=1, P=2, w[0,0]=1.
Dane są: liczba n, wielomian w[a,b] = abP+aA+bB o współczynnikach P, A i B
Wybieramy liczbę pierwszą p, dbając, by p nie dzieliło n.
Dzielimy z resztą n = p*q+r, 0<r<p.
Szukamy rozwiązań równania kongruencyjnego jako par liczb (a,b)
w[a,b] = r (mod p)
po prostu wstawiając kolejno a=0, ..., a=p-1. Uzyskamy wartość b.
Tworzymy równanie zmiennych C=A, D=B
q = (A+aP)*C+(B+bP)*D
które traktujemy modulo min(A,B) uzyskując równanie rekurencyjne jednej zmiennej. Jeżeli rozwiązanie tego równania istnieje i jest (w kolejności) mniejsze niż p, równe 0, wielokrotnością p, w takim razie mamy dużą szansę, że para (a,b) stanowi lepsze przybliżenie dzielników.
Modyfikujemy w kolejności: A=A+aP, B=B+bP, P=p*P, n = (n-w[a,b]):p
Powtarzamy tak długo, dopóki n jest dodatnie.
Wartość n=0 oznacza, że udało nam się dokładnie wyznaczyć ciągi współczynników dzielników, czyli A, B są dzielnikami N.
A oto pułapka. Może się zdarzyć, że któryś ze współczynników A+aP, B=bP jest dzielnikiem oraz n>0, wtedy algorytm się nie zatrzymuje, lecz szuka dalej.
Przykład numeryczny, 8934053
n = (8934053-1):2 = 4467026
p = 3, w[a,b] = 2ab+a+b
rozwiązaniem równania w[a,b]=2 są pary (2,0), (0,2), zaś (1,1) daje równanie sprzeczne. Przyjmujemy a=(2,0) dzieki równaniu (n-2):3 = 5A+B
nowy wielomian w[a,b] = 6ab+a+5b
Dalej p=5, równanie dla a=4 jest 4b+4=3 (5)
Podaję w[a,b] oraz p
w[a,b] = 30ab+7a+29b, p=7, rozwiązanie (0,1) z D=4
w[a,b] = 210ab+37a+29b, p=3, rozwiązanie (0,2) z C=21=7*3
w[a,b] = 1470ab+457a+29b, p=7, rozwiązanie (0,6) z D=7
i tu następuje ominięcie, gdyż właściwym rozwiązaniem jest (6,1), wtedy współczynnik przy D jest dzielnikiem 1087 liczby 8934053 oraz następne n jest równe n=0.
Biorąc liczbę pierwszą większą od 11, też uzyskamy rozkład na dzielniki po tych sześciu iteracjach. A=8219, B=1087.
N = q1*q2
gdzie wszystkie wartości są traktowane schematem Hornera. Oznacza to dla wartości zapisanych trzycyfrowo
q1=(a2*p1+a1)*p0+a0,
dla dowolnych liczb p1, p2 i ciągach współczynników (a), (b) zapis:
N = ((a2*p1+a1)*p0+a0) * ((b2*p1+b1)*p0+b0) =
(( a2*b2*p1*p0 + p0*(a2*b1+a1*b2) + (a2*b0+a0*b2) )*p1 +
(a1*b1*p0 + (a1*b0+a0*b1) )*p0 + a0*b0 =
n*p1*p0 + w[a1,b1]*p0 + a0*b0
Jeśliby udało się jednoznacznie dopasowywać kolejne współczynniki ciągów 'cyfr' a, b, mamy do czynienia z algorytmem działającym na coraz mniejszych wartościach, w którym wyszukiwanie zostałoby przerzucone na znajdywanie odpowiednich par (a,b) stanowiących kolejne wyrazy.
Owszem, wielomian w[a,b] pojawiający się we wzorze jest łatwo wyznaczany, współczynniki znajdujemy z równań kongruencyjnych, których liczność zależy od wartości współczynnika danego poziomu p. Jest jednak duży kłopot z jednoznacznością.
Można temu częściowo zaradzić, biorąc za p niewielką liczbę pierwszą. Wtedy mamy do rozwiązania p-1 kongruencji, z których jedna jest sprzeczna, jak to pisałem w poprzednim poście.
Nauczyłem się także wskazywać, które z par (a,b) nie warto badać przy szukaniu dzielników. Zatem ciągi (a), (b) powinny z pomocą (p) skonwertować się na dzielniki. A wtedy pojawia się pułapka.
Oto heurystyka, która pozwala przybliżać orbity dzielników zapisanych jako a*P+A,
gdzie P jest iloczynem wartosci p0*p1*p2*..., A jest liczbą uzyskaną z obciętego ciągu (a0, a1, ..., am) wzorem A=a0+p0*(a1+p1*(...(am)...)).
Inicjujemy n=(N-1):2, A=1, B=1, P=2, w[0,0]=1.
Dane są: liczba n, wielomian w[a,b] = abP+aA+bB o współczynnikach P, A i B
Wybieramy liczbę pierwszą p, dbając, by p nie dzieliło n.
Dzielimy z resztą n = p*q+r, 0<r<p.
Szukamy rozwiązań równania kongruencyjnego jako par liczb (a,b)
w[a,b] = r (mod p)
po prostu wstawiając kolejno a=0, ..., a=p-1. Uzyskamy wartość b.
Tworzymy równanie zmiennych C=A, D=B
q = (A+aP)*C+(B+bP)*D
które traktujemy modulo min(A,B) uzyskując równanie rekurencyjne jednej zmiennej. Jeżeli rozwiązanie tego równania istnieje i jest (w kolejności) mniejsze niż p, równe 0, wielokrotnością p, w takim razie mamy dużą szansę, że para (a,b) stanowi lepsze przybliżenie dzielników.
Modyfikujemy w kolejności: A=A+aP, B=B+bP, P=p*P, n = (n-w[a,b]):p
Powtarzamy tak długo, dopóki n jest dodatnie.
Wartość n=0 oznacza, że udało nam się dokładnie wyznaczyć ciągi współczynników dzielników, czyli A, B są dzielnikami N.
A oto pułapka. Może się zdarzyć, że któryś ze współczynników A+aP, B=bP jest dzielnikiem oraz n>0, wtedy algorytm się nie zatrzymuje, lecz szuka dalej.
Przykład numeryczny, 8934053
n = (8934053-1):2 = 4467026
p = 3, w[a,b] = 2ab+a+b
rozwiązaniem równania w[a,b]=2 są pary (2,0), (0,2), zaś (1,1) daje równanie sprzeczne. Przyjmujemy a=(2,0) dzieki równaniu (n-2):3 = 5A+B
nowy wielomian w[a,b] = 6ab+a+5b
Dalej p=5, równanie dla a=4 jest 4b+4=3 (5)
Podaję w[a,b] oraz p
w[a,b] = 30ab+7a+29b, p=7, rozwiązanie (0,1) z D=4
w[a,b] = 210ab+37a+29b, p=3, rozwiązanie (0,2) z C=21=7*3
w[a,b] = 1470ab+457a+29b, p=7, rozwiązanie (0,6) z D=7
i tu następuje ominięcie, gdyż właściwym rozwiązaniem jest (6,1), wtedy współczynnik przy D jest dzielnikiem 1087 liczby 8934053 oraz następne n jest równe n=0.
Biorąc liczbę pierwszą większą od 11, też uzyskamy rozkład na dzielniki po tych sześciu iteracjach. A=8219, B=1087.
28 lutego 2014
Próba przyspieszenia znajdowania dzielników
Każdą liczbę naturalną można
jednoznacznie przedstawić jako sumę uporządkowanego ciągu iloczynów
współczynnika i potęgi podstawy systemu
\sum_{k>0} a_k*p^k
(budowa systemu pozycyjnego, dla p=10 dziesiątkowego).
Nic nie stoi na przeszkodzie, by potęgi podstawy systemu (p_k)_k zastąpić
losowym, lecz uporządkowanym ciągiem liczb, np. kolejnych liczb
pierwszych.
Wtedy liczbę naturalną N jednoznacznie przedstawiamy wyrażeniem, np.127 = ((1*7+1)*5+2)*3+1
za pomocą ciągów [1,1,2,1] oraz ciągu liczb pierwszych [7,5,3].
Rozpatrując to wyrażenie w ogólnej postaci, kiedy N jest liczbą złożoną, uzyskujemy dla liczb p i q pierwszych oraz szukanych współczynników a,b,c,d,e,f następującą postać:
N = ((ap+c)q+e) * ((bp+d)q+f) =
[[abpq + (ad+cb)q + (af+eb)]p +
(cdq+(cf+ed))]q + ef
Zauważmy, że w każdym bardziej
zewnętrznym nawiasie prostokątnym potrzeba coraz mniej
współczynnikow. Zatem znalezienie pary (e,f), z jej pomocą pary (c,d)
oraz ostatecznie pary (a,b) wyznaczy nam ciągi zapisu dzielników N.
Wyrażenie jest przedstawione jednoznacznie (własności działań po utożsamieniu wyrażeń przemiennych ab = ba).
Kiedy wylosujemy dowolną liczbę q, wyznaczamy parę (e,f). Wstawiamy
n_1 = (N- ef)/q
oraz kontynuujemy z mniejszą wartościa n_2 = n_1*p+(cdq+cf+ed) = n_1*p+w(), w której nieznane są tylko c i d. Z niej wyznaczamy parę (c,d)
kongruencjami
w() = (n mod p) modulo p.
Wielomiany w() dla kolejnych poziomów modyfikują się w bardzo specyficzny sposób. Są to równania kongruencyjne nieliniowe dla największego indeksu k, kiedy to mamy iloczyn obu współczynników oraz wszystkich liczb pierwszych ciągu. Do tego składnika dodawane są liniowe fragmenty (a_k*b_i + a_i*b_k) mnożone przez iloczyn liczb pierwszych na pozycjach nie przekraczajacych i<k. Współczynniki a_i oraz b_i zostały już wyznaczone wcześniej.
Jeśli liczba p jest pierwsza, znajdziemy p-1 par (c,d) takich, że c*d = n%p (mod p) dla kolejnych wartości c. Jedna z kongruencji jest w praktyce sprzeczna.
Przemienność może nam popsuć wybór, dlatego reszta z dzielenia przez pierwszą z ciągu liczb pierwszych nie powinna być kwadratem a^2. W przeciwnym przypadku mamy dwie możliwości: (a,a) oraz (-a,-a). Zapewni to jednoznaczność wyboru par (c,d) dla różnych c, z których w przypadku liczby złożonej dokładnie jedna staje się częścią ciągów dzielników N.
W ogólności tworzymy drzewo z krotnością gałęzi danego poziomu nie wiekszą niż wybrana liczba pierwsza. Drzewo to przyda się przycinać przez wybór odpowiedniej pary.
Aby zapobiec kłopotom z wyborem, warto uzywać liczb pierwszych jeszcze nie występujacych w ciągu. W przeciwnym przypadku kongruencje będą się upraszczać kosztem większych kłopotów z doborem pary. Tak działał jeden z wcześniejszych algorytmów faktoryzacji.
Zauważyłem pewne kryterium, działające dla początkowego wyboru par.
Wartość n modulo kwadrat wybranej liczby pierwszej jest najbliższe reszcie z dzielenia różnicy n-w() przez iloczyn początkowych liczb pierwszych ciagu z niedomiarem, lecz to kryterium w pewnym momencie staje się fałszywe. Nie udaje mi się znaleźć lepszego, wszystie w pewnym momencie zaczynają wskazywać inne pary niż właściwa.
Przykład numeryczny N = 8934053.
Ciąg podstaw w kolejności od cyfr
bardziej znaczących: [13, 11, 7, 5, 3].
Ponieważ N = 2 (mod 3), przyjmujemy
A=[2], B=[1], czyli para (a,b) = (2,1).
Do dalszego przebiegu używamy n =
(8934053-2*1)/3.
Pobieramy nową liczbę pierwszą q=5,
modyfikujemy wielomian
w() = 3ab+(a+2b)
i szukamy kolejnych par w() =
n%5 (mod 5). Wybieramy tę, która jest najlepsza, tu okazała się nią para (4,2). Wstawiamy wartości (a,b) = (4,2) wyznaczając wartość wielomianu w(), którą odejmujemy od n. Kontynuujemy.
Podsumowując cały rozkład.
Tablice dla dzielników są postaci:
A = 8219 = [7,
1, 1, 4, 2];
B = 1087 = [0, 10, 2, 2, 1].
Ostatnia para (a,b) = (7,0) (zapisana jako pierwsze elementy w tablicach ciągów) jako
jedyna wsród kandydatek byla równa n przez wstawienie do wielomianu tego poziomu
w() = 11*7*5*3ab + 105(10a+b) +
15(2a+b) + 3(2a+4b) + (a+2b)
Zakończyła ona proces znajdowania ciągów przez wyzerowanie n.
Liczby pierwsze nie wyzerują n, oraz dla kongruencji nie znajdziemy rozwiązań w małych liczbach naturalnych, ograniczonych liczbami pierwszymi.
Nie potrafię policzyć złożoności, czy proces jest szybszy niż standardowe dzielenia. Czy przycinanie drzew jest wystarczająco dobre.
Jest to po prostu kolejny sposób spojrzenia na faktoryzację.
30 stycznia 2014
Metoda 'kolejnych dzieleń' od końca
Przyglądanie się zmianom wyglądu liczb zapisanych jako iloczyny małych wartości oraz liczby w systemie Fibonacciego spowodowały powstanie kolejnego algorytmu na poszukiwania dzielników.
Złożoność jest podobna jak przy przeskakiwaniu na kolejny system minimalizujący resztę, kiedy liczba jest dwucyfrowa. W odróżnieniu jednak od tamtego algorytmu można przeskakiwać część wartości, oraz sprawdzanych jest około połowa dopuszczalnych wartości - liczb od 1 do pierwiastka z n.
Założenie. Niech liczba nieparzysta n będzie przedstawiona wzorem
n = p*q+r,
gdzie p jest liczbą nieparzystą ograniczającą pierwiastek kwadratowy n z góry, analogicznie q jest liczbą bliską pierwiastkowi z n.
Algorytm jest następujący:
while (1<p) {
p-=2;
r+=2*q;
q+= r/p;
r = r%p;
if( 0=r ) return dzielniki(p,q);
}
return liczba pierwsza;
Wartość p maleje liniowo do 1, q rośnie skokami w zależności od części całkowitej ilorazu r/p. Ponieważ startujemy z wartości bliskich sobie, iloraz r/p początkowo jest równy 2 i bardzo powoli rośnie (niemonotonicznie). Dopiero kiedy sprawdzimy wszystkie dzielniki tej samej krotności cyfr dziesiętnych, iloraz ten przekroczy 20 (dzielenia rzędu 19/1). Reszta równa 0 oznacza, że mamy do czynienia z dzielnikami n = p*q. Przechodzimy tylko po wartościach p nieparzystych, gdyż dla parzystych nie ma dzielników przy nieparzystym n.
Dla bardzo dużych wartości można dostrzec dodatkową własność, która pozwoli pominąć część przypadków. Z reszty ubywa elementów według ciągu rozbieżnego, który lokalnie daje się wyznaczać.
Zatem sposób ten nadaje się jako konkurent rozkładu metodą dzielenia z resztą przez kolejne liczby nieparzyste (ta sama krotność przypadków), przy czym startuje od bardzo dużych wartości. Dodatkowo te ogromne są nieco prostsze w liczeniu.
Przykład numeryczny, 8934053 = 2989*2988+2921
p = 2989, q = 2988, r = 2921, (r+2q)/p = 2
na początku nowego przebiegu pętli mamy wartości
p = 2987, q = 2990, r = 2923, (r+2q)/p = 2
blisko pierwiastków
p = 1091, q = 8188, r = 945, (r+2q)/p = 15
p = 1089, q = 8203, r = 986, (r+2q)/p = 16
p = 1087, q = 8219, r = 0
rozkład 8934053 = 1087*8219.
W tym przypadku trzeba sprawdzać 952 przypadki zamiast 543 jak przy metodzie kolejnych dzieleń ze względu na położenie dzielników.
Wyjaśnienie (próba): liczba zapisana systemem Fibonacciego np. q=100100101
jest przekształcana do iloczynu (hybryda)
p*q = p00p00p0p.
Odjęcie 2 na każdej niezerowej pozycji powoduje zmniejszenie iloczynu o 2q, który jest przedstawiany jako nowa hybryda liczby Fibonacciego z taką samą wartością (p-2) na wartościach niezerowych oraz resztą r trzymaną oddzielnie.
Przechodząc z systemu Fibonacciego na dziesiątkowy powstaje ww algorytm.
Złożoność jest podobna jak przy przeskakiwaniu na kolejny system minimalizujący resztę, kiedy liczba jest dwucyfrowa. W odróżnieniu jednak od tamtego algorytmu można przeskakiwać część wartości, oraz sprawdzanych jest około połowa dopuszczalnych wartości - liczb od 1 do pierwiastka z n.
Założenie. Niech liczba nieparzysta n będzie przedstawiona wzorem
n = p*q+r,
gdzie p jest liczbą nieparzystą ograniczającą pierwiastek kwadratowy n z góry, analogicznie q jest liczbą bliską pierwiastkowi z n.
Algorytm jest następujący:
while (1<p) {
p-=2;
r+=2*q;
q+= r/p;
r = r%p;
if( 0=r ) return dzielniki(p,q);
}
return liczba pierwsza;
Wartość p maleje liniowo do 1, q rośnie skokami w zależności od części całkowitej ilorazu r/p. Ponieważ startujemy z wartości bliskich sobie, iloraz r/p początkowo jest równy 2 i bardzo powoli rośnie (niemonotonicznie). Dopiero kiedy sprawdzimy wszystkie dzielniki tej samej krotności cyfr dziesiętnych, iloraz ten przekroczy 20 (dzielenia rzędu 19/1). Reszta równa 0 oznacza, że mamy do czynienia z dzielnikami n = p*q. Przechodzimy tylko po wartościach p nieparzystych, gdyż dla parzystych nie ma dzielników przy nieparzystym n.
Dla bardzo dużych wartości można dostrzec dodatkową własność, która pozwoli pominąć część przypadków. Z reszty ubywa elementów według ciągu rozbieżnego, który lokalnie daje się wyznaczać.
Zatem sposób ten nadaje się jako konkurent rozkładu metodą dzielenia z resztą przez kolejne liczby nieparzyste (ta sama krotność przypadków), przy czym startuje od bardzo dużych wartości. Dodatkowo te ogromne są nieco prostsze w liczeniu.
Przykład numeryczny, 8934053 = 2989*2988+2921
p = 2989, q = 2988, r = 2921, (r+2q)/p = 2
na początku nowego przebiegu pętli mamy wartości
p = 2987, q = 2990, r = 2923, (r+2q)/p = 2
blisko pierwiastków
p = 1091, q = 8188, r = 945, (r+2q)/p = 15
p = 1089, q = 8203, r = 986, (r+2q)/p = 16
p = 1087, q = 8219, r = 0
rozkład 8934053 = 1087*8219.
W tym przypadku trzeba sprawdzać 952 przypadki zamiast 543 jak przy metodzie kolejnych dzieleń ze względu na położenie dzielników.
Wyjaśnienie (próba): liczba zapisana systemem Fibonacciego np. q=100100101
jest przekształcana do iloczynu (hybryda)
p*q = p00p00p0p.
Odjęcie 2 na każdej niezerowej pozycji powoduje zmniejszenie iloczynu o 2q, który jest przedstawiany jako nowa hybryda liczby Fibonacciego z taką samą wartością (p-2) na wartościach niezerowych oraz resztą r trzymaną oddzielnie.
Przechodząc z systemu Fibonacciego na dziesiątkowy powstaje ww algorytm.
07 stycznia 2014
Systemy niedziesiątkowe, wygląd liczby złożonej wskazującej dzielniki
Końcówka roku upłynęła na oglądaniu liczby złożonej przedstawianej w różnych systemach liczbowych, celem poszukiwania możliwości zmniejszenia krotności rozpatrywanych przypadków.
Kilka z algorytmów publikowanych na tym blogu korzysta z przedstawienia liczby w systemie o podstawach nieparzystych. Wtedy liczba złożona n=p*q w systemie o podstawie q ma cyfrę jedności równą 0.
Należy sprawdzać do połowy z pierwiastka kwadratowego z n przypadków.
Kiedy podstawy są parzyste, można oprócz konwersji na systemy p-2, p+2 łatwo przechodzić na systemy 2p, p/2. Miałem nadzieję, że manipulując tymi czteroma przekształceniami uda się dotrzeć do dzielników.
Przekształcenia działają, pojawiają się dodatkowe postacie błyskawicznie zwracające dzielnik, np:
n = a*(p-1)+a, tzn. naprzemienna suma cyfr dzieli p;
n = a*(p+1)+(p-a), tzn. suma cyfr dzieli p;
Wzory działają także dla liczb większych niż dwucyfrowe. Trudno je jednak znaleźć bez przeszukiwania.
Jeżeli mamy szczęście, to dodatkowo suma (naprzemienna) cyfr spełnia ten warunek także dla podstaw będących wielokrotnościami o potęgi całkowite z 2.
I na tym się kończy. Wiadomo jeszcze, że w pierwszej z tych postaci a jest nieparzyste (bo p-1 jest parzyste oraz wartość sumy ma być nieparzysta).
Myślałem, że dla kolejnych podstaw (p-1)*2^k wystarczy brać najmniejsze z możliwych a (nieparzyste), ale przeskakiwana jest wtedy szukana postać.
Drugie przypuszczenie: sumy i sumy naprzemienne majace w dzielnikach określoną postać w sąsiedztwie będą zbiegać do właściwych. Lecz nie do końca. Pojawiają się inne małe wartości dla sum bliskich szukanych z daleka od dzielników. Także przy samych dzielnikach wartości nie są monotoniczne.
Liczność przypadków może być nieco mniejsza, ale i tak wymagane jest sprawdzanie wielu przypadków (podstawy parzyste). Czasem widać dzielnik po przeliczeniu zaledwie połowy potrzebnych przypadków.
Nasuwa mi się możliwość zastosowania w tym przypadku specyficznego dodawania + mającego wiele argumentów. Jeden z argumentów jest dany relacją jednoargumentową '< r()', która zachowuje sumę tak długo, dopóki jest mniejsza od r, jeśli jest większa, blokuje dodawanie zwracajac wartość pustą NIL. Np. +(1, 2, 3, 4, 5, <7) zwraca NIL, bo 1+2+3+4>7. Z kolei +(1,2,3, <16, 4, 5) = 1+2+3+4+5=15<16 zwraca 15.
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).
Subskrybuj:
Posty (Atom)