Faktoryzacja z układu niezmienników

Kilka lat temu zauważyłem, że faktoryzację można przeprowadzić rozwiązując układ równań. Rozwiązywanie jest dwuetapowe, i wymaga znalezienia specjalnej postaci liczby zapisanej jako wielomian nad systemem pozycyjnym.

Generalna uwaga: na pozycjach cyfr występują liczby całkowite, nie tylko dodatnie, ale i ujemne. Teraz wzmocniłem nieco tamto podejście, wtedy potrzebowałem liczby 'pięciocyfrowej', teraz chcę pokazać niezmienniki pozwalające uzyskać rozkład z liczby 'czterocyfrowej'.

Odwracam zwykłe pisemne liczb trójcyfrowej oraz dwucyfrowej:
  [a,b,c]_{p} * [d,e]_{p} = [a*d, a*e+b*d, b*e+c*d, c*e]_{p} = [A,B,C,D]_{p}   (1)
Liczba jest w systemie o podstawie p ma postać [La, Lb, Lc, Ld]_{p}

Pierwsze niezmienniki dotyczą przeniesień:
a*d + k1 = A+k1 = La
(a*e+b*d) + k2 - p*k1 = Lb
(b*e+c*d) + k3 - p*k2 = Lc
c*e - p*k3 = D - p*k3 = Ld

Kolejne przekształcenia kierują do specjalnej budowy liczby, w której wyrazy skrajne [A, B, C, D]_{p} oddziaływują na wyrazy środkowe według formuły (1), czyli  wyrazy B, C zawierają jako składniki wielokrotności dzielnika A, oraz wielokrotności dzielnika D; dodatkowo składniki te mają występować w tych samych proporcjach, związanych z szukanym czynnikiem [d,e]_{p}. Podczas uzgadniania wyrazów skrajnych, wyrazy środkowe mają dodatkowe składniki w proporcjach takich samych jak kolejne pary. 

Jeśli wartości są dobrane poprawnie, konwersja na system o podstawie {p+1} postaci:
[a,b,c]_{p} = [a, b-2*a, a-b+c]_{p+1}
[d,e]_{p} = [d, e-d]_{p+1}
pozostaje niezmiennicza względem wartosci liczby. Nie stosuję tu konwersji całej liczby postaci
[a,b,c,d]_{p} = [a, b-3*a, c-2*b+3*a, d-c+b-a] _{p+1}
co pomaga lepiej szacować wartości wyrazów A, B, C, D.

Chcę pokazać budowę tych specjalnych postaci liczb dla dwu wartości. Dodatkowo przedstawię je w takim systemie, by maksymalnie uprościć współczynnik b, dokładniej on się tu zeruje. W tym celu użyję podstawy, przy której oba dzielniki są bliskie pierwiastka kwadratowego, czyli prawidłowy zapis liczby w systemie pozycyjnym jako najmniejszej liczby 'pięciocyfrowej' skracam do czterech cyfr. 

7387 = [1,1,1,1,7]_{9} zapisane jako [10, 1, 1, 7]_{9}.
Tu dodatkowo dochodzą własności, że a=1 oraz d jest duże. Docelową postacią jest
[9, 8, 18, 16] _{9}
Z tej postaci łatwo odczytujemy dzielniki [9,8]_{9} * [1, 0, 2]_{9} = 89*83

7031 = [9, 5, 7, 2]_{9}
Postaciami są [9, 8, -18, -16]_{9}  oraz  [8, 7, 64, 56]_{9}.
Dzielnikami są: 7031 = [9, 8]_{9} * [1, 0, -2]_{9} = [1, 0, 8]_{9} * [8, 7]_{9} =  89*79

Podstawowym celem jest dopuszczenie w myśleniu takiego zapisu liczby, zaś same przekształcenia są dodatkiem kierujacym się na wyłapywaniu wspólnych składników na różnych pozycjach. Np. moim pierwszym krokiem było zwiększenie wartości, by można było przenosić całe wielokrotności p. Co okazywało się fałszywym tropem, gdyż zaraz trzeba było zmniejszać dla uzgadniania proporcji, już o p.

Do pełnej postaci potrzeba znać dzielniki. Lecz już cząstkowe uzgodnienia, zwłaszcza 'cyfry jedności' sugeruje postać dzielnika. Dla liczb pierwszych uzgodnienie jest zawsze psute jakąś drobną różnicą. Dla liczb złożonych już często sprowadzenie dwu 'cyfr' (setek i jedności) do tej samej większej niż podstawa wartości powoduje, że największy wspólny dzielnik pozostałej pary 'cyfr' (tysięcy i dziesiątek) jest większy niż 1. Przypadek czy test pierwszości? 

Algorytmy faktoryzacji, kod w Pythonie

Dziś przedstawię nieco kodów. Python, liczba na początku oznacza poziom wcięcia.

Pierwszy kod to przechodzenie po kolejnych wielokrotnościach kolejnych liczb nieparzystych będących większymi niż liczba rozkładana. Można je podejrzeć wypisując zawartość linii assert. Kod zwraca tablicę wszystkich dzielników pierwszych rozkładu. 

0def factorize(n):
1 divisors = []
1 while 0==(n%2):
2    divisors.append(2)
2    n//=2
1  while 0==(n%3):
2    divisors.append(3)
2    n//=3
1  p = 3
1  d = 3-(n%3)
1  r = (n+d)//p
1  while 1<n:
2    a = d+(r<<1)
2    p += 2
2    d = a%p
2    r = r-(a-d)//p
2     assert(d,". ",p,"*",r,"=",n+d)
2    if 0==d:
3      while 0==(n%p):
4        divisors.append(p)     
4        n//=p
3      if p*p>n: # no prime divisors less than sqrt n
4        break
2    if p>=r: # finish due p*r = p*p > n
3      break
1  if 1<n: divisors.append(n)
1  return divisors

Obserwując działanie (assert = print), początkowo następuje szybkie zbliżanie się wartości p oraz r. Z czasem jednak następuje tłumienie, pierwsza różnica ukryta w wyrażeniu (a-d)/p maleje geometrycznie do 2. Dla dużych wartości odległość między p oraz r potrafiła być gigantyczna, zaś ich zbliżanie się już było rzędu dziesiątek, czyli sprawdzanie, sprawdzanie, sprawdzanie... Tylko, operacje asymptotycznie dążą do działań na liczbach rzędu pierwiastka kwadratowego z liczby rozkładanej n.

Algorytm cały czas sprawdza kolejne wartości nieparzyste, nawet gdy wiem z własności zachowania się (lokalna monotoniczność przedziałami), że dla tych wartości w otoczeniu nie ma dzielników.

Zatem lepiej jest podmienić ten algorytm dla r<p^2 na inny, w którym można wykorzystać informacje o braku dzielników w sąsiedztwie. Wybór padł na konwersje systemów niedziesiątkowych. Konwersja poniższa jest uogólnieniem konwersji z pracy "Conversion of Number System and Factorization" dla liczby trójcyfrowej [a,b,c]_p. 

#conversion number [a,b,c]_{p} into [x,y,z]_{p+k}
0 def conversion(a,b,c,p,k):
1  b = b-k*a
1  c = c-k*b
1  b = b-k*a
1  p += k
1  d=0
1  if 0>c:
2    d = 1+(-c//p)
2    c += d*p
1  elif c>p-1:
2    d = -(c//p)
2    c = (c%p)
1  if c==p:
2    d=-1
2    c=0
1  b -= d
1  d=0
1  if 0>b:
2    d = 1+(-b//p)
2    b += d*p
1  a -= d
1  return [a,b,c,p]

Jest tu jedna instrukcja if, która wydaje się być nie na miejscu, gdyż warunek jest spełniany if wcześniej. W praktyce u mnie bez tego ostatniego if c==p, wartość c okazywała się być równa p wbrew wcześniejszemu c = (c%p).Kod nie zawsze działa poprawnie dla większych k - przyczyną może być ignorowanie "naprawy" przeniesieniami wartości b po pierwszej instrukcji, kiedy b staje się ujemne.

I okazuje się, że dla bardzo, bardzo dużych wartości można zmniejszyć krotność rozpatrywanych iteracji, i to nawet z dwu stron. Trzeba robić to jednak bardzo delikatnie. W pierwszym przypadku reszty z dzielenia przez kolejne wartości tworzą ciąg lokalnie monotonicznie rosnący, w drugim lokalnie monotonicznie malejący. Zaś możliwość wystąpienia dzielnika jest tylko przy tych wartościach, w których następuje przeniesienie.
Sprowadza się to do uważnego śledzenia, kiedy warunek wyrażenia b-2*a*k oraz wartość bk po konwersji o k rozjeżdzają się. Dokładniej, jeśli e = (b-2*a*k)-bk przy spełnieniu pewnych warunków (nie zawsze można wziąć b z marszu), to mamy dla zapisu liczby n w przedziale podstaw systemu [p, p+k] wystąpienie dokładnie e przeniesień. Czyli dla końcówki algorytmu -- z tysięcy wartości wystarczy zlokalizować i sprawdzić kilka. I to w obszarze bardzo obciążającym obliczeniowo procesory działające innymi metodami.