Na początku grudnia opublikowałem zarys algorytmu, w którym podstawową ideą jest utworzenie płaszczyzny z liczby, oraz szukanie śladu dzielników dla dużych wartości.
Na razie coś mnie blokuje przed napisaniem tego w C++, ale jest już wersja w Pythonie.
Miałem kłopoty, jak szybko ominąć część kodu po znalezieniu małego dzielnika. Zdecydowałem się na pętlę.
[code python]
def nwd(a,b,s):
if( 0==b ):
return a
if( a<b ): #swap(a,b)
c=a
a=b
b=c
while( 1 ): # zmniejszamy b, bo wiemy, ze sa one wzglednie pierwsze
if( 0==b%2 ):
b/=2
continue
if( 0==b%3 ):
b/=3
continue
if( 0==b%5 ):
b/=5
continue
if( 0==b%7 ):
b/=7
continue
break
if( 0==a%b ):
return b
a = a%b
if( a+a>b ) :
a = b-a
if( s>a ): #zatrzyma: wiemy ze nie ma tak malych dzielnikow
return 1
return nwd(b,a,s)
def modstop(s):
# zwieksza s do najblizszej liczby nieparzystej nie podzielnej przez 3, 5, 7
while( 1 ):
s+=2
if( 0==s%3 ):
continue
if( 0==s%5 ):
continue
if( 0==s%7 ):
continue
#break
return s
s=7 # lub wiecej, jesli jestesmy pewni
n= duza liczba
a=int (sqrt n) # jak to policzyc na liczbach calkowitych?
#a=int (n**(0.5)) # blad w liczeniu pozniejszych operacji
b= 1+(n-a*a)//a
c= n-(a+b)*a # w algorytmie c<0, trzymamy je jako ujemna wartosc
print a,b,c, (a+b)*a+c # sprawdzenie niezmiennika (a+b)*a+c = n
while( 1 ):
if( 0==n%2 ): # male dzielniki
n /=2
print "Dzielnik 2"
break
if( 0==n%3 ):
n/=3
print "Dzielnik 3"
break
if( 0==n%5 ):
n/=5
print "Dzielnik 5"
break
if( 0==n%7 ):
n/=7
print "Dzielnik 7"
break
while ( 1 ):
f = a%s # obliczenie (a+b)*a+c modulo s czyli spr. male dzielniki
e = ((f+b)*f+c)%s
if( 0==e ):
print "Dzielnik z modulo ",s
break
s = modstop(s) # zwiekszenie s
e = nwd(a,-c,s) # nwd wyrazu wolnego i podstawy systemu
if( 1<e ):
print "Dzielnik ze sladu ",e
break;
a-=1 # pelni role podstawy systemu, zmniejszamy o 1
b+=2 # wyrugowanie wplywu podstawy na wyraz wolny
c+=b-1 # modyfikacja wyrazu wolnego
if( 0==c ) :
print "Dzielnik z wyrazu wolnego ",a
break
if( 0<c ) : // wyraz wolny wiekszy niz podstawa a
b+=1
c-=a
#mozna wypisac - dla testow
#print "liczba (1*",a,"+",b,")*",a,c," stop = ",s
if( s>a ) : # duze lub male dzielniki nie znalezione, brak mozliwosci innych
print "Liczba pierwsza ",n
break
break
[\code python]
Dla malych wartosci dziala. Dla liczby 90-bitowej po dobie liczenia iloraz wartosci do sprawdzenia do oszacowania dzielnika z dołu jest mniejszy niż 2000, następnego dnia poniżej 1000. Sprawdzane są równocześnie małe i duże kandydaci, dlatego zakres możliwości się zmniejsza z obu stron.
Drugiego dnia zostało sprawdzone około 0,1% wszystkich możliwości, zastosowany wzór: iloraz różnicy podstawy a i minimalnej możliwej liczby pierwszej s przez liczność początkową (a-s)/(sqrt n). Można zatem spróbować oszacować czas przebiegu na tej maszynie na kilkaset dni. Choć podejrzewam istnienie przedziałów, w których po wejściu do nwd nastąpi od razu wyjście ( bo |c|<s ).
Mam jeszcze jeden pomysł związany z tą płaszczyzną, choć będzie wymagał sporo pamieci i moze być nieefektywny. Przechowywanie podstawy w postaci docelowo kanonicznej... Jako skutek uboczny pojawi się tablica liczb pierwszych.
18 grudnia 2019
09 grudnia 2019
Obliczanie silni większych liczb
Spróbowałem kilkoma sposobami obliczyć silnię ze 100.
W tym celu należy dobrze przygotować strukturę. Liczbę trzymam na liście dwukierunkowej o węzłach link postaci:
[int cyfra, link prev, link next, int flaga]
Okazało się, że same wskaźniki na sąsiednie cyfry prev oraz next nie wystarczą. Aby dobrze zaalokować liczbę na liście potrzebna była flaga, gdyż w przeciwnym przypadku były kłopoty z jednocyfrową liczbą 0. Nie była ona odróżniana od węzła pustego. Ponieważ węzły należały do dużego kontenera, szukanie w tym kontenerze wolnego miejsca na kolejną cyfrę wskazywało taką liczbę jako puste miejsce i nadpisywało...
Sama liczba ma dwa pola: link value oraz link baze. Składa się z cyfr dziesiętnych (docelowo), chwyconych za cyfrę najmniej znaczącą oraz podstawę systemu liczbowego.
Oprócz tego występują funkcje mnożenia przez 10 "mul0()" oraz przez małą wartość "mula( int a)". Dołączyłem także konwersje między systemami niedziesiątkowymi "convert (int new_baze)".
Wykorzystałem najpierw standardowe mnożenie dla liczenia silni z liczby n dla sprawdzenia wyników:
Liczba a;
for( i=n; 0<i; i-- )
a.mula(i);
a.show();
Następnie konwersje, by uniknąć kłopotów rozbiłem je na trzy fazy:
- liczenie 9! małymi iloczynami - by uzyskać odpowiednią wartość początkową;
- przekształcenie wartości do liczby w systemie o podstawie n;
- reszta obliczeń, by wynik wyszedł w dziesiątkowym:
Liczba a = 9!;
a.baze( n );
a.mul0();
for( i=n-1; 9<i; i-- ) {
a.convert(i);
a.mul0();
};
a.show();
Wreszcie rozklad kanoniczny czyli mieszanka konwersji oraz mul0() w odpowiednich systemach. I tu się sypało, uzyskana wartość była znacznie za mała.
W pozostałych też były kłopoty. Okazało się, że zarówno przy mnożeniu mula(int), jak i przy konwersjach nie mogę sobie pozwolić na przechowywanie reszt w dodatkowej zmiennej przechodząc między węzłami cyfr. Należy je przekazać sąsiedniej cyfrze zaraz po ich zaobserwowaniu. Inaczej czasem miewałem błędy oraz uzyskana wartość była zawyżana.
A oto szczególy funkcji:
mul0() { // mnożenie przez 10 (dokładniej: przez baze)
link cyfr = t.add(); // nowy węzeł z kontenera;
cyfr.cyfra = 0;
value->next = cyfr;
cyfr->prev = value;
value = cyfr;
};
czyli zwykłe doczepienie węzła z zerem na końcu liczby i uchwycenie go jako liczby.
mula( int a ) {}
// pierwszy przebieg: mnożenie każdej cyfry liczby przez a
link b = value;
while( b ) { b->item *= a; b=b->prev; }
// przeniesienia wartości
b = value;
int s=0;
while( 1 ) {
s = b->cyfra;
b->cyfra = s% 10; // system dziesiątkowy
s /= 10;
if( s ) {
if( b-> prev ) b->prev->cyfra += s;
else {
link c=t.add();
c->cyfra = s;
c->next = b; b->prev = c; // przypięcie węzła cyfry najbardziej znaczącej
}
};
s=0;
if( b->prev ) b = b->prev; else break;
};
Konwersja też jest podzielona na dwa moduły, w pierwszym wymagana jest liczba dwucyfrowa oraz ustawienie węzła blokującego b na cyfrę poprzedzajacą najbardziej znaczącą. Drugi węzeł link c przebiega cyfry bardziej znaczące od b ze wzorem s = cyfra*r, której część całkowita zasila bieżącą cyfrę wskazywaną przez c, zaś reszta jest odkładana do cyfry mniej znaczącej c->next->cyfra. Podstawa baze maleje, zatem wartość liczby rośnie.
int r = baze->cyfra - new-baze;
int s=0;
// modyfikacja pól baze
link b = value;
while( b->prev ) b = b->prev;
while( 1 ) {
if( b->next ) b = b->next; else break;
c = b->prev;
while( 1 ) {
s = c->cyfra * r;
c->cyfra += s/new_baze;
c->next->cyfra += s % new_baze;
if( c->prev ) c = c->prev; else break;
};
c=b; ... // druga część
};
W drugiej części zmieniamy kierunek c=b oraz przenosimy nadmiary podobnie jak przy mnozeniu mula(int).
Potem przesuwamy b = b->next; aż przekształcane będą wszystkie cyfry liczby.
Obydwoma sposobami w przeciągu sekund wypisało mi wartość wszystkich 158 cyfr dziesiętnych silni 100!
W tym celu należy dobrze przygotować strukturę. Liczbę trzymam na liście dwukierunkowej o węzłach link postaci:
[int cyfra, link prev, link next, int flaga]
Okazało się, że same wskaźniki na sąsiednie cyfry prev oraz next nie wystarczą. Aby dobrze zaalokować liczbę na liście potrzebna była flaga, gdyż w przeciwnym przypadku były kłopoty z jednocyfrową liczbą 0. Nie była ona odróżniana od węzła pustego. Ponieważ węzły należały do dużego kontenera, szukanie w tym kontenerze wolnego miejsca na kolejną cyfrę wskazywało taką liczbę jako puste miejsce i nadpisywało...
Sama liczba ma dwa pola: link value oraz link baze. Składa się z cyfr dziesiętnych (docelowo), chwyconych za cyfrę najmniej znaczącą oraz podstawę systemu liczbowego.
Oprócz tego występują funkcje mnożenia przez 10 "mul0()" oraz przez małą wartość "mula( int a)". Dołączyłem także konwersje między systemami niedziesiątkowymi "convert (int new_baze)".
Wykorzystałem najpierw standardowe mnożenie dla liczenia silni z liczby n dla sprawdzenia wyników:
Liczba a;
for( i=n; 0<i; i-- )
a.mula(i);
a.show();
Następnie konwersje, by uniknąć kłopotów rozbiłem je na trzy fazy:
- liczenie 9! małymi iloczynami - by uzyskać odpowiednią wartość początkową;
- przekształcenie wartości do liczby w systemie o podstawie n;
- reszta obliczeń, by wynik wyszedł w dziesiątkowym:
Liczba a = 9!;
a.baze( n );
a.mul0();
for( i=n-1; 9<i; i-- ) {
a.convert(i);
a.mul0();
};
a.show();
Wreszcie rozklad kanoniczny czyli mieszanka konwersji oraz mul0() w odpowiednich systemach. I tu się sypało, uzyskana wartość była znacznie za mała.
W pozostałych też były kłopoty. Okazało się, że zarówno przy mnożeniu mula(int), jak i przy konwersjach nie mogę sobie pozwolić na przechowywanie reszt w dodatkowej zmiennej przechodząc między węzłami cyfr. Należy je przekazać sąsiedniej cyfrze zaraz po ich zaobserwowaniu. Inaczej czasem miewałem błędy oraz uzyskana wartość była zawyżana.
A oto szczególy funkcji:
mul0() { // mnożenie przez 10 (dokładniej: przez baze)
link cyfr = t.add(); // nowy węzeł z kontenera;
cyfr.cyfra = 0;
value->next = cyfr;
cyfr->prev = value;
value = cyfr;
};
czyli zwykłe doczepienie węzła z zerem na końcu liczby i uchwycenie go jako liczby.
mula( int a ) {}
// pierwszy przebieg: mnożenie każdej cyfry liczby przez a
link b = value;
while( b ) { b->item *= a; b=b->prev; }
// przeniesienia wartości
b = value;
int s=0;
while( 1 ) {
s = b->cyfra;
b->cyfra = s% 10; // system dziesiątkowy
s /= 10;
if( s ) {
if( b-> prev ) b->prev->cyfra += s;
else {
link c=t.add();
c->cyfra = s;
c->next = b; b->prev = c; // przypięcie węzła cyfry najbardziej znaczącej
}
};
s=0;
if( b->prev ) b = b->prev; else break;
};
Konwersja też jest podzielona na dwa moduły, w pierwszym wymagana jest liczba dwucyfrowa oraz ustawienie węzła blokującego b na cyfrę poprzedzajacą najbardziej znaczącą. Drugi węzeł link c przebiega cyfry bardziej znaczące od b ze wzorem s = cyfra*r, której część całkowita zasila bieżącą cyfrę wskazywaną przez c, zaś reszta jest odkładana do cyfry mniej znaczącej c->next->cyfra. Podstawa baze maleje, zatem wartość liczby rośnie.
int r = baze->cyfra - new-baze;
int s=0;
// modyfikacja pól baze
link b = value;
while( b->prev ) b = b->prev;
while( 1 ) {
if( b->next ) b = b->next; else break;
c = b->prev;
while( 1 ) {
s = c->cyfra * r;
c->cyfra += s/new_baze;
c->next->cyfra += s % new_baze;
if( c->prev ) c = c->prev; else break;
};
c=b; ... // druga część
};
W drugiej części zmieniamy kierunek c=b oraz przenosimy nadmiary podobnie jak przy mnozeniu mula(int).
Potem przesuwamy b = b->next; aż przekształcane będą wszystkie cyfry liczby.
Obydwoma sposobami w przeciągu sekund wypisało mi wartość wszystkich 158 cyfr dziesiętnych silni 100!
03 grudnia 2019
Algorytm największego wspólnego dzielnika
Klasa największego wspólnego dzielnika.
Wiemy skąd inąd, że jeśli wartości zmniejszą się poniżej wartości zewnętrznej stop, liczby są względnie pierwsze.
Jak modna hermetyzacja, to proszę :) Bardziej zahermetyzować się już chyba nie da.
class Gcd {
private:
Liczba a, b;
Liczba( Liczba an, Liczba bn ) : a(an) b(bn) {};
Liczba ngcd();
friend Liczba nwd(Liczba a, Liczba b);
};
Liczba Gcd :: ngcd() {
if( a<b ) swap(a,b);
a = a%b;
if( !a ) return b;
if( a+a>b ) a = b-a;
if( a<stop ) return Liczba(1);
return ngcd();
};
Funkcja wywoławcza:
Liczba gcd( Liczba an, Liczba bn ) {
Nwd c(an,bn);
return c.ngcd();
};
Algorytm jest nieco szybszy niż klasyczny. Warunek a+a>b pozwala jeszcze raz zabrać b uzyskując wartość ujemną w a, do której w następnym kroku byśmy dodawaii. Mnożenie przez (-1) pozwala przywrócić wartość dodatnią, co odbywa się niejawnie.
Wartości funkcji są zewnętrzne do wywołania funkcji rekursywnej, dzięki czemu jest mniejsze zużycie stosu pamięci. Można pójść dalej, w funkcji gcd utworzyć petlę nieskończoną wywołującą c.ngcd(), która zwraca wartość np. -1 jak ma liczyć dalej, 0 jak względnie piersze lub 1 gdy zwracana Liczba jest już w b. Czy wtedy będzie jeszcze widać, że ta funkcja jest rekursywna?
Wiemy skąd inąd, że jeśli wartości zmniejszą się poniżej wartości zewnętrznej stop, liczby są względnie pierwsze.
Jak modna hermetyzacja, to proszę :) Bardziej zahermetyzować się już chyba nie da.
class Gcd {
private:
Liczba a, b;
Liczba( Liczba an, Liczba bn ) : a(an) b(bn) {};
Liczba ngcd();
friend Liczba nwd(Liczba a, Liczba b);
};
Liczba Gcd :: ngcd() {
if( a<b ) swap(a,b);
a = a%b;
if( !a ) return b;
if( a+a>b ) a = b-a;
if( a<stop ) return Liczba(1);
return ngcd();
};
Funkcja wywoławcza:
Liczba gcd( Liczba an, Liczba bn ) {
Nwd c(an,bn);
return c.ngcd();
};
Algorytm jest nieco szybszy niż klasyczny. Warunek a+a>b pozwala jeszcze raz zabrać b uzyskując wartość ujemną w a, do której w następnym kroku byśmy dodawaii. Mnożenie przez (-1) pozwala przywrócić wartość dodatnią, co odbywa się niejawnie.
Wartości funkcji są zewnętrzne do wywołania funkcji rekursywnej, dzięki czemu jest mniejsze zużycie stosu pamięci. Można pójść dalej, w funkcji gcd utworzyć petlę nieskończoną wywołującą c.ngcd(), która zwraca wartość np. -1 jak ma liczyć dalej, 0 jak względnie piersze lub 1 gdy zwracana Liczba jest już w b. Czy wtedy będzie jeszcze widać, że ta funkcja jest rekursywna?
Algorytm rozkładu z płaszczyzny uzyskanej z liczby
W połowie listopada podałem, jak z liczby n zrobić płat F(b,d)=c o niezmienniku
n = d^2+bd+c = d(b+d)+c
Okazało się, że zarówno dzielniki, jak i ich wielokrotności tworzą łatwo wykrywalne cięcia przy stałych d, z których łatwo znajdziemy dzielniki. Im mniejszy dzielnik, tym tych specyficznych cięć jest więcej, regularnie rozłożonych. Oznacza to, że możemy ograniczyć się tylko do badania liczb bliskich pierwiastowi kwadratowego z liczby rozkładanej.
Algorytm wygląda zatem następująco:
szukamy pierwiastka liczby sqrt n; jego wyraz wolny równy 0 wyznacza rozkład z pierwiastka;
obliczamy wartość modulo mała liczba stop z wyrażenia d(b+d)+c; możemy też teraz zwiększyć stop;
obliczamy największy wspólny dzielnik wyrazu wolnego c oraz podstawy d, co sprawdza przynależność do tego specyficznego cięcia;
przejście na sąsiednie cięcie płata F:
d--;
b+=2;
c+=b-1;
czasem tak uzyskana wartość c jest większa niż podstawa d, zmniejszamy ją w pętli o operacjach przypominających przeniesienie między cyframi liczby w systemie pozycyjnym
{ b++; c-=d; }
jeśli uzyskamy 0=c, poznamy dzielnik.
Kryterium stopu jest
stop > d
po spełnieniu którego wiemy, że liczba jest pierwsza.
Dochodzi kilka drobnych, lecz istotnych szczegółów. Możemy przyspieszyć wzrost wartości stop do kolejnych liczb niepodzielnych przez 3, 5, 7 (oczywiście większych niż 7, by ewentualny tak mały dzielnik nam nie umknął). Oznacza to, że stop staje się minimum z badanych wartości, co możemy uwzględnić w liczeniu największej wspólnej wielokrotności. Nachylenie tej krzywej dyskretnej stop(d) jest między 2d a 10d. Od dużych wartości jest ograniczone przez d, Zmniejsza to obszar poszukiwań pierwiastka.
Jeżeli iloraz dzielników liczby n=p*q, czyli p/q>1, napotkamy po drodze więcej niż jedno cięcie wskazujące dzielnik. Zatrzymujemy się na pierwszym. Na ogół wtedy znaleziony wspólny dzielnik jest iloczynem małej wartosci i dzielnika n. Może być parzysty, dlatego nie możemy wykluczyć parzystych wartości d. Potrzebne będzie wtedy jeszcze jedno liczenie znalezionego wspólnego dzielnika z liczbą rozkładaną n.
Kod zwiększania wartości stop jest następujący:
if( nie_podzielne(...) ) // test: stop | d(b+d)+c
else while(1) {
stop+=2;
if( 7<stop && (!(stop%3) || !(stop%5) || !(stop%7)) ) continue;
else break;
};
Kod obliczania dzielnika podam w kolejnym poście. Jest zbyt specyficzny i pozwoli zerknąć na nowe możliwości obliczania funkcji rekursywnych.
n = d^2+bd+c = d(b+d)+c
Okazało się, że zarówno dzielniki, jak i ich wielokrotności tworzą łatwo wykrywalne cięcia przy stałych d, z których łatwo znajdziemy dzielniki. Im mniejszy dzielnik, tym tych specyficznych cięć jest więcej, regularnie rozłożonych. Oznacza to, że możemy ograniczyć się tylko do badania liczb bliskich pierwiastowi kwadratowego z liczby rozkładanej.
Algorytm wygląda zatem następująco:
szukamy pierwiastka liczby sqrt n; jego wyraz wolny równy 0 wyznacza rozkład z pierwiastka;
obliczamy wartość modulo mała liczba stop z wyrażenia d(b+d)+c; możemy też teraz zwiększyć stop;
obliczamy największy wspólny dzielnik wyrazu wolnego c oraz podstawy d, co sprawdza przynależność do tego specyficznego cięcia;
przejście na sąsiednie cięcie płata F:
d--;
b+=2;
c+=b-1;
czasem tak uzyskana wartość c jest większa niż podstawa d, zmniejszamy ją w pętli o operacjach przypominających przeniesienie między cyframi liczby w systemie pozycyjnym
{ b++; c-=d; }
jeśli uzyskamy 0=c, poznamy dzielnik.
Kryterium stopu jest
stop > d
po spełnieniu którego wiemy, że liczba jest pierwsza.
Dochodzi kilka drobnych, lecz istotnych szczegółów. Możemy przyspieszyć wzrost wartości stop do kolejnych liczb niepodzielnych przez 3, 5, 7 (oczywiście większych niż 7, by ewentualny tak mały dzielnik nam nie umknął). Oznacza to, że stop staje się minimum z badanych wartości, co możemy uwzględnić w liczeniu największej wspólnej wielokrotności. Nachylenie tej krzywej dyskretnej stop(d) jest między 2d a 10d. Od dużych wartości jest ograniczone przez d, Zmniejsza to obszar poszukiwań pierwiastka.
Jeżeli iloraz dzielników liczby n=p*q, czyli p/q>1, napotkamy po drodze więcej niż jedno cięcie wskazujące dzielnik. Zatrzymujemy się na pierwszym. Na ogół wtedy znaleziony wspólny dzielnik jest iloczynem małej wartosci i dzielnika n. Może być parzysty, dlatego nie możemy wykluczyć parzystych wartości d. Potrzebne będzie wtedy jeszcze jedno liczenie znalezionego wspólnego dzielnika z liczbą rozkładaną n.
Kod zwiększania wartości stop jest następujący:
if( nie_podzielne(...) ) // test: stop | d(b+d)+c
else while(1) {
stop+=2;
if( 7<stop && (!(stop%3) || !(stop%5) || !(stop%7)) ) continue;
else break;
};
Kod obliczania dzielnika podam w kolejnym poście. Jest zbyt specyficzny i pozwoli zerknąć na nowe możliwości obliczania funkcji rekursywnych.
Subskrybuj:
Posty (Atom)