share/sym/elem.lisp

   1 ; Fichier elem.lsp
   2
   3 ;       ***************************************************************
   4 ;       *                    MODULE SYM                               *
   5 ;       *       MANIPULATIONS DE FONCTIONS SYMETRIQUES                *
   6 ;       *        (version01: Commonlisp pour Maxima)                 *
   7 ;       *                                                             *
   8 ;       *                ----------------------                       *
   9 ;       *                  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 ;               DECOMPOSITION D'UN POLYNOME SYMETRIQUE
  21 ;                   PAR LES SYMETRIQUES ELEMENTAIRES
  22
  23 ; appelle avec elem([card,e1, e2, ...],sym(x,y,..,z),[x,y,...,z])
  24 ; ou multi_elem pour des polyn\^omes multisym\'etriques
  25 ;=============================================================================
  26 (in-package :maxima)
  27 (macsyma-module elem macros)
  28
  29
  30
  31 (mdefprop $elem
  32     ((lambda ()) ((mlist) $valei $sym $lvar)
  33      ((mprog) (($operation)) (($elem_init) $valei $sym $lvar)))
  34     mexpr)
  35
  36 ;; IT APPEARS ARGS WAS A MACRO. THERE IS NO ARGS MACRO AT PRESENT.
  37 ;; DUNNO IF THE ABSENCE OF ARGS CAUSES ANY INCORRECT BEHAVIOR IN SYM
  38 ;; (args $elem '(3 . 3))
  39
  40 (add2lnc '(($elem) $valei $sym $lvar) $functions)
  41
  42 (mdefprop $multi_elem
  43     ((lambda ()) ((mlist) $lvalei $pc $llvar)
  44      ((mprog) (($operation)) (($multi_elem_init) $lvalei $pc $llvar)))
  45     mexpr)
  46
  47 (add2lnc '(($multi_elem) $lvalei $pc $llvar) $functions)
  48
  49 ;================================================================
  50 ; fonction bidon de chargement pour eviter de construire pour detruire
  51 ; lorsque l'on appelle une fonction de elem a partir d'un autre
  52 ; fichier du module sym
  53 (defun $bidon ())
  54 ;---------------------------------------------------------------------------
  55 ;           VARIABLES DECLAREES SPECIALES PAR LE COMPILATEUR
  56 (progn
  57   (defvar listei)
  58   (defvar $elem)
  59   (defvar nb1)
  60   (defvar lgI)
  61   (defvar coei)
  62   (defvar nblib))
  63
  64 ;***************************************************************************
  65 ;          MISE SOUS FORME INTERNE DU POLYNOME SYMETRIQUE
  66 ;                SUIVANT LES FORMES EXTERNES DONNEES
  67 ; Donnees :
  68 ; valei = ((mlist) card e1 e2 ...)
  69 ; sym est un polynome symetrique pouvant etre represente
  70 ; de plusieurs manieres en entree .
  71 ; lvar = ((mlist) x1 x2 ...) les variables de sym.
  72 ; Representation interne : REP([pol]) = [lppart](2)
  73 ;                          listei=(card e1 e2 ...)
  74
  75 ;----------------------------------------------------------------------------
  76 ;                MULTIDECOMPOSITION
  77 ; Le polynome donne est multi-symetrique sous forme contractee
  78 ;----------------------------------------------------------------------------
  79 (defun $multi_elem_init ($multi_lelem $multi_pc $llvar)
  80   (multi_elem (mapcar 'cdr (cdr $multi_lelem)) $multi_pc
  81               (cdr $llvar)))
  82
  83 ; cf. e_red1 plus loin
  84
  85 (defun multi_elem (multi_lelem $multi_pc l$lvar)
  86   (cond
  87     ((meval (list '($is) (list '(mequal) $multi_pc 0))) 0)
  88     ((null l$lvar) $multi_pc)
  89     (t (multi_elem (cdr multi_lelem)
  90               (if (meval (list '($is) (list '(mequal) $multi_pc 0))) 0
  91                   (e_red1 (car multi_lelem)
  92                           (lgparts (ch2repol
  93                                  (mapcar 'cdr
  94                                          (cdr (meval
  95                          (list '($cont2part) $multi_pc
  96                                          (car l$lvar)))))))))
  97            (cdr l$lvar)))))
  98
  99 ;---------------------------------------------------------------------------
 100
 101
 102 (defun $elem_init (valei sym $lvar)
 103   (let ((sauvlistei
 104             (cdr (flet ((franz.boundp (name)
 105                             "equivalent to Franz Lisp 'boundp'."
 106                             (and (boundp name)
 107                                  (cons nil (symbol-value name)))))
 108                    (franz.boundp 'listei)))))
 109     (prog1 (case $elem
 110              (1 ; sym = polynome contracte
 111               (if (meval (list '($is) (list '(mequal) sym 0))) 0
 112                   (e_red1 (cdr valei)
 113                       (lgparts (ch2repol
 114                                   (mac2lisp (meval
 115                         (list '($cont2part) sym $lvar))))))))
 116              (2 ;le polynome symetrique en entier ou en partie
 117               (if (meval (list '($is) (list '(mequal) sym 0))) 0
 118                   (e_red1 (cdr valei)
 119                       (lgparts (ch2repol
 120                                  (mac2lisp (meval
 121                            (list '($partpol) sym $lvar))))))))
 122              (3 ; sym=REP([pol])(1) mais pas forcement ordonne'
 123                 ; mais les monomes sont tous suppose's distincts
 124               (e_red1 (cdr valei)
 125                       (lgparts (ch2repol (mapcar 'cdr (cdr sym))))))
 126              (4 ; sym est le polynome symetrique
 127                 ; on test egalement sa symetrie
 128               (let ((pol (lgparts (ch2repol
 129                                     (mac2lisp (meval
 130                       (list '($tpartpol) sym $lvar)))))))
 131                 (e_red2 ($degrep pol) pol (cdr valei) )))
 132              (5 ; sym = (REP([pol])(2) + longueurs) retirer les "mlist"
 133                (e_red1 (cdr valei) (mapcar 'cdr (cdr sym))))
 134              (6 ; sym = REP([pol])(2)
 135                (e_red1 (cdr valei) (lgparts (mapcar 'cdr (cdr sym)))))
 136              (t "erreur $elem n'a pas de valeur"))
 137       (setq listei sauvlistei))))
 138
 139 (defun e_red1 (l ppart)
 140   (e_red2 ($degrep ppart)
 141           (sort ppart '$e_lexinv) l))
 142
 143
 144 (defun e_red2 (degpol ppart l)
 145   (cond
 146     ((eql 0 (lgi ppart)) (coei ppart)) ; on n'a qu'une constante
 147     (t (setq listei
 148              (rangei l
 149                      (if (and l (numberp (car l)))
 150                          (min (car l) degpol) ; le cardinal est impose
 151                          degpol)
 152                      (list-length l)))
 153          ; autant que l'inf du cardinal  et du degre du polynome
 154        ($reduit (if (numberp (car l)) (car l) degpol) ppart))))
 155
 156 ;---------------------------------------------------------------------------
 157 ;           CREATION DE LA LISTE listei DES VALEURS DES ELEMENTAIRES
 158 ;l=(card e1 e2 ... e(lg))  card est le cardinal de l'alphabet.
 159 ; avec ki < k(i+1)
 160 ;----------------------------------------------------------------------------
 161 ; on range les plus grand en premier
 162
 163 (defun rangei (l n lg)
 164   (if (eql (1+ n) lg)
 165       l (append l (rangei2 nil lg n))))
 166
 167 (defun rangei2 (lesei i n)
 168   (if (< n i)
 169       (nreverse lesei)
 170       (rangei2 (cons (flet ((franz.concat (&rest args)
 171                                 "equivalent to Franz Lisp 'concat'."
 172                                 (values (intern
 173                                          (format nil "~{~A~}" args)))))
 174                        (franz.concat '$e i))
 175                      lesei)
 176                (1+ i)
 177                n)))
 178
 179 ;--------------------------------------------------------------------------
 180 ;                   LA BOUCLE PRINCIPALE
 181 ; sym = [lppart](2) ordonnee dans l'ordre lexicographique decroissant.
 182 ;-------------------------------------------------------------------------
 183
 184 (defun $reduit (card sym)
 185   (let ((I (moni sym)))
 186     (if (or (null sym) (syele I)) (e_ecrit sym)
 187         ($reduit card
 188                  (somme (cdr sym)
 189                         (devel1 (factpart I)
 190                                 (coei sym) (lgi sym) card)
 191                         '$e_lexinv)))))
 192 ;-------------------------------------------------------------------------
 193 ;                        FACTORISATION DE I
 194 ;--------------------------------------------------------------------------
 195 (defun factpart (i)
 196   (let ((test nil) (alt nil))
 197     (let ((j (mapcar #'(lambda (puiounb)
 198                         (setq alt (null alt))
 199                         (if alt
 200                             (if (eql 1 puiounb)
 201                                 (and (setq test 't) nil)
 202                                 (1- puiounb))
 203                             puiounb))
 204                      i)))
 205       (cond
 206         (test
 207          (setq nb1 (car (last j)))
 208          (nbutlast (nbutlast j)))
 209         (t
 210          (setq nb1 0) j)))))
 211 ;---------------------------------------------------------------------------
 212 ;                             REECRITURE DE I
 213 ;                  Developpement de ei*J ou i= lgI = nb1 + lgJ
 214 ;                J=(pui1 n1 pui2 n2 .....) avec puik > pui(k-1)
 215 ;----------------------------------------------------------------------------
 216
 217 (defun devel1 (J coeI lgI card)
 218    (let ((coeJ ($mult_sym coeI (nth lgI listei)))
 219                (nblib (- card lgI)))
 220          (nconc (and (plusp nblib)
 221                      (devel2 J nblib (cons nil nil)))
 222                 (and (or (not (numberp coeJ))
 223                          (null (zerop coeJ)) )
 224                      (list (list* (- lgI nb1) coeJ J))))))
 225
 226
 227 (defun devel2 (J nblib pilesol)
 228    (devel3 pilesol J 0 (cadr J) (cons -1 nil) nil)
 229    (cddr pilesol)) ; pilesol=(nil I . listparti)
 230
 231 ;----------------------------------------------------------------------------
 232 ;   r le nombre d'elements passant a la meme puissance superieure, pui1 + 1.
 233 ; r vaut n1 au depart et decroit jusqu'a une valeur inf non negative.
 234 ; Ou inf est choisie afin que la forme monomiale representee
 235 ; par la partition ramenee soit non nulle relativement au cardinal, card, de
 236 ; l'alphabet concidere. En fait il faut que la longueur de cette partition
 237 ; qui est de nbpui1 + lgI soit inferieure ou egal a card.
 238 ;    solu est la partition partielle d'une partition solution en construction
 239 ;    pile contient les partitions en construction mais mise en instance
 240 ;    pilesol contient les partition solutions deja construites
 241 ;-----------------------------------------------------------------------
 242
 243 (defun devel3 (pilesol J nbpui1 r solu pile)
 244
 245    (if (null  J)
 246        (progn (rplacd pilesol
 247                       (list (ramsolfin nbpui1 (+ nbpui1 nb1) solu)))
 248               (and pile
 249                   (devel3 (cdr pilesol); pas apply pour recursivite terminale
 250                            (car pile)
 251                            (cadr pile)
 252                            (caddr pile)
 253                            (car (cdddr pile))
 254                            (car (last pile)))))
 255        (let* ((reste (- (cadr J) r))
 256               (nnbpui1 (+ nbpui1 reste)))
 257 ; on met le cas r --> r-1 en instance (si nnbpui1 + lgI < card) en empilant,
 258 ; et on passe tout de suite a r --> n2 pour continuer a construire une
 259 ; partition solution.
 260              (devel3 pilesol
 261                      (cddr J)                         ; (pui2 n2 .....)
 262                      nnbpui1                          ; lg(M) >= nbpui1 + lgI
 263                      (cadr (cddr J))                  ; n2
 264                      (ramsol (car J) reste r solu)
 265                      (if (and (plusp r)
 266                               (> nblib nnbpui1))      ; **
 267                          (list J nbpui1 (1- r) solu pile)
 268                           pile) ))))
 269
 270 ; ** pour ne pas depasser le cardinal de l'alphabet
 271
 272 (defun ramsol (pj nbj r solu)
 273   (if (eql 0 r) (list* (car solu) nbj pj (cdr solu))
 274       (let ((solu (ramsol2 pj r (car solu) (cdr solu))))
 275         (if (eql 0 nbj) solu (list* (car solu) nbj pj (cdr solu))))))
 276
 277 (defun ramsol2 (pj r coe solu)
 278   (if (eql (cadr solu)
 279              (1+ pj))
 280       (list* (calcoe coe (car solu) r)
 281              (+ (car solu) r)
 282              (cdr solu))
 283       (list* coe r
 284              (1+ pj)
 285              solu)))
 286 ; tnb1=0 si sol=I et que nb1=0
 287 (defun ramsolfin (nbpui1 tnb1 solu)
 288   (if (eql 1 (caddr solu))
 289       (list* (+ lgI nbpui1)
 290              ($mult_sym coei (calcoe (car solu) tnb1 (cadr solu)))
 291              (reverse (cons (+ tnb1 (cadr solu))
 292                             (cddr solu))))
 293       (list* (+ lgI nbpui1)
 294              ($mult_sym coei (car solu))
 295              (reverse (list* tnb1 1 (cdr solu))))))
 296 ;-------------------------------------------------------------------------
 297 ;         CALCUL DU COEFFICIENT BINOMIAL C(n+c,c)
 298 ;--------------------------------------------------------------------------
 299 (defun calcoe (coe c n)
 300   (let ((nc (+ n c)))
 301     (* coe (calcoe2 (inferieur n c) nc
 302                         (1- nc)
 303                         2))))
 304
 305 (defun calcoe2 (n res nc i)
 306   (if (eql (1+ n)
 307            i)
 308       res
 309       (calcoe2 n
 310                (div (* res nc)
 311                     i)
 312                (1- nc)
 313                (1+ i))))
 314 ;---------------------------------------------------------------------------
 315 ; syele teste si une partition est celle d'une fonction
 316 ;  symetrique elementaire
 317 (defun syele (mon)
 318   (and mon (and (eql 1 (car mon)) (null (cddr mon)))))
 319 (defun inferieur (a i) (and a (min a i)))
 320 ; ---------------------------------------------------------------------------
 321 ;                       L'ECRIVAIN
 322 ;----------------------------------------------------------------------------
 323 ; une constante
 324 (defun e_ecrit (solu)
 325   (let ((solu (nreverse solu)))
 326     (cond
 327       ((null solu) 0)
 328       ((eql 0 (lgi solu))
 329        (e_ecrit2 (cdr solu) (cdr listei) (coei solu) 1))
 330       (t (e_ecrit2 solu (cdr listei) 0 1)))))
 331 (defun e_ecrit2 (solu listei mpol i_init)
 332   (let ((i (lgi solu)))
 333     (cond
 334       ((null solu) mpol)
 335       ((eql i i_init)
 336        (e_ecrit2 (cdr solu) listei
 337            ($add_sym mpol ($mult_sym (coei solu) (car listei))) i_init))
 338       (t (setq listei
 339                (flet ((franz.nthcdr (ind lis)
 340                           "equivalent to Franz Lisp 'nthcdr'."
 341                           (let ((evalind (eval ind)))
 342                             (if (minusp evalind) (cons nil lis)
 343                                 (nthcdr evalind lis)))))
 344                  (franz.nthcdr
 345                      (- i i_init)
 346                      listei)))
 347          (e_ecrit2 (cdr solu) listei
 348              ($add_sym mpol ($mult_sym (coei solu) (car listei))) i)))))
 349
 350
 351
 352
 353