home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / clisp / polytop1.lsp < prev    next >
Lisp/Scheme  |  1993-10-23  |  20KB  |  534 lines
  1. ;; Berechnung des Volumens konvexer Polytope
  2. ;; Bruno Haible 5.9.1990
  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. ;; Wir bewegen uns im R^d, d>=0. Ein Punkt ist ein Vektor aus seinen d
  173. ;; Koordinaten (rationale Zahlen).
  174.  
  175. ; Gleichheit von Punkten:
  176. (defun punkt= (p1 p2)
  177.   (equalp p1 p2) ; Vergleich aller Koordinaten
  178. )
  179.  
  180. ;; Ein Simplex ist eine (ungeordnete) Liste von k Punkten, k>=1, deren
  181. ;; konvexe Hülle die Dimension k-1 hat.
  182.  
  183. ; Dimension eines Simplex:
  184. (defun simplex-dim (simplex)
  185.   (1- (length simplex))
  186. )
  187.  
  188. ; Gleichheit von Simplizes:
  189. (defun simplex= (simplex1 simplex2)
  190.   (and (= (simplex-dim simplex1) (simplex-dim simplex2)) ; müssen gleiche Dimension haben
  191.        (dolist (punkt1 simplex1 t) ; und jede Ecke von simplex1
  192.          (unless (position punkt1 simplex2 :test #'punkt=) ; muß auch Ecke von simplex2 sein
  193.            (return nil)
  194.        ) )
  195. ) )
  196.  
  197. ; Bildet ein (k-1)-dimensionales Simplex aus k Punkten:
  198. (defun make-simplex (punkte)
  199.   (if (endp punkte)
  200.     (error "Zu wenig Punkte für ein Simplex.")
  201.     (if ; Löse das LGS x1*p1 + ... + xk*pk = 0, x1+...+xk=0
  202.       (let* ((d (length (first punkte)))
  203.              (k (length punkte))
  204.              (A (make-array (list (+ d 1) k))))
  205.         ; LGS aufstellen:
  206.         (dotimes (j k)
  207.           (let ((punkt (elt punkte j)))
  208.             (assert (= (length punkt) d))
  209.             (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  210.             (setf (aref A d j) 1)
  211.         ) )
  212.         ; LGS lösen. Der Kern muß leer sein, also (p1,1), ..., (pk,1) linear
  213.         ; unabhängig, also p1,...,pk affin unabhängig.
  214.         (null (solve-homLGS A))
  215.       )
  216.       punkte
  217.       (error "Punkte für ein Simplex sind affin abhängig: ~S" punkte)
  218. ) ) )
  219.  
  220. ; Stellt fest, ob ein Punkt in der affinen Hülle eines Simplex liegt.
  221. (defun on-affine-p (punkt simplex)
  222.   ; Löse das LGS x1*p1 + ... + xk*pk + x*p = 0, x1+...+xk+x=0
  223.   (let* ((d (length punkt))
  224.          (k (length simplex))
  225.          (A (make-array (list (+ d 1) (+ k 1)))))
  226.     ; LGS aufstellen:
  227.     (dotimes (j k)
  228.       (let ((punkt (elt simplex j)))
  229.         (assert (= (length punkt) d))
  230.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  231.         (setf (aref A d j) 1)
  232.     ) )
  233.     (dotimes (i d) (setf (aref A i k) (aref punkt i)))
  234.     (setf (aref A d k) 1)
  235.     ; LGS lösen:
  236.     (let ((kern (solve-homLGS A)))
  237.       (assert (<= (length kern) 1)) ; LGS durfte nicht mehrdeutig lösbar sein
  238.       (not (endp kern)) ; falls es unlösbar war, liegt der Punkt nicht drin.
  239. ) ) )
  240.  
  241. ; Stellt fest, ob zwei Punkte auf derselben Seite der affinen Hülle eines
  242. ; Simplex liegen.
  243. (defun on-same-side-p (punkt1 punkt2 simplex)
  244.   ; Löse das LGS x1*p1 + ... + xk*pk + y1*punkt1 + y2*punkt2 = 0, x1+...+xk+y1+y2=0
  245.   (let* ((d (length punkt1))
  246.          (k (length simplex))
  247.          (A (make-array (list (+ d 1) (+ k 2)))))
  248.     ; LGS aufstellen:
  249.     (dotimes (j k)
  250.       (let ((punkt (elt simplex j)))
  251.         (assert (= (length punkt) d))
  252.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  253.         (setf (aref A d j) 1)
  254.     ) )
  255.     (let ((j k))
  256.       (assert (= (length punkt1) d))
  257.       (dotimes (i d) (setf (aref A i j) (aref punkt1 i)))
  258.       (setf (aref A d j) 1)
  259.     )
  260.     (let ((j (+ k 1)))
  261.       (assert (= (length punkt2) d))
  262.       (dotimes (i d) (setf (aref A i j) (aref punkt2 i)))
  263.       (setf (aref A d j) 1)
  264.     )
  265.     ; LGS lösen:
  266.     (let ((kern (solve-homLGS A)))
  267.       (ecase (length kern)
  268.         (0 (error "Simplex und zwei Punkte sind windschief."))
  269.         (1 (let ((lsg (first kern))) ; Lösungvektor #(x1 ... xk y1 y2)
  270.              ; y1*y2>0 -> eine echte Konvexkombination von punkt1 und punkt2
  271.              ;            liegt in der Ebene, also liegen die Punkte in
  272.              ;            verschiedenen Halbebenen der Ebene.
  273.              (<= (* (signum (aref lsg k)) (signum (aref lsg (+ k 1)))) 0)
  274.         )  )
  275.         (2 t) ; beide Punkte liegen in der affinen Hülle des Simplex
  276. ) ) ) )
  277.  
  278. ; Liefert den Mittelpunkt eines Simplex:
  279. (defun simplex-mittelpunkt (simplex)
  280.   (let ((k (length simplex))) ; Anzahl der Punkte
  281.     (map 'vector #'(lambda (x) (/ x k)) ; durch sie dividieren, nachdem man
  282.       (apply #'map 'vector #'+ simplex) ; alle Koordinaten addiert hat
  283. ) ) )
  284.  
  285. ; Liefert das Volumen eines Simplex der Dimension d:
  286. (defun simplex-vol (simplex)
  287.   (let* ((d (length (first simplex)))
  288.          (d+1 (+ d 1)))
  289.     (if (= (length simplex) d+1)
  290.       ; volldimensionaler Simplex conv({p0,...,pd}) hat Volumen
  291.       ;  1          | p00 ... p0d 1 |
  292.       ; --- abs det | ............. |
  293.       ; d !         | pd0 ... pdd 1 |
  294.       (let ((A (make-array (list d+1 d+1))))
  295.         (dotimes (i d+1)
  296.           (let ((punkt (elt simplex i)))
  297.             (dotimes (j d) (setf (aref A i j) (aref punkt j)))
  298.             (setf (aref A i d) 1)
  299.         ) )
  300.         (/ (abs (det A)) (! d))
  301.       )
  302.       ; niederdimensionaler Simplex hat d-Volumen 0
  303.       0
  304. ) ) )
  305.  
  306. ; Stellt fest, ob ein Punkt in einem Simplex liegt.
  307. (defun on-simplex-p (punkt simplex)
  308.   ; Löse das LGS x1*p1 + ... + xk*pk + x*p = 0, x1+...+xk+x=0
  309.   (let* ((d (length punkt))
  310.          (k (length simplex))
  311.          (A (make-array (list (+ d 1) (+ k 1)))))
  312.     ; LGS aufstellen:
  313.     (dotimes (j k)
  314.       (let ((punkt (elt simplex j)))
  315.         (assert (= (length punkt) d))
  316.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  317.         (setf (aref A d j) 1)
  318.     ) )
  319.     (dotimes (i d) (setf (aref A i k) (aref punkt i)))
  320.     (setf (aref A d k) 1)
  321.     ; LGS lösen:
  322.     (let ((kern (solve-homLGS A)))
  323.       (assert (<= (length kern) 1)) ; LGS durfte nicht mehrdeutig lösbar sein
  324.       (if (endp kern) ; falls es unlösbar war, liegt der Punkt nicht drin.
  325.         nil
  326.         (let ((lsg (first kern)))
  327.           ; x auf -1 normieren:
  328.           (let ((-x (- (aref lsg k))))
  329.             (dotimes (j k) (setf (aref lsg j) (/ (aref lsg j) -x)))
  330.           )
  331.           ; Nun ist x1*p1 + ... + xk*pk = p .
  332.           ; Falls hierbei alle xi>=0 sind, liegt p im Simplex.
  333.           (dotimes (j k t)
  334.             (unless (>= (aref lsg j) 0) (return nil))
  335.           )
  336. ) ) ) ) )
  337.  
  338. ; Stellt fest, ob ein Punkt im Innern eines Simplex liegt.
  339. (defun in-simplex-p (punkt simplex)
  340.   ; Löse das LGS x1*p1 + ... + xk*pk + x*p = 0, x1+...+xk+x=0
  341.   (let* ((d (length punkt))
  342.          (k (length simplex))
  343.          (A (make-array (list (+ d 1) (+ k 1)))))
  344.     ; LGS aufstellen:
  345.     (dotimes (j k)
  346.       (let ((punkt (elt simplex j)))
  347.         (assert (= (length punkt) d))
  348.         (dotimes (i d) (setf (aref A i j) (aref punkt i)))
  349.         (setf (aref A d j) 1)
  350.     ) )
  351.     (dotimes (i d) (setf (aref A i k) (aref punkt i)))
  352.     (setf (aref A d k) 1)
  353.     ; LGS lösen:
  354.     (let ((kern (solve-homLGS A)))
  355.       (assert (<= (length kern) 1)) ; LGS durfte nicht mehrdeutig lösbar sein
  356.       (if (endp kern) ; falls es unlösbar war, liegt der Punkt nicht drin.
  357.         nil
  358.         (let ((lsg (first kern)))
  359.           ; x auf -1 normieren:
  360.           (let ((-x (- (aref lsg k))))
  361.             (dotimes (j k) (setf (aref lsg j) (/ (aref lsg j) -x)))
  362.           )
  363.           ; Nun ist x1*p1 + ... + xk*pk = p .
  364.           ; Falls hierbei alle xi>0 sind, liegt p im Innern des Simplex.
  365.           (dotimes (j k t)
  366.             (unless (> (aref lsg j) 0) (return nil))
  367.           )
  368. ) ) ) ) )
  369.  
  370. ;; Eine Triangulation eines Polytops P der Dimension k ist eine Liste von
  371. ;; Simplizes der Dimension k, deren Vereinigung = P und deren paarweiser
  372. ;; Durchschnitt jeweils höchstens (k-1)-dimensional ist.
  373. ;; Ein allgemeines Polytop wird als Triangulation dargestellt.
  374.  
  375. ; Dimension eines Polytops:
  376. (defun polytop-dim (polytop)
  377.   (if (endp polytop)
  378.     -1 ; leeres Polytop hat Dimension 0
  379.     (simplex-dim (first polytop)) ; sonst die Dimension aller Simplizes
  380. ) )
  381.  
  382. ; Simplex als Polytop:
  383. (defun simplex-polytop (simplex)
  384.   (list simplex) ; Triangulation besteht als dem Simplex selbst
  385. )
  386.  
  387. ; Stellt fest, ob ein Punkt in einem Polytop liegt.
  388. (defun on-polytop-p (punkt polytop)
  389.   ; Liegt der Punkt in einer der Simplizes der Triangulation?
  390.   (dolist (simplex polytop nil)
  391.     (when (on-simplex-p punkt simplex) (return t))
  392. ) )
  393.  
  394. ; Liefert zu einem Polytop der Dimension k eine Liste aller Randsimplizes,
  395. ; alle derselben Dimension k-1.
  396. (defun simplex-rand (simplex)
  397.   (if (zerop (simplex-dim simplex)) ; Hat Dimension 0 ?
  398.     nil ; Rand ist leer
  399.     ; P = conv{p0,...,pk} mit k>0 -> Die Randsimplizes sind die Simplizes
  400.     ; conv{p0,...,pj-1,pj+1,...,pk} (pj fehlt) für j=0,...,k.
  401.     (let ((rand nil))
  402.       (do ((L1 nil (cons (car L2) L1))
  403.            (L2 simplex (cdr L2)))
  404.           ((endp L2))
  405.         (push (revappend L1 (cdr L2)) rand)
  406.       )
  407.       (nreverse rand)
  408. ) ) )
  409.  
  410. ; Liefert zu einem Polytop eine Liste aller Randsimplizes. Sie haben alle
  411. ; eine um 1 geringere Dimension als das Polytop selbst.
  412. (defun polytop-rand (polytop)
  413.   ; Der Rand setzt sich zusammen aus den Randsimplizes der Triangulierung.
  414.   ; Solche Randsimplizes, die dabei doppelt vorkommen, liegen im Innern des
  415.   ; Polytops und werden daher gestrichen.
  416.   (let ((rand nil))
  417.     (dolist (simplex polytop)
  418.       (let ((rand-neu (simplex-rand simplex)))
  419.         (dolist (rand-simplex rand-neu)
  420.           (if (member rand-simplex rand :test #'simplex=)
  421.             (setq rand (delete rand-simplex rand :test #'simplex=))
  422.             (setq rand (cons rand-simplex rand))
  423.     ) ) ) )
  424.     (nreverse rand)
  425. ) )
  426.  
  427. ; Stellt fest, ob ein Punkt im Innern eines Polytops liegt.
  428. (defun in-polytop-p (punkt polytop)
  429.   ; Liegt der Punkt im Polytop, aber nicht auf dem Rand?
  430.   (and (on-polytop-p punkt polytop)
  431.        (not (some #'(lambda (randsimplex) (on-simplex-p punkt randsimplex))
  432.                   (polytop-rand polytop)
  433. ) )    )    )
  434.  
  435. ; Liefert die konvexe Hülle eines Polytops und eines Punktes.
  436. (defun polytop+punkt (polytop punkt)
  437.   (if (endp polytop) ; Polytop leer?
  438.     (simplex-polytop (make-simplex (list punkt))) ; gibt einpunktiges Polytop
  439.     ; liegt der Punkt in der affinen Hülle des Polytops (= der affinen
  440.     ; Hülle eines beliebigen Simplex der Triangulation) ?
  441.     (if (on-affine-p punkt (first polytop))
  442.       ; Dimension vergrößert sich nicht: Bilde die Randsimplizes des Polytops,
  443.       ; die den Punkt vom Polytop trennen, und verbinde sie jeweils mit dem
  444.       ; Punkt. Diese neuen Simplizes kommen zur Triangulation dazu.
  445.       (let ((einpunkt (simplex-mittelpunkt (first polytop))) ; Ein Punkt im Polytop
  446.             (neue-simplizes nil))
  447.         (dolist (rand (polytop-rand polytop))
  448.           (unless (on-same-side-p punkt einpunkt rand)
  449.             (push (make-simplex (cons punkt rand)) neue-simplizes)
  450.         ) )
  451.         (append polytop (nreverse neue-simplizes))
  452.       )
  453.       ; Dimension vergrößert sich um 1: Verbinde den Punkt mit allen
  454.       ; Simplizes des Polytops. Dies liefert die neue Triangulation.
  455.       (mapcar #'(lambda (simplex) (make-simplex (cons punkt simplex))) polytop)
  456. ) ) )
  457.  
  458. ; Liefert die konvexe Hülle einer Liste von Punkten als Polytop.
  459. (defun conv (punkte)
  460.   (let ((polytop nil)) ; leeres Polytop
  461.     (dolist (punkt punkte) ; schrittweise um je einen der Punkte erweitern
  462.       (setq polytop (polytop+punkt polytop punkt))
  463.     )
  464.     polytop
  465. ) )
  466.  
  467. ; Liefert das Volumen (der Dimension d) eines Polytops.
  468. (defun polytop-vol (polytop)
  469.   (apply #'+ (mapcar #'simplex-vol polytop))
  470. )
  471.  
  472. ; Liefert die Liste der Ecken der konvexen Hülle von endlich vielen Punkten.
  473. (defun conv-ecken (punkte)
  474.   (setq punkte (remove-duplicates punkte :test #'punkt=))
  475.   (remove-if
  476.     #'(lambda (punkt)
  477.         (on-polytop-p punkt (conv (remove punkt punkte :test #'punkt=)))
  478.       )
  479.     punkte
  480. ) )
  481.  
  482.  
  483. ;; Beispiel: Hn(r) bei festem n.
  484. ;; Ecken sind die n! Permutationsmatrizen, aufgefaßt im R^(n^2), durch
  485. ;; Weglassen der rechten und unteren Spalte projiziert in den R^((n-1)^2).
  486.  
  487. ; Permutation von {0,...,n-1} in Punkt umwandeln:
  488. (defun perm-als-punkt (perm)
  489.   (let* ((n (length perm))
  490.          (n-1 (- n 1))
  491.          (A (make-array `(,n-1 ,n-1) :initial-element 0))
  492.          (p (make-array (* n-1 n-1) :displaced-to A)))
  493.     (dotimes (i n)
  494.       (let ((j (elt perm i)))
  495.         (when (array-in-bounds-p A i j)
  496.           (setf (aref A i j) 1)
  497.     ) ) )
  498.     p
  499. ) )
  500.  
  501. ; Liste aller Permutationen von {0,...,n-1} produzieren:
  502. (defun alle-perm (n)
  503.   (if (zerop n)
  504.     (list '())
  505.     (let ((perm_n-1 (alle-perm (- n 1)))
  506.           (perm_n nil))
  507.       ; produziere die Permutationen von {0,...,n-1} aus denen von {0,...,n-2}:
  508.       ; (n-1)-Liste unverändert. Vorne n-1 anfügen.
  509.       ; In der (n-1)-Liste alle n-2 in n-1 umwandeln. Vorne n-2 anfügen.
  510.       ; In der (n-1)-Liste alle n-3 in n-2 umwandeln. Vorne n-3 anfügen.
  511.       ; ...
  512.       ; In der (n-1)-Liste alle 0 in 1 umwandeln. Vorne 0 anfügen.
  513.       (let ((i (- n 1)))
  514.         (loop
  515.           (setq perm_n
  516.             (append (mapcar #'(lambda (p_n-1) (cons i p_n-1)) perm_n-1)
  517.                     perm_n
  518.           ) )
  519.           (when (zerop i) (return))
  520.           (decf i)
  521.           (setq perm_n-1 (subst (+ i 1) i perm_n-1))
  522.       ) )
  523.       perm_n
  524. ) ) )
  525.  
  526. (defun H-ecken (n) (mapcar #'perm-als-punkt (alle-perm n)))
  527. ; n=1
  528. (setq H1 (conv (H-ecken 1)))
  529. ; n=2
  530. (setq H2 (conv (H-ecken 2)))
  531. ; n=3
  532. (setq H3 (conv (H-ecken 3)))
  533.  
  534.