Kod faktoryzacji przeglądem zupełnym, zakres podstaw większy niż sqrt n

Kod w poprzednim poście można przedłużyć, i teraz zaprezentuję to przedłużenie. To też jest algorytm faktoryzacji, tutaj "odwracający" kierunek sprawdzeń w tym sensie, że współczynnik a maleje od pierwiastka kwadratowego z liczby rozkładanej n do 1. 

Kod ma 'tylko' pierwiastek kwadratowy z n iteracji. Jednak w testach okazał się wolny. Zaprezentuję go z funkcją użytą w poprzednim kodzie, conversion1(n,p), która wyraża liczbę w postaci c*p*p+a*p+b = n dla c=0. Chociaż dla tych wartości lepiej użyć modyfikacji pobranych z rozwinięcia (a-1)*(b+1) = a*b+b-a-1. Dzielenie będzie zastąpione dodawaniem i odejmowaniem.
Jeszcze jedna uwaga, użyta funkcja floorSqrt(n) zwraca całkowito-liczbowy pierwiastek kwadratowy oraz jego różnicę od wartości, tzn. jeśli f,g = floorSqrt(n), to n = f*f+g, 0 <= g < 2*f+1. 

Numery linii oznaczają wcięcia, jak przy bardziej dopracowanych, publikowanych przeze mnie kodach w Pythonie 2.

0def factor4(n):
1  a,b = floorSqrt(n)
1  if 0==b and n==a*a: return [a,a]
1  p = a
1  k=1
1  while 1<a:
2    if b<a:
    #  if 1==p%2: #dzielnik bywa pomijany np. 55
    #    p+=2
    #    k=2
    #  else:
3      p+=1
3      k=1
3      c,a,b = conversion1(n,p)
2    else:
3      k = b/a
3      if 0==b%a:
4        p += k
4        q = n/p
4        print "dzielnik p",p, q
4        return [p,q]       
3      p = p+k+1
3      c,a,b = conversion1(n,p)
2    print "[",a,",",b,"]_",p,"   ",k  
1  return [n]
 

Zakomentowana instrukcja if miała zmiejszyć krotność iteracji dla wartości nieco większych niż sqrt n, by sprawdzać tylko nieparzyste, ale okazało się, że niektóre iloczyny nie są wtedy rozkładane. Taką liczbą jest 55 = 5*11.

Kod ma kolejną wadę opóźniajacą faktoryzację: zwraca wartość będącą pierwszą [n], lub dwie [p,q], ale nie wiadomo, czy p oraz q są pierwsze czy złożone. Zatem funkcja wywołująca factor4(n) musi ją wywołać ponownie dla tych wartości. Przy programowaniu współbieżnym wyniki można odkładać do tablic, i rzut wskaźnikiem (czytaj: rzut okiem) na długości tych tablic wskazuje, czy dzielniki są znalezione, jakie są, ale to nie będzie pełny rozkład. 

Tak też można...

Kod faktoryzacji szukaniem przeniesień podczas konwersji

Udało mi się zbudować kolejny kod faktoryzacji, w którym lokalizuję i sprawdzam przeniesienia występujące podczas konwersji na coraz większe podstawy. Dla małych wartości ma narzut obliczeniowy, który jednak jest kompensowany dla liczb kilkunastocyfrowych. Liczba, którą innym przeglądem zupełnym rozkładałem (i zatrzymywałem po kilkunastu godzinach gdy podstawa systemów przekroczyła milion), tą wersją uległa rozłożeniu na czynniki w 2h 13 minut. 

Program nie ma już pętli wykonującej się (sqrt N) /2 razy, przeskakuje część iteracji, gdy wiadomo, że nie ma przeniesienia podczas konwersji, a zatem nie ma kandydata na dzielnik. Ile razy wykona się pętla? Nie wiem. Są obszary, kiedy w sekundy obrabia setki tysięcy podstaw w sekundę (liczba 16-cyfrowa), a są inne, w którym przeniesienia podczas konwersji są co kolejną podstawę. 

Kiedy zażyczyłem sobie, żeby program starał się jednak zwiększać podstawę o co najmniej 2, mógł przeskoczyć dzielnik, ale kiedy dołączyłem warunek parzystości podstawy, ta sama liczba co wspomniana już została rozłożona na czynniki w godzinę 8 minut... Prawie 50 procent szybciej.

Podaję kod 'jak powstał', w którym conversion1(n,p) polega na przedstawieniu liczby w systemie o podstawie p, zaś dalej(...,p,k) to konwersja na system p+k bez wykonania przeniesień ([a,b,c] => [a, b-2k*a, c-k*b+k*k*a]), by je lokalizować w tym przedziale. Na podstawie tych postaci przekształconej liczby dobieramy skok k. Można go jeszcze usprawnić... 

def factorize(n):
  divisors = []
  while 0==n%2:
    divisors.append(2)
    n/=2
  p=p1=3
  a1,b1,c1 = conversion1(n,3)
  k=2
  a=0 #,b1,c1,p1 = conversion(a,b,c,p,k)
  while a1!=a:
    a,b,c,p = a1,b1,c1,p1
    if 0==c:
      if 0==n%p1:
        divisors.append(p1)
        n/=p1
        a1,b1,c1 = conversion1(n,p1)
        continue
      c = n%p1
      b = ((n-c)/p)%p
      a = (n-b*p*a)/(p*p)
    p1 = p+k
    a1,b1,c1 = conversion1(n,p1)
    #print "[",a,",",b,",",c,"]_",p,"   ",
    print "[",a1,",",b1,",",c1,"]_",p1         
 
  print "dalej teraz a maleje do 1"
  print "[",a1,",",b1,",",c1,"]_",p1,"  druga czesc"
  if b1<a1:
    a1-=1
    b1+=p1
  while 0<a1: #szukanie binarne k
    k3 = b1/a1/2
    k1,k2 = 0,k3     
    a3,b3,c3,p3 = a1,b1,c1,p1
    a,b,c,p = a1,b1,c1,p1
    while k1<k2-1:
      k = k1+(k2-k1)/2
      a,b,c,p = dalej(a1,b1,c1,p1,k)
      p3 = p1+k
      a3,b3,c3 = conversion1(n,p3)
      #print "[",a,",",b,",",c,"]_",p,"   ",k1,k,k2,"  ",
      #print "[",a3,",",b3,",",c3,"]_",p3,"   "
      if b3>a3 and 0<c:
        k1=k
      else:
        k2=k
    print "[",a1,",",b1,",",c1,"]_",p1,"   ",
    print "[",a,",",b,",",c,"]_",p,"   ",k1,k,k2,"  ",
    print "[",a3,",",b3,",",c3,"]_",p3,"   "
    if b3<a3:
      a1,b1,c1,p1 = a3-1,b3+p3,c3,p3
    else:
      if 2>k and 1==p1%2:
        k=2
      p1 = p1+k
      a1,b1,c1 = conversion1(n,p1)
    print "po [",a1,",",b1,",",c1,"]_",p1
    if 0==c1:
      while 0==n%p1:
        divisors.append(p1)
        n/=p1
      a1,b1,c1 = conversion1(n,p1)
     
  if 1<n:
    divisors.append(n)
  return divisors

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.