Konwersja do systemu Fibonacciego

W poście z 12 czerwca

Faktoryzacja oraz reszty z dzielenia w systemie Fibonacciego
wspomniałem o konieczności zamieniania liczby dowolnego systemu na system Fibonacciego. Ale pokazałem też przebieg który zwiększa cyfry, wstawiając na odpowiednich pozycjach potęgi 2. Ten schemat można odwrócić!

Dołączmy jeszcze jeden schemat do podanych w tamtym poście:
[0 0 2 0 0] = [0 1 0 0 1];
[0 0 3 0 0] = [1 0 0 0 1];
[0 0 4 0 0] = [1 0 1 0 1];
pominę przypadki brzegowe.
Teraz algorytm rekursywny w bardzo szybkiej pętli zmniejsza N aż 4 razy. Wewnątrz ma wolniejszy przebieg, gdyż wykorzystując schematy 'naprawia' liczbę do postaci Fibonacciego. 
Oto ten przebieg, wykorzystujący postać N w systemie czwórkowym:
IN N_4
i = 2^k, 2^k < N < 2^(k+2), czyli i<N<4i;
F_0+R = F_0+N%i, liczba Fibonacciego F_0 z resztą R (reszta zawiera cyfry N poza najbardziej znaczącą, te najbardziej znaczące dodajemy do postaci F_i)
powtarzaj dopóki i>1
  i /= 4;
  F_i = 4*F_i +HIGH(R); // dodajemy najbardziej znaczącą cyfrę R
  napraw F_i do systemu Fibonacciego;
OUT liczba w systemie Fibonacciego

Pomijając przekształcenia naprawiające liczbę, jest tu bardzo szybkie zchodzenie.
Przykład: N = 8934053d = 202011022211_4
Mamy N = 2^23 + 2011022211_4 = 2*2^22 + 545445, i = 2^22
F_0 = 10F z resztą 545445
i/=4, czyli naprawiamy liczbę [4 0+0]F = 10000F
następna iteracji zaczyna się od 40002 i kończy na
100 00000F z resztą 11022211_4 = 21157
Teraz jeden schemat działający na 400 00000 dostarcza
10101 00000F z tą samą resztą.
Kolejna naprawa 40404 00001 przekształca się do
101 00010 00100F z resztą 1022211_4 = 4773
I dalej:
1 00101  01000 10000F z resztą 22211_4 = 677
1001 00010 10100 01001F z tą samą resztą
10 00101 00100 01010 10010F z resztą 2211_4=165
10001 00000 00100 10010 00100F z resztą 211_4=37
100 00010 10100 00010 10100 10001F z resztą 11_4 = 5
1 00000 00100 01001 00100 10000 10000F z resztą 1_4 = 1
i ostatnie, gdyż mamy już i=2^2=4
101 01001 00001 00101 01010 01000 10101F, co jest wartością Fibonacciego liczby N.

Mix dzielenia pisemnego i przez zmianę podstaw

Swego czasu opublikowałem dzielenie przez zmianę podstaw. Przypomnę przebieg na przykładzie

N = 8934053 : 8, 10-8=2, iloraz d[] z resztą r, val oznacza wyliczoną wartość

d_5 = n_6 = 8; 

d_{n-1} = 2*val a_{n+1} + a_n, n=5..0,

d_{-1} = 2*d_0+a_0 = 8*k+r

d_0 = d_0+k

czyli 

8   9    3      4    0      5     3 

8  25  53 110 220 445 893 = 8*111+5

Wynik uzyskujemy przez dziesiętne dodawanie ww wartości i przesunięcie o jedną pozycję, czyli 

cyfra jedności (445+111)%10 = 556%10 = 6

cyfra dziesiątek (220+55)%10 = 275%10 = 5 itd.

cyfra 'najbardziej znacząca' (8+3)%10 = 11%10, pojawiła się dodatkowa cyfra.

Duże wartości, wydaje się, że dzielenie pisemne jest prostsze, gdyż 8:8=1, podobnie 9. 

Liczbę N można rozbić na dwie: 8832048 + 102005, oraz porachować osobno: 1104006 + 102005:8. Użyjemy tego pomysłu do mixu, dla kolejnych cyfr.

Mix dzieleń polega na połączeniu obu metod przy przebiegu. Dzielimy w cyfrze przez 8, dopiero resztę mnożoną przez 2 dodajemy do cyfry mniej znaczącej. Cyfry są przesunięte o jedną pozycję.  Poprawkę stanowią współczynniki przy 8.

A oto dokładny przebieg: 

d_6 = 8 = 1*8+0 => 0

d_5 = (2*0+9) = 1*8+1 => 0

d_4 = (2*1+3) = 5 = 0*8+5 => 1

d_3 = (2*5+4) = 14 = 1*8+6 => 5 

d_2 = (2*6+0) = 12 = 1*8+4 => 6

d_1 = (2*4+5) = 13 = 1*8+5 => 4

d_0 = (2*5+3) = 13 = 1*8+5 => 5

odczytujemy wynik

0 0 1 5 6 4 5 , 5   cyfry i reszta

1 1 0 1 1 1 1        poprawki

sumujemy poprawki 

1 1 1 6 7 5 6, reszta 5

Pamiętając o poszczególnych cyfrach i poprawkach, uzyskujemy bardzo małe wartości jak przy dzieleniu pisemnym, zachowując prostotę obliczeń dzielenia przez zmianę systemów. 

Dla większych dzielników możemy zastosować dzielenie chłopów rosyjskich, dzięki czemu wartości obrabiane nie będą nigdy większe niż kwadrat dzielnika.

Faktoryzacja oraz reszty z dzielenia w systemie Fibonacciego

System Fibonacciego zawiera na pozycjach cyfr zamiast potęg np. 10 liczby Fibonacciego  przesunięte o 3 pozycje. Pozwala to uzyskać system zawierający tylko 0 i 1 oraz spełniający schemat
  1)    [1 0 0] = [0 1 1].
Liczba jest sumą dwu poprzednich, z których pierwsze są 0 i 1. W tym systemie pozycja jedności F_3=1.
Np. 9 = [1 0 0 0 1] = 1*8+0*5+0*3+0*2+1*1.
Na potrzeby faktoryzacji ignoruję fakt, że na każdej pozycji jest tylko 1 albo 0, dopuszczam dowolną naturalną wartość z zerem włącznie. Wtedy liczba w systemie Fibonacciego jest podzielna przez a, gdy jej zapis w systemie Fibonacciego M składa się tylko z a oraz zer. Gdy liczba nie jest podzielna przez a, resztę można umieścić w cyfrze najmniej znaczącej, np.
[3 0 3 0 4] = 3*8+0*5+3*3+0*2+4 = 3(8+3+1)+1 = 37

Dodawanie w systemie Fibonacciego sprowadza się do orowania cyfr, zaś gdy obie są ustawione, do mnożenia przez 2.
Szczególne przypadki cyfr najmniej znaczących (z definicji cyfr systemu):
[. 1 0 1] = [. 0 2 0] = 4
[. 0 1 0] = [. 0 0 2] = 2
Dla pozostałych pozycji mamy schematy, które się dodaje do aktualnych cyfr
  2)  [1 0 0 1] = [0 2 0 0]
  3)  [1 0 0 0 1] = [0 0 3 0 0]
Dowody
[1 0 0 1] = [0 1 1 1] = [0 1+1 0 0] = [0 2 0 0];
[1 0 0 0 1] = [0 1 1 0 1] = [0 0 2 1 1] = [0 0 2+1 0 0] = [0 0 3 0 0].

Pierwszy pomysł na faktoryzację z systemu Fibonacciego. Zapisuję liczbę N w tym systemie jako M_{}, następnie za pomocą uogólnień schematów 1) i 2)
  [a 0 0] = [0 a a];
  4)  [a 0 0 a] = [0 2a 0 0]
przekształcam by uzyskać odpowiednią ustaloną wartość zamiast jedynek. np.
z liczby
[1 0 0 0 1 0 1 0 0 1]
można zrobić wartość podzielną przez 2
[0 0 3 0 0 0 0 2 0 0] = [3 0 0 0 0 2 0 0], przeniesienie reszty do cyfry jedności
[3 0 0 0 0 2 0 0] = [2 1 1 0 0 2 0 0] = [2 0 2 1 0 2 0 0] = [2 0 2 0 1 3 0 0] =
[2 0 2 0 0 4 1 0] = [2 0 2 0 0 2 3 2] = [2 0 2 0 0+2 0 1 2] = [2 0 2 0 2 0 0 4] =
[2 0 2 0 2 0 2 0]
lub przez 3:
[3 0 0 0 0 2 0 0] = [3 0 0 0 0 0 2 2] = [3 0 0 0 0 0 3 0]
Ale jak widać w przykładzie, występują nawroty. Zwłaszcza kłopotliwe są fragmenty [. a b 0 a .]

Zatem drugi pomysł. Chcąc uzyskać resztę z przedziału (2^(i-1), 2^i], przekształcamy liczbę do postaci Fibonacciego M_i tak, by zamiast jedynek były liczby 2^i, podwajając cyfry z pomocą 4)
[1 0 0 0 1 0 1 0 0 1] = [0 1 1 0 1 0 1 0 0 1] = [0 0 2+1 0 0 0 1 0 0 1] =
[2 1 1 0 1 0 0 1] = [2 0 2+1 0 0 0 0 1] = [2 0 2 1 1 0 0 1] = [2 0 2 0 2 1 0 1] =
[2 0 2 0 2 0 1 2] = [2 0 2 0 2 0 2 0] =    (M_1)
[2 0 0 2 4 0 2 0] = [0 4 0 0 4 0 2 0] = [4 0 0 4 0 0 4] =    (M_2)
[0 8 0 0 0 0 4]     (M_3)
Jest to liniowe i nie ma nawrotów.
Jeśli teraz chcemy znaleźć resztę przez d<2^i, jest ona równa iloczynowi postaci
M_i *(2^i - d) + r ,
gdzie r jest resztą (cyfra jedności M_i lub M_i-2^i, gdy M_i ma ustawioną tę cyfrę ).
Ten drugi sposób jest dowodem na cechę podzielności przez d w systemie.

I znowu, działania na dużej liczbie zostały zastąpione działaniami na znacznie mniejszych, gdyż M_{i+1} ma co najmniej dwie cyfry mniej niż M_i, 2^i-d<2^{i-1), zaś reszta też jest mniejsza niż 2^{i-1}.

Pozostaje sposób konwersji N na postać Fibonacciego M_0. Chyba najbardziej opłacalny jest klasyczny szukania największej liczby Fibonacciego F_k nie większej niż N i odejmowanie. Jeśli chcemy zmniejszać N wcześniej, możemy odejmować kolejne liczby Fibonacciego, a kiedy F_k >N-\sum{F_{i<k}}, odejmujemy jeszcze 2, razem zostało odjęte F_{k+2}.
Wynika to z przekształceń
[1 1 . . . 1 1] = [1 0 1 . . . 0 1 0 0] lub [1 0 1 . . . 1 0 0 1]
w zależności od liczności jedynek. Uzyskane ciągi mają dokładnie jedną pozycję więcej. Brakuje dokładnie 2, by cały ten ciąg używając 1) stał się kolejną liczbą Fibonacciego, zyskując dodatkową pozycję. Z programowania dynamicznego lub rekursji ustawiamy kolejne jedynki M_0.

Reszta z dzielenia wstecz - poprawka

Opisałem przekształcenie przekształcające postacie binarne, a okazało się, że zastosowanie zwykłego odejmowania jest równie proste i skuteczne.
Chciałem uzyskać binarny wynik z dzielenia przez liczbę. W tym celu potrzebowałem binarnej postaci dzielenia przez pobliską liczbę - nie jest ona potrzebna.

Aby uzyskać binarną postać ilorazu liczb nieparzystych przy dzieleniu wstecz a/p wystarczy odjąć dzielnik i podzielić przez 2. Wpisujemy 1. Jeśli uzyskana wartość  jest parzysta wpisujemy 0 i znów dzielimy przez 2. Kończymy gdy wartość
b= (a-p)/2 or a/2
jest mniejsza niż dzielnik. Postać ta sama: b*2^m, gdzie m jest bliska krotności cyfr binarnych dzielnika
Ciąg zer i jedynek w odwrotnej kolejności jest ilorazem z dzielenia wstecz, o ile zadziała. W przeciwnym razie daje mniejszą wartość. 

Na razie nie widzę, czy reszta przy dzieleniu wstecz zachowuje się w sposób równie przewidywalny jak reszta z dzielenia zwykłego, na pewno nie dla małych i średnich wartości, chociaż występuje tendencja do zmniejszania wartości binarnych, niestety, ciągle zasilana z reszty odbudowującej się blisko potęg 2.

Jeśli chcę uzyskać wynik w innym systemie liczenia, z tabliczki mnożenia szukam cyfry c, by
a%p = c*p (p)
oraz od liczby a odejmuję c*p wypisując cyfrę c, np. w dziesiątkowym 323 : 17, c=9, bo 9*7 = 3 (10), odejmujemy 9*17 = 153.
Najlepiej to robić w systemach z jednoznacznie wyznaczoną cyfrą, np. jedenastkowym. Zwykła rzecz przy dzieleniu wstecz...

Reszta z dzielenia wstecz

Ostatnie kilka miesięcy poświęciłem na sposoby usuwania nieoznaczoności przy znajdywaniu dzielników. Zmiany systemów pozwalają znajdować dokładnie cyfry, o ile uda nam się w jakiś sposób znaleźć krotność cyfr dzielników. W przeciwnym razie pojawia się nieoznaczoność. Sprawdzanie zapisów liczb w różnych systemach do pewnego miejsca działa dobrze, ale nie wiadomo, jak długo.
Zwłaszcza występujące przeniesienia wartości między kolejnymi cyframi bardzo przeszkadzają.

Wróciłem do trial division, tym razem jednak zdecydowałem, że reszta z dzielenia przez p będzie przy cyfrze najbardziej znaczącej. Budowa reszty w systemie binarnym jest następująca:
r = a*2^m
gdzie m+1 jest krotnością cyfr kandydata na dzielnik b, czyli
N = r + b*p,  log_2 b = m+1
Używam dzielenia wstecz, rzadko działającego, w którym iloraz znajdujemy od cyfr najmniej znaczących. 

Przekształcenie badane przechodzi między dzielnikami nieparzystymi p i q  za pomocą przeniesień. Dodatkowo, znaleziony sposób jest liniowy, w którym występują dodawania, odejmowania i czasem mnożenie przez 2 lub 3. Zatem przekształcenie jest liniowe. Co więcej, w odróżnieniu od opublikowanego już algorytmu zmiany podstaw, przekształcenie pozostaje liniowe dla q-p<p.
Czynnością generującą najwięcej pomyłek obliczeniowych dla człowieka jest dbanie o poprawność przeniesień. Reszta to proste działania.

Przykład N = 8934053, reszta z dzielenia przez 59 jest (hexadecymalnie)
N = 38*2^16 + 59*0x1AA9F
Znalezione przekształcenie modyfikuje 0x1AA9F próbując pobrać przy ustawionym bicie i*59 wartość 61, a następnie zależnie od parzystości wyniku albo zeruje, albo ustawia bit i w ilorazie N / 61. Możemy to traktować jako liczbę binarną, stosując dodatkową zmienną interpretującą ustawione bity jako 59, które parzyste przenoszą wcześniej połowę wartości, zaś nieparzyste wykorzystując pożyczki odejmują 61 od wartości też przenosząc wcześniej połowę parzystej różnicy. Reszta nie może być traktowana jako iloczyn 59*i, ale dla przebiegu to i tak ma niewielkie znaczenie. Uzyskujemy nową resztę 69*2^16, zaś postać liczby ulega zmianie do
69*2^16 + 61*0x11A89
Powstała postać binarna jest najlepiej dopasowana do nowej wartości, co dzieje się automatycznie.
Jak to przebiega szczegółowo. Końcówka
1 1 1 1 oznacza 59 59 59 59.
Chcemy zabrać 61 od ostatniego 59, ale nie można. Zatem pożyczamy 1 1 = 0 3 uzyskując postać
1 1 0 3  czyli 59 59 0 183
Od 183 odejmujemy 61 i dzielimy przez 2 uzyskując 58. Wartość parzysta zeruje bit w tablicy ? ? 0 1. Przenosząc na kolejną wcześniejszą pozycję mamy
58/2+1*59 = 88.
Kolejna pozycja to 88/2+59 = 103, nieparzysta, stąd dla tabeli przy 61 mamy ostatnie bity 1 0 0 1.
Jeśli na danej pozycji wartość nieparzysta jest mniejsza niż q, pożyczamy.
Końcówka, na pozycji 2^16 mamy 59+38, z przeniesień dodajemy 33 uzyskując 61+69. Wartość 61 to bit w tablicy ewentualnego dzielnika, zaś 69 reszta.
Z dzielnikiem mamy do czynienia, gdy reszta jest równa 0.

Jeśli poprawnie uda się pobierać pożyczki z binarnych cyfr wcześniejszych, możemy zapamiętując tylko wielokrotność a, tablicę (liczba binarna) ewentualnego ilorazu i dzielnej p, przekształcając je liniowo.
Wtedy algorytm faktoryzacji szacuję na O() = n log n. 

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.

'trial division' algorytm w jeszcze innej odsłonie

Tym razem, po długiej przerwie, kolejny algorytm 'trial division' na rozkład liczby N na czynniki. Jak moje najnowsze podejścia, znów jest prosty dla bardzo dużych wartości, zaś daje w kość dla małych dzielników.

Przedstawiam liczbę N w postaci iloczynu dwu wartości p, q i reszty r:
N = p*q + r,
gdzie p>q są liczbami nieparzystymi, r jest liczbą nieujemną, zmniejszaną by nie przekraczała podwojonego q.
Podstawowym przekształceniem jest
  (p+k)*(q-k) = pq + k(q-p) - k^2      (1) ,
albo w używanej postaci szczególnej
  (p+2)*(q-2) = pq - 2(p-q) - 4      (2) .

Inicjacja następuje dla nieparzystych p = q < sqrt(N)  lub p < sqrt N < q = p+2 gdy sqrt N jest parzyste. Liczba
  r = N - pq
jest parzysta.
W kolejnych iteracjach wykorzystujemy stary iloczyn pq by dostać nowe r następująco:
p = p+2;
q = q-2;
s = r + 2(p-q) + 4;
r = s (mod 2q);  // s można zapamiętać by mniej liczyć
p = p+ floor( s/ 2q );
jeśli r =0, dzielnikami są p i q.

niezmienniki: p, q nieparzyste, z r zabierana parzysta wielokrotność q, która jest dodawana do p.
W ten sposób p rośnie wykładniczo, zwłaszcza dla małych q, q maleje liniowo. Zmiany są bardzo małe, zwłaszcza dla q mniejszych niż pierwiastek sześcienny z 2N. Z kolei dla q mniejszych efekt wykładniczy modyfikuje p bardzo silnie.
W praktyce przeważająca krotność przypadków przypada na małe zmiany.

Do tego stopnia zmiany są małe, że dla dużych q możemy korzystać ze cech podzielności dla q i korzystając z (1) przeskakiwać nie po kolejnych nieparzystych, ale wybierać te, które nie mają dzielników 3, 5. Podobna zakusa istnieje też dla p, ale tylko dla tak dużych wartości, by p/q <2.

Przykład: fragment
sqrt 8934053  jest bliskie 2988
p = 2989, q = 2987, r = 5910
różnica 2(p-q)+4 = 8, zatem
p = 2991, q = 2985, r = 5918
jeszcze dwa razy i mamy
p = 2997, q = 2979, r = 5990 = 2*2979+32,
za duże r, dochodzi do głosu dzielenie r/ 2q , zatem uzyskujemy zwiększenie p:
p = 2999, q = 2979, r = 32

przy kolejnym zmniejszaniu q przenoszenia do p występują coraz częściej, aż od
p = 4169+2, q = 2141, r = 2*2141+3942
stają się powszechne
dochodzimy w pobliże dzielnika,
p = 8189+14 = 8203, q = 1089, r = 14*1089+986
p = 8205+14 = 8219, q = 1087, r = 14*1087
p = 8221+12 = 8233, q = 1085, r = 12*1085+1248

poczynając od wartości
p = 33971+258, q = 261, r = 258*261+284
lepiej zmienić algorytm. Tu już q jest stosunkowo bardzo małe względem p.

Czy istnieje logarytmiczny sposób rozkładu liczb

Czy istnieje sposób rozkładu liczb, w którym kolejne obliczenia wyznaczają kolejne cyfry dzielników?
Wydaje się to niemożliwe dzięki przeniesieniom.

Lecz testuję pewien sposób dzielenia, który w dziesiątkowym działa w znikomej krotności przypadków, i on jest w stanie to ominąć.
Według jednego z pomysłów należałoby zbadać wykładniczą liczność przypadków dopasowania cyfr. Badane dzielenie wyłuskuje prostą zależność między cyframi dzielników, które ogranicza tę liczność. Dodatkowy test wyznacza konkretne cyfry, zwłaszcza kiedy w trakcie przebiegu zostanie złamana symetria.

Występuje nieoznaczoność, gdyż teoretycznie dostaję 2 wartości dla liczb pierwszych, co najmniej 4 dla złożonych. Są one symetryczne, i ta symetria nie pozwala przewidywać dokładnie. Jednak większość przypadków po złamaniu symetrii jestem w stanie wyznaczać jednoznacznie. Łamanie symetrii następuje tuż po pierwszym wystąpieniu pary różnych cyfr.

Dla wygody rachuję w binarnym. Wtedy dzielenie zwraca wielomian postaci np.
a+b+[0]^*+[1]^* = 1
z dwoma niewiadomymi. Zapis [0]^* oznacza, że dodawane jest kilka zer, ale nie muszą wystąpić. Po uproszczeniu mamy a+b=1 lub a+b=0, gdzie a i b są cyframi binarnymi, dopisywanymi do kandydatów na dzielniki.
Drugi z testów jest związany z parzystością różnicy liczby rozkładanej i iloczynu cząstkowego znajdowanych kandydatów na dzielniki.

Zmniejszenie krotności przypadków dla bardzo dużych liczb, poprawka

Wspomniany w styczniu sposób zmniejszania krotności nie jest tak dobry, jak początkowo uważałem.
Przetestowałem kilka sposobów dla postaci
N = ap^2 + bp +c

1) trzymam jak najwięcej w a, dopasowuję skok k by c-bk > p
mamy pierwiastek sześcienny przypadków, a początkowo błyskawicznie maleje, po przekroczeniu p > root[3] N wartość a zmniejsza się o 1 na każde obliczenie.
Jednak istnieje duże ryzyko chybienia dzielników. Podczas testu wykryłem  dzielnik tylko dlatego, że go znałem, i w odpowiednim miejscu zastosowałem cechę podzielności. Było bardzo blisko.

2) Jak największe b, a=1
Nie różni się licznością przypadków od 'trial division' przez kolejne liczby nieparzyste.

3) c jest kwadratem dużych wartości, bez nawrotów przy zmniejszaniu a (by zastosować (sqrt a *p+c)^2 +/- bp )
Szybkie zmniejszanie a i c, ale wartość w c za szybko stała się na tyle mała, że sposób przekształcił się w standardowe konwersje o 2.
Myślałem, że jest możliwość sprawdzania równocześnie małych i dużych dzielników, ale do tego trzeba nawrotów, tzn. przenoszenia wartości z a do c i obliczanie kolejny raz dużych kwadratów.

Wnioski. Dla dzielników mniejszych niż pierwiastek sześcienny root[3] {6N} i tak się stosuje konwersje o 2. Dla dzielników większych, a jest na tyle małe, że b zaczyna dominować nad zmianami c. Duże b pociąga znów konwersję o 2.
Zmiany c przy ustalonych a, b i wolno rosnącym p to funkcja kwadratowa z wierzchołkiem w środku przedziału, zatem trudno jest dopasować skoki poza zboczem rosnącym.
Zatem nie uda się zmniejszyć liczności przypadków sprawdzanych do pierwiastka sześciennego.

Zmniejszenie krotności przypadków rozkładu dla bardzo dużych liczb

Ten pomysł nie jest najlepszy obliczeniowo, ale w pobliżu pierwiastka z liczby rzadko natrafia się na dzielniki. Znalazłem sposób zmniejszenia krotności przypadków, szacuję, że nawet do pierwiastka sześciennego z rozkładanej liczby. Ale koszt jest duży.
Stosując konwersję odwrotną dla liczby N=(a*p+b)*p+c, gdy p^2<N<(p+1)^2, a=1, b jest zerem lub jedynką, normalnie przenosiłem nadmiar z c do b.
Można tego zaniechać. Wtedy współczynniki dla kolejnych iteracji p-k są funkcjami rosnącymi, a właściwie szeregami
c[0] = N-(p+b)*p;
c[k+1] = b[k]+1 + c[k];
b[0] = 0 | 1; (zależy od liczby N)
b[k+1] = 2*k+b[0];

Dosyć szybko c[k]>p-k. Z dzielnikiem mamy do czynienia wtedy, gdy p-k dzieli c[k]. Wyrazy ciągu c[k] można łatwo wyznaczyć, interesują nas tylko te, kiedy c[k] ma resztę bliską p-k. 
Najciekawszy efekt jest wtedy, gdy zapisujemy c[k] w systemie o podstawie p-k.
Mamy wtedy przekształcenia kroku iteracji k:
dec p;
konwersja c[k] na system o podstawie 1 mniejszą;
c[k] + b[k]_p; (b[k]_p oznacza liczbę b[k] w systemie o podstawie p)
sprawdzenie, czy cyfra jedności c[k] jest zerem.

Korzystając z monotoniczności c[k] można wyznaczyć sumę ciągu (b[i]+1)  tak, by wyznaczyć kolejną wartość przy występującym skoku cyfry jedności c[k], oraz wyznaczyć c[k], p. 

W pewnym momencie c[k] zaczyna być liczbą trzycyfrową. Wtedy lepiej zmienić algorytm. Gdyż kontynuacja spowoduje, że c[k] będzie się rozrastać do wielkości rzędu liczby N.

Przykładowo dla N=8934053 mamy jak niżej. Zapis [1,2]_7 oznacza liczbę zapisaną siódemkowo 1*7+2
b[0] = 1, c[0] = 2921, p = 2988
b[7] = 15, c[7] = 2977; p = 2981
b[8] = 17, c[8] = 2993 = [1,13]_2980; p = 2980
b[9] = 19, c[9] = 3011 = [1,32]_2979; p = 2979
następne można pominąć aż do b[53], gdyż jest skok cyfry jedności c[k]_p
b[54] = 109; c[54] = 5891 = [2,23]_2934; p = 2934
b[1901] = 3803; c[1901] = 3618623 = [3,68,0]_1087; p=1087  dzielnik
b[2981] = 5963; c[2981] = 8892263; p=7

Podnoszenie do kwadratu lub mnożenie

Podobny sposób, który był użyty do znajdowania pierwiastka kwadratowego w ostatnim poście, może też być użyty do podnoszenia liczby do kwadratu.
Są jednak pewne komplikacje - albo pracuje się na dużych liczbach, albo odcina się reszty, które są wyznaczonymi już cyframi wyniku.

Podczas badania odcinania reszt, wpadłem na pewną metodę. Każdy z czynników mnożenia zawiera 'schemat' dojścia konwersjami do potęgi 10, drugi czynnik podążając za tym schematem modyfikuje się w bardzo specyficzny sposób. Kiedy pierwszy stanie się już potęgą 10, drugi staje się wartością iloczynu.
Najciekawsze jest to, że krotność konwersji jest nie większa niż podwojona liczność cyfr liczby, z której pobieramy schemat. 
Rozpracowałem już mnożenie, ale zanim je opublikuję, chcę je skonsultować z naukowcami. Powód jest prosty, szacując złożoność, uzyskuję że jest ona logarytmiczna. Zaś mnożenie nie powinno być mniej pracochłonne niż dodawanie czy odejmowanie. Tym bardziej, że przy mnożeniu korzystam z kilku dzieleń.

Pierwiastkowanie złożoności logarytmicznej

Ostatni algorytm konwersji wskazał możliwość uzyskania pierwiastka kwadratowego stosunkowo małym kosztem. Jako pierwszy znalazłem sposób dla licz zapisanych w postaci binarnej. Istnieje też sposób w systemie dziesiętnym, który jest nieco bardziej rozbudowany.

Liczbę będącą kwadratem p^2 zapiszemy jako czwórkę [1,0,0,p], w której pierwsze pozycje to cyfry rozwinięcia w systemie o podstawie p, czyli
N = 1*p^2 + 0*p+0 .
Oto wszystkie używane manipulacje:
1) [1, a+2, b+a+1, p] = [1, a, b, p+1]
2) [1, a, 2*p+c, p ] = [1, a+2, c, p]
3) [100, 10*a, b, p] = [1, a, b, 10*p] . 
Pierwsza to konwersja na system o 1 większy.

Bierzemy liczbę N i odcinamy grupy po 2 cyfry. Następnie wyznaczamy  przybliżoną wartość pierwiastka. Mając ją, mnożymy wszystkie wartości przez 100 i stosujemy konwersję 3). Uzyskujemy w ten sposób przybliżenie pierwiastka uciętej liczby z niedomiarem, który poprawiamy za pomocą 1) lub 2). Powtórzą się one nie więcej niż 9 razy.
Staramy się doprowadzać do tego, by wszystkie wyrazy były nieujemne, oraz na pozycji drugiej była wartość 0 lub 1 ([1, 0|1, b, p]. Bardzo rzadko pojawia się przypadek, kiedy uzyskamy [1, 2, 0, p] - o 1 mniej niż kwadrat. Zamieniamy go na niepoprawne w systemie [1, 1, p, p]. Uzyskamy w ten sposób dokładną kolejną cyfrę pierwiastka.
Zatem złożoność to logarytm z liczby znajdowanych cyfr plus stały narzut.

Jeśli kto chce znaleźć właściwą wartość pierwiastka i nadmiar, wartość pierwiastka jest przybliżana przez p, zaś nadmiar z [1, a, b, p] jest liczony jako a*p+b większy od p^2.

Przykład: 8 93 40 53
kwadrat z 8 to [1, 2, 0, 2]=[1,1,2,2] (2^2+2*2) - jest to ten wyjątkowy przypadek z drugą wartością większą niż 1.
Mnożymy przez 100 i dodajemy kolejne 2 cyfry uzyskując 893
[100, 100*1/10, 100*2+93, 2] = [1, 10, 293, 20]  z 3)
konwersje 1) i poprawki 2) przybliżają p do pierwiastka:
[1, 10, 293, 20] = [1, 8, 284, 21] = [1, 10, 242, 21] = ...
= [1,1, 23, 29] jest kolejnym przybliżeniem do pierwiastka z przyciętej wartości.
Dokładamy kolejne cyfry 40 uzyskując po konwersji 3):
[1, 10, 2340, 290] = ... = [1, 1, 238, 298]
Ostatnie cyfry 53 i mamy
[1, 10, 23853, 2980] = ... = [1, 1, 2921, 2988]
pierwiastkiem całkowitym jest 2988 z resztą 1*2988+2921.
Możemy kontynuować dokładając '00', uzyskamy w ten sposób kolejne cyfry części ułamkowej z niedomiarem. 

Konwersje a rozkład liczb o dużych podstawach

W ostatnim poście wspominałem o sposobie 'trial division' dla liczb postaci
n = a*p+b   (1)
gdzie a,b<p oraz p^2>n.
Przy zmieniającym się a szukamy dzielników kiedy a=b (cecha podzielności, bo a*p+a = a*(p+1) = n ).
Zastanawiałem się nad przechodzeniem między różnymi p. Mamy tu aż n-sqrt(n) możliwości, przy czym dla dzielników wystarczy znajomość sqrt(n) wartości p. Albo wystarczy umieć znajdować p dla zmienianego a. Jest to obarczone błędem, jednak łatwym do skorygowania.

Niech r oznacza różnicę między p-q, gdzie
n = (a+1)*q+c = a*p+b .
Różnica ta jest bardzo duża dla małych a, zaś dla większych maleje, w szczególności dąży do 1 gdy a dąży do sqrt(n).
Kiedy konwertujemy a*p+b o k-tą wielokrotność r, wykonujemy przekształcenia:
a*p+b = a*(p-rk)+(b+ark)  (2)
Dobierając r bliskie a*r=p, możemy w ten sposób zwiększać a o prawie dowolną wartość. Oczywiście, dla większych a błąd jest mniejszy. Pozwala to pomijać podczas faktoryzacji kolejne przypadki, kiedy a nie jest dzielnikiem.

Występuje błąd, który należy skorygować. Jeśli uzyskamy wartość c>a+k, dokonujemy dzielenia c/a, co jest znacznie prostsze i szybsze niż dzielenie p/a. Uzyskaną część całkowitą odejmujemy od c uzyskując właściwą resztę.
W przypadku niedomiaru, dopasowujemy współczynnik konwersji
1+(q-c)/(a+k),
by uzyskać c z przedziału [0,a+k]. Brzegowe wartości oznaczają dzielnik. Przyrost r(a) należy modyfikować, np. odejmując iloraz nadmiaru przez przyrost k .

Zobaczmy na przykładzie.
Liczba 8934053 ma przedstawienie a*p+b = 91 * 98176 + 37,
r bliskie 98176 / 91 jest równe np. 1078.
Przejdziemy z a na 97, przyrost 6. Zatem najpierw wyznaczamy przybliżenie podstawy systemu q = p-6r = 98176 - 6*1078 = 91708.
Teraz dokonujemy właściwej konwersji, co sprowadza się do obliczenia
c = 37 + 91*1078*6 - 6*91708 = 38377
W tym miejscu mamy już a zwiększone o 6 (odejmowanie 6*q oznacza a+6),
dokonujemy korekty, w czasie której uzyskujemy nadmiar dzieleniem
38377 / 97 = 395. 
Nasza konwersja przeskoczyła o 395 za daleko. Modyfikujemy r o 395/6 oraz wracamy by c było nie większe niż a+k:
38377 - 97*395 = 62
wtedy q = 91708 + 395 = 92103.
Mamy przedstawienie liczby n = (a+6)*q+c = 97 * 92103 + 62

W ten sposób pamiętając r można przebiegać po różnych a testując kandydatów na dzielniki - liczby pierwsze a.  Po każdym przekształceniu należy jednak zmodyfikować r, tu r przechodzi na 1078 - 395/6 = 1012, co tylko przybliża p/a = 949, lecz jest wystarczające do dalszych obliczeń.

Jeśli uzyskamy wartość c>q, możemy ją zignorować. Korekta to naprawi pozostawiając właściwą wartość.

Nowa metoda faktoryzacji 'trial division'

Cechą podzielności n przez liczbę a jest możliwość zapisania n = a*p+a.
Zastanawiało mnie, czy można łatwo przechodzić z liczby postaci a*p+b do (a+1)*q+c, gdzie q jest bliskie p zaś b<2a, c<2(a+1).
Można, co doprowadziło do uzyskania kolejnego siłowego algorytmu faktoryzacji typu 'trial division'.
W postaci n=a*p+b największą trudność wydaje się zachowywać b w przedziale (0,2a), lecz gorsze jest znajdowanie kolejnej podstawy systemu q. Znajdujemy ją za pomocą dzielenia n/(a+1), co nie zawsze trafia we właściwą wartość.

Oto jedna z wersji tego algorytmu rozkładu liczby n:
a = b = 1; p=n-1;
while( a<p ) {
  a++;
  b = b-p%a;
  p = p-p/a;
  if( b<a ) { p--; b = b+a; }
  if( a==b ) return dzielniki a, p+1
}
return liczba pierwsza;

Sposób ten niezwykle szybko pozwala znajdować kolejne reszty b, ale tylko dla sąsiedniej wartości a.
Oto początek rozkładu liczby 8934053:
a  b  p
1  1  8934052 
2  3  4467025  (4467025:3=1489008, 4467025-1489008-1 = 2978016)
3  5  2978016  (2978016:4=744504, 2978016-744504 = 2233512)
4  5  2233512  (2233512:5=446702, 2233512-446702-1 = 1786809)
5  8  1786809
6 11 1489007
7  9  1276292
...

Wartości p początkowo szybko maleją dążąc do pierwiastka kwadratowego z n, potem różnica między p na wejściu oraz wyjściu z pętli zmniejsza się do jedynki. Nie jest znana dokładnie, lecz błąd jest rzędu a, czyli stosunkowo mały. Dzielenia są prostsze niż w standardowym algorytmie dzielenia przez kolejne liczby.

Jeśli chcemy mieć mniejsze reszty b, należy zapamiętać wartość c=b z początku pętli oraz zmienić warunek pod if na tę kopię if(c<a)... Wtedy zamiast
5 8 1786809
mamy
5 3 1786810
Przebieg ten sam. Algorytm rekursywny (nie równoważny powyższemu, robiący to samo)
factd(a,p,b) {
  if( a==b ) return dzielniki a, p+1;
  if( a>p ) return liczba pierwsza;
  q=p-p/(a+1);
  c=p%b;
  if( b<a+1 ) {q--; c=c+a+1; }
  return factd(a+1, q, c);
}

Liczby złożone nieparzyste charakteryzują się nieparzystymi a i b=a oraz parzystymi p. Pojawia się chęć przeskoków tylko po nieparzystych a, jak w algorytmie standardowym. Jest to możliwe, lecz wtedy należy użyć konwersji. Tracimy łatwość wyznaczania reszt b. Pojawia się dodatkowe mnożenie, oraz konieczność korekty, gdyż b ma często zbyt dużą wartość. Wyliczane kolejne ilorazy nie są równe. Tworzą zaszumiony ciąg statystycznie nierosnący (fragmenty ciągu kolejnych ilorazów mogą być postaci: ..., 273, 270, 268, 265,..., 2, 1, 1, 1, 2, 1, 1, ).