lion137
2018-11-28 00:08

TDD in C++

Ciąg dalszy, po tam jakiś teoretycznych próbkach, napisałem coś bardziej praktycznego:): ADT geneyczny graf w C++. Niestety założenia uczyniły interfejs tłustym, jak tłusty bankowy kocur:); przytaczam klasę w całości:

/*
 * template class Graph ADT (undirected)
 */

template <class T>
class Graph {

    int V;
    int E;  
    std::unordered_map<T, std::forward_list<T>> adj;
    std::vector<T> vert;

    void construct(int _V){
        V = _V;
        E = 0;
        count = 0;
        std::forward_list<T> a;
        T el;
        for (int i = 0; i < V; i++){
            vert.push_back(el);
        }
    }

    public:

    std::unordered_map<T, bool> marked;
    int count;  // #connected components

    Graph(int _V) {
        construct(_V);
    }

    Graph(std::string fname){
        std::ifstream in(fname);
        int a, e;
        in >> a;
        in >> e;
        construct(a);
        T ver;
        std::forward_list<T> a_list;
        for (int i = 0; i < V; i++){
            in >> ver;
            vert[i] = ver;
            marked[ver] = false;
            adj[ver] = a_list;
        }
        T v, w;
        for (int i = 0; i < e; i++){
            in >> v >> w;
            add_edge(v, w); 
        }
    }

    ~Graph() {}

    int vertices() {return V;}  // num of vertices
    int edges() {return E;}     // num of edges
    T const& vertex(int v) {return vert[v];} // vertex v

    void add_edge(T _elem_v, T _elem_w){
        adj[_elem_v].push_front(_elem_w);
        adj[_elem_w].push_front(_elem_v);
        E++;
    }

    std::forward_list<T> const& adjacent(T v){
        return adj[v];
    } 

    // helper functions
    bool is_marked(T v){
        return marked[v];
    }

    void print_graph(){
        std::cout << V << " vertices, " << E << " edges"<<"\n";
            for (auto x: adj){
                std::cout << x.first << ": ";
                for (auto it = x.second.begin(); it != x.second.end(); ++it){
                    std::cout << *it << " ";
                } 
                std::cout<< "\n";
            } 
            std::cout << '\n';
    }       
};

Jak widać sporo dodatkowych struktur, jak marked, czy vert, ale za to upraszczamy kod klienta, może być umieszczony w oddzielnym module czy przestrzeni nazw (nie potrzebuje dodatkowych klas, ale coś za coś, paskudna mutowalność grafu:/). Np., sprawdzanie czy spójny(przy założeniu, że mamy działający dfs), może wyglądać tak:

template <class T>
bool is_connected(Graph<T>& g, T v){
    dfs(g, v);
    if (g.count == g.vertices()) return true;
    else return false;
}

Będa updaty:)
Podziękowania dla:

Pijak

Gnwno kod, ale zawsze można się pochwalić.

lion137

Bardzo dziękujemy za merytoryczny komentarz:-)

lion137
2018-11-22 12:16

TDD w C++

Tym razem:). Potestowałem sobie (przyda się) googletest. Powiem, że całkiem przyzwoite extreme programming wychodzi:)
I powstała, a jakże immutable list w stylu funkcyjnym, w C++ (jeśli w ogóle mając wskaźniki coś może być immutable:)). Kod, jeszcze nie uczesany:), tutaj:
https://github.com/lion137/C-Immutable-List/tree/master
#theory #c++ #functional

slsy

@lion137 nie znam się dobrze na scali, ale wydaje mi sie, to o czym napisałeś to po prostu hack pamiętający czasy Javy pozwalający na osiągnięcie sum type wykorzystując abstrakcyjny trait i jego rózne implementacje. Zresztą Dotty wprowadza bardziej elegancki mechanizm: https://dotty.epfl.ch/docs/reference/enums/adts.html . Niestety w C++ nie możesz ograniczyć możliwości dziedziczenia po interfejsie ImmutableList tak jak w Scali, więc musisz się spodziewać, że ktoś będzie mógł stworzyć własną implementację tego typu.

lion137

@jarekr000000: Under development, may change:)

lion137
2018-11-05 22:44

Kto mówi, że

W C nie da się używać TDD:) Oto taki sobie side project, toy project w C: generic linked list (via void pointers). Napisany właśnie w stylu TDD:
https://github.com/lion137/C_[...]ee/master/generic_linked_list
A po co to? Palce mi się zastały od tego pisania w Pythongu, musiałem coś w normalnym języku stworzyć:-P
#C #theory

vpiotr

Czulem ze te intuicyjne inaczej nazwy maja dluga historie.

lion137

Czyli sposób jaki został wybrany w Common Lisp. Osobiście taki design mi się podoba, mamy swoje is_empty które się z tym nie kłóci. Natomiast, rzeczywiście, w dzisiejszych czasach, programiści spodziewają się, że jak dojdzie do pop z pustej kolekcji, to nie dzieje się dobrze.

lion137
2018-11-02 20:41

Pakiety w Pythonie

Dawno nie było nic o Pythonie, to jest. Znowu dość praktycznie: pakiety. Jak tworzymy, jak są izolowane, jak wygląda w nich praca.
Jeśli chcemy napisać coś większego w Pythonie, dobrze jest zacząć od stworzenia środowiska wirtualnego:

➜  test git:(master) ✗ python3.7 -m virtualenv --no-site-packages venv
Using base prefix '/usr/local'
New python executable in /home/lion/code/workspace/python/tests/test/venv/bin/python3.7
Not overwriting existing python script /home/lion/Code/workspace/python/tests/test/venv/bin/python (you must use /home/lion/Code/workspace/python/tests/test/venv/bin/python3.7)
Installing setuptools, pip, wheel...
done.

Aktywujemy je:

➜  test git:(master) ✗ source venv/bin/activate 

I jak teraz sprawdzimy, jakie mamy zainstalowane pakiety, lista jest pusta (pip, jako nakładka na setuptools służy do instalacji, listowania, itp...):

(venv) ➜  venv git:(master) ✗ pip freeze
(venv) ➜  venv git:(master) ✗ 

Lista jest pusta. Kilka uwag, Komenda --no-site-packages słuźy do kompletnego odizolowania środowiska (tak dla pewności, Python domyślnie może tego nie robić). W naszym mini unixowym drzewie katalogów mamy plik wykonawczy Pythona (i kilka innych), jest to dokładnie ta sama wersja, którą użyliśmy do stworzenia środowiska(i będzie do niego przypisana na zawsze). Jesteśmy odseparowani od reszty, cokolwiek nie robimy jesteśmy w naszym venv.
Instalujemy coś:

(venv) ➜  venv git:(master) ✗ pip install pytest
Collecting pytest
 Successfully installed atomicwrites-1.2.1 attrs-18.2.0 more-itertools-4.3.0 pluggy-0.8.0 py-1.7.0 pytest-3.9.3 six-1.11.0
(venv) ➜  venv git:(master) ✗ python
Python 3.7.0 (default, Sep 15 2018, 00:51:17) 
[GCC 5.4.0 20160609] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pytest
>>> pytest.__file__
'/home/lion/code/workspace/python/tests/test/venv/lib/python3.7/site-packages/pytest.py'
>>> 

Jak widać pytest jest w naszym lokalnym (nie globalnie) venv.
Teraz pip freeze:

(venv) ➜  venv git:(master) ✗ pip freeze
atomicwrites==1.2.1
attrs==18.2.0
more-itertools==4.3.0
pluggy==0.8.0
py==1.7.0
pytest==3.9.3
six==1.11.0
(venv) ➜  venv git:(master) ✗ 

Zależności, jeżeli przekierujemy to wyjście do pliku:

(venv) ➜  venv git:(master) ✗ pip freeze > requirements.txt && cat requirements.txt
atomicwrites==1.2.1
attrs==18.2.0
more-itertools==4.3.0
pluggy==0.8.0
py==1.7.0
pytest==3.9.3
six==1.11.0
(venv) ➜  venv git:(master) ✗ 

To mamy gotową listę zależności, teraz możemy ją spakować (plus cokolwiek co mamy w projekcie) i w nowym środowisku wirtualnym możemy my (lub ktokolwiek inny), zainstalować nasz projekt:

(venv) ➜  venv git:(master) ✗ deactivate 
➜  venv git:(master) ✗ python3.7 -m virtualenv --no-site-packages venv2     
Using base prefix '/usr/local'
New python executable in /home/lion/code/workspace/python/tests/test/venv/venv2/bin/python3.7
Also creating executable in /home/lion/code/workspace/python/tests/test/venv/venv2/bin/python
Installing setuptools, pip, wheel...
done.
➜  venv git:(master) ✗ ls
bin  include  lib  requirements.txt  venv2
➜  venv git:(master) ✗ source venv2/bin/activate
(venv2) ➜  venv git:(master) ✗ pip install -r requirements.txt
Successfully installed atomicwrites-1.2.1 attrs-18.2.0 more-itertools-4.3.0 pluggy-0.8.0 py-1.7.0 pytest-3.9.3 six-1.11.0
(venv2) ➜  venv git:(master) ✗ 

Ciekawostka: jeśli zedytujemy requirements.txt, zmienimy wersję jakiegoś pakietu i ponownie odpalimy pip install -r req..., to pip zainstaluje nową i usunie poprzednią, czyli nie mamy opcji zrobienia chaosu z różnymi wersjami pakietu w projekcie:)
To tyle, jak wszystko poszło dobrze:), to powinno śmigać:). Pochwalcie się czymś na githubie:)
#Python #python #theory

KageYoshi

W windowsie ręczna instalacja jest tragedią :/ setki dziwnych błędów, u każdego spowodowane czym innym. Dla własnego zdrowia psychicznego używam już tylko condy.

lion137
2018-09-27 11:29

Typowany Python?

Witam, po przerwie, jak w tytule, co tam z tzw. "type annotations" w Pythonie, otóz mają się nieźle; mypy coraz lepiej zintegrowane z trójką (najlepiej działa 3.5>=). Praca polega (przynajmnije na Linux) na zainstalowaniu mypy:

$ python -m pip install mypy

Co ono robi? Sprawdza typy zgodnie z naszymi annotacjami, które [te annotacje] nie mają wpływu na wykonanie programu, tzn.:

$ python my_file.py

będzie działać, chociaż:

$ mypy my_file.py

zwróci błędy. Czyli do istniejącego kodu, można je dodawać stopniowo i sprawdzać, testować, etc... Kilka przykładów, jak ktoś chce to znajdzie resztę:)

Deklaracja wygląda tak:

from typing import List
x: List[int] = []
x.append(1)

mypy zaprotestuje, jeśli dodamy do listy cokolwiek innego niż int.

Mamy generyki i type variables :

from typing import Generic, TypeVar

A = TypeVar('A')

class Tree(Generic[A]):
    def __init__(self, value: A, left: 'Tree[A]', right: 'Tree[A]') -> None:
        self.payload = value
        self.left_tree = left
        self.right_tree = right

Żeby interpreter się nie wysypał, parametry typów w __init__ umieszczamy w apostrofach. Funkcja zwracająca None, to odpowiednik void z Javy, czy C/C++.

Kolejna funkcja to mapa (w sensie jak w programowaniu funkcyjnym), zwracająca stream - w Pytonie są to generatory. Callable[[A], A] - przyjmuje argument typu A i zwracająca takiż element : (A) -> A. Drugi argument Iterable[A] jest supertypem dla listy (może to być, również Dict, Tuple, Range):

from typing import Callable, Iterable, Generator

def my_map_generator(f: Callable[[A], A], xs: Iterable[A]) -> Generator:
    for x in xs:
        yield f(x)

Jak akurat nie chcemy strumienia, to proszę bardzo:

def my_map(f: Callable[[A], A], xs: List[A]) -> List[A]:
    for ind, x in enumerate(xs):
        xs[ind] = f(x)
    return xs

Ta wersja bierze taką samą funkcję, jak poprzednia, i listę typu A, a zwraca listę typu A.
To tyle, opcja warta rozważenia w dużych projektach, przy skryptach, to typowy "overkill". Pozdrawiam.
#Python #python #theory #programming

Leroy

Pracowałem w korpo pythonowym projekcie na Python'ie 2 gdzie w sumie wdrazylem MyPy.
Samo narzedzie jest fajne, ale dalej, Python byl, jest i bedzie dynamicznie typowany.
Spielismy to z CI, napisalismy wlasne toole, ktore pozwalaly na ustawianie roznych poziomow skanowania i walidacji w roznych podprojektach i pakietach.
Ogolnie wg mnie MyPy dziala, ale to tylko proteza, ktora mimo ze pomaga, ma swoje wady i czesto natrafialismy na sciane gdzie musielismy pisac wlasne stub'y.
Najwiekszym plusem bylo to, ze PyCharm po prostu zaczal dzialac.
O wiele lepszym rozwiazaniem jest po prostu statycznie typowany jezyk.

lion137
2018-07-24 12:32

Python, Dwie Sztuczki i Algorytm

Dawno nic nie pisałem, tak to już jest, czasem człowiek jest busy, albo lazy:). Do rzeczy, dwa małe, praktyczne, triki (jakby ktoś nie znał). Pierwszy, otwieranie pliku przy pomocy iteratora(załóżmy, że mamy otwarty plik, na przykład z logami, jako f):

def good_lines(f):
    for line in f:
        line = line.strip("\n ")
        if line.startswith('#'):
            continue
        if not line:
            continue
        yield line

Utworzony iterator, omija linie komentarzy, puste linie, obcina spacje i znaki końca linii, teraz można już robić coś z interesującymi nas partiami pliku:

for line in good_lines(f):
    do_stuff(line)

Argumentem do funkcji good_lines(f) niekoniecznie musi być plik, może to być lista, zbiór, generalnie każda struktura wspierająca iterację.

Dalej,
W Pythonie nie ma łatwej możliwości przerwania podwójnej pętli, trzeba tworzyć dodatkowe zmienne kontrolne, warunki, co raczej zaciemnia kod, utrudnia testowanie, itd... Np.:

for y in range(height):
    for x in range(width):
        if sth_found(x, y):
            break # nie przerwie obu pętli

Ale gdy tylko się da możemy zrobić tak:

def 2d_range(width, height):
    for y in range(height):
        for x in range(width):
            yield x, y

for col, row in 2d_range(width, height):
    if sth_found(col, row):
        break # Cool, works!:)

Na koniec Fibonacci, jak można znaleźć w necie jest kilka sposobów liczenia tej funkcji, od najgorszego:

  • naiwny rekurencyjny, O(n) i pamięć również O(n), więc po chwili rozsadza nam stos;
  • dynamic programming, dalej O(n), ale w stałej pamięci, więc można jako tako używać;
  • matrix exponentiation, lepiej O(logn), stała pamięć, bardzo szybki;
  • double fibonacci identities, ulepszenie poprzedniego, złożność ta sama, ale mniej operacji stałych.

Wszystko to, dobrze opisane, np., tutaj. Właśnie ten przedostatni algorytm poniżej, w zasadzie jest to zapisanie identyczności z tytułu, trzeba tylko pamiętać o exponentiation by squaring, które można też, jak widać, stosować dla macierzy. Może się okazać, że ten tip jest nawet praktyczny (interview), przed kliknięciem zachęcam do samodzielnego zakodowania, link tutaj:
https://nbviewer.jupyter.org/[...]trix_Exponetiation_Fibo.ipynb

#Python #python #algorithms #theory

lion137

Pewnie, że się da, ale czasem mogą być problemy z wysłaniem i mutowaniem zmiennych (nie można sobie wyslać adresu jak w C++ :)), więcej klepania, a tak jest prościej i przejrzyściej.

Afish

No ale jakie wysyłanie zmiennych? Masz domknięcie, nie musisz nic przekazywać przez parametry. Pokaż może przykład, gdzie to nie zadziała.

lion137
2018-06-06 02:27

Różniczkowanie Symboliczne

Ha, ha, dawno już nie było teorii, wszyscy stęsknieni:P; to do roboty: Symbolic differentiation!
Temat jest bliski moim zainteresowaniom, tzn. AI; więc, studiując słynne PAIG plus (obowiązkowo!) SICP natkąłem się na jak w temacie.
Jako że ciężko mi było nadążyć z oryginalnym kodem w Lispie:

(defun deriv-poly (p x)
"Return the derivative, dp/dx, of the polynomial p."
" If P is a number or a polynomial with main-var > x,
" then p is free of x, and the derivative is zero;
" otherwise do real work.
" But first, make sure X is a simple variable,
" of the form #(X 0 1).
(assert (and (typep x 'polynomial) (= (degree x) 1)
(eql (coef x 0) 0) (eql (coef x 1) 1)))
(cond
« numberp p) 0)
«var> (main-var p) (main-var x)) 0)
(var= (main-var p) (main-var x))
d(a + bx + cx2 + dx3)/dx = b + 2cx + 3dx A 2
;; So. shift the sequence p over by 1. then
;; put x back in. and multiply by the exponents
(let «r (subseq pI)))
(setf (main-var r) (main-var x))
(loop for i from 1 to (degree r) do
(setf (coef r i) (poly*poly (+ i 1) (coef r i))))
(normalize-poly r)))
(t ;; Otherwise some coefficient may contain x. Ex:
d(z + 3x + 3zx2 + z2x3)/dz
;; = 1 + 0 + 3x2 + 2zx A 3
;; So copy p. and differentiate the coefficients.
(let «r (copy-poly p)))
(loop for i from 0 to (degree p) do
(setf (coef r i) (deriv-poly (coef r i) x)))
(normalize-poly r)))))

Plus pomoc z SICP (Scheme):

(define (deriv exp var)
(cond ((number? exp) 0)
((variable? exp) (if (same-variable? exp var) 1 0))
((sum? exp) (make-sum (deriv (addend exp) var)
(deriv (augend exp) var)))
((product? exp)
(make-sum
(make-product (multiplier exp)
(deriv (multiplicand exp) var))
(make-product (deriv (multiplier exp) var)
(multiplicand exp))))
(else
(error "unknown expression type: DERIV" exp))))

(Nawiasem pisząc, formatka 4programmers ssie (nie formatuje Lispu)), to wysmażyłem swój w Pythonie (kompilowalny pseudokod :-D). Idea jest prosta: Na najniższym poziomie - pochodna liczby jest zerem, pochodna zmiennej jest jedynką, jeśli to ta sama zmienna, zerem jeśli inna (jak stała). Problemem jest, jak to odwzorować w notacji infiksowej (jak Python i inne języki głównonurtowe). Okazuje się, że można wykorzystać algorytm parsowania formuły do (binarnego na razie!!) drzewa składniowego.
Idea: Pochodna sumy jest sumą pochodnych, a pochodna iloczynu, jest sumą iloczynu pierwszego członu razy pochodna drugiego plus pochodna pierwszego razy drugi; doskonale wpisuje się to w drzewo składniowe:

def evaluate(tree):
    opers = {'+': sym_add, '-': op.sub, '*': sym_mul, '/': op.truediv}

    leftT = tree.getLeftChild()
    rightT = tree.getRightChild()

    if leftT and rightT:
        fn = opers[tree.getRootVal()]
        return fn(evaluate(leftT), evaluate(rightT))
    else:
        return tree.getRootVal()

def evaluate_parse_tree(tree, var):
    opers = {'+': sym_add, '-': op.sub, '*': sym_mul, '/': op.truediv}

    leftT = tree.getLeftChild()
    rightT = tree.getRightChild()

    if leftT and rightT:
        fn = opers[tree.getRootVal()]
        # changes start:
        if fn is sym_mul:
            fn = sym_add
            return fn(sym_mul(evaluate(leftT), evaluate_parse_tree(rightT, var)), sym_mul(evaluate_parse_tree(leftT, var),
                    evaluate(rightT)))
        else:
            return fn(evaluate_parse_tree(leftT, var), evaluate_parse_tree(rightT, var))
    else:
        return derrive(tree.getRootVal(), var)

def derrive(exp, var):
    if is_number(exp): return 0
    elif is_variable(exp):
        return 1 if same_variable(exp, var) else 0

Dwie funkcje evaluate_parse_tree i evaluate działają zgodnie z powyższym schematem; chociaż widać tu dużo braków: Trzeba upraszczać końcowe wyrażenia (albo poprawić sym-add i sym_mul). Lecz pomysł działa. Kod (już nie GitHubie, bo go nie ma:-D), z funkcjami sym_add i tak dalej (w razie nieścisłości zapraszam do dyskusji) tutaj: https://bitbucket.org/lion137[...];fileviewer=file-view-default
Działa?

 form = "((1 + (2 * x)) * (1 + x))"
   p_tree = BinaryTree('')
   p_tree = build_parse_tree(form)
   print(evaluate_parse_tree(p_tree, "x"))# -> 1(1 + 2x) + (0 + 2 + 0x)*(1 + x)
   form = "(x + 1)"
   p_tree = BinaryTree('')
   p_tree = build_parse_tree(form)
   print(evaluate_parse_tree(p_tree, "x")) # -> 1

Jak widać tak (chociaż): Wyrażenie musi być poprawnie znawiasowane i forma wyniku przedstawia sporo do życzenia.
Ale droga jest właściwa! :-D)
#theory #Python

#theory #Python

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.