diff --git a/misc/peano.pl b/misc/peano.pl
new file mode 100644
index 0000000..4e3ef56
--- /dev/null
+++ b/misc/peano.pl
@@ -0,0 +1,82 @@
+% Peano axioms in Prolog
+%% Definition of a natural number
+nat(o).
+nat(succ(N)) :- nat(N).
+
+%% Operations
+
+add(X, o, X) :- nat(X), !.					% X + o = X
+add(X, succ(Y), succ(R)) :- add(X, Y, R). 	% X + succ(Y) = X + Y
+
+mul(X, o, o) :- nat(X), !. 					% X * o = o
+mul(X, succ(Y), R) :- 						% X * succ(Y) = X + (X * Y)
+	mul(X, Y, Tmp),
+	add(X, Tmp, R).
+
+% https://math.stackexchange.com/a/461959
+sub(X, o, X) :- nat(X), !.					% X - o = X
+sub(succ(X), succ(Y), R) :- sub(X, Y, R).	% succ(X) - succ(Y) = X - Y
+% sub(succ(succ(o)), succ(o), R)
+% sub(succ(o), o, R)
+%    nat(succ(o))
+% R = succ(o)
+
+fact(N, R) :- fact(N, succ(o), R).
+fact(o, R, R) :- !.
+fact(N, Acc, R) :-
+	sub(N, succ(o), NN),
+	mul(Acc, N, AccN),
+	fact(NN, AccN, R).
+
+peano2dec(N, R) :- peano2dec(N, 0, R).
+peano2dec(o, R, R) :- !.
+peano2dec(N, C, R) :-			% "list length"
+    sub(N, succ(o), NN),		% cdr
+    CC is C + 1,
+    peano2dec(NN, CC, R).
+
+dec2peano(N, R) :- dec2peano(N, o, R).
+dec2peano(0, R, R) :- !.
+dec2peano(N, Acc, R) :-			% "list of N elements"
+    add(Acc, succ(o), AccX), 	% cons
+    NN is N - 1,
+    dec2peano(NN, AccX, R).
+
+:- begin_tests(peano).
+
+test(addition) :-
+    add(o, o, R1), R1 = o,
+    add(succ(o), o, R2), R2 = succ(o),
+    add(succ(succ(o)), succ(succ(o)), R3), R3 = succ(succ(succ(succ(o)))).
+
+test(multiplication) :-
+    mul(o, o, R1), R1 = o,
+    mul(succ(succ(o)), o, R2), R2 = o,
+    Nine = succ(succ(succ(succ(succ(succ(succ(succ(succ(o))))))))),
+    mul(succ(succ(succ(o))), succ(succ(succ(o))), R3), R3 = Nine.
+
+test(substraction) :-
+    sub(o, o, R1), R1 = o,
+    sub(succ(succ(succ(o))), succ(succ(o)), R2), R2 = succ(o),
+    sub(succ(succ(succ(o))), succ(succ(succ(o))), R3), R3 = o.
+
+test(peano_to_decimal) :-
+    peano2dec(succ(succ(o)), R1), R1 =:= 2,
+    peano2dec(o, R2), R2 =:= 0,
+	peano2dec(succ(succ(succ(succ(o)))), R3), R3 =:= 4.
+
+test(decimal_to_peano) :-
+    dec2peano(3, R1), R1 = succ(succ(succ(o))),
+    dec2peano(0, R2), R2 = o,
+    dec2peano(4, R3), R3 = succ(succ(succ(succ(o)))).
+    
+test(factorial_p) :-
+    fact(succ(succ(succ(succ(succ(o))))), PR),
+    peano2dec(PR, DR),
+    DR =:= 120.
+
+test(extra) :-
+    Num is 5,
+    dec2peano(Num, PN), fact(PN, PR), peano2dec(PR, DN), DN =:= 120.
+
+:- end_tests(peano).