home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / examples / polytop2.lsp < prev    next >
Encoding:
Lisp/Scheme  |  1993-10-23  |  22.0 KB  |  592 lines
  1. ;; Berechnung des Volumens konvexer Polytope
  2. ;; Bruno Haible 5.9.1990, 31.5.1991
  3.  
  4. ;; Ein Polytop sei als konvexe Hülle endlich vieler Punkte (z.B. seiner Ecken)
  5. ;; gegeben. Es wird eine Triangulation und das Volumen berechnet.
  6.  
  7. #|
  8. (require 'lgs_Q) ; Funktion solve-homLGS
  9. |#
  10. ;; Lösen linearer Gleichungssysteme über Q
  11. ;; Bruno Haible 25.4.1989
  12.  
  13. ; Elimination:
  14. ; Löst ein homogenes LGS, d.h. bestimmt den Kern einer Matrix.
  15. ; > A: eine mxn-Matrix rationaler Zahlen
  16. ; < eine Liste der Basisvektoren des Kernes von A,
  17. ;   das sind jeweils Vektoren der Länge n aus teilerfremden ganzen Zahlen
  18. ; < 2.Wert: Anzahl der eliminierten Zeilen (>=0, <=m)
  19. (defun solve-homLGS (A)
  20.   (let ((m (array-dimension A 0))
  21.         (n (array-dimension A 1)))
  22.     ; A kopieren:
  23.     (setq A
  24.       (let ((B (make-array `(,m ,n))))
  25.         (dotimes (i m) (dotimes (j n) (setf (aref B i j) (aref A i j))))
  26.         B
  27.     ) )
  28.     (let ((spaltenabsätze nil) 2.Wert)
  29.       ; A auf Treppennormalform bringen:
  30.       (let ((i 0)) ; i = Zeile, ab der noch ausgeräumt werden muß
  31.         (dotimes (j n) ; j = untersuchte Spalte
  32.           ; Suche betragsgrößtes Element unterhalb (i,j):
  33.           (let ((i0 nil))
  34.             (do ((i1 i (1+ i1))
  35.                  (v 0))
  36.                 ((= i1 m))
  37.               (if (> (abs (aref A i1 j)) v)
  38.                 (setq i0 i1 v (abs (aref A i1 j)))
  39.             ) )
  40.             ; i0 = Zeile des betragsgrößten Elements, oder
  41.             ; i0 = nil falls alle fraglichen Elemente =0 sind.
  42.             (if i0
  43.               (progn
  44.                 ; Zeilen i und i0 vertauschen
  45.                 (unless (= i i0)
  46.                   (dotimes (j0 n) (rotatef (aref A i j0) (aref A i0 j0)))
  47.                 )
  48.                 ; Zeile i normalisieren
  49.                 (let ((x (/ (aref A i j))))
  50.                   (dotimes (j0 n) (setf (aref A i j0) (* (aref A i j0) x)))
  51.                 )
  52.                 ; Zeilen unterhalb von Zeile i ausräumen durch Subtraktion
  53.                 ; eines Vielfaches von Zeile i:
  54.                 (do ((i1 (1+ i) (1+ i1)))
  55.                     ((= i1 m))
  56.                   (let ((x (aref A i1 j)))
  57.                     (dotimes (j1 n)
  58.                       (setf (aref A i1 j1) (- (aref A i1 j1) (* x (aref A i j1))))
  59.                 ) ) )
  60.                 ; Spaltenabsatz
  61.                 (push i spaltenabsätze)
  62.                 (incf i)
  63.               )
  64.               (push nil spaltenabsätze)
  65.         ) ) )
  66.         (setq 2.Wert (- m i))
  67.       )
  68.       ; A hat Treppennormalform, z.B.:
  69.       ;         1 .............
  70.       ;         0 1 ...........
  71.       ;         0 0 0 0 1 .....
  72.       ;         0 0 0 0 0 0 0 0
  73.       ; Dann enthält die Liste der Spaltenabsätze, in umgekehrter Reihenfolge:
  74.       ;         0 1 n n 2 n n n
  75.       ; (also die Zeilennummer beim Absatzanfang, n=NIL sonst)
  76.       ; Nun den Kern aufbauen:
  77.       (let ((Kern nil))
  78.         (do ((j (1- n) (1- j))
  79.              (L spaltenabsätze (cdr L)))
  80.             ((< j 0))
  81.           (if (car L)
  82.             ; Spaltenabsatz in Spalte j
  83.             (dolist (B Kern)
  84.               (do ((i (car L))
  85.                    (s 0)
  86.                    (j1 (1+ j) (1+ j1)))
  87.                   ((= j1 n) (setf (aref B j) (- s)))
  88.                 (setq s (+ s (* (aref A i j1) (aref B j1))))
  89.             ) )
  90.             ; kein Spaltenabsatz in Spalte j
  91.             (progn
  92.               (dolist (B Kern) (setf (aref B j) 0)) ; Komponente j frei wählbar
  93.               (let ((B (make-array n)))
  94.                 (setf (aref B j) 1)
  95.                 (do ((j1 (1+ j) (1+ j1)))
  96.                     ((= j1 n))
  97.                   (setf (aref B j1) 0)
  98.                 )
  99.                 (push B Kern) ; neues Basiselement #(..... 1 0 ... 0) im Kern
  100.         ) ) ) )
  101.         ; alle Vektoren des Kerns ganzzahlig und teilerfremd machen:
  102.         (dolist (B Kern)
  103.           (let ((nenner (do ((j 0 (1+ j))
  104.                              (x 1 (lcm x (denominator (aref B j)))))
  105.                             ((= j n) x)
  106.                ))       )
  107.             (dotimes (j n) (setf (aref B j) (* (aref B j) nenner)))
  108.           )
  109.           (let ((faktor (do ((j 0 (1+ j))
  110.                              (x 0 (gcd x (aref B j))))
  111.                             ((= j n) x)
  112.                ))       )
  113.             ; faktor /=0, da sonst B der Nullvektor wäre.
  114.             (dotimes (j n) (setf (aref B j) (/ (aref B j) faktor)))
  115.         ) )
  116.         (values Kern 2.Wert)
  117. ) ) ) )
  118.  
  119. ; Determinante:
  120. ; Bestimmt die Determinante einer quadratischen Matrix.
  121. ; > A: eine nxn-Matrix rationaler Zahlen
  122. ; < deren Determinante
  123. (defun det (A)
  124.   (let ((n (array-dimension A 0)))
  125.     (assert (= n (array-dimension A 1)))
  126.     ; A kopieren:
  127.     (setq A
  128.       (let ((B (make-array `(,n ,n))))
  129.         (dotimes (i n) (dotimes (j n) (setf (aref B i j) (aref A i j))))
  130.         B
  131.     ) )
  132.     ; A auf Treppennormalform bringen:
  133.     (let ((det 1)) ; det = bisherige Faktoren der Determinante
  134.       (dotimes (i n) ; i = Zeile, ab der noch ausgeräumt werden muß
  135.         (let ((j i)) ; j = untersuchte Spalte
  136.           ; Suche betragsgrößtes Element unterhalb (i,j):
  137.           (let ((i0 nil))
  138.             (do ((i1 i (1+ i1))
  139.                  (v 0))
  140.                 ((= i1 n))
  141.               (if (> (abs (aref A i1 j)) v)
  142.                 (setq i0 i1 v (abs (aref A i1 j)))
  143.             ) )
  144.             ; i0 = Zeile des betragsgrößten Elements, oder
  145.             ; i0 = nil falls alle fraglichen Elemente =0 sind.
  146.             (if i0
  147.               (progn
  148.                 ; Zeilen i und i0 vertauschen
  149.                 (unless (= i i0)
  150.                   (dotimes (j0 n) (rotatef (aref A i j0) (aref A i0 j0)))
  151.                   (setq det (- det))
  152.                 )
  153.                 ; Zeile i normalisieren
  154.                 (let ((x (aref A i j)))
  155.                   (setq det (* det x))
  156.                   (setq x (/ x))
  157.                   (dotimes (j0 n) (setf (aref A i j0) (* (aref A i j0) x)))
  158.                 )
  159.                 ; Zeilen unterhalb von Zeile i ausräumen durch Subtraktion
  160.                 ; eines Vielfaches von Zeile i:
  161.                 (do ((i1 (1+ i) (1+ i1)))
  162.                     ((= i1 n))
  163.                   (let ((x (aref A i1 j)))
  164.                     (dotimes (j1 n)
  165.                       (setf (aref A i1 j1) (- (aref A i1 j1) (* x (aref A i j1))))
  166.               ) ) ) )
  167.               (return-from det 0)
  168.       ) ) ) )
  169.       det
  170. ) ) )
  171.  
  172. (require 'simplex)
  173. ; Simplex-Algorithmus:
  174. ; Input: A, eine mxn-Matrix von Zahlen,
  175. ;        b, ein m-Vektor von Zahlen,
  176. ;        c, ein n-Vektor von Zahlen.
  177. ; Löst eine kanonische Optimierungsaufgabe :
  178. ; ( C ) Finde unter den y ∈ R^n mit A*y=b, y≥0 diejenigen mit <c,y> = Min
  179. ; ( C* ) Finde unter den z ∈ R^m mit x=A`*z+c≥0 diejenigen mit <b,z> = Min
  180. ; Siehe Koulikov, Algèbre et théorie des nombres, Chap. 6.1.
  181. ; Output: 1. Wert: ob lösbar, NIL oder T,
  182. ;         falls lösbar:
  183. ;           2. Wert: Lösung von ( C ), dem LP: Liste, bestehend aus
  184. ;                    - der Lösung y
  185. ;                    - dem optimalen Wert von <c,y>
  186. ;           3. Wert: Lösung von ( C* ), dem Dualproblem DP: Liste, bestehend aus
  187. ;                    - der Lösung x
  188. ;                    - den konstanten Gliedern Z für die z's
  189. ;                    - den Abhängigkeiten ZZ zwischen den z's
  190. ;                    - dem optimalen Wert von <b,z>
  191. ;                    Die allgemeine Lösung ist so aufgebaut:
  192. ;                    Falls ZZ[i,i]=T, ist die Variable z_i frei wählbar.
  193. ;                    Sonst ist ZZ[i,i]=NIL, und die Variable z_i errechnet sich als
  194. ;                    z_i = Z[i] + Summe(j=1,...,m mit Z[j,j]=T , ZZ[i,j]*z_j) .
  195. ;           4. Wert: ob die angegebenen Lösungen die gesamte Lösungsmenge
  196. ;                    ausmachen (T) oder ob möglicherweise noch eine andere
  197. ;                    Lösung existiert.
  198. ;                    In jedem Fall ist die Lösungsmenge konvex.
  199.  
  200. ;; Wir bewegen uns im R^d, d>=0. Ein Punkt ist ein Vektor aus seinen d
  201. ;; Koordinaten (rationale Zahlen).
  202.  
  203. ; Gleichheit von Punkten:
  204. (defun punkt= (p1 p2)
  205.   (equalp p1 p2) ; Vergleich aller Koordinaten
  206. )
  207.  
  208. ;; Ein Simplex ist eine (ungeordnete) Liste von k Punkten, k>=1, deren
  209. ;; konvexe Hülle die Dimension k-1 hat.
  210.  
  211. ; Dimension eines Simplex:
  212. (defun simplex-dim (simplex)
  213.   (1- (length simplex))
  214. )
  215.  
  216. ; Gleichheit von Simplizes:
  217. (defun simplex= (simplex1 simplex2)
  218.   (and (= (simplex-dim simplex1) (simplex-dim simplex2)) ; müssen gleiche Dimension haben
  219.        (dolist (punkt1 simplex1 t) ; und jede Ecke von simplex1
  220.          (unless (position punkt1 simplex2 :test #'punkt=) ; muß auch Ecke von simplex2 sein
  221.            (return nil)
  222.        ) )
  223. ) )
  224.  
  225. ; Bildet ein (k-1)-dimensionales Simplex aus k Punkten:
  226. (defun make-simplex (punkte)
  227.   (if (endp punkte)
  228.     (error "Zu wenig Punkte für ein Simplex.")
  229.     (if ; Löse das LGS x1*p1 + ... + xk*pk = 0, x1+...+xk=0
  230.       (let* ((d (length (first punkte)))
  231.              (k (length punkte))
  232.              (A (make-array (list (+ d 1) k))))
  233.         ; LGS aufstellen:
  234.         (dotimes (j k)
  235.           (let ((punkt (elt punkte j)))
  236.             (assert (= (length punkt) d))
  237.             (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  238.             (setf (aref A d j) 1)
  239.         ) )
  240.         ; LGS lösen. Der Kern muß leer sein, also (p1,1), ..., (pk,1) linear
  241.         ; unabhängig, also p1,...,pk affin unabhängig.
  242.         (null (solve-homLGS A))
  243.       )
  244.       punkte
  245.       (error "Punkte für ein Simplex sind affin abhängig: ~S" punkte)
  246. ) ) )
  247.  
  248. ; Stellt fest, ob ein Punkt in der affinen Hülle eines Simplex liegt.
  249. (defun on-affine-p (punkt simplex)
  250.   ; Löse das LGS x1*p1 + ... + xk*pk + x*p = 0, x1+...+xk+x=0
  251.   (let* ((d (length punkt))
  252.          (k (length simplex))
  253.          (A (make-array (list (+ d 1) (+ k 1)))))
  254.     ; LGS aufstellen:
  255.     (dotimes (j k)
  256.       (let ((punkt (elt simplex j)))
  257.         (assert (= (length punkt) d))
  258.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  259.         (setf (aref A d j) 1)
  260.     ) )
  261.     (dotimes (i d) (setf (aref A i k) (aref punkt i)))
  262.     (setf (aref A d k) 1)
  263.     ; LGS lösen:
  264.     (let ((kern (solve-homLGS A)))
  265.       (assert (<= (length kern) 1)) ; LGS durfte nicht mehrdeutig lösbar sein
  266.       (not (endp kern)) ; falls es unlösbar war, liegt der Punkt nicht drin.
  267. ) ) )
  268.  
  269. ; Stellt fest, ob zwei Punkte auf derselben Seite der affinen Hülle eines
  270. ; Simplex liegen.
  271. (defun on-same-side-p (punkt1 punkt2 simplex)
  272.   ; Löse das LGS x1*p1 + ... + xk*pk + y1*punkt1 + y2*punkt2 = 0, x1+...+xk+y1+y2=0
  273.   (let* ((d (length punkt1))
  274.          (k (length simplex))
  275.          (A (make-array (list (+ d 1) (+ k 2)))))
  276.     ; LGS aufstellen:
  277.     (dotimes (j k)
  278.       (let ((punkt (elt simplex j)))
  279.         (assert (= (length punkt) d))
  280.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  281.         (setf (aref A d j) 1)
  282.     ) )
  283.     (let ((j k))
  284.       (assert (= (length punkt1) d))
  285.       (dotimes (i d) (setf (aref A i j) (aref punkt1 i)))
  286.       (setf (aref A d j) 1)
  287.     )
  288.     (let ((j (+ k 1)))
  289.       (assert (= (length punkt2) d))
  290.       (dotimes (i d) (setf (aref A i j) (aref punkt2 i)))
  291.       (setf (aref A d j) 1)
  292.     )
  293.     ; LGS lösen:
  294.     (let ((kern (solve-homLGS A)))
  295.       (ecase (length kern)
  296.         (0 (error "Simplex und zwei Punkte sind windschief."))
  297.         (1 (let ((lsg (first kern))) ; Lösungvektor #(x1 ... xk y1 y2)
  298.              ; y1*y2>0 -> eine echte Konvexkombination von punkt1 und punkt2
  299.              ;            liegt in der Ebene, also liegen die Punkte in
  300.              ;            verschiedenen Halbebenen der Ebene.
  301.              (<= (* (signum (aref lsg k)) (signum (aref lsg (+ k 1)))) 0)
  302.         )  )
  303.         (2 t) ; beide Punkte liegen in der affinen Hülle des Simplex
  304. ) ) ) )
  305.  
  306. ; Liefert den Mittelpunkt eines Simplex:
  307. (defun simplex-mittelpunkt (simplex)
  308.   (let ((k (length simplex))) ; Anzahl der Punkte
  309.     (map 'vector #'(lambda (x) (/ x k)) ; durch sie dividieren, nachdem man
  310.       (apply #'map 'vector #'+ simplex) ; alle Koordinaten addiert hat
  311. ) ) )
  312.  
  313. ; Liefert das Volumen eines Simplex der Dimension d:
  314. (defun simplex-vol (simplex)
  315.   (let* ((d (length (first simplex)))
  316.          (d+1 (+ d 1)))
  317.     (if (= (length simplex) d+1)
  318.       ; volldimensionaler Simplex conv({p0,...,pd}) hat Volumen
  319.       ;  1          | p00 ... p0d 1 |
  320.       ; --- abs det | ............. |
  321.       ; d !         | pd0 ... pdd 1 |
  322.       (let ((A (make-array (list d+1 d+1))))
  323.         (dotimes (i d+1)
  324.           (let ((punkt (elt simplex i)))
  325.             (dotimes (j d) (setf (aref A i j) (aref punkt j)))
  326.             (setf (aref A i d) 1)
  327.         ) )
  328.         (/ (abs (det A)) (! d))
  329.       )
  330.       ; niederdimensionaler Simplex hat d-Volumen 0
  331.       0
  332. ) ) )
  333.  
  334. ; Stellt fest, ob ein Punkt in einem Simplex liegt.
  335. (defun on-simplex-p (punkt simplex)
  336.   ; Löse das LGS x1*p1 + ... + xk*pk + x*p = 0, x1+...+xk+x=0
  337.   (let* ((d (length punkt))
  338.          (k (length simplex))
  339.          (A (make-array (list (+ d 1) (+ k 1)))))
  340.     ; LGS aufstellen:
  341.     (dotimes (j k)
  342.       (let ((punkt (elt simplex j)))
  343.         (assert (= (length punkt) d))
  344.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  345.         (setf (aref A d j) 1)
  346.     ) )
  347.     (dotimes (i d) (setf (aref A i k) (aref punkt i)))
  348.     (setf (aref A d k) 1)
  349.     ; LGS lösen:
  350.     (let ((kern (solve-homLGS A)))
  351.       (assert (<= (length kern) 1)) ; LGS durfte nicht mehrdeutig lösbar sein
  352.       (if (endp kern) ; falls es unlösbar war, liegt der Punkt nicht drin.
  353.         nil
  354.         (let ((lsg (first kern)))
  355.           ; x auf -1 normieren:
  356.           (let ((-x (- (aref lsg k))))
  357.             (dotimes (j k) (setf (aref lsg j) (/ (aref lsg j) -x)))
  358.           )
  359.           ; Nun ist x1*p1 + ... + xk*pk = p .
  360.           ; Falls hierbei alle xi>=0 sind, liegt p im Simplex.
  361.           (dotimes (j k t)
  362.             (unless (>= (aref lsg j) 0) (return nil))
  363.           )
  364. ) ) ) ) )
  365.  
  366. ; Stellt fest, ob ein Punkt im Innern eines Simplex liegt.
  367. (defun in-simplex-p (punkt simplex)
  368.   ; Löse das LGS x1*p1 + ... + xk*pk + x*p = 0, x1+...+xk+x=0
  369.   (let* ((d (length punkt))
  370.          (k (length simplex))
  371.          (A (make-array (list (+ d 1) (+ k 1)))))
  372.     ; LGS aufstellen:
  373.     (dotimes (j k)
  374.       (let ((punkt (elt simplex j)))
  375.         (assert (= (length punkt) d))
  376.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  377.         (setf (aref A d j) 1)
  378.     ) )
  379.     (dotimes (i d) (setf (aref A i k) (aref punkt i)))
  380.     (setf (aref A d k) 1)
  381.     ; LGS lösen:
  382.     (let ((kern (solve-homLGS A)))
  383.       (assert (<= (length kern) 1)) ; LGS durfte nicht mehrdeutig lösbar sein
  384.       (if (endp kern) ; falls es unlösbar war, liegt der Punkt nicht drin.
  385.         nil
  386.         (let ((lsg (first kern)))
  387.           ; x auf -1 normieren:
  388.           (let ((-x (- (aref lsg k))))
  389.             (dotimes (j k) (setf (aref lsg j) (/ (aref lsg j) -x)))
  390.           )
  391.           ; Nun ist x1*p1 + ... + xk*pk = p .
  392.           ; Falls hierbei alle xi>0 sind, liegt p im Innern des Simplex.
  393.           (dotimes (j k t)
  394.             (unless (> (aref lsg j) 0) (return nil))
  395.           )
  396. ) ) ) ) )
  397.  
  398. ;; Eine Triangulation eines Polytops P der Dimension k ist eine Liste von
  399. ;; Simplizes der Dimension k, deren Vereinigung = P und deren paarweiser
  400. ;; Durchschnitt jeweils höchstens (k-1)-dimensional ist.
  401. ;; Ein allgemeines Polytop wird als Eckenliste und als Triangulation
  402. ;; dargestellt.
  403.  
  404. (defstruct (polytop (:type list))
  405.   (ecken nil :type list)
  406.   (triang nil :type list)
  407. )
  408.  
  409. ; Dimension eines Polytops:
  410. (defun polytop-dim (polytop)
  411.   (let ((triang (polytop-triang polytop)))
  412.     (if (endp triang)
  413.       -1 ; leeres Polytop hat Dimension -1
  414.       (simplex-dim (first triang)) ; sonst die Dimension aller Simplizes
  415. ) ) )
  416.  
  417. ; Simplex als Polytop:
  418. (defun simplex-polytop (simplex)
  419.   (make-polytop :ecken simplex
  420.                 :triang (list simplex) ; Triangulation besteht als dem Simplex selbst
  421. ) )
  422.  
  423. ; Stellt fest, ob ein Punkt in der konvexen Hülle einer Punktmenge liegt.
  424. (defun on-conv-p (punkt punkte)
  425.   ; Läßt sich der Punkt aus den "Ecken" konvex kombinieren?
  426.   ; Simplex-Algorithmus anwenden:
  427.   (let* ((d (length punkt))
  428.          (m (1+ d))
  429.          (n (length punkte))
  430.          (A (make-array `(,m ,n)))
  431.          (b (make-array m))
  432.          (c (make-array n :initial-element 0)))
  433.     (dotimes (j n)
  434.       (let ((punkt (pop punkte)))
  435.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  436.         (setf (aref A d j) 1)
  437.     ) )
  438.     (dotimes (i d) (setf (aref b i) (aref punkt i)))
  439.     (setf (aref b d) 1)
  440.     (values (simplex A b c))
  441. ) )
  442.  
  443. ; Stellt fest, ob ein Punkt in einem Polytop liegt.
  444. (defun on-polytop-p (punkt polytop)
  445.   ; Läßt sich der Punkt aus den Ecken konvex kombinieren?
  446.   (on-conv-p punkt (polytop-ecken polytop))
  447. )
  448.  
  449. ; Liefert zu einem Polytop der Dimension k eine Liste aller Randsimplizes,
  450. ; alle derselben Dimension k-1.
  451. (defun simplex-rand (simplex)
  452.   (if (zerop (simplex-dim simplex)) ; Hat Dimension 0 ?
  453.     nil ; Rand ist leer
  454.     ; P = conv{p0,...,pk} mit k>0 -> Die Randsimplizes sind die Simplizes
  455.     ; conv{p0,...,pj-1,pj+1,...,pk} (pj fehlt) für j=0,...,k.
  456.     (let ((rand nil))
  457.       (do ((L1 nil (cons (car L2) L1))
  458.            (L2 simplex (cdr L2)))
  459.           ((endp L2))
  460.         (push (revappend L1 (cdr L2)) rand)
  461.       )
  462.       (nreverse rand)
  463. ) ) )
  464.  
  465. ; Liefert zu einer Triangulation eines Polytops eine Liste der Randsimplizes.
  466. ; Sie haben alle eine um 1 geringere Dimension als das Polytop selbst.
  467. (defun triang-rand (triang)
  468.   ; Der Rand setzt sich zusammen aus den Randsimplizes der Triangulierung.
  469.   ; Solche Randsimplizes, die dabei doppelt vorkommen, liegen im Innern des
  470.   ; Polytops und werden daher gestrichen.
  471.   (let ((rand nil))
  472.     (dolist (simplex triang)
  473.       (let ((rand-neu (simplex-rand simplex)))
  474.         (dolist (rand-simplex rand-neu)
  475.           (if (member rand-simplex rand :test #'simplex=)
  476.             (setq rand (delete rand-simplex rand :test #'simplex=))
  477.             (setq rand (cons rand-simplex rand))
  478.     ) ) ) )
  479.     (nreverse rand)
  480. ) )
  481.  
  482. ; Stellt fest, ob ein Punkt im Innern eines Polytops liegt.
  483. (defun in-polytop-p (punkt polytop)
  484.   ; Liegt der Punkt im Polytop, aber nicht auf dem Rand?
  485.   (and (on-polytop-p punkt polytop)
  486.        (not (some #'(lambda (randsimplex) (on-simplex-p punkt randsimplex))
  487.                   (triang-rand (polytop-triang polytop))
  488. ) )    )    )
  489.  
  490. ; Liefert die Liste der Ecken der konvexen Hülle von endlich vielen Punkten.
  491. (defun conv-ecken (punkte)
  492.   (setq punkte (remove-duplicates punkte :test #'punkt=))
  493.   (remove-if
  494.     #'(lambda (punkt)
  495.         (on-conv-p punkt (remove punkt punkte :test #'punkt=))
  496.       )
  497.     punkte
  498. ) )
  499.  
  500. ; Liefert die konvexe Hülle einer Liste von Punkten als Polytop.
  501. (defun conv (punkte)
  502.   (setq punkte (conv-ecken punkte)) ; unnötige Punkte streichen
  503.   (make-polytop
  504.     :ecken punkte
  505.     :triang
  506.       ; Triangulation durch schrittweise Erweiterung um je einen der Punkte
  507.       ; gewinnen:
  508.       (let ((triang nil)) ; leere Triangulation
  509.         (dolist (punkt punkte)
  510.           (setq triang
  511.             ; punkt zu triang dazunehmen:
  512.             (if (endp triang) ; Polytop leer?
  513.               (list (make-simplex (list punkt))) ; gibt einpunktiges Polytop
  514.               ; liegt der Punkt in der affinen Hülle des Polytops (= der affinen
  515.               ; Hülle eines beliebigen Simplex der Triangulation) ?
  516.               (if (on-affine-p punkt (first triang))
  517.                 ; Dimension vergrößert sich nicht: Bilde die Randsimplizes des Polytops,
  518.                 ; die den Punkt vom Polytop trennen, und verbinde sie jeweils mit dem
  519.                 ; Punkt. Diese neuen Simplizes kommen zur Triangulation dazu.
  520.                 (let ((einpunkt (simplex-mittelpunkt (first triang))) ; Ein Punkt im Polytop
  521.                       (neue-simplizes nil))
  522.                   (dolist (rand (triang-rand triang))
  523.                     (unless (on-same-side-p punkt einpunkt rand)
  524.                       (push (make-simplex (cons punkt rand)) neue-simplizes)
  525.                   ) )
  526.                   (append triang (nreverse neue-simplizes))
  527.                 )
  528.                 ; Dimension vergrößert sich um 1: Verbinde den Punkt mit allen
  529.                 ; Simplizes des Polytops. Dies liefert die neue Triangulation.
  530.                 (mapcar #'(lambda (simplex) (make-simplex (cons punkt simplex))) triang)
  531.         ) ) ) )
  532.         triang
  533. ) )   )
  534.  
  535. ; Liefert das Volumen (der Dimension d) eines Polytops.
  536. (defun polytop-vol (polytop)
  537.   (apply #'+ (mapcar #'simplex-vol (polytop-triang polytop)))
  538. )
  539.  
  540.  
  541. ;; Beispiel: Hn(r) bei festem n.
  542. ;; Ecken sind die n! Permutationsmatrizen, aufgefaßt im R^(n^2), durch
  543. ;; Weglassen der rechten und unteren Spalte projiziert in den R^((n-1)^2).
  544.  
  545. ; Permutation von {0,...,n-1} in Punkt umwandeln:
  546. (defun perm-als-punkt (perm)
  547.   (let* ((n (length perm))
  548.          (n-1 (- n 1))
  549.          (A (make-array `(,n-1 ,n-1) :initial-element 0))
  550.          (p (make-array (* n-1 n-1) :displaced-to A)))
  551.     (dotimes (i n)
  552.       (let ((j (elt perm i)))
  553.         (when (array-in-bounds-p A i j)
  554.           (setf (aref A i j) 1)
  555.     ) ) )
  556.     p
  557. ) )
  558.  
  559. ; Liste aller Permutationen von {0,...,n-1} produzieren:
  560. (defun alle-perm (n)
  561.   (if (zerop n)
  562.     (list '())
  563.     (let ((perm_n-1 (alle-perm (- n 1)))
  564.           (perm_n nil))
  565.       ; produziere die Permutationen von {0,...,n-1} aus denen von {0,...,n-2}:
  566.       ; (n-1)-Liste unverändert. Vorne n-1 anfügen.
  567.       ; In der (n-1)-Liste alle n-2 in n-1 umwandeln. Vorne n-2 anfügen.
  568.       ; In der (n-1)-Liste alle n-3 in n-2 umwandeln. Vorne n-3 anfügen.
  569.       ; ...
  570.       ; In der (n-1)-Liste alle 0 in 1 umwandeln. Vorne 0 anfügen.
  571.       (let ((i (- n 1)))
  572.         (loop
  573.           (setq perm_n
  574.             (append (mapcar #'(lambda (p_n-1) (cons i p_n-1)) perm_n-1)
  575.                     perm_n
  576.           ) )
  577.           (when (zerop i) (return))
  578.           (decf i)
  579.           (setq perm_n-1 (subst (+ i 1) i perm_n-1))
  580.       ) )
  581.       perm_n
  582. ) ) )
  583.  
  584. (defun H-ecken (n) (mapcar #'perm-als-punkt (alle-perm n)))
  585. ; n=1
  586. (setq H1 (conv (H-ecken 1)))
  587. ; n=2
  588. (setq H2 (conv (H-ecken 2)))
  589. ; n=3
  590. (setq H3 (conv (H-ecken 3)))
  591.  
  592.