クライン群の極限集合を描画する

draw-limiting-set.lisp

(defparameter *inf* (* most-positive-fixnum #c(1 1)))

(defclass matrix ()
  ((a :accessor mat-a :initform 1 :initarg :a)
   (b :accessor mat-b :initform 0 :initarg :b)
   (c :accessor mat-c :initform 0 :initarg :c)
   (d :accessor mat-d :initform 1 :initarg :d)))

(defgeneric mat* (a b))
(defgeneric inverse (m))
(defgeneric mat-trace (m))

(defmethod mat* ((m1 matrix) (m2 matrix))
  (with-accessors ((m1-a mat-a) (m1-b mat-b) (m1-c mat-c) (m1-d mat-d)) m1
    (with-accessors ((m2-a mat-a) (m2-b mat-b) (m2-c mat-c) (m2-d mat-d)) m2
      (make-instance 'matrix
                     :a (+ (* m1-a m2-a) (* m1-b m2-c))
                     :b (+ (* m1-a m2-b) (* m1-b m2-d))
                     :c (+ (* m1-c m2-a) (* m1-d m2-c))
                     :d (+ (* m1-c m2-b) (* m1-d m2-d))))))

(defmethod mat* ((coefficient complex) (m1 matrix))
  (with-accessors ((m1-a mat-a) (m1-b mat-b) (m1-c mat-c) (m1-d mat-d)) m1
    (make-instance 'matrix
                   :a (* coefficient m1-a)
                   :b (* coefficient m1-b)
                   :c (* coefficient m1-c)
                   :d (* coefficient m1-d))))


(defmethod mat* ((coefficient integer) (m1 matrix))
  (with-accessors ((m1-a mat-a) (m1-b mat-b) (m1-c mat-c) (m1-d mat-d)) m1
    (make-instance 'matrix
                   :a (* coefficient m1-a)
                   :b (* coefficient m1-b)
                   :c (* coefficient m1-c)
                   :d (* coefficient m1-d))))

(defmethod inverse ((m matrix))
  (let ((determinant (- (* (mat-a m) (mat-d m))
			(* (mat-b m) (mat-c m)))))
    (if (= determinant 0)
	(make-instance 'matrix
		       :a *inf*
		       :b *inf*
		       :c *inf*
		       :d *inf*)
	(make-instance 'matrix
		       :a (/ (mat-d m) determinant)
		       :b (/ (* -1 (mat-b m)) determinant)
		       :c (/ (* -1 (mat-c m)) determinant)
		       :d (/ (mat-a m) determinant)))))

(defmethod mat-trace ((m matrix))
  (+ (mat-a m) (mat-d m)))

(defmethod print-object ((m matrix) s)
  (write-char #\{ s)
  (princ (mat-a m) s)
  (write-char #\, s)
  (princ (mat-b m) s)
  (write-char #\LineFeed s)
  (write-char #\space s)
  (princ (mat-c m) s)
  (write-char #\, s)
  (princ (mat-d m) s)
  (write-char #\} s))

(defun fix-plus ((m matrix))
  (with-accessors ((a mat-a) (b mat-b) (c mat-c) (d mat-d)) m
    (/ (+ (- a d) (sqrt (- (expt (mat-trace m) 2) 4)))
       (* c 2))))

(defun mobius-on-point ((m matrix) p)
  (with-accessors ((a mat-a) (b mat-b) (c mat-c) (d mat-d)) m
    (cond ((= *inf* c) (cond ((= c 0) *inf*)
			     (t (/ a c))))
          (t (let ((numerix (+ (* a p) b))
                   (denominator (+ (* c p) d)))
               (cond ((= denominator 0) *inf*)
                     (t (/ numerix denominator))))))))

(defun init (max-level epsilon)
  (defparameter *fix-point* (make-array '(5 4)))
  (defparameter *gens* (make-array '(5)))
  (defparameter *max-level* max-level)
  (defparameter *word* (make-array (list (+ *max-level* 1))))
  (defparameter *tags* (make-array (list (+ *max-level* 1))))
  (defparameter *epsilon* epsilon)
  (defparameter *line-points-list* ())
  (defparameter *level* 1))

(defun print-tags ()
  (loop for i from 1 to *level* do
       (princ (aref *tags* i)))
  (print ""))

(defun set-gens ((t_a complex) (t_b complex))
  (let* ((t_ab (/ (- (* t_a t_b)
		     (sqrt (- (* (expt t_a 2) (expt t_b 2))
			      (* 4 (+ (expt t_a 2) (expt t_b 2))))))
		  2))
        (z0 (/ (* (- t_ab 2) t_b)
	       (+ (- (* t_b t_ab) (* t_a 2))
		  (* t_ab #c(0 2)))))
        (a (make-instance 'matrix
                               :a (/ t_a 2)
                               :b (/ (+ (- (* t_a t_ab) (* t_b 2)) #c(0 4))
				     (* z0 (+ (* t_ab 2) 4)))
                               :c (/ (* z0
					(- (* t_a t_ab) (* t_b 2) #c(0 4)))
				     (- (* t_ab 2) 4))
                               :d (/ t_a 2)))
        (b (make-instance 'matrix
                               :a (/ (- t_b #c(0 2)) 2)
                               :b (/ t_b 2)
                               :c (/ t_b 2)
                               :d (/ (+ t_b #c(0 2)) 2)))
	(_A (inverse a))
	(_B (inverse b)))

    (setf (aref *fix-point* 1 1) (fix-plus (reduce #'mat* (list b _A _B a))))
    (setf (aref *fix-point* 1 2) (fix-plus a))
    (setf (aref *fix-point* 1 3) (fix-plus (reduce #'mat* (list _B _A b a))))

    (setf (aref *fix-point* 2 1) (fix-plus (reduce #'mat* (list _A _B a b))))
    (setf (aref *fix-point* 2 2) (fix-plus b))
    (setf (aref *fix-point* 2 3) (fix-plus (reduce #'mat* (list  a _B _A b))))

    (setf (aref *fix-point* 3 1) (fix-plus (reduce #'mat* (list _B a b _A))))
    (setf (aref *fix-point* 3 2) (fix-plus _A))
    (setf (aref *fix-point* 3 3) (fix-plus (reduce #'mat* (list b a _B _A))))

    (setf (aref *fix-point* 4 1) (fix-plus (reduce #'mat* (list a b _A _B))))
    (setf (aref *fix-point* 4 2) (fix-plus _B))
    (setf (aref *fix-point* 4 3) (fix-plus (reduce #'mat* (list _A b a _B))))

    (setf (aref *gens* 1) a)
    (setf (aref *gens* 2) b)
    (setf (aref *gens* 3) _A)
    (setf (aref *gens* 4) _B)

    (setf (aref *word* 0) (make-instance 'matrix))
    (setf (aref *word* 1) (aref *gens* 1))
    (setf (aref *tags* 1) 1)))

(defun go-forward ()
  (if (eq (branch-termination-p) nil)
      (progn
	(incf *level*)
	(setf (aref *tags* *level*)
	      (abs (mod (+ (aref *tags* (- *level* 1)) 1) 4)))
	(if (= (aref *tags* *level*) 0)
	    (setf (aref *tags* *level*) 4))
	(setf (aref *word* *level*)
	      (mat* (aref *word* (- *level* 1)) (aref *gens* (aref *tags* *level*))))
	(go-forward))))
(compile 'go-forward)


(defun go-backward ()
  (decf *level*)
  (if (and (/= *level* 0) (not (available-turn-p)))
      (go-backward)))
(compile 'go-backward)

(defun available-turn-p ()
  (let ((next-next-tag (abs (mod (+ (aref *tags* *level*) 2) 4)))
	(pre-tag (- (aref *tags* (+ *level* 1)) 1)))
    (if (= next-next-tag 0) (setf next-next-tag 4))
    (if (= pre-tag 0) (setf pre-tag 4))
    (if (= next-next-tag pre-tag)
	nil
	t)))

(defun turn-and-go-forward ()
  (setf (aref *tags* (+ *level* 1))
	(abs (mod (- (aref *tags* (+ *level* 1)) 1) 4)))
  (if (= (aref *tags* (+ *level* 1)) 0)
      (setf (aref *tags* (+ *level* 1)) 4))
  (if (= *level* 0)
      (setf (aref *word* 1) (aref *gens* (aref *tags* 1)))
      (setf (aref *word* (+ *level* 1))
	    (mat* (aref *word* *level*)
		  (aref *gens* (aref *tags* (+ *level* 1))))))
  (incf *level*))

(defun branch-termination-p ()
  (let ((z (make-array '(4))))
    (loop for i from 1 to 3 do
	 (setf (aref z i)
	       (mobius-on-point (aref *word* *level*)
				(aref *fix-point* (aref *tags* *level*) i))))
    (if (or (= *level* *max-level*)
	    (and (<= (abs (- (aref z 2) (aref z 1))) *epsilon*)
		 (<= (abs (- (aref z 3) (aref z 2))) *epsilon*)))
	(setf *line-points-list*
	      (append *line-points-list* (list (aref z 1) (aref z 2)  (aref z 2) (aref z 3))))
	nil)))

(defun explore-limiting-set-with-dfs ()
  (labels ((explore ()
	     (go-forward)
	     (go-backward)
	     (turn-and-go-forward)
	     (if (or (/= *level* 1) (/= (aref *tags* 1) 1))
		 (explore))))
    (explore)))
(compile 'explore-limiting-set-with-dfs)

(defun get-limiting-set (t_a t_b max-level epsilon)
  (init max-level epsilon)
  (set-gens t_a t_b)
  (explore-limiting-set-with-dfs))

(ql:quickload "cl-glut")
(defparameter *width* 500)
(defparameter *height* 500)
(defparameter *magnification* 8)
(defclass my-window (glut:window)
  ()
  (:default-initargs :title "limiting set" :width *width* :height *height*
		     :mode '(:single :rgb :depth)))

(defmethod glut:display-window :before ((w my-window))
  (gl:clear-color 1 1 1 0)
  (gl:matrix-mode :projection)
  (gl:load-identity)
  (gl:ortho 0 *width* *height* 0 -1 1))

(defun set-line-vertexes (lines)
  (mapcar #'(lambda (point)
	      (gl:vertex (* *magnification* (realpart point))
			 (* *magnification* (imagpart point))
			 0))
	  lines))

(defmethod glut:display ((window my-window))
  (gl:clear :color-buffer-bit)
  (%gl:color-3f 0 0 0)
  (gl:push-matrix)
  (gl:translate (/ *width* 2) (/ *height* 2) 0)
  (gl:begin :lines)
  (set-line-vertexes *line-points-list*)
  (gl:end)
  (gl:pop-matrix)
  (gl:flush))

(defun draw-limiting-set ()
  (glut:display-window (make-instance 'my-window)))


;; (defun go-forward ()
;; 	(incf *level*)
;; 	(setf (aref *tags* *level*)
;; 	      (abs (mod (+ (aref *tags* (- *level* 1)) 1) 4)))
;; 	(if (= (aref *tags* *level*) 0)
;; 	    (setf (aref *tags* *level*) 4))
;; 	(setf (aref *word* *level*)
;; 	      (mat* (aref *word* (- *level* 1)) (aref *gens* (aref *tags* *level*)))))

;; (defun go-backward ()
;;   (decf *level*))

;; (defun explore-limiting-set-with-dfs ()
;;   (labels ((explore ()
;; 	     (loop while (eq (branch-termination-p) nil) do
;; 		  (go-forward))
;; 	     (go-backward)
;; 	     (loop while (and (/= *level* 0) (not (available-turn-p))) do
;; 		  (go-backward))
;; 	     (turn-and-go-forward)))
;;     (explore)
;;     (loop while (or (/= *level* 1) (/= (aref *tags* 1) 1)) do
;; 	 (explore))))