lion137
2018-12-13 17:50

Sito Erastotenesa

Było coś obiecane z nowej książk i jest. Najpierw kod:

#include <iostream>
#include <vector>
#include <iterator>
#include <cctype>

#define RandomAccessIterator typename
#define Integer typename

template <RandomAccessIterator I, Integer N>
void mark_array(I first, I last, N factor) {
    first = false;
    while (last - first > factor) {
        first = first + factor;
        *first = false;
    }
}

template <RandomAccessIterator I, Integer N>
void erastosthenes_sieve(I first, N n) {
    I last = first + n;
    std::fill(first, last, true);
    N i(0);
    N index_square(3);
    N factor(3);
    while (index_square < n) {
        if (first[i]) {
            mark_sieve(first + index_square, last, factor);
        }
        ++i;
        index_square += factor;
        factor += N(2);
        index_square += factor;
    }
}

Szkoda, że nie można sobie już definiować konceptów / pojęć? Tylko takie haki preprocesorem, ale każdy się chyba domyśli, jaki typ dać w miejsce Integer:). Jak widać nie ma nigdzie tablicy z wartościami liczb pierwszych, nie potrzeba, można je sobie wyciągnąć z mutowanej(void) w sicie tablicy, generalnie czegokolwiek, na co będzie wskazywał Iterator(true jest na 2i + 3 miejscu).

std::vector<int> primes;
primes.push_back(2); // dodanie dwójki
    for (int i = 0; i < vec1.size(); i++){ // vec1 - wektor "wysłany" do sita 
        if (it1[i]){ // wskazuje na vec1
            primes.push_back(2 * i + 3);
        }
    }

W miarę szybkie (milion liczb pierwszych w 0.1 sek.), do matematyki rekreacyjnej i na uczelnię wystarczy:)
#theory #computation #primes

lion137

Myślałem o tym, ale chciałem zanęcić trochę concepts, które mają się pojawić niedługo:)

lion137
2018-02-04 01:24

Faktoryzacja 2

Wypadałoby, poza naiwnym dzieleniem, jak tutaj dodać, jakieś efektwniejsze (ale jeszce w miarę proste) metody faktoryzacji. Na początek Pollard's Rho - artykuł jest jasny, jak ktoś chce poeksperymentować, to kod w Pythonie:

def rho(n):
    if miller_rabin(n):
        return n
    def g(x, m):
        return (x*x - 1) % m
    x = 2
    y = 2
    d = 1
    while d == 1:
        x = g(x, n)
        y = g(g(y, n), n)
        d = gcd(abs(x - y), n)
    if d == n:
        return False
    else:
        return d

Funkcja, nie zawsze musi zwrócić dzielnik pierwszy, w praktyce trzeba by ograniczyć liczę prób, a jest skuteczny gdy czynniki nie są duże; złożoność: O(n^1/4).

Do większych czynników Fermat, prosta obserwacja, że każda liczba nieparzysta jest równa różnicy dwóch kwadratów, prowadzi bezpośrednio do algorytmu, Python:

def fermat_factorization(n):
    def is_square(a):
        if a < 0:
            a = -a
        return int(sqrt(a)) * int(sqrt(a)) == a
    a = int(sqrt(n))
    b = a * a - n
    while not is_square(b):
        a += 1
        b = a * a - n
    return [a - int(sqrt(b)), a + int(sqrt(b))]

Ta zaś funkcja, przeciwnie do rho, zwraca najszybciej czynniki bliskie pierwiatkowi z n. Algorytmy się ładnie uzupełniają, mogłyby współpracować w faktoryzacji (najlepiej równolegle,).

Na koniec jeden z algorytmów faktoryzujących duże liczby(znowu Pollard) Pollard p -1 algorytm - artykuł tłumaczy jasno, kod(Python):

def pollard_p_minus_1(n, max_iter, b):
    """Performing Pollard p- 1 factorization,
    with odd number input"""
    for _ in range(max_iter):
        primes = prime_sieve(b)
        for i in range(len(primes)):
            primes[i] = primes[i] ** (int(log(n, primes[i])))
        m = reduce(op.mul, primes)
        g = gcd(modularexponenation(2, m, n) - 1, n)
        if 1 < g < n:
            return g
        if g == 1:
            b += (b // 5)
        if g == n:
            b -= (b // 5)

Początkowe b, a także ilość iteracji, trzeba dobrać, rekordy algorytmu tutaj: https://members.loria.fr/PZimmermann/records/Pminus1.html.
Dodam tylko, że jest używany przez Great Internet Mersenne Prime Research, i w moich paru testach spisywał się bardzo dobrze, złożoność: O(b × log b × log2 n), jako podstwa znowu pojawia się Małe Twierdzenie Fermata. To tyle, pozostała tylko notka o generowaniu liczb pierwszych, ale to już natępnym razem, do przeczytania! #Primes #Numerical #Python #python

lion137
2018-01-31 18:45

GMP C++

Była ostatnio dyskusja na forum o nakładkach w C++ na GMP, postanowiłem sprawdzić interfejs od twórców. Prosta funkcja, Test Lucasa Lehmera, sprawdza pierwszość liczb Mersenne'a, używany do znajdowania wielkich liczb pierwszych: https://www.mersenne.org/. Kod:

#include <iostream>
#include <gmpxx.h>
#include <gmp.h>

int lucas_lehmer(unsigned long int p){
    mpz_class s = 4;
    mpz_class base = 2;
    mpz_class m = 0;
    mpz_pow_ui (m.get_mpz_t(), base.get_mpz_t(), p);
    m -= 1;
    for (unsigned long i = 0; i < p - 2; ++i){
        //s = s * s - 2;
        //mpz_mod (s.get_mpz_t(), s.get_mpz_t(), m.get_mpz_t());
        s = ((s * s) - 2) % m;
    }
    if (s == 0) return 1;
    else return 0; 
} 

int main(int argc, char **argv)
{   
    unsigned long int p = 44497;
    cout << lucas_lehmer(p) << endl; // -> True
    return 0;
}

Używanie i instalacja są proste, jak już mamy GMP, to trzeba tylko doinstalować libgmpxxv4-4, na linuxie: sudo apt-get install libgmpxxv4-4, do pliku linkujemy nagłówki, jak widać na listingu, a compilacja: g++ -Wall -std=c++11 main.cxx -lgmpxx -lgmp.
Wiekszośc operatorów arytmetycznych jest przeładowana (jak również IO streamy i konstruktory) z elementami typu mpz_class, jeśli funkcji nie ma w przestrzeni mpz_class, to mamy metodę get_mpz_t() "wyciagającą" zmienną mpz z klasy i taki element możemy wysłać do odpowiedniej funkcji w "czystym" GMP. Gdyby mod nie był przeładowany to w ten właśnie sposób (wersja ukomentowana) trzeba by napisać wewnętrzną pętlę.Jak działają funkcje mpz_pow_ui i mpz_mod można poczytać w dokumentacji To tyle, do następnego. #Numerical #Primes #C++

lion137
2018-01-30 13:34

Faktoryzacja Do Liczb Pierwszych

Jak już były testy(hashtag Primes), to teraz faktoryzacja, artykuł na Wikipedii szczegółowo opisuje dostępne metody. Wiele z nich (tych do dużych liczb) jest skomplikowana i raczej dla nas (jeśli nie pracujemy w kryptografii czy coś w tym stylu) mało użyteczna.

Ale w matematyce rekreacyjnej, jak adventofcode, hackerrank, projecteuler czy inne, wystarczy, w miarę szybko faktoryzować do, powiedzmy 10^16. A to można łatwo zrobić mając obliczoną tablicę liczb pierwszych do pierwiastka z n.

def factorize_trial(n, primes):
    """Factorizing n"""
    if miller_rabin(n):
        return n
    result = []
    for p in primes:
        while n % p == 0:
            result.append(p)
            n //= p
        if miller_rabin(n):
            result.append(n)
            break
    return result

primes można wygenerować używając, na przykład tego: simple_sieve.
Użyłem też, trochę zoptymalizowanego miller_rabin:

def miller_rabin(n):
    d = n - 1
    s = 0
    while d % 2 == 0:
        d >>= 1
        s += 1
    for repeat in range(20):
        a = 0
        while a == 0:
            a = random.randrange(n)
        if not miller_rabin_pass(a, s, d, n):
            return False
    return True

def miller_rabin_pass(a, s, d, n):
    a_to_power = pow(a, d, n)
    if a_to_power == 1:
        return True
    for i in range(s - 1):
        if a_to_power == n - 1:
            return True
        a_to_power = (a_to_power * a_to_power) % n
    return a_to_power == n - 1

Algorytm, jak widać iteruje po primes, próbując dzielić i dołączać do tablicy wynikowej, na końcu pętli jest sprawdzenie pierwszości, żeby ewentualnie zakończyć, gdy po podzieleniu n stanie się pierwsze.
Złożoność będzie primes(sqrt(n)) - liczba liczb pierwszych do pierwiastka z n, czyli będzie działał najwolniej, gdy liczba będzie mieć dwa duże czynniki pierwsze, ale dla wiekszości jest w miarę szybki. Kod. #Primes #Theory #Python #Numerical

lion137
2018-01-29 09:20

Lucas - Lehmer Test

Żeby zakończyć przegląd efektywnych Prime Tests, część druga, część pierwsza dodam test Lucasa - Lehmera, sprawdza on tylko liczby pierwsze Mersenne'a - Zastosowania w kryptografii - jest szybszy niż Miller - Rabin. Artykuł na Wikipedii wyjaśnia dokładnie ideę i implementację. Dodam tylko kod w Pythonie:

def lucas_lehmer(p):
    s = 4
    m = exp(2, p) - 1
    for i in range(p - 2):
        s = ((s * s) - 2) % m
    return True if s == 0 else False

W sumie jest w tej samej klasie złożoności co Miller Rabin ale wykonuje mniej operacji. Program sprawdził 27 liczbę Mersenne'a w 79 sekund, gdzie Miller Rabin wyłączyłem po kilkunastu minutach.
Złożoność algorytmu O(p^2 log p log log p) - zakładając liniowy algorytm mnożenia (jak np. w GMP)

Czyli generalne program do sprawdzania dowolnych liczb: najlepszy wybór Miller - Rabin.
Do liczb Mersenne'a: Lucas - Lehmer.

Z probablistycznych jeszcze warto wspomnieć o teście Frobeniusa, jego jedna runda jest dłuższa, ale osiąga szybsze ograniczenie prawdopodobieństwa.
Są jeszcze algorytmy deterministyczne, jak na przykład AKS, ale jest niewydajny O((logn)^6).
Poza tym, AKS jest deterministyczny, ale w praktyce, jeśli Miller Rabin I AKS zwrócą różne wyniki na tej samej liczbie, to jest większa szansa, że to promieniowanie kosmiczne(!) przypadkiem zmmieniło bit i wynik w komputerze, na którym był zapuszczony AKS (gdyż liczył dłużej!). #Primes #Python #Theory #Numerical

Aryman1983

Takie wpisy to ja rozumiem, a nie płacze nad swoją niedolą :-) Dzięki!

lion137
2018-01-26 23:38

Miller Rabin Test

Obiecany Miler Rabin, część pierwsza. Test Fermata, był prawie dobry, tylko dla pewnej klasy liczb (niewielkiej, ale jednak) zawsze się mylił (były złożone, ale przechodziły). Okazuje się, że można to usunąć, (Miller Rabin).

Najpierw Twierdzenie Fermata trochę inaczej zapisane:
Jeśli n jest liczbą pierwszą, a a dowolną liczbą dodatnią mniejszą od n, to a podniesione do potęgi n - 1przystaje do 1 modulo n.

Czyli(opis algorytmu):
Bierzemy losową liczbę 0 < a < n i podnosimy ją do potęgi n - 1 modulo n;
Podnosząc do potęgi, w każdym momencie podnoszenia do kwadratu szukamy "nietrywialnego pierwiastka kwadratowego jedynki modulo n", czyli sprawdzamy czy kolejna liczba, nie równa jeden lub n - 1 i podniesiona do kwadratu nie jest równa jeden modulo n (jeśli taką znajdziemy to liczba jest złożona). Dla liczby nieparzystej i złożonej, dla co najmniej połowy liczb a < n -1, sprawdzanie tym sposobem odkryje nietrywialny pierwiastek. Dlatego procedura Miller Rabin nie może być oszukana i można w granicy zredukowac błąd do zera gdy liczba prób dąży do nieskończoności. Pseudokod:

test_nontrivial(a, n):
    # t and u, t >=1, u odd and 2^tu = n - 1
    x[0] = modular-exponenation(a, u, n);
    for i -> 1 to t:
        x[i] = x[i - 1]^2 mod n;
        if x[i] == 1 and x[i - 1] != 1 x[i - 1] != n - 1:
            return True (composite)
        if x[t] != 1:
            return True
    return False

Dzięki sztuczce z dobraniem t i u, na początku i ustaleniu x[0], potem, aby podnieść a do n -1 modulo n, wystarczy podnieść x do kwadratu t razy, co umożliwia sprawdzanie warunku. Po iteracji, gdy nie mamy jedynki if x[t] != 1 to zwracamy True i finalnie False (nie znaleźliśmy nietrywialnych pierwiastków).
Teraz wystarczy tylko wykonać odpowiednia ilość testów w pętli, aby uzyskać zadaną dokładność, pseudokod:

miller-rabin(n, s):
    for j <- 1 to s:
        a = random(1, n - 1);
        if nontrivial_root(a, n):
            return False (composite)
        else:
            return True

I kod w Pythonie (nie wiele się różniący:)):

# 2. Miler - Rabin

# procedure nontrivial_root, which looking for a non trivial square root of one mod n

def nontrivial_root(a, n):
    """checking Fermat Theorem, and non trivial square root"""

    # find t and u, such that u odd, t >= 1 and n - 1 = 2^tu:
    t, u = 0, n - 1
    while u % 2 == 0:
        t += 1
        u //= 2

    x0 = modularexponenation(a, u, n)
    for i in range(1, t + 1):
        x1 = x0 ** 2 % n
        if x1 == 1 and x0 != 1 and x0 != n - 1:
            return True
    if x1 != 1:
        return True
    return False

def miller_rabin(n, s):
    """Miler Rabin Test"""

    # Return True if n = 2
    if n == 2:
        return True
    # Return False if even
    if n % 2 == 0:
        return False

    for i in range(s):
        a = randint(1, n - 1)
        if nontrivial_root(a, n):
            return False
    return True

Dokładność, można udowodnić, że błąd (gdy miller_rabin zwróci True) wynosi najwyżej 1 / 2^s, czyli, zgodnie z tą odpowiedzią na stacoverflow: https://stackoverflow.com/a/6330138, s = 40 wystarcza w zupełności. I na koniec złożoność obliczeniowa: O(k * (logn)^3) - jak daleko odeszliśmy od wykładniczej (co prawda w bitach:)). Kod na githubie. Mam nadzieję, że wywód był jasny, to tyle, dziękuję za cierpliwość. #Theory #Primes #Python #python #Numerical

enedil

Tbh, pominąłeś wyjaśnienie najbardziej kluczowej kwestii, dlaczego właściwie można się ograniczyć do sprawdzania tych właśnie liczb, jako kandydatów na nietrywialne pierwiastki.

lion137

@enedil: Zobacz też w części pierwszej, nic nie pominałęm, tak jest definicja liczb Carmichael'a, i procedura nontrivial_root te liczby wykrywa, co powoduje, że miller-rabin jest ulepszonym testem Fermata - Przeprowadza go dla s liczb, przy okazji sprawdzając czy poza tym, że spełniają one warunki twierdzenia Fermata, nie są relatywnie pierwsze do n (czyli Carmichael numbers). Finalnie wracamy do warunków z pierwszego testu - mamy liczbę pierwszą z prawdobodobieństwem.

lion137
2018-01-26 01:49

Primality Tests, Test Fermata

Jak już powstał wpis o algorytmach numerycznych, to w temacie nie można nie wspomnieć o liczbach pierwszych, a zacząć wypada od testów pierwszości. Wypada również pominąć naiwne testy, oparte na szukaniu dzielników, bezużyteczne dla dużych liczb (w teorii liczb bilion to mała liczba) - mają złożoność wykładniczą w bitach.

Jako pierwszy, chyba najprostszy, opierający się na Małym Twierdzeniu Fermata:
Jeśli n jest liczbą pierwszą, a jest jakąkolwiek dodatnią liczbą mniejszą od n, to a podniesione do n - tej potęgi przystaje do a modulo n
Dwie liczby przystają do siebie modulo n, jeśli mają tą samą reszte przy dzieleniu przez n.

Idea, Jeśli n jest złożone, to większość liczb nie będzie spełniać relacji z twierdzenia, co prowadzi do algorytmu:
Mając dane n, weź liczbę 0 < a < n i oblicz a do potęgi n modulo n, jeśli wynik nie równa się a, to liczba jest (na 100%) złożona, jeśli się równa, to są spore szanse na to, że n jest pierwsza (skoro w przypadku złożonej większość liczb a nie spełnia zależności); teraz weź kolejne losowe 0 < a < n i jeśli spełnia, to szanse na pierwszość n rosną jeszcze bardziej, i każde kolejne sprawdzenie zwieksza naszą pewność. Taki algorytm, to Test Fermata.

Potrzebne do algorytmu potęgowanie modulo, znajduje się, na przykład tutaj. Kod (Python):

# Primality Tests
# Copyleft 2018 lion137
# 1. Fermat Test

from Euler import modularexponenation
import random

def fermat_test_single(n):
    """performs single Fermat test"""
    def test(a):
        return modularexponenation(a, n, n) == a
    return test(random.randint(1, n - 1))

Dość oczywista zgodność kodu ze specyfikacją. Następna funkcja wykonuje test zadaną ilośc razy, dając zaskakująco skuteczny test:

def fermat_extended(n, cnt):
    """performing fermat tests in a loop,
    for #cnt of numbers (0, n) """
    for i in range(cnt):
        if not fermat_test_single(n):
            return False
    return True

Można ją potestować, ale najpierw o ograniczeniach. Po pierwsze, jako algorytm zrandomizowany nie zawsze jest poprawny, tzn., istnieje bardzo duża szansa, że daje dobry wynik, ale tylko szansa, precyzyjniej, jeśli test nie przejdzie, to liczba, na pewno jest złożona, ale jeśli przejdzie, to jest pierwsza z prawdopodobieństwem. Ale i to nie jest prawdą, są liczby, które są złożone, ale przechodzą test (mają własność, że dla każdego a, a do n przystaje do a modulo n), są to Carmichael Numbers. Liczby te są niezwykle rzadkie - poniżej stu milionów jest ich tylko 255, ale jednak. Pora na testy, tutaj liczby:

    print(fermat_extended(2 ** 2203 - 1, 5)) # -> True (a ** b w Pythonie to a do potęgi b)
    print(fermat_extended(2 ** 4423 - 1, 5)) # -> True
    print(fermat_extended(137, 6))           # -> True
    print(fermat_extended(6601, 6))          # -> True, złożona (Carmichael)
    print(fermat_extended(2821, 6))          # -> True, złożona (Carmichael)
    print(fermat_extended(999999999999999999, 6)) # -> False (oczywiście złożona)

Jesli chcemy przetestować losową liczbę pierwszą, to ta procedura, z praktycznych względów będzie najlepsza, prawdopodobieństwo błędu na losowym m bitowym numerze dąży do zera, gdy m dąży do nieskończoności (szansa błędu na losowej 1024 bitowej liczbie jest jeden do 10^41 całkiem nieźle:)). Ale gdy liczby nie są losowe to już tak nie jest - na liczbach Carmichael'a, błędu wyeliminować się nie da. To zrobi dopiero (ale później:) ) algorytm, bazujący na teście Fermata, Miller - Rabin. Jeszcze tylko złożoność :O(k (logn)^2 log log n * log log log n),całkiem przyzwoita - dowód raczej do pominięcia:-D. To tyle.
Żródła: Wikipedia, CLRS. #Theory #Primes #python #Python #Numerical