share/sym/schur.lisp

   1 ; schur.lsp
   2
   3 ;       ***************************************************************
   4 ;       *                    MODULE SYM                               *
   5 ;       *       MANIPULATIONS DE FONCTIONS SYMETRIQUES                *
   6 ;       *        (version01: Commonlisp pour Maxima)                 *
   7 ;       *                                                             *
   8 ;       *                ----------------------                       *
   9 ;       *          Philippe ESPERET Annick VALIBOUZE                  *
  10 ;       *                    GDR MEDICIS                              *
  11 ;       *  (Mathe'matiques Effectives, De'veloppements Informatiques, *
  12 ;       *           Calculs et Ingenierie, Syste`mes)                 *
  13 ;       *             LITP (Equipe Calcul Formel)                     *
  14 ;       *                 Universite' Paris 6,                        *
  15 ;       *        4 place Jussieu, 75252 Paris cedex 05.               *
  16 ;       *              e-mail : avb@sysal.ibp.fr                      *
  17 ;       ***************************************************************
  18
  19 ;----------------------------------------------------------------------------
  20 ;                 FONCTIONS DE SCHUR et CHANGEMENTS DE BASES
  21 ;----------------------------------------------------------------------------
  22 ; DIFFERENCE ENTRE COMPILE ET INTERPRETE : voir fonction schur2comp_pol
  23 ;----------------------------------------------------------------------------
  24 ;----------------------------------------------------------------------------
  25 ;            PASSAGE DES FORMES MONOMIALES AUX FONCTIONS DE SCHUR
  26
  27 (in-package :maxima)
  28 (macsyma-module schur)
  29
  30 (mdefprop $mon2schur
  31     ((lambda ()) ((mlist) $part)
  32      ((mprog) (($operation)) (($mon2schur_init) $part)))
  33     mexpr)
  34 (add2lnc '(($mon2schur) $part) $functions)
  35 (mdefprop $kostka
  36     ((lambda ()) ((mlist) $part1 $part2)
  37      ((mprog) (($operation)) (($kostka_init) $part1 $part2)))
  38     mexpr)
  39 (add2lnc '(($kostka) $part1 $part2) $functions)
  40 (mdefprop $treinat
  41     ((lambda ()) ((mlist) $part)
  42      ((mprog) (($operation)) (($treinat_init) $part)))
  43     mexpr)
  44 (add2lnc '(($treinat) $part) $functions)
  45 ;            PASSAGE DES FONCTIONS DE  SCHUR AUX COMPLETES
  46 (mdefprop $schur2comp
  47     ((lambda ()) ((mlist) $comp $listofvars)
  48      ((mprog) (($operation)) (($schur2comp_init) $comp $listofvars)))
  49     mexpr)
  50 (add2lnc '(($schur2comp) $comp $listofvars) $functions)
  51
  52 ;========================================================================
  53 ;                 PASSAGE DES FONCTIONS DE SCHUR
  54 ;                     AUX FONCTIONS COMPLETES
  55 ; On se donne un polyno^me en h1, h2, ... repre'sentant un polyno^me
  56 ; en les fonctions comple`tes.
  57 ; on recupere une liste dont chaque element est une liste dont
  58 ; le premier terme est un entier et le reste une partition renversee
  59 ; representant la fonction de Schur associ\'ee.
  60 ; REP(partition) = [partition](1)
  61 ;========================================================================
  62 ; l'entree est un polynome en les hi
  63 ; l'entree est une liste que l'on ordonne
  64
  65 (defun $schur2comp_init ($comp $listofvars)
  66   (cond
  67     ((eql '$pol $schur2comp) (schur2comp_pol $comp (cdr $listofvars)))
  68     (t (cons '(mlist)
  69              (schur2comp
  70                  (ch2rep (sort (cdr $comp) '>)))))))
  71 ;.........................................................................
  72 ; Si la donnee est un polynome en les fonctions completes.
  73 ; schur2comp rend '(mlist) en tete de chaque terme partitionne
  74 ; le coefficient est donc en cadr
  75
  76 (defun schur2comp_pol ($pol listofvars)
  77   (do ((polpart (pol2part $pol listofvars) (cdr polpart)) (sol 0))
  78       ((null polpart) sol)
  79     (setq sol
  80           ($add_sym
  81               ($fadd_sym
  82                         (mapcar #'(lambda (l)
  83                                      (let ((coef (caar polpart)))
  84                                            ($fmult_sym
  85                                                 (list (cadr l) coef
  86                                                 (cons '($S array)
  87                                                       (cddr l))))))
  88                                  (schur2comp (cdar polpart))))
  89               sol))))
  90
  91 ; Ordre Lexicographique pour des polynomes partitionnes de type 1.
  92 ; l'egalite ne nous importe pas.
  93
  94 (defun lexinv_type1 (terme1 terme2)
  95   (2lexinv_type1 (cddr terme1) (cddr terme2)))
  96
  97 (defun 2lexinv_type1 (1part 2part)
  98   (cond
  99     ((null (car 1part)) nil)
 100     ((null (car 2part)) t)
 101     ((< (car 1part) (car 2part))
 102      nil)
 103     ((> (car 1part) (car 2part))
 104      t)
 105     (t (2lexinv_type1 (cdr 1part) (cdr 2part)))))
 106 ;........................................................................
 107 ; pol2part permet a` partir d'un polynome en les monomes
 108 ; h^a = h1^a1 ... hn^an
 109 ; de recuperer la partition [1,a1,2,a2,...,n,an] sous type 2 (cf. structure)
 110 ; en representation creuse ,i.e. si hi=0 on ne retrouve par le couple (i,ai).
 111 ;........................................................................
 112 (defun pol2part ($pol listofvar)
 113   (let ((i 0) (varetdegre (chvaretdegre listofvar)))
 114     (mapcar #'(lambda (l)
 115                (setq i (1+ (cdr varetdegre)))
 116                (cons (car l)
 117                      (mapcan #'(lambda (nb)
 118                                 (setq i (1- i))
 119                                 (and (not (eql nb 0)) (list i nb)))
 120                              (nreverse (cdr l)))))
 121             (lect $pol
 122                   (cons 'aa (lvar_lettre (cdr varetdegre) nil
 123                       (flet ((franz.concat (&rest args)
 124                                  "equivalent to Franz Lisp 'concat'."
 125                                  (values (intern
 126                                           (format nil "~{~A~}" args)))))
 127                         (franz.concat '$ (car varetdegre)))))))))
 128
 129 (defun chvaretdegre (listofvar)
 130   (let ((hj (cdr (flet ((franz.exploden (arg)
 131                             "equivalent to Franz Lisp 'exploden'."
 132                             (map 'list #'char-code
 133                                  (prin1-to-string arg))))
 134                    (franz.exploden
 135                        (car (last (sort listofvar
 136                                      'string-lessp)))))))
 137         (i 1))
 138      (cons (flet ((franz.ascii (charcode)
 139                      "equivalent to Franz Lisp 'ascii'."
 140                      (intern (string (code-char charcode)))))
 141             (franz.ascii (car hj)))
 142           (apply '+
 143                  (mapcar #'(lambda (nbascii)
 144                             (prog1 (* i (- nbascii 48))
 145                               (setq i (* i 10))))
 146                          (reverse (cdr hj)))))))
 147 ;........................................................................
 148 ; si au depart on a :
 149 ;REP[part]=[partition](2) au depart
 150 ;REP[part]=[partition](1) en sortie mais sous forme :
 151 ; ( ((mlist).termpart) ...) qui permet d'utiliser la fonction : somme
 152 ; du fichier util.l qui s'attend a trouver la longueur des partition en tete
 153 ; de chaque terme partitionne. On le remplace donc par '(mlist)
 154 ; qui n'est pas inutile
 155 ; [partition](2)
 156
 157 (defun schur2comp (part)
 158   (let ((part1 (ch1rep part)))
 159     (mapcar #'(lambda (2part)
 160                (cons '(mlist)
 161                      (cons (kostka 2part part1) (reverse 2part))))
 162             (mapcar 'ch1rep (treinat part)))))
 163 ;========================================================================
 164 ;                 PASSAGE DES FORMES MONOMIALES
 165 ;                     AUX FONCTIONS DE SCHUR
 166 ; On donne une partition renversee repre'sentant une fonction de Schur
 167 ; on recupere un polynome symetrique contracte.
 168 ; REP(partition) = [partition](1)
 169 ;========================================================================
 170 ; dans util.l ==> ($x1 ... $xn)
 171 ; pour $fadd_sym
 172
 173 (defun $mon2schur_init ($rpart)
 174   (let ((part (reverse (cdr $rpart)))
 175         (lvar (lvar (apply '+ (cdr $rpart)) nil)))
 176              ($fadd_sym
 177                  (cons 0
 178                       (mapcar #'(lambda (2part)
 179                                  (ecrit_mon 2part lvar (kostka part 2part)))
 180                               (mapcar 'duale21 (treinat (duale12 part))))))))
 181
 182 ; etant donne un partition de representation [partition](1), on
 183 ; recupere sa forme duale sous forme [partition](2)
 184 (defun duale12 (partition)
 185   (let ((ai 0))
 186     (nreverse
 187         (mapcon #'(lambda (restpart)
 188                    (setq ai
 189                          (1+ ai))
 190                    (cond
 191                      ((null (cdr restpart))
 192                       (append restpart (list ai)))
 193                      (t (let ((mi (- (car restpart)
 194                                      (cadr restpart))))
 195                           (and (< 0 mi)
 196                                (list mi ai))))))
 197                 partition))))
 198 ;etant donne un partition de representation [partition](2), on
 199 ; recupere sa forme duale sous forme [partition](1)
 200 ;(defun duale21 (partition)
 201 ;  (let ((m1 (cadr partition)))
 202 ; (2duale21 (cddr partition) (list m1)
 203 ;          (* m1 (car partition))
 204 ;         m1)))
 205 ;(defun 2duale21 (part2 part1 p1 p2)
 206 ;  (cond ((null part2) (nconc part1 (make-list (- p1 p2) :initial-element 1)))
 207 ;       (t (let ((nxpart (+ (cadr part2) (car part1))))
 208 ;         (2duale21 (cddr part2)
 209 ;                  (cons nxpart
 210 ;                        part1)
 211 ;                      (+ p1 (* (car part2) (cadr part2)))
 212 ;                     (+ p2 nxpart))))))
 213 (defun duale21 (partition)
 214   (let ((lmultiplicites_lparts
 215             (chmultiplicites_parts partition nil nil)))
 216     (2duale21 (car lmultiplicites_lparts)
 217         (cons 0 (cdr lmultiplicites_lparts)) nil)))
 218 (defun 2duale21 (lmulti lpart partition1_duale)
 219   (cond
 220     ((null (cdr lmulti))
 221      (nconc partition1_duale
 222             (make-list
 223                 (- (cadr lpart) (car lpart))
 224                  :initial-element  (car lmulti))))
 225     (t (2duale21 (cdr lmulti) (cdr lpart)
 226            (nconc partition1_duale
 227                   (make-list
 228                         (- (cadr lpart) (car lpart))
 229                        :initial-element  ($fadd_sym lmulti)))))))
 230
 231 (defun chmultiplicites_parts (partition lmulti lpart)
 232   (cond
 233     ((null partition) (cons lmulti lpart))
 234     (t (chmultiplicites_parts (cddr partition)
 235            (cons (cadr partition) lmulti) (cons (car partition) lpart)))))
 236
 237 ;========================================================================
 238 ;                  NOMBRES DE KOSTKA
 239 ;            (Algorithme : Philippe ESPERET)
 240 ;========================================================================
 241 ; REP(partition) = [part](1)
 242
 243 (defun $kostka_init ($1part $2part)
 244   (kostka (cdr $1part) (cdr $2part)))
 245
 246 (defun kostka (l m)
 247  (list-length (good_tab0 l (make-list (apply '+ l) :initial-element 0) m)))
 248
 249 (defun schur-firstn (n l)
 250   (cond
 251     ((null l) nil)
 252     ((plusp n)
 253      (cons (car l)
 254            (schur-firstn (1- n)
 255                    (cdr l))))
 256     (t nil)))
 257 ; normalement cette fonction existe en Common  sous le nom de "last"
 258 (defun lastn (l n)
 259    (nreverse (schur-firstn n (reverse l))))
 260
 261 (defun good_tab0 (l lcont ltas)
 262   (let ((l1 nil) (rep nil) (relais nil))
 263        (cond
 264          ((eql 1 (list-length l))
 265           (mapcar 'list (good_line (car l) lcont ltas)))
 266          (t
 267           (setq l1 (good_line (car l) lcont ltas))
 268           ;; (print "tete des tableaux possibles " L1)
 269           (do nil
 270               ((null l1))
 271             (setq relais
 272                   (good_tab0 (cdr l) (car l1) (new_tas0 (car  l1) ltas)))
 273             ;; (print " car L1 future tete "(car L1) " et relais "relais)
 274             (if (not relais) (setq l1 (cdr l1))
 275                 (setq rep (nconc rep (insert_tete (car l1) relais))
 276                       l1 (cdr l1))))
 277           rep))))
 278
 279 ;L liste de listes : retourne la meme liste ou les listes ont ete modifiees
 280 ; par insertion de i en tete
 281 (defun insert_tete (i l)
 282   (if (null l) (list (list i))
 283       (append (mapcar #'(lambda (z) (cons i z)) l))))
 284 ;ote du tas Ltas les elts de L1,avec les not ci-dessus Ltas=(3 2 1) cad
 285 ; 3 "1", 2 "2" et 1 "3"
 286 (defun new_tas0 (l1 ltas)
 287   (if (null l1) ltas
 288       (new_tas0 (cdr l1)
 289           (append (schur-firstn (1- (car l1))
 290                           ltas)
 291                   (list (1- (nth (1- (car l1)) ltas)))
 292                   (lastn ltas (- (list-length ltas) (car l1))
 293                          )))))
 294 (defun good_line (taille lcontrainte ltas)
 295   (good_length taille (good_line0 taille lcontrainte ltas)))
 296
 297
 298 (defun good_line0 (taille lcontrainte ltas)
 299   (let ((i 0) (lotas (list-length ltas)) (avanti nil) (rep nil))
 300           ; (print "taille = "taille "  Ltas" Ltas "GREP "rep)
 301        (unless (or (null lcontrainte) (zerop taille))
 302          (setq i (1+ (car lcontrainte)))
 303          (do nil
 304              ((< lotas i))
 305            (if (zerop (nth (1- i) ltas))
 306                (setq i (1+ i))
 307                (setq rep
 308                      (append rep
 309                              (insert_tete
 310                               i
 311                               (good_line0 (1- taille)
 312                                           (cdr lcontrainte)
 313                                           (append
 314                                            (make-list (1- i)
 315                                                       :initial-element 0)
 316                                            (list (1- (nth (1- i) ltas)))
 317                                            (lastn ltas (- lotas i))
 318                                            ))))
 319                      i (1+ i)
 320                      avanti t)))
 321          (if avanti rep nil))))
 322
 323 (defun good_length (taille l)
 324   (if (null l) nil
 325       (if (eql taille (list-length (car l)))
 326           (cons (car l) (good_length taille (cdr l)))
 327           (good_length taille (cdr l)))))
 328
 329 ;=======================================================================
 330 ;               TREILLIS DES PARTITIONS DANS L'ORDRE NATUREL
 331 ; ETANT DONNE UNE PARTITION I ON RAMENE LES PARTITIONS DE MEME
 332 ; POIDS INFERIEURES A I.
 333 ;=======================================================================
 334 ; REP(partition) = [part](1) en entree et en sortie
 335
 336 (defun $treinat_init ($partition1)
 337   (cons '(mlist)
 338         (mapcar #'(lambda (part) (cons '(mlist) (ch1rep  part)))
 339                 (treinat (ch2rep (cdr $partition1))))))
 340 ; REP(partition) = [part](2) en entree et en sortie
 341 (defun treinat (part2) (soustreillisnat (list part2) nil))
 342 ; prendre a chaque fois la plus petite dans l'ordre lexicographique
 343 (defun soustreillisnat (lpartnt lpartt)
 344   (cond
 345     ((null lpartnt) lpartt)
 346     (t (soustreillisnat
 347            (union_sym (cdr lpartnt)
 348                      (sort (crefils_init (car lpartnt)) '$lex))
 349            (cons (car lpartnt) lpartt)))))
 350 ; deux listes identiquement ordonnees par lex
 351 (defun union_sym (l1 l2)
 352   (cond
 353     ((null l2) l1)
 354     ((null l1) l2)
 355     ((equal (car l1) (car l2)) (union_sym l1 (cdr l2)))
 356     (($lex (car l1) (car l2)) (union2 l1 l2) l1)
 357     (t (union2 l2 l1) l2)))
 358 (defun union2 (l1 l2)
 359   (and l2
 360        (cond
 361          ((null (cdr l1)) (rplacd l1 l2))
 362          (t (let ((part1 (cadr l1)) (part2 (car l2)))
 363               (cond
 364                 ((equal part1 part2) (union2 l1 (cdr l2)))
 365                 (($lex part1 part2) (union2 (cdr l1) l2))
 366                 (t (let ((ll1 (cdr l1)))
 367                      (union2 (cdr (rplacd l1 l2)) ll1)))))))))
 368 (defun crefils_init (part) (crefils (reverse part) nil nil))
 369 ; part = (an mn ... a2 m2 a1 m1) an > ... > a2 > a1
 370 ; debut = (a(i-1) m(i-1) ... a1 m1)
 371 ; rfin = (mi ai m(i+1) a(i+1) ...)
 372 ; evite l'identite
 373 (defun crefils (rfin debut lfils)
 374   (cond
 375     ((null rfin) lfils)
 376     (t (let ((ai (cadr rfin)) (mi (car rfin)) (rfin (cddr rfin)))
 377          (cond
 378            ((and (null rfin) (eql 1 mi)) lfils)
 379            (t (crefils rfin (cons ai (cons mi debut))
 380                        (cons (tombecube rfin ai mi debut) lfils))))))))
 381 ; ai --> ai-1 et mi reste si ai > 1
 382 ; disparition ai --> ai-1 = 0
 383 ; disparition ai --> ai-1 = 0
 384 ; ai --> ai + 1 = 2
 385 ; il en reste mi-2 egales a 1
 386 (defun tombecube (rfin ai mi debut)
 387   (cond
 388     ((eql 1 mi)
 389      (cond
 390        ((eql 1 ai) (reverse (arrivecube rfin)))
 391        (t (nconc (reverse (arrivecube rfin))
 392                  (schur-met (1- ai)
 393                       debut)))))
 394     (t (cond
 395          ((eql 1 ai)
 396           (cond
 397             ((eql 2 mi) (reverse (rmet 2 rfin)))
 398             (t (reverse (cons (- mi 2)
 399                               (cons 1 (rmet 2 rfin)))))))
 400          (t (cond
 401               ((eql 2 mi)
 402                (nconc (reverse (rmet (1+ ai)
 403                                      rfin))
 404                       (schur-met (1- ai)
 405                            debut)))
 406               (t (nconc (reverse (cons (- mi 2)
 407                                        (cons ai
 408                                         (rmet
 409                                          (1+ ai)
 410                                          rfin))))
 411                         (schur-met (1- ai)
 412                              debut)))))))))
 413 ; rpart = (m a ...)
 414 ; aj = a(i-1) ==> m(i-1) --> m(i-1) +1
 415 (defun rmet (aj rpart)
 416   (cond
 417     ((null rpart) (list 1 aj))
 418     ((eql aj (cadr rpart))
 419      (cons (1+ (car rpart))
 420            (cdr rpart)))
 421     (t (cons 1 (cons aj rpart)))))
 422 ; part = part2 sens croissant des parts
 423 (defun schur-met (aj part)
 424   (cond
 425     ((null part) (list aj 1))
 426     ((eql aj (car part))
 427      (cons aj
 428            (cons (1+ (cadr part))
 429                  (cddr part))))
 430     (t (cons aj (cons 1 part)))))
 431 ; part = (mj aj ...) un aj passe a  aj+1
 432 ; mj = 1
 433 (defun arrivecube (rpart)
 434   (cond
 435     ((eql 1 (car rpart))
 436      (rmet (1+ (cadr rpart))
 437            (cddr rpart)))
 438     (t (cons (1- (car rpart))
 439              (cons (cadr rpart)
 440                    (rmet (1+ (cadr rpart))
 441                          (cddr rpart)))))))
 442
 443
 444
 445
 446
 447