Problem Goldbacha, spostrzeżenia

Goldbach pyta Eulera w liście: Czy każda liczba parzysta jest sumą dwu liczb pierwszych? 

Euler nie potrafił odpowiedzieć. Problem trafił na sumę problemów do rozwiązania. 

A oto kilka moich spostrzeżeń: suma liczb pierwszych większych niż 2 jest zawsze parzysta. A czy liczba parzysta jest sumą liczb pierwszych? Wiele prób, choć już pojawił się dowód słabszej wersji, wspomniany w Wikipedii pod hasłem 'hipoteza Goldbacha'. 

Każdą liczbę nieparzystą, w tym liczby pierwsze różne od 2 także, można przedstawić jako różnicę kolejnych kwadratów:
f(p) = a*a - (a-1)*(a-1)
I tak f(1) = 1*1 - 0*0 = 1 (dla Golbacha -- u niego 1 była liczbą pierwszą)
f(2) = 2*2 - 1*1 = 4-1 = 3
f(3) = 3*3 - 2*2 = 5
f(4) = 4*4 - 3*3 = 7
f(5) = 5*5 - 4*4 = 9
i tak dalej... 

Dodajmy do siebie:
f(a) + f(b) = a*a - (a-1)*(a-1) + b*b - (b-1)*(b-1) =
(a-b+1)*(a+b-1) + (b-a+1)*(a+b-1) =
2*(a+b-1)
Parzysta... Według oczekiwań, lecz jest postać tej liczby parzystej.

Poszukajmy okresu, czyli orbit jednego dzielnika:
jeśli k dzieli f(a), to k dzieli f(a+k):
f(a+k) + f(a) = 2*(a+k+a-1) \equiv 2*(2*a-1) = f(a)+f(a) \equiv 0 modulo k. 

Kiedy f(a) jest liczbą pierwszą; może lepiej spytać: kiedy f(a) NIE JEST liczbą pierwszą:
dla 3 = f(2): f(5) = 9, f(8) = 15, f(11) = 21, ...
Pierwsza wartość jest, powstała orbita f(5+3N) to liczby złożone.
Dla 5 = f(3): f(8), f(13), f(18), ogółem f(8+5N) to liczby złożone.
Szybko: f(4+7+7N) = f(11+7N) to liczby złożone..., kwadrat 7 to f(25).
Podobnie dla kolejnych liczb pierwszych, kwadrat jakiejś liczby generuje orbitę liczb złożonych, nawet z kilkoma wartościami wcześniej, te dodatkowo wpadają w orbitę mniejszej liczby pierwszej. 

Niezwykle przypomina to sito Eratostenesa, chociaż mamy więcej liczb pierwszych niż w sicie, zatem własności sita Eratostanesa są ograniczeniem górnym. Liczb pierwszych wyznaczonych przez f() jest nieskończenie wiele, nawet więcej niż oszacowanie górne n/(lg n), ale w granicy... Interesują nas liczby spoza orbit, czyli pierwsze.

Przejdżmy do pytania Golbacha w mocnej wersji. Mam liczbę parzystą 2n, czy jest ona sumą dwu liczb pierwszych?
Wystarczy dopasować takie a oraz b, by a+b=n+1 oraz by ani f(a), ani f(b) nie znalazło się na żadnej z orbit wyznaczonych przez liczby pierwsze przez f(). Orbity są na tyle 'dziurawe', że raczej znajdzie się miejsce na parkę, a może nawet więcej możliwości. 

Jakie jest prawdopodobieństwo, że jedna z ostatnich wartości, np. n-1, jest liczbą pierwszą? Oszacujmy prawdopodobieństwo, że nie jest pierwszą, że jest liczbą złożoną.
Po pierwsze, przesunięcie dla liczby pierwszej p. Orbity zaczynają się od pozycji f^{-1}(p)+p, czyli około p/2+p = 3/2 p. W szukaniu prawdopodobieństwa możemy pominąć 1/3 największych liczb pierwszych mniejszych niż n.
Po drugie, liczby pierwsze a, b są względnie pierwsze, oraz sumę prawdopodobieństw z ich orbit można łatwo wyznaczyć: 1/a + 1/b - 1/(a*b) = (b+a-1)/(a*b) dla wartości podzielnych przez obie. Dotyczy to także wartości względnie pierwszych. Orbity od liczb złożonych zawierają się w orbitach generowanych przez dzielniki liczb -- są już uwzględnione. Tak począwszy od k=8 mamy 1/3 + 1/5 - 1/15 = 7/15, że liczba n-1 jest podzielna przez 3 lub 5. Dołączając orbitę z 7 mamy 7/15 + 1/7 - 1/105 = 63/105. Dołączymy orbitę z kolejnej liczby pierwszej, i znowu... i kolejną...
Ciąg prawdopodobieństw uzyskiwany z kolejnych orbit jest rosnący, lecz ograniczony, bo liczb pierwszych w granicy jest nieskończenie wiele, i musi się dla nich znaleźć pozycja. Prawdopodobieństwo trafienia w liczbę pierwszą pozostaje dodatnie.

Suma odwrotności liczb pierwszych jest szeregiem rozbieżnym (Delta wrzesień 1999), chociaż tu mamy ograniczenie wielkości składników do mniej niż 2/3 n oraz dodatkowe zmniejszanie. Zatem wszystko się może zdarzyć... I dlatego problem stawia opór.

Np. 2n = 156: a+b=78+1, suma liczb pierwszych 156 = 17+139 czyli f(9)+f(70), zaś f(a) dla a<9 niemal zawsze wpada we wspomniane orbity. Wyjątkowe jest 13: dzielnik 156, które mi namieszało w poszukiwaniach sugerując siebie.

Wskazówka do zmniejszania liczności sprawdzanych przypadków podczas faktoryzacji

Dotychczas moje algorytmy przekształcały rozkładaną liczbę. Aby podliniowo faktoryzować, należy skorzystać z właściwości dzielników. Można je uzyskiwać nawet nie znając tych dzielników. Zainteresowych odsyłam, by zapisać kolejne cyfry liczby w innych systemach, oraz sposób odwrócenia kolejności dwu ostatnich podstaw, co już daje nieco informacji o dzielnikach...

Pomocą może być także funkcja wskazująca rozmieszczenie liczb pierwszych, bliźniaczych, czworaczych, opublikowana już na tym blogu. Rozszerzenie nie musi być o kolejną liczbę pierwszą, lecz o dowolną liczbę pierwszą. 

Liczba jest pierwsza, gdy wspomniana funkcja przyjmuje wartość 1 dla każdego rozszerzenia mniejszego niż dana liczba. Pojawiającą się drobną nieoznaczoność przy nie kolejnych liczbach pierwszych można łatwo namierzyć, np. dla k=6 i rozszerzaniu o liczbę p będzie sugerowana własność pierwszości dodatkowo tylko dla p oraz 5*p, gdy nie uwzględniamy liczby pierwszej 5. Dla p pierwszej wartość funkcji tych miejsc ulega modyfikacji z 1 do p. Funkcja z liczby ma wartość a>1, gdy liczba ta jest podzielna przez a.

A kiedy weźmiemy za mało danych, do gry wchodzi chińskie twierdzenie o resztach, dzięki któremu obszar poszukiwań dzielników drastycznie się kurczy. 

Wystarczy funkcja dla p \in {2, 3, 5, 7}, a znajdziemy rozkład 7387, mając reszty modulo: 6, 30=6*5, 42=6*7. Te same dane (już bez chińskiego twierdzenia o resztach) wystarczyły na wykazanie pierwszości 127. I tylko warto przyjmować rozkład liczb ujemnych, np. testując -5 = (-1)*5 lub -5 = 1*(-5), i jeszcze parę innych związanych zależnościami, by dopasować iloczyny reszt dzielników. Unikamy w ten sposób testowania wielu wartości. 

Nie chcę podawać tu szczegółów - to jest za szybkie na obecne zabezpieczenia cyberbezpieczeństwa. Przy rozkładzie 7387 użyte tylko cztery liczby pierwsze spośród 23 kandydatów już wskazały rozkład... Wystarczyło uzgodnić 30*a+r1 = 42*b+r2, gdzie a, b małe liczby; r1, r2 pozycje funkcji z wartością 1 dla parametrów {2,3, 5 | 7}. w której r1,r2 przystają do wyznaczonej reszty, np. -5 = 37-42, zaś 37 jest resztą iloczynu a*b (mod 42). Reszta będąca liczbą pierwszą to nie jest przypadek...

 

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

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.

Wprawka ćwiczeniowa o większych możliwościach, test pierwszości

Geeks for geeks opublikowało programy do przedstawienia liczby jako sumy kolejnych liczb 

Count ways to express a number as sum of consecutive numbers, 

Jest pewien szczegół techniczny. Ta metoda dla liczb nie będących potęgami 2 stanowi nie tylko test pierwszości, ale i pozwala znaleźć dzielnik liczby. 

Liczby pierwsze mogą być przedstawione w postaci takiej sumy na dokładnie dwa sposoby. W zakresie sprawdzania, sposób jest tylko 1. Liczby złożone można przedstawić na wiecej sposobów w postaci takich sum. Zatem, jeśli wynik programu da
0, liczba jest potęgą 2,
1, liczba jest pierwsza,
2 lub więcej, liczba jest złożona. 

Ostatni algorytm o złożoności O(sqrt n) jest 'szybkim' testem pierwszości. Algorytm o złożoności O(n) korzystający z dwóch wskaźników left, right, wskaże dzielnik jako jakąś podwielokrotność różnicy
right-left.

I tylko szybkość algorytmu dla gigantycznych liczb może nieco nużyć. Ostatni można zakończyć szybciej niż po znalezieniu wszystkich sposobów,  

Jest to praktyczne dopracowanie problemu 'własność liczb złożonych' z 25 kwietnia 2018 roku

Przegląd zupełny z użyciem systemu Fibonacciego

 Spróbowałem napisać faktoryzację przeglądem zupełnym z użyciem systemu Fibonacciego. Udało się, choć mam jakąś dziwną usterkę.

Przypomnę zasadę tego typu: liczbie n zapisanej pierwotnie w systemie Fibonacciego zwiększam wartości 'cyfr' usiłując uzyskać tylko dwie wartości: liczbę p albo 0. Jeśli się uda, p jest dzielnikiem liczby n, jeśli nie, zwiększam p oraz za pomocą przekształceń niezmienniczych wartości liczby modyfikuję wartości stojące na miejscu 'cyfr'. n = [p,0,0,p,0,....,p,0, 0 albo p] = p *coś

Mam cztery przekształcenia niezmiennicze:
1) [a,0,0] = [0,a,a]
2) [0,a,a] = [a,0,0]
3) [0,2*a,0,0] = [a,0,0,a]
4) [0,0,3*a,0,0] = [a,0,0,0,a]
uważając na przypadki szczególne pozycji najmniej znaczących. 

Program główny. Najpierw zapisuje liczbę w systemie Fibonacciego jako ciąg jedynek oraz zer, po czym dla danego p porównuje 'cyfrę' z p.
Pętla zewnętrzna zwiększa p o 1, oraz od 'cyfr' najbardziej znaczących modyfikuje wartości podanymi schematami. Pięknie działa... dla małych wartości n. 

A dla dużych? Pojawia się dosyć specyficzny błąd wyprowadzenia liczby z postaci, z której można odczytać dzielnik, Co nie znaczy, że wartość uległa zmianie. Z niejasnych dla mnie powodów nie można zastosować zwykłej pętli typu
for c in range(len(liczba), 0, -1) {...}
gdyż czasem pojawia się nieco inna wartość na pozycji którejś z cyfr. Zamiana pętli na równoważną:
c = len(liczba)
while( c>0 ) {
c-=1
{...}
}
powoduje, że wyświetlany podgląd 'cyfr' liczby może, choć nie musi ominąć ten błąd prowadząc do innego wyniku. wbrew definicji algorytmu jako 'jednoznacznego sposobu uzyskania wyniku z danych wejściowych'. Zawsze błąd pojawia się na pozycji 3 najmniej znaczącej (odpowiednik 'tysięcy'), kiedy reszta jest całkowitą wielokrotnością 3. 

Praktyczny sposób postępowania w pętli wewnętrznej, kolejność dla a:
jeśli cyfra to 0, następne przekraczają p, stosuję 2) z a=p;
jeśli cyfra większa od p, do nadmiaru a stosuję 1), nawrót do wcześniejszej cyfry
cyfrę mniejszą od p uzupełniam, najchętniej do p, stosując 2) do minimum z dopełnienia do p, kolejnych cyfr; po wyzerowaniu którejś z kolejnych zadziała 3) albo 4); 
jeśli da rady uzyskać na danej pozycji wartości p, tylko mniejszą, przenoszę całość korzystając z 1);
przechodzę do kolejnej, mniej znaczącej cyfry;
kiedy docieram do najmniej znaczącej, powinna przyjąć jedną z dwu wartości: 0 albo coś nie większego od p. Jeśli to jest 0 albo p, pojawia się sugestia: 'p jest dzielnikiem'.

W programie w przypadku sugestii znalezienia dzielnika jako 'cyfry najmniej znaczącej' należy ją sprawdzić, np. sprawdzając krotność róźnych wartości 'cyfr'. Jeśli są trzy, ta różniąca się nie musi występować tam gdzie chcę ją sprawdzać, i iteracja musi się powtórzyć dla tego samego p by ją (już zmienioną) wstawić w odpowiednie miejsce dla odczytania wyniku. 

Istnieje bardzo powolna tendencja zmniejszania długości ciągu potrzebnego na zapis 'cyfr'. 

Kryterium stopu jest pierwiastek kwadratowy z liczby, bo okazało się, że dla liczb pierwszych jest co prawda szybsze kryterium stopu, spełniane jednak często przez liczby złożone. Zatem nie do użycia.

Wielomianowy algorytm faktoryzacji, rodzina odwracania działania

Mam już napisany program do faktoryzacji z rodziny rozwiązań odwracających działania. Nie wydaje się mi się bardzo szybki, ale to dlatego, że każę mu wypisywać w każdej iteracji drzewo możliwości, aby móc obserwować jego działanie. To wypisywane jest właśnie bardzo wolne. 

Algorytm do działania potrzebuje listy początkowych liczb pierwszych primes = [3,5,7,...] oraz iloczynów początkowych liczb pierwszych product = [6,30,210,...] = [2*3, 2*3*5, 2*3*5*7, ...]. Liczba pierwsza 2 nie jest uzwględniona w tablicy primes, ale w product już tak.
Tablice te w razie potrzeby trzeba rozszerzać, choć dla testowanych przypadków kilkunastucyfrowych wystarczy naledwie kilka liczb pierwszych oraz ich iloczynów. Rozszrzerzanie tych tablic w razie potrzeby o kolejną liczbę pierwszą, jej iloczyn z pozostałymi, stanowi część funkcji.

Zasada iteracji jest następująca: na drzewie mam już odłożone iloczyny a*b mające tę samą resztę co product[i-1]. Pobieram jeden. Rozszerzenie do iloczynu product[i] o liczbę pierwszą p ma specyficzną budowę:
a = a+k*product[i-1],  k=0..p-1
b = b+j*product[i-1],  j=0..p-1
Jest ich p^2. Spośród nich wybieram te a*b, które mają tę samą resztę co liczba rozkładaną n modulo product[i], oraz nie mogą być zbyt duże. Te odkładam do drzewa (a,b,i), chyba, że iloczyn pokrywa sie z n, wtedy następuje koniec funkcji, zwrócenie wyniku n = a*b. W drzewie trzymam dodatkowo indeks i, będący zarazem poziomem drzewa. Przycinanie przepuszcza maksymalnie do p wartości, i to tylko na niższych poziomach. Inne wersje algorytmów tej rodziny mogą mieć mocniejsze przycinanie, kosztem poziomów. Po przekroczeniu pewnej wartości pojawia się bardzo mało nowych węzłów wstawianych do drzewa. 

Inicjacja, gdy n modulo 6 jest 1, odkładam dwie pary (1,1, indeks 1), (5,5, 1), gdy n modulo 6 jest 5 odkładam (1,5, 1).
Dla liczb kilkunastocyfrowych poziom drzewa nie przekracza 10, a kiedy pisałem ten post, zanim zapisałem pierwsze dwa akapity, już uzyskałem rozkład:
9198491200007 = 197 * 109 * 428374759
poziom drzewa czyli długość tablic pomocniczych nie przekraczał 8. 

Podstawowa trudność spowalaniająca przy odkładaniu a*b była tą, że te same obliczenia były wykonywane także dla równoważnej pary b*a. Niepotrzebnie. Po wyłapaniu tych przypadków program zyskał zauważalnie na szybkości. Oraz osobno sprawdzam podzielność przez 2, 3.
Program lepiej działa dla większych wartości. Liczby gładkie mają zbyt dużo równoważnych iloczynów, co powoduje mocny rozrost drzewa. Z kolei liczby pierwsze i półpierwsze są rozpoznawane znacznie szybciej. 

Nie wiem dlaczego, ale intuicja podpowiada mi, że dopasowywanie par (a,b,) warto zacząć od większych wartości, a zmniejszać w kierunku wartości z poprzedniej iteracji odejmując product[i]. Elementy product[] mają wzrost wykładniczy, 

Odwracanie systemów o podstawach będących liczbami złożonymi, jak np. dziesiątkowego, mają zawsze tendencję do powolnego rozrostu drzewa wydającego się wykładniczym. Jednak ten wzrost wykładniczy jest dla podstawy systemu, czyli nie przekracza m^(log _m n) = n, dla binarnego 2^(lg n), gdzie lg oznacza logarytm binarny liczby. Najmocniejsze przycinanie będące porównaniem reszt po konwersji do wielokrotności podstawy praktycznie blokuje rozrost wykładniczy w wielu przypadkach, zostaje kilka ścieżek prowadzących do różnych czynników. Z tego powodu ośmielam się nazywać tę rodzinę algorytmów wielomianowymi.
Praktycznie, kiedy ktoś mi podsuwa klucz publiczny RSA, jeżeli ostatnią 'cyfrą' jest 0xF, to nie ma innych możliwości modulo 8 niż iloczyny 0x3*0x5 albo 0x1*0x7.

Logarytm dyskretny, obliczanie

Logarytmy dyskretne wciąż są traktowane jako funkcje jednokierunkowe. A z tą jednokierunkowością... Poniżej przedstawiam program liczący potęgi na dwie różne permutacje modulo n. Przypominam, bo wiele lat temu to publikowałem, teraz mam gotowy kod. 

Moją słabością jest brak zapamiętywania - zawsze na nowo wyprowadzam, i okazywało się, że jestem w stanie łatwo pomylić ze sobą potęgę oraz podstawę potęgi w kodzie. Dlatego to co piszę o kodzie nie zawsze daje właściwy obraz. Dlatego wolę opisywać kod niż go prezentować. Tutaj prawie kopiuj - wklej gotowego.

g jest podstawą potegi, s = '1000..._g' jest wartością potegi,modulo n
Numer w linii oznacza poziom wcięcia kodu w Pythonie, standardowo ładuje mi się wersja 2.

0 def dysLog(g,s,n):
1  if s==g: return 1
1  r = n-g
1  c = 1
1  a = g
1  while 1:
2    a = -a*r
2    c+=1
2    #print a,r, dla zobaczenia co sie dzieje z danymi, funkcja % może dawać zły wynik
2    if 0>a:
3      b=(-a)//n
3      a+=b*n
3      a+=n
2    #print a koncowa wartosc
2    if a==s: return c
2    if a==g and c>n-1: return -c #nie ma rozwiazań, cykl z małego twierdzenia Fermata, osłabienie and do or wyznaczy długość cyklu
 
Dane np.:
g=11
n=17
0 for a in range(n):
1   print g,"^", dysLog(g,a,n), "=",a,
0 print
0 for i in range(1,n):
1  print i,":",(g**i)%n,"     ",
0 print

Rekursywna konwersja potęgi na system o podstawie n, w którym wystarczy obserwować wartość na pozycji najmniej znaczącej modulo n.

Dzielenie, gdy iloczyn jest obcięty modulo

 Wśród operacji arytmetycznych występuje taka, że iloczyn dwu liczb jest obcinany do reszty z dzielenia  przez jakąś wartość. Czyli brany modulo. 

Czy mając dane taki iloczyn oraz jeden z czynników, możemy wyznaczyć drugi? Odpowiedź jest twierdząca, a nawet łatwa do wyznaczenia, chociaż w ogólności moze być niejednoznaczna - może być więcej wartości, które przycięte dadzą tę samą wartość iloczynu. 

Przykład liczbowy, Iloczyn 65*x = 175 (mod 1000). Czyli należy podzielić 175 / 65 w liczbach całkowitych modulo 1000.
Zapiszmy 1000 w systemie o podstawie 65: 1000 = 15*65+25.
Teraz, aby uzyskać wartość podzielną przez 65, należy dodać tyle tysięcy, by uzyskana wartość była podzielna przez 65. Szukamy k, by
  175 + (25*k)
było wielokrotnością 65.
Szybko znajdujemy, że współczynnikiem k jest 6, czyli 6*1000+175 = 6175 jest wielokrotnością 65, dokładniej 6175 = 65*95. 

Znaleźliśmy: 175 / 65 = 95 (mod 1000). Dzielenie przeszło na dodawanie i sprawdzanie podzielności, co zapisane w systemie znanego dzielnika sprowadza się tylko do dodawań i przeniesień...

Pierwiastek kwadratowy z dużych liczb, różne podejścia algorytmów

Problem obliczania pierwiastka kwadratowego oraz zwrotu wartości całkowitej został zaproponowany na portalu GeeksForGeeks 12 sierpnia 2025 roku. Zaproponowano 4 algorytmy: naiwny z mnożeniem, wyszukiwanie binarne, oraz dwa korzystające z biblioteki wbudowanej, np. Python.math. 

Dołączyłem jeszcze dwa swoje: naiwny z dodawaniem +=(2p+1), wspomniany miesiąc temu przy podawaniu złożoności wielomianowego algorytmu faktoryzacji, oraz jego rozwinięcie wykorzystujące fakt, że każdy kwadrat mozna zapisać jako 1 0 0 w systemie o pewnej podstawie p. Kod ostatniego niemiły, zaczyna działać dopiero dla wartości większych niż 100 w systemie dziesiętnym.
10*p + k, gdzie k wyznaczane z krotności możliwości konwersji przycięcia liczby N_p -> N_{p+k} do odpowiedniej długości względem p, jedna cyfra p na dwie cyfry z N = N[n]...N[1] N[0] w dziesiętnym
2p+1 < 100*reszta+10*N[b+1]+N[b]

Użyłem bardzo dużej liczby: 27997833911221327870829467638722601621070446786955428537560009929
26128400107609345671052955360856061822351910951365788637105954482
006576775098580557613579098734950144178863178946295187237869221823983

Algorytmy miały działać w parę sekund - tego testu naiwe nie przeszły.
Te, które w teorii były najszybsze, złożoności stałej O(1) zarówno pamięci jak i czasu, czyli
int(math.sqrt(n))
int(math.exp(0.5*math.log(n)))
podały błędne wyniki -- wystąpiła propagacja błędu obliczeń liczb zmiennoprzecinkowych, uzyskałem np. wartość
52912979420196447 944286770683838014
89465575494900455707262162396795715834972575453540520946581897216
Na portalu wspomniano, że wynik może być niezbyt dokładny, dokładniej nieco zaniżony. A okazał się znacznie zawyżony.

Poprawny wynik dały zaledwie dwa algorytmy z sześciu bez czekania:
52912979420196447 553235400217339415 00571661886712719898054010160087776775983652376577898717296799172

Zgodność z 'najszybszym' tylko na siedemnastu cyfrach dziesiętnych. 

Jak to trzeba sprawdzać zakres stosowania bibliotek... 

Kodowanie Huffmana

 Po angielsku Huffman Encoding, jest to sposób zapisu symboli informacji w taki sposób, by zawatość zajęła mniej mniejsca niż początkowo.
Wydaje mi sie, że kodowanie to moze mieć jeszcze jedno zastosowanie w grach RPG bez przemocy. Mianowicie częstostliwość występowania symbolu może być uznana za koszt użytego... powiedzmy sobie -- czaru, ataku... którym nokautuje, pomaga się oponentowi, by przejść dalej.

Portal geeksforgeeks podaje algorytm zachłanny wyznaczania tego kodu. Proponuje kopce lub kolejki priorytetowe. Spróbowałem 'po swojemu', i teraz podam, że wystarczą listy oraz drzewa zapisane jako tablice. Najgorszy moment - w ostatniej fazie chciałem przenumerować węzły, ale okazało się to niepotrzebne. Odpowiednie drzewo było niemal gotowe. 

Opis algorytmu, który jest taki sam we wszystkich podejsciach.
Mamy tablice symboli oraz ich częstotliwość występowania. Sortujemy je rosnąco. Wybieramy dwa najmniejsze, i tworzymy z nich nowy węzeł, w którym częstotliwość jest sumą składowych. Powtarzamy tak długo, aż zostanie jeden węzeł, w którym częstotliwość jest równa sumie częstotliwości wszystkich symboli.
W drugiej fazie, dane te umieszczamy na drzewie. Przechodząc po tym drzewie prefiksowym po dotarciu do liścia uzyskujemy symbol zapisany tym kodem. 

W moim kodzie od razu generuję węzły drzewa, choć częstotliwość jako ostatecznie zbędna wprowadzam do drugiej pomocniczej tablicy węzłów wsp.
Zatem tworzę węzły z pary (symbol, częstotliwość, pozycja), gdzie pozycja to indeks węzła w tablicy, niezbędny, by dobrze je łączyć. Umieszczam je na liście tree. Kopiuję je do pomocniczej tablicy wsp. Warunek, nie mogę mieć symbolu 0 służącego jako flaga separująca liście od węzłów drzewa binarnego.

Wyszukuję wskaźnikami dwa węzły a, b o dwu najmniejszych częstotliwościach w tablicy wsp, tworzę z nich nowy węzeł o budowie (0, suma_częstotliwości, indeks), który dołączam do tablicy wsp, zaś równocześnie do tree dołączam węzeł o budowie (0, a, b), gdzie a i b są wyszukanymi indeksami węzłów. Z tablicy wsp usuwam węzły o indeksach a oraz b. Tak długość wsp w każdej iteracji maleje o 1.
Kontynuuję, dopóki we wsp nie zostanie tylko jeden węzeł, dołączany na końcu. Węzeł ten jest korzeniem drzewa, które już mam zaszyte w tablicy tree!

Druga faza to u mnie kosmetyka. We wszystkich węzłach z symbolami zeruję pola inne niż nazwa symbolu, przez co węzeł staje się liściem (symbol, 0,0), i drzewo gotowe! 

Korzeń jest w ostatnim wężle wsp postaci np. (0,x,y). Na pozycji x jest węzeł tree.left, na pozycji y węzeł tree.right, liście mają wyzerowane te odsyłacze. Odczyt preOrder(tree, len(tree)-1, s="") schodzi rekursywnie s+'0' dla left, s+'1' dla right. 

Wielomianowy algorytm faktoryzacji, jeden z wolniejszych

 W rodzinie faktoryzacji zwiazanych z rozwinięciem liczby N w szereg odkryłem sposób, który ma bardzo specyficzne wartości, i jest łatwy do wyznaczenia złożoności. 

Niech p[k] oznacza k-tą liczbę pierwszą. Współczynniki rozwinięcia w szereg są iloczynami początkowych kolejnych liczb pierwszych oraz 
2*3*5*...*p[k] <= N < 2*3*5*...*p[k]*p[k+1] .
Rozwinięcie w szereg polega na odejmowaniu kolejnych wielokrotności w[k] iloczynów początkowych liczb pierwszych, z których każdy kolejny powstaje przez mnożenie przez coraz mniej wartości, kolejne iloczyny powstają przy zmniejszaniu się k, przedostatni byłby 2*3=6, ostatni przez 2. 

Każdy taki iloczyn zabiera jak największą wartość nieujemną z N generując współczynniki w[i], i<k. Wtedy N miałoby przedstawienie:
N = 2*...*p[k]*w[k] + 2*...*p[k-1]*w[k-1] + ... + 2*3*w[2] + 2*w[1] ,
tylko - możemy zacząć od liczby pierwszej d=p[k+1], sprawdzamy, który z iloczynów 2*...*p[i] > d, oraz pozostałe łączymy w jedną wartość r zapisaną w systemie o podstawie p[k+1]. Podobnie iloczyn 2*...*p[i] zapisujemy w systemie o podstawie d=p[k+1], a pozostałe -- wystarczą tylko współczynniki w[i], ..., w[k] wyznaczane w tej kolejności. Wartości te będziemy wyliczali po drodze, każdy mnożony przez małą liczbę pierwszą p[k-j] dla j=k-i+1 malejące do 0. 

Szereg ma tę własność, że jeśli modyfikowane d będzie liczbą p[i]-gładką, to często począwszy od wyrazu i+1 uzyskamy 0 lub wielokrotność liczby pierwszej, które będzie sie propagowało w kierunku w[k]. Podzielność przez początkowe liczby pierwsze sprawdzimy np. bezpośrednio licząc największy wspólny dzielnik iloczynu 2*...*p[k] z N. 

Zatem pamiętamy: współczynniki w[k],...,w[i] rozwinięcia w szereg z przypisanymi im liczbami pierwszymi, liczby p[i] oraz r, zapisane w systemie o podstawie d, samo d.
Szukamy wartości d, dla której wartość szeregu przystaje do 0 modulo d. Wszystkie współczynniki są małymi liczbami, tylko p[i], d oraz r są większe. Aby znaleźć wartość reszty N modulo d, reszta z dzielenia przez 2*...*p[i] jest dana jawnie, zaś by uzyskać resztę z dzielenia przez 2*...*p[i]*p[i+1] mnożymy ją przez p[i+1] modulo d. Kontynujemy aż do p[k]. Wyznaczone reszty w kombinacji liniowej z w[j] liczone modulo d dają resztę N%d.
Jeśli iloczyn 2*...*p[i] = 1*d+1, na końcu iteracji przemieszczamy go do reszty r = r+w[i]*(d+1), oraz zaczynamy z wcześniejszym iloczynem
w[i+1] * (2*...*p[i] * p[i+1]) = w[i+1] * (p[i+1]*d + p[i+1])
pamiętając od teraz tę wartość. 

Przebieg iteracji: zwiększamy d+2, konwersje p[i] oraz r postaci
a*(d-2)+b = a*d+(b-2a)
z możliwym przeniesieniem a*d+b = (a-1)*d+(d+b), wyznaczamy współczynniki u[j]=2*...*p[j] modulo d dla i<j<=k, obliczamy kombinację liniową sum _{k>=j>=i} w[j]*u[j] i przyrównujemy ją do 0. wszystko na stosunkowo małych liczbach ograniczonych przez d, zatem złożoności maksymalnej obliczeń O(log n).
Wyrazów nie jest zbyt dużo. Dla RSA 200 okazało się, że wystarczy mniej niż 95 kolejnych liczb pierwszych, które powoli przemieszczają sie do r, czyli jakieś O(log log n) dla wyznaczenia wartości kombinacji liniowej. Tłumienie wartosci w r jest jak na modyfikacje wyrazów jest szybsze choć liniowe - staje się stałą wartością mniejszą od d.

Jeśli zauważymy, że musimy z d dla liczby N pierwszej dotrzeć do pierwiastka kwadratowego sqrt N, nasuwa się złozoność wykładnicza. Policzymy inaczej, kryterium stopu nie jest tutaj widoczne, powstaje, gdy suma szeregu rozbieżnego zaczynajacego się w p[i+1]^2 o wyrazie 4*d+4 staje sie większa niż N. Takich sumowań jest mniej niż N, zaś, jeśli N ma dzielniki (wszystkie większe niż p[i]), krotność iteracji jest zmniejszonym o p[i+1] rzędem najmniejszego z dzielników właściwych w grupie addytywnej, zatem O(n).
Podobnie jest w teście pierwszości AKS, gdzie w analogicznej grupie muliplikatywnej (małe twierdzenie Fermata) można czasem otrzymać rząd dzielnika równy N. Tu nawet nie docieramy to tak duzej wartości.

Łącznie szacuję złożoność na O(n * log n * log log n). Wielomianowe. 

Lecz algorytm jest jednym z wolniejszych, mimo że ma automatyczne ograniczenie stosowanych wielkości. Niewidoczne jest kryterium stopu. Propagowanie reszt powoduje, że niejako 'przy okazji' rozpoznajemy część wartości d jako złożonych, przy których kombinacja liniowa i tak nie będzie zerem choćby ze względu na dodatnie r.
Inną ciekawostką jest bardzo mało 'pamiętania', współczynniki w[] są wielkości liczb pierwszych, tylko dwie 'większe' wartości 'dwucyfrowe', z których używamy 'cyfr jedności' oraz podstawa systemu odpowiadająca dziesiątce d może osiągać wartość sqrt N -- wtedy lepiej jest stosować inny algorytm tej rodziny. Uzyska mniejszą złożoność dla tak dużych d, dokładniej już od d approx root[3] N inne algorytmy mają mniej prostszych obliczeń.  

Praktyczny zapis szeregu przy d=3:
7387= 3*[210*11] + 2*[30*7] + 1*[6*5] + 1*[2*3+0] + 1
w[] = [1, 1, 2, 3], r=1.