Faktoryzacja z podstawy systemu

Wspomniałem 18 grudnia o nowym pomyśle na faktoryzację, związanym z przedstawianiem podstawy systemu w postaci kanonicznej.
Pomysł jest tak specyficzny, że działa także dla normalnych liczb, choć potrzeba pamiętania ogromnych wartości z ich dzielnikami powoduje, że sposób jest nieopłacalny, i nikt o nim nie wspomina.

Podobnie jest i teraz, choć wymagam tylko znajomości tablicy liczb pierwszych nie przekraczajacych pierwiastka czwartego stopnia z rozkładanej liczby. Nie mamy? Nic straconego, delikatna rekursywna modyfikacja i już taką uzyskujemy, niejako przy okazji.Ale teraz podaję jakby tablica ta była znana.

Zacznijmy od nieformalnego podziału liczb pierwszych w zależności od rozkładanej liczby N:
- małe, mniejsze niż pierwiastek czwartego stopnia z N;
- średnie, większe niż pierwiastek czwartego stopnia, ale pojawią się w co najmniej dwu warstwach badanych;
- duże, wystąpią raz, rozpoznajemy je po tym, że podstawa systemu jest albo taką liczbą pierwszą, albo jej podwojeniem.

W pętli przeszukiwania zupełnego duże liczby pierwsze są sprawdzane i ignorowane. Potrzebujemy struktury kopca trzymajacego strukturę (liczba, lista jej dzielników), sortowaną malejąco po liczbie. Wkładając do kopca liczbę z dzielnikiem, sprawdzamy jej obecność. Jeśli już jest, tylko dopisujemy kolejny dzielnik do listy. Kopiec osiąga rozmiar równy krotności wszystkich małych i średnich liczb pierwszych.
Zaczynamy od liczby n będącej pierwiastkiem kwadratowym z N, w kopcu umieszczamy wartości (n- (n%k), k), gdzie k są małymi liczbami pierwszymi.  Wiemy, że są one już dzielnikami podstawy systemu, trzeba jeszcze sprawdzić, czy są dzielnikami N.Dlatego specjalnie je zaznaczamy do sprawdzenia. Każdą liczbę pierwszą sprawdzamy dokładnie raz.

Mamy niezmiennik
N = (a+b)*a+c
Jeśli liczba k dzieli a oraz c, dzieli też N. 
W pętli algorytmu przechodzimy na kolejną warstwę operacjami:
a-=1;
b+=2;
c+=b-1;
stosując przemieszczenie nadmiaru (z liczby w systemie o podstawie a), by c było ograniczone przez a. Może być też ujemne. Wartosć b tylko zwiększamy, może przekroczyć a. Dzielnik znajdziemy jako wspólny dzielnik gcd(c,a). Wystarczy tylko warunek 0= c%a, zamiast szukać wspólnego dzielnika gcd(c,a). Jest to pierwszy krok szukania tegoż dzielnika, i właśnie tyle wystarczy przy założeniu znajomości małych liczb pierwszych. 

Algorytm przebiega następująco:
- pobieramy liczbę o 1 mniejszą, czyli n:=n-1
jeżeli z kopca, dzielimy ją przez jej dzielniki (z listy) by uzyskać jak najmniejszą wartość d, w przeciwnym przypadku generujemy;
uzyskujemy w ten sposób liczbę pierwszą d albo 1;
- dzielimy c%d i sprawdzamy reszty, tak samo z małymi liczbami pierwszymi z pobranej listy, reszta zero oznacza dzielnik, który zwracamy. Wyjście z algorytmu. Przy reszcie niezerowej te liczby pierwsze już nie muszą być sprawdzane, zaznaczają tylko kolejne wartości dla wygodniejszego szukania d;
- dla wszystkich uzyskanych średnich i małych liczb pierwszych k wprowadzamy do kopca pary (n-k, k);- przechodzimy na kolejną warstwę a-1;
- powtarzamy, aż osiągniemy połowę pierwiastka kwadratowego z N, zwracamy że N jest liczb pierwszą;

Załóżmy N=893453
pierwiastek kwadratowy sqrt N = 2988, którego kolejny pierwiastek to 54.
Bierzemy listę liczb pierwszych mniejszych niż 55 i wprowadzamy do kopca w postaci (c2 oznacza konieczność sprawdzenia podzielności 2 przez N):
(2988, [c2,c3])
(2987, [c29])
(2985, [c5])
(2983, [c19])
(2982, [c7])
(2981, [c11])
(2977, [c13])
(2976, [c31])
(2975, [c17])
(2968, [c53])
(2967, [c23, c43])
(2961, [c47])
(2960, [c37])
(2952, [c41])

Przystępujemy do algorytmu:
pobieramy n=a=2988, które ma dwa dzielniki: 2 i 3. Po podzieleniu przez nie n (czasem kilkakrotnie uzyskujemy nową średnią liczbę pierwszą d=83.
Niezmiennik dostarcza a=2988, b=1, c=2921
Sprawdzamy, czy N jest podzielne przez 2, 3, 83, dzieląc c przez te wartości. Zawsze występuje niezerowa reszta, zatem odkładamy do kopca nowe pozycje (2988-83, [83]) = (2905, [83]), (2985, [3]), (2986, [2])
Pozycja wyznaczona przez 2985 przyjmuje teraz wartość (2985, [3, c5])

Następna pozycja (2987, [29]), wskazuje kolejną liczbę pierwszą 103. Ale poznieważ c było duże, możemy przenieść 2987 z c jako jedynkę do b uzyskując
nową postać niezmiennika a=2897, b=4, c=-64
Sprawdzamy obydwie liczby 'dzieląc' 64%29, 64%103, odkładając na kopiec kolejne dwie wartości (2884, [103]) oraz (2958, [29]).

I tak początkowo co iterację dostajemy nową liczbę pierwszą, którą sprawdzamy. Ale wkrótce mamy wziąć 2978, na kopcu jest (2978, [2]), i mamy rozkład 2978 = 2*1489. Jest to duża liczba pierwsza w myśl podziału. Wartości (1489, [1489]) nie odkładamy do kopca.
Wartość (2976, [2,3,c31]) ma rozkład wykorzystujący tylko te cyfry, d=1. Z niezmiennika, c urosło do 101, sprawdzamy 101%31. Przy następnej takiej wartości 2970 nie ma już nic do sprawdzenia - wszystkie użyte dzielniki zostały sprawdzone.
Nieco wcześniej mamy 2971, której to wartości nie ma w kopcu. Jest to duża liczba pierwsza d=2971, która jest za duża, by być dzielnikiem. Niezmiennik ma wtedy znacznie mniejsze c=256. 
W ten sposób znajdujemy kolejne liczby pierwsze; zeszlibyśmy z a do połowy 2988 czyli 1494, gdyby nie dzielnik, który znajdziemy dla pozycji o niezmienniku a=2174, b=1935 lub b=1936, c=|1087|, d=1087 = |c|.

Teraz modyfikacje.
Jeśli nie znamy wszystkich małych liczb pierwszych, możemy wywołać procedurę dla pierwiastka rekursywnie. Jednak wtedy każda znaleziona wartość pierwsza (także duza) generuje nowy obiekt do kopca do sprawdzenia. Przy wywołaniu rekursywnym niezmiennik jest nieistotny, nie porównujemy ani nie sprawdzamy. Uzyskamy w ten sposób tablicę liczb pierwszych.
Może być, że w petli głównej wartość jest iloczynem mniejszych liczb pierwszych, dla średnich liczb pierwszych wkrótce napotkamy kolejną warstwę z tym dzielnikiem, który powienien być lepiej widoczny (bo np. zmieni się parzystość). Ale dla dużych liczb pierwszych obawiam się przeoczenia tego dzielnika.
Na szczęście kierunek przechodzenia warstw można odwrócić.Proponowanym punktem startowym jest połowa pierwiastka z N, czyli w tym przykładzie a=1494, b=3a=4482 (choć lepiej 4486, gdyż wtedy c=-67).

Szczególny przypadek konwersji o większą wartość

Kiedy opisywałem faktoryzację korzystając z podstaw parzystych, wspomniałem o konwersji o większą wartość. Przykładowo przechodzenie na system o 8 więcej, o 16 więcej itp.

Kiedyś proponowałem coś takiego: mając liczbę [a, b] konwersja o 8 polega na przekształcceniu [a, b-8a], po czym naprawa powstałych cyfr za pomocą przeniesień.Wartości przechodzą na cyfrę mniej znaczącą, potem wracają na początkowe miejsce.
Okazuje się, że istnieje możliwość takiego przekształcenia, że do cyfry mniej znaczącej trafia tylko reszta. Nawet mocniej. Nie zawsze musimy mnożyć przez różnicę podstaw docelowej i początkowej.

Pokażę to na przykładzie konwersji z podstawy 20 na 28. Różnica podstaw 8. Ale największy wspólny dzielnik podstaw to 4. Możemy niejawnie 'przejść' na system 20/4 =5, przeliczyć nową wartość różnicy jak dla 8/4=2, oraz 'wrócić z powstałego' siódemkowego. Te przechodzenia zawierają przekształcenia przeciwne, czyli wzajemnie się niwelują. Pozostaje mnożenie przez 2 cyfry w dwudziestkowym, oraz reszta dla cyfry mniej znaczącej jest mnożona przez 4.
Więcej, możemy od razu oszacować wartość w nowym systemie, dobierając resztę tak, by nie mieć nadmiaru bądź niedomiaru w cyfrze.

Przykład: przeliczyć [17 14] _20 na system o podstawie 28
największy wspólny dzielnik systemów to 4, oraz 28/4=7. Jest to wielokrotność, która pozostaje przy cyfrze. 
2*17/7 daje modyfikatory 4 i resztę 6, ale możemy być 'sprytniejsi' i weżmiemy 5 oraz resztę -1.  Wtedy zamiast odejmować 6*4 od 14 dodamy doń 1*4.
Teraz odejmiemy znaleziony modyfikator od cyfry dziesiątek, resztę od cyfry jedności:
[17-5, 14-(-1)*4] = [12, 18] _ 28
Cyfry najmniej znaczącej nie przekształcamy. Otrzymujemy poprawną wartość cyfr w nowym systemie!
Modyfikatory odejmujemy dopiero po użyciu cyfr. Obliczenia są tak łatwe, że przy odrobinie wprawy wiemy kiedy brać nieco więcej, oraz rzadko kiedy potrzebne są dodatkowe przeniesienia między cyframi.

Przypominam, że konwersja polega na iteracyjnym dołączeniu cyfry od liczby, przekształceniu wszystkich cyfr uzyskanego podciągu.

Fragment takich obliczeń (dołączamy kolejną cyfrę 21):
z systemu 28 na system 36 (wspólny dzielnik 4), zmodyfikowana podstawa to  36/4=9
[4, 13, 11, 21] _28  podciąg cyfr podczas konwersji, ((4*36+13)*36+11)*28+21
[8, 26, 22, ]    podwojone cyfry
[1*9-1, 3*9-1, 2*9+4, ]  wyłuskanie modyfikatorów - wielokrotności 9 i reszt
[4-1, 13-3, 11-2, 21]   uwzględnienie modyfikatorów
[4-1, 13-3+4, 11-2+4, 21-16]   uwzględnienie poczwórnych reszt
[3, 14, 13, 5] _ 36  nowa wartość podciągu ((3*36+14)*36+13)*36+5

Jeśli największy wspólny dzielnik podstaw to 1, nie możemy skorzystać z tego sposobu.

Przekształcenie np. z systemu binarnego na dziesiątkowy i odwrotnie

Konwersje z systemu o podstawie p na p*q generują proste przekształcenia, które można zastosować np. na przekształcenie z systemu binarnego na dziesiątkowy i odwrotnie. Wykorzystywana jest konwersja
a_n(p^n q^n) + ... + a_1 (p*q) + a_0 = q^n a_n * p^n + ... + q a_1 *p + a_0 .
Czyli potęgi q są w podstawie, a po drugiej stronie w 'cyfrach'.
Korzystamy cały czas z założenia, że tymczasowo na miejscu cyfr stoją liczby. Następnie przeniesieniami 'naprawiamy' liczbę.

Mając liczbę [a_n ... a_0] w systemie o podstawie p, z każdej cyfry a_i wyłączamy wielokrotność q
a_i = b_i*q+c
wstawiając w miejsce a_i wartość b_i (dzieląc po prostu zmodyfikowaną wartość a_i przez q). Resztę c przy i dodatnim przenosimy do 'cyfry' mniej znaczącej
a_{i-1} = a_{i-1}+p*c.
Dla i=0 zostawiamy - jest to cyfra w systemie o podstawie p*q.
Do 'liczby' [b_n ... b_1] stosujemy rekursję. Cyfra najmniej znacząca nie jest brana pod uwagę aż do powrotu z rekursji.

Przykład. Przedstawić 1000 0010b w dziesiątkowym.
Mamy p=2, q=5, zatem przy przenoszeniu pozostawiamy piątki. Bardzo szczegółowo:
2 0 0   0 0 1 0
4 0   0 0 1 0
5+3   0 0 1 0
5   5+1 0 1 0
5   5 2 1 0
5   5 0 5 0
5   5 0 5 0
cyfra najmniej znacząca 0, liczbę b = [1 1 0 1] poddajemy rekursji:
1 1 0 1
      5 3
Następną cyfrą jest 3, i uzyskujemy ciąg jednoelementowy, który daje najbardziej znaczącą cyfrę 1.
Wracamy uzyskując 1000 0010b = 130.

Jeśli chcemy liczbę z systemu o podstawie p*q przedstawić w systemie o podstawie p, odcinamy cyfrę najmniej znaczącą, po czym mnożymy przez q. Następnie stosujemy rekursję do wyniku póki liczba nie 'zniknie' i 'naprawiamy'.

Przykład. przedstawić dziesiątkowe 128 w piątkowym.
Ponieważ 10=2*5, mnożymy cyfry przez 2 po odcięciu cyfry najmniej znaczącej 8:
2 4
Rekursja po odcięciu 4:
4
Powrót: 4 4 8
Naprawa: 4 4 8 => 1 0 0 3, czyli 128 = 1003_5. Zgadza się, 128 = 5^3+3.

Przekształcenia te mogą być pomocne przy moim ostatnim pomyśle na faktoryzację, w którym znajduję własności dzielników liczby nie znając tychże. Na razie go nie upubliczniam, w metodzie porównuję zbiory rozwiązań bardzo prostych kongruencji nieliniowych.

Równoczesne sprawdzanie dwu kandydatów na dzielniki

Dokończyłem standardowy przykład, w którym wyrażałem n w systemach o podstawach s*t, s<t oraz s i t są liczbami nieparzystymi sumującymi się do stałej wartości parzystej równej sufitowi z pierwiastka z n.

Wykorzystałem
  n = a*(s*t)+b,
gdzie s=3+2k, t = 2(1/2 int sqrt n)-1-2k (nieparzysta wartość naturalna ograniczajaca pierwiastek z n z dołu).

Liczba n jest równaniem kwadratowym zmiennej k, zatem druga różnica jest stałą równą 8. Ewentualny dzielnik liczny n dzieli także resztę b, przy czym należy sprawdzić obydwa przypadki dla s oraz dla t.

Startujemy z s=3, t odpowiednio dopasowane do n. Przekształcamy n zwiększając równocześnie s o 2 oraz zmniejszając t o 2. Oznacza to przechodzenie na system o podstawie (s+2)(t-2), który różni się od danego o stałą (drugą różnicę, czyli 8a). Wartość a dla małych s bardzo szybko maleje. Przy małych a możemy zapisywać współczynnik b w systemach o podstawach s oraz t. Te przekształcenia są liniowe.
Mając iloczyn st, iloczyn (s+2)(t-2) uzyskuje się przez dodanie 8a do zapamiętanej wartości, ta zaś jest dodawana do podstawy. Ze wzrostem s (suma s+t jest stała) wartość ta maleje.

Jest to nadal trial division. Lecz wydaje mi się, że pomijając skrajne przypadki małych s uzyskamy algorytm mający mniej operacji aniżeli algorytm konwersji na kolejne systemy liczenia. Zwłaszcza jeśli nie zmieni się współczynnik a, gdyż po zmianie s, t zmniejszamy współczynnik b o pewną stale malejącą o 8 wartość (modulo st).
Niestety, występuje dzielenie b/s oraz b/t. Brak części ułamkowej tych ilorazów wskazuje dzielnik. Nie można przewidzieć jego wystąpienia.

Algorytm kończy się, gdy s stanie się większe od t. Praktycznie wyglądało to tak, że zapowiadał się koniec dla współczynnika a=4, dla którego jest kilkadziesiąt procent przypadków. Jednak w przykładzie pojawiły się wartości iloczynów st, dla których a było równe 3.
Obliczenia pomijając początkowe przypadki są bardzo proste. Jeszcze je badam.

Czy można pominąć część dzielników podczas faktoryzacji

Przed burzami miewam czasem ciekawsze pomysły. Połączyłem metodę faktoryzacji, która sprawdza od razu duże i małe dzielniki opublikowaną 7 sierpnia 2013 roku
'Faktoryzacja, przepis z prostymi przekształceniami'
z jedną z pierwszych faktoryzacji przez konwersję na inny system.
Okazało się wtedy, że nie muszę sprawdzać nie tylko liczb parzystych jako kandydatów na dzielniki (po wyłączeniu 2 na pewno nie są dzielnikami), ale mogę nawet pomijać wszystkich kandydatów mniejszych niż około 1/3 pierwiastka kwadratowego kosztem około log_3 n obliczeń nwd.

Łącznie można pominąć n/6 wartości szukając dzielników liczby n. Dzielnik znajdziemy przeliczając nie więcej niż trzykrotność jego wartości.
W poniższym sposobie dzielnik d pojawia się co d wartości, na kolejny ślad obecności trafimy nie później niż zmiejszając podstawę systemu o 3d.
Dodatkowo, dla bardzo dużych liczb jesteśmy w stanie łatwo oszacować na podstawie reszty, że nie ma blisko dzielników, lecz nie można tego łatwo sformalizować. Praktycznie można przyspieszyć stosując w niektórych przypadkach konwersję nie o 2, lecz o 6. Tutaj ignorowane.

Lemat 1. Jeśli k jest względnie pierwsze z a, to
  nwd(a,b) = nwd(a, k*b).
Jeśli jest wspólny dzielnik d = nwd(a,b), to d dzieli a oraz b, zatem dzieli także k*b.
Jeśli d = nwd(a, k*b), to d dzieli równocześnie a oraz k*b. Skoro a oraz k są względnie pierwsze, zatem d dzieli b. Stąd nwd(a,b) = d.

Lemat 2. Jeśli n = ap^2+bp+c, to dzielnika n można szukać jako nwd(p,c) dla pewnego p.
Zapiszmy n nieco inaczej n = (ap+b)p+c. Liczba się dzieli przez d, gdy c się dzieli przez d oraz iloczyn (ap+b)p się dzieli przez d. Ponieważ p jest modyfikowane, można się ograniczyć do szukania tylko nwd(p,c).

Ciekawe rezultaty są wtedy, gdy d jest małą liczbą, np. 3. Wtedy co trzecia podstawa jest podzielna przez 3, co piąta jest podzielna przez 5. Po wyeliminowaniu 2 jako dzielnika podstawy systemów są nieparzyste.

Konwersja na system p, kiedy podstawa systemu jest równa p+2 dla liczb 'trójcyfrowych' n = ap^2+bp+c wygląda następująco:
a  b+2a
a  b+4a  c+(b+2a)
co można zapisać równaniami:
a'=a;
b'=b+4a;
c'=c+(b+2a);
Jeśli którakolwiek z b', c' jest większa niż p, wykonujemy operacje (c'-=p; b'++;) oraz (b'-=p; a'++;) odpowiednio. Nie więcej niż dwukrotnie każda.

W praktyce, zaczynając od podstawy będącej liczbą nieparzystą mniejszą od pierwiastka z n, mamy tak długo do czynienia z trzema 'cyframi', że trzecia część pierwiastka kwadratowego z n jest jeszcze liczbą 'trójcyfrową'.  

             Algorytm:
Inicjacja: niech p będzie największą liczbą nieparzystą mniejszą niż pierwiastek kwadratowy z liczby nieparzystej n. Przedstawiamy liczbę n = p^2+bp+c, gdzie b jest mniejsze niż 4 oraz -1<c<p. Zapamiętujemy podłogę z p/3 jako x.
Powtarzamy w pętli dopóki nie znajdziemy dzielnika:
Jeżeli c jest zerem, dzielnikiem jest p.
Jeśli p jest podzielne przez 3, szukamy liczby q takiej że p=q*3^i, gdzie i jest całkowite oraz 3 nie dzieli q. Obliczamy d = nwd(c,q). Jeśli d>1, d jest dzielnikiem.
Zmniejszamy p o 2 i porównujemy z x. Jeśli p>x wracamy do pętli, w przeciwnym przypadku n jest pierwsza.

Ponieważ co szóste p jest nieparzyste podzielne przez 3, w algorytmie przyjmujemy za p tylko wartości nieparzyste, zatem co trzecie p jest podzielne, a na dodatek niech t=p/3. Wtedy co trzecią iterację mamy wartość t-2, zapamiętując je nie musimy dzielić przez 3 więcej niż raz, a następnie zastosować funkcję rekursywną:  t[0]=p/3, t[j+1] = 1/3 t[j],
zmienna wewnętrzna i odmierza co 3 wejście, inicjowana jest jako reszta i=2 (t[j] = 3r+1), i=1 (t[j] = 3r+2), i=0 (t[j] = 3r), r naturalne;
new f( t[j] ): { t[j+1] = floor t[j]/3; i=3-(t[j]%3)); }, t[j]>1
tworzy t[j+1].ty element. Funkcja zwraca dodatnią wartość nwd albo ujemną wiadomość 'jest reszta z dzielenia przez 3':
f (t[j])(i) = {
  i--;
  if( 0<i) return -i; else {
    i = 3;
    t[j] -= 2;
    k = f(t[j+1]);
    if(0<k) return k; else return nwd(t[j],c);
  }
};

Złożoność.
Konwersje mają złożoność arytmetyczną stałą, pamięciowo liniową.
Obliczanie nwd ma najgorszą złożoność, i to ono decyduje o złożoności algorytmu. Przy okazji, wartości dla nwd są ograniczone, jedna przez p, druga przez p/3. Obydwie mniejsze niż pierwiastek z n.


Kontrprzykład numeryczny, liczba pierwsza n = 1447 = 38^2+3
a  b  c   p   t[j]  nwd(d,c)
1  0  3   38
1  2  4   37    // algorytm
1  6  12 35
1 10 28 33 t[0]=11 nwd(11,28)=1
1 15 21 31
1 20 26 29
1 26 16 27 t[0]=9 t[1]=3 t[2]=1 nwd(1,16)=1
2  7  22 25
2 16 11 23
3  5  19 21 t[0]=7  nwd(7,19)=1
4  0  3  19
5  0  2  17
6  6  7  15  t[0]=5  nwd(5,7)=1
8  7  4  13
Wszystkie możliwe liczby nieparzyste uwzględnione, od 13 do 37 przy p, mniejsze przy liczeniu nwd. Jest przeliczonych 13 spośród 19 wartości nieparzystych (2/3 oraz nadmiar). Pozostałe uwzględnione są w czterech nwd. 

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. 

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.

Algorytm faktoryzacji, jak on działa?

Kiedy nie wiadomo co robić, zaczynam szukać innych form, to było granie w gry na sucho, w edytorze tekstu. Gdyż myślałem, że z teorii liczb nic już nie wycisnę. Myliłem się.

Powstał nowy algorytm. Jest bardzo dobry na maszyny wieloprocesorowe, ma stosunkowo mało przypadków do sprawdzenia, i do końca nie wiem jak działa.

Ustalmy 3 procesory. Pierwszy delikatnie przekształca liczbę faktoryzowaną n zapisaną w postaci trójki
n = a*b+c ,
przy czym a rośnie od parzystej liczby ograniczajacej pierwiastek z n z góry do podwojonej tej wartości po liczbach parzystych. Ma zatem 2^i przypadków, gdzie 2^(2i) < n < 2^(2i+1).
Przekształcenie polega na konwersji na system o 2 wyższy, czyli dany tożsamościami:
a*b+c = (a+2)*d+e, gdzie d in {b-2, b-1, b} w zależności od znaku (c-2b).
Dokładniej:
a*b+c = (a+2)*b + (c-2b), -1 < c-2b;
a*b+c = (a+2)*(b-1) + (c-2b+a+2), -a-3 < c-2b <0;
a*b+c = (a+2)*(b-2) + (c-2b+2(a+2)), c-2b < -a-2.
Np. jeśli n jest postaci n = 12441^2+987, to a=12442, b=12440, c=987-12441+12442 = 988; następna wartość jest równa a=12444, zaś b i c są znajdowane według wzoru (tu jest b=12438).
Dla n postaci n=12442^2+c przyjmujemy a=12444, pozostałe są odpowiednio dopasowane, by b było bliskie a, zaś c jest resztą.
Wartości (a,b,c) są przekazywane drugiemu procesorowi.

Drugi z procesorów przyjmujący (a,b,c), gdzie a jest parzyste, usiłuje dokonać konwersji, w której a jest połowione, b jest podwajane, zaś c w zależności od znaku (c<a) jest albo pozostawiane, albo zmniejszane o a, z równoczesnym zwiększeniem b o 1. Wzorami
(2a)*b+c = a*(2b)+c, c<a;
(2a)*b+c = a*(2b+1)+(c-a), c>a-1.
Pętlę kończymy, gdy wartość zmniejszana a jest nieparzysta. Mamy do czynienia co najwyżej z logarytmem binarnym z a przypadków (b jest nie większe niż podwojony sufit pierwiastka kwadratowego z n). Ciąg wartości c jest ciągiem nierosnącym.
Watości (a,c) oraz (b,c) są przesyłane do trzeciego z procesorów. Jedna z tych wartości jest stosunkowo mała.

Trzeci procesor ma najwięcej roboty. Sprawdza największy wspólny dzielnik przesłanych wartości, jeśli jest on większy od 1, znaleziona wartość jest zarazem dzielnikiem liczby n.
Bardzo podobnie wygląda cecha podzielności przez 9 lub 11, która stała się prekursorem algorytmu. Jeśli cyfry liczby dwucyfrowej zapisanej w systemie niedziesiątkowej są równe, to podstawa systemu zwiększona o 1 jest dzielnikiem liczby.

W praktyce nie zawsze musimy przesyłać dane do trzeciego procesora. Ale musimy je przesłać, gdy następuje zmiana wartości c, które nie jest zbyt małe. Przy c=1 możemy zignorować ten krok, przy c pierwszych przyspieszyć sprawdzanie testując podzielność przez c.

Przykładowy fragment rozkładu znajdujący dzielnik: 
n = a*b+c = (2990 * 2987 + 2923) = 8934053
Procesor pierwszy: przelicza od a=2990 do a=5980:
(2990 , 2987 , 2923 ) = (2992 , 2985 , 2933) = ...
= (4348 , 2054 , 3261) = ... 
Procesor drugi (4348, 2054, 3261):
(4348 , 2054 , 3261) = (2174 , 4109 , 1087) = (1087 , 8219 , 0)
Procesor trzeci:
nwd(4348, 3261) = 1;
nwd(2054, 3261) = 1;
nwd(2174, 1087) = 1087, znaleziony dzielnik 1087
Zaś 8934053 = 1087 * 8219.

Pętla działania procesora drugiego tworzy 'grzebień'. Tzn. dla połowy przypadków jest jedna iteracja, każda dłuższa pętla sąsiaduje z dwoma najkrótszymi. Krotność iteracji tworzy podciągi (1, 2, 1, 4, 1, 2, 1, 8, 1, 2, 1, 4, ... ).

Algorytm powstał, kiedy połączyłem cechy podzielności przez liczby sąsiadujące z podstawą systemu, konwersje o potęgi dwójki oraz konwersje o 2.

Pierwiastek kwadratowy przy systemach niedziesiątkowych

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

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

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

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

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


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

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

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

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

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

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

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

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

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

Szukanie najwiekszego wspolnego dzielnika

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

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

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

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

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