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.

Problem komiwojażera z obwiedni wypukłych, szczególny przypadek

 W październiku 2022 roku napisałem tu, że korzystając z obwiedni wypukłych można wyznaczyć trasę komiwojażera (ang. Salesman problem). 

Rzecz dzieje się na płaszczyźnie, mamy n miast, należy wyznaczyć trasę między wszystkimi miastami o najmniejszym koszcie, Tutaj potrzebuję współrzędnych miast, nie tylko symetrycznej macierzy kosztów podróży. Te dodatkowe informacje są potrzebne do wyznaczenia obwiedni wypukłych.

Heurystyka wyznacza najpierw obwiednię wypukłą A, która finalnie przekształci się w trasę komiwojażera K. Następnie wyznacza z pozostałych miast obwiednię wypukłą B. Miasta z B będą przemieszczane do A i równocześnie do K jako nadzbioru A.
Drobny kłopot związany z pierwszą krawędzią, oraz posiłkując się kolejnymi miastami obwiedni, dwa z A oraz najbliższe kolejne dwa z obwiedni B siłowo wyznaczam (dodatkowo ustalone już jest miasto startowe), kończę na obwiedni na której przesuwam się o miasto według orientacji.
Cały czas korzystam z pary miast, dwa z obwiedni A oraz miast doń przyległych wyznaczonej trasy; drugie z obwiedni B. Zatem mam do czynienia z trzech do czterech miast ułożonych liniowo dla wyznaczenia kosztu podróży, by przesunąć miasto z B do A. Permutacja jest tylko między sposobami zazębiania grafów. Np. a,c in A; b,d in B, trasy: acbd, abdc, abcd, także przed/za miasto zca, cde gdzie z, e in K są sąsiadami odpowiednio a oraz c trasy K.
Kiedy zamknę cykl, obwiednie przekształcą się w cykliczny graf skierowany przechodzące trasą przez wszystkie uwzględnione miasta, myślałem, że o najmniejszej długości.
W październiku brałem do nowo generowanego w kolejnych iteracjach A miasta obwiedni B, oraz wyznaczałem nową obwiednię B z miast jeszcze nie odwiedzonych. I tu pojawia się kłopot. Ta trasa niekoniecznie jest optymalną. 

Potrafi wystąpić przypadek niezbyt intuicyjny, kiedy trasa biegnie miastami obwiedni zewnętrznej, oraz suma długości trasy jest MNIEJSZA od odległości sąsiednich miast obwiedni wewnętrznej. Oszacowałem trasę, i wydaje mi się, że w tym przypadku miasta takie nie mogą być usuwane z A podczas przechodzenia do kolejnej iteracji, nowa obwiednia A powinna składać się z miast obwiedni B oraz dodatkowo tych miast. 

Przy tej modyfikacji naprawia się przypadek, kiedy jakieś miasto z B leży daleko od miast październikowej obwiedni A, ale wystarczajaco blisko trasy K, by je połączyć w jedną trasę od miast K. 

Oczywiście, podobnie jak w pażdzierniku, kiedy zbiór nieodwiedzonych jeszcze miast staje sie pusty, wyznaczony graf cykliczny skierowany staje się bardziej optymalną trasą komiwojażera.

Labirynty, rodzaje i ich generowanie

W tym poście dla odmiany zaprezentuję kilka pseudokodów na tworzenie labiryntow. W założeniu wszystkie pola w prostokącie mają być dostępne, co najwyżej przedzielone ścianami. Zatem rozmiar tablicy na labirynt jest iloczynem dwu liczb dodatnich nieparzystych, a do każdego z pól o obu współrzędnych parzystych można dotrzeć -- indeksowanie pól od zera. Nie ma większych sal, oraz labirynt gęsto pokrywa dostępny obszar.

Klasycznie podaje się labirynt generowany za pomocą stosu -- jako przykład zastosowania stosu. Taki labirynt jest drzewem. I to jest najczęściej spotykana przeze mnie implementacja w wielu źródłach.
W takim labiryncie trudno o sekretne przejścia czy cykle. Trzeba je specjalnie tworzyć. 

Labirynt oparty na liście dostępnych pól ze ścianami już podczas generowania może wygenerować ściany mogące być sekretnymi przejsciami, a nawet cykle - korytarze po których obchodzimy w kółko jakiś filar. Listę inicjujemy tylko częścią pól -- tylko tam gdzie są ściany do wyburzenia, a przy odwiedzaniu usuwamy pole z listy. Ten sposób wymaga resetów przy tworzeniu. 

A pseudokody? Okazuje się, że modyfikacje są tak nieznaczne, że można opisywać je równocześnie... 

 

Założenia: Wypełniamy obszar labiryntu ścianami, które będziemy burzyć podczas odwiedzania. Potrzebna jest dodatkowa funkcja zwracająca wektor lub pole bitowe pól odległych o dwa, czyli dwa pola obok w pionie, poziomie jako następne do zwiedzania. Wektor zwraca, czy ściany występują (pola nieodwiedzone), czy też już wyburzone (pole odwiedzone).  Dlaczego nie sąsiednie? Ponieważ potrzebujemy miejsce na ewentualną ścianę, która też jest polem prostokąta. Inaczej z labiryntu utworzy się plac wypełniajacy cały obszar. 

Obszar to tablica wypełniona liczbami o interpretacji np. 0 ściana widoczna, 1 ściana niewidoczna ukryta w mgle, 2 przejście widoczne, 3 przeście niewidoczne, 4 sekret nieodkryty, 5 sekret niewidoczny nieodkryty, 6 sekret odkryty widoczny...
Przy rysowaniu mapy wyświetlane tylko pola z wartościami parzystymi, a zmniejszaniem o 1 zajmie się osoba wędrująca po labiryncie odkrywajac go równocześnie.

1) Start. dla stosu oraz podczas pierwszego resetu dla list wchodzimy na pole początkowe, które odwiedzamy

2) Jesteśmy na odwiedzonym polu, odkładamy je na stos albo, w przypadku listy tuż po resecie, szukamy na liście pola, które przylega (krokiem 2) do już odwiedzonego pola. Łączymy je korytarzem, któremu można nadać status ukrytego przejścia

3) Wywołujemy funkcję zwracajacą zabudowę okolicznych pól dystansu 2

4) Losujemy kierunek, w którym pójdziemy. W przypadku stosu ograniczamy się tylko do nieodwiedzonych pól, o ile istnieją, wtedy czasem wyboru nie ma -- ślepy zaułek;
W przypadku listy z możliwymi cyklami przebijamy korytarz w wylosowanym kierunku zanim sprawdzimy wartość funkcji z 3). W końcu programistyczny kilof w garści jest labirynt[pozycja]=dostęp... 
 

5) Gdy wylosowany kierunek wskazuje nieodwiedzone pole, odwiedzamy je niszcząc także ścianę je przedzielajacą;
gdy jednak jesteśmy na już odwiedzonym polu, w przypadku stosu zabieramy element ze szczytu by wrócić pole wcześniej; zaś w przypadku list następuje reset.
Sprawdzamy także warunek końca algorytmu: stos ma być opróżniony albo lista pól do odwiedzenia pusta, czyli albo punkt 6), albo 2) gdy coś jeszcze zostało do wyburzenia / sprawdzenia

6) Stos pusty albo lista pusta -- cały obszar jest wypełniony labiryntem.
Można dołożyć drobiazgi czyli mgła na polach (np. wartość+1), skarby, NPC, oponenci. 

Kiedy implementowałem labirynty listą, okazało się, że są generowane nieco szybciej niż stosem, gdyż ten po wypełnieniu całego obszaru zwykle zawiera dużo pól. Nieodwiedzone pola wyznaczane z list były podłączane znacznie wcześniej. I nawet resetów nie było zbyt dużo.
Zaś w jednym wylosowanym przypadku napotkałem korytarz, whodzę -- sekret, mijam, kolejny sekret, i tak powstał długi korytarz wypełniony sekretnymi przejściami, widocznych już wcześniej na mapie labiryntu, ale dostępnych dopiero po znalezieniu przejścia daleko w przeciwnej części obszaru. W sam raz na znalezienie NPCa lub ukrywajacego się bandyty. 

Miłego zwiedzania...

To nie heurystyka, działa bardziej ogólnie korzystając z cech podzielności

Wspomniałem o szybkiej heurystyce rozkładu liczb. 

 

Liczbę szacuję za pomocą potęg 4, oraz zapisuję N jako parabolę o równaniu
a*p*p + b*p + c = N = q*r,

gdzie początkowo a jest jedną z trzech wartości power(4,k). 2*power(4,k), 3*power(4,k), a<p, Podstawa p też poczatkowo jest potęgą 4 bliską a, współczynniki b i c nieujemne są bliskie p, najlepiej z przedziału (-p, 3p).

 

Będziemy to modyfikować pilnując, by p było jak największe. Przy b i c stosunkowo niewielkich małe dzielniki mogą być niewykrywalne - iloraz dzielników r/q, r>q jest zbyt dużą wartośc ią, zaś algorytm czynimy zbieżnym do tego ilorazu. W tym przypadku należy sprawdzić, czy można doprowadzić do istnienia wspólnego dzielnika wszystkich współczynników.

 

Przekształcenie niezmiennicze dla wartości N przesuwa podstawę

[a, b, c, p] = [a, b-2k, c-2k*b+k*k*a, p+k]  (jak np. (1, 4, 4, 8] = [1, 0, 0, 10]

dla dowolnego k całkowitego. Dla tego przekształcenia c ze wzrostem p jest malejace, oraz szukamy jednej z dwu wartości, gdzie a+b+c  albo a-b+c sa całkowitymi wielokrotnościami p-1, p+1. Warto poszukać zmian znaku c dla jednego, lub niezależnie obu funkcji liniowych. Nie zwracają jednej wartości.

Kolejne przekształcenia niezmiennicze wartości N:

[a, b, c, p] = [a, b+1, c-p, p]  np. [5, 3, -1, 10] = [5, 2, 9, 10] = 529

[4a, 2b, c, p] = [a, b, c, 2p] np. [12, 2, 3, 10] = [3, 1, 3, 20]

płaszczę nieco parabolę przesuwając znalezioną wartość. Wymaga znów szukania miejsca, gdzie c ze zmianą p zmienia znak.

Kiedy spłaszczymy parabolę na tyle, że a+b+c=(p-1)*const, dzielnikiem liczby N jest p-1; Gdy a-b+c=0 + p*const, dzielnikiem N jest p+1. Przypomina to cechy podzielności przez 9 oraz 11 w dziesiątkowym - wszak to jest uogólnienie tych przypadków. 

Jęsli dopasujemy wartości a+b+c=p-1 oraz a-b+c=0 odległe dokładnie o 2, między nimi jest postać [a,b,0,p], skąd odczytujemy rozkład N = p * (a*p+b).

Teoretycznie jest to bardzo szybkie. Dla wyboru a oraz jednej funkcji mamy nie więcej niż 3*220 płaszczeń i wyszukiwań liniowych w przypadku RSA 200.

 

Przykładowo 259 zwraca swój dzielnik dla [4, 0, 3, 8] (suma 4+0+3=8-1=7, suma naprzemienna o wspólnym dzielniku [7, 35, 7, 4]).
Niektóre dzielniki można wychwycić dla dwu różnych wartości a - funkcje wykażą je od 'przeciwnych' stron, mniejszej albo większej.