Sunday 27 January 2008

Gödel's beta function

I read a most incredible proof in the book, Lectures On The Curry Howard Isomorphism that I decided to created a program from it.


(defun factorial (n)
(if (= n 0) 1
(* n (factorial (- n 1)))))

(defun extended-euclidean-algorithm (a b)
;; function extended_gcd(a, b)
;; if (a mod b = 0)
;; return {0, 1}
;; else {x, y} := extended_gcd(b, a mod b)
;; return {y, x-y*(a div b)}
(if (= (mod a b) 0) (values 0 1)
(multiple-value-bind (x y)
(extended-euclidean-algorithm b (mod a b))
(values y (- x (* y (floor a b)))))))

(defun extended-euclidean-algorithm-prime (a b)
;; extended-euclidean-algorithm solves:
;; u*a + v*b = gcd
;; This one solves:
;; u*a - v*b = gcd
(multiple-value-bind (u v)
(extended-euclidean-algorithm a b)
(if (= b 1)
(values 1 (- a 1))
(values u (- v)))))

;; For every finite sequence k_0, k_1, ..., k_r there exist two numbers m, n, such that (β m n j) = k_j, for j = 0..r.
(defun β (m n j)
(mod m (+ (* n (+ j 1)) 1)))

;; Produce n m, such that (β m n j) = k_j
(defun β-solve (k)
(let* ((n (factorial (apply #'max (- (length k) 1) k)))
(a (loop for j below (length k) collect (+ (* n (+ j 1)) 1))))
(labels ((next-m (l m)
(if (= l (- (length k) 1)) m
(let ((a* (apply #'* (subseq a 0 (min (+ l 1) (length k))))))
(multiple-value-bind (u v)
(extended-euclidean-algorithm-prime (elt a (+ l 1)) a*)
(declare (ignore u))
(next-m (+ l 1)
(+ m (* (* (- m (elt k (+ l 1))) v) a*))))))))
(values (next-m 0 (elt k 0)) n))))


;; CL-USER> (multiple-value-bind (n m)
;; (β-solve (map 'list #'char-code "You win!"))
;; (map 'string #'code-char (loop for j below 8 collect (β n m j))))
;; "You win!"

No comments: