share/sym/partpol.lisp

   1 ; partpol.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 ;       CHANGEMENTS DE REPRESENTATIONS SUR k[y1,..., yn][x1,..., xn]
  21 ;=========================================================================
  22 ; explose
  23 ; contract (avec ou sans test)
  24 ; partpol   ..................
  25 ; part2cont et cont2part
  26 ; orbit
  27 ;=============================================================================
  28 ;                         INTERFACE
  29
  30 (in-package :maxima)
  31 (macsyma-module partpol)
  32
  33 (meval '((mdefine) (($tpartpol) $mpol $lvar)
  34                    ((mprog) (($operation)) (($tpartpol_init) $mpol $lvar))))
  35
  36 (meval '((mdefine) (($p_tpartpol) $mpol $lvar)
  37                    ((mprog) (($operation)) (($p_tpartpol_init) $mpol $lvar))))
  38
  39 ; Passage d'un polynome symetrique sous la forme rat a ses partitions
  40 (meval '((mdefine) (($partpol) $pol $lvar)
  41                    ((mprog) (($operation)) (($partpol_init) $pol $lvar))))
  42
  43 ; CONTRACTION D'UN POLYNOME SYMETRIQUE
  44
  45 (meval '((mdefine) (($tcontract) $psym $lvar)
  46                    ((mprog) (($operation)) (($tcontract_init) $psym $lvar))))
  47
  48 (meval '((mdefine) (($contract) $psym $lvar)
  49                    ((mprog) (($operation)) (($contract_init) $psym $lvar))))
  50
  51 ; PASSAGE D'UN POLYNOME SYMETRIQUE CONTRACTE A LA LISTE
  52 ; ET INVERSEMENT
  53
  54 (meval '((mdefine) (($cont2part) $pcont $lvar)
  55                    ((mprog) (($operation)) (($cont2part_init) $pcont $lvar))))
  56
  57 (meval '((mdefine) (($part2cont) $ppart $lvar)
  58                    ((mprog) (($operation)) (($part2cont_init) $ppart $lvar))))
  59
  60 (meval '((mdefine) (($explose) $pc $lvar)
  61                    ((mprog) (($operation)) (($explose_init) $pc $lvar))))
  62 ;*****************************************************************************
  63
  64
  65
  66
  67 ;                   DECLARATIONS AU COMPILATEUR
  68 ; $tcontract_init $contract_init
  69 ; $p_tpartpol_init p_tpartpol
  70 ; $tpartpol_init tpartpol
  71 ; $partpol_init partpol  appellees par apply
  72 ; $cont2part_init cont2part
  73 ; $partipol partipol
  74 ;** FTOC. WARNING:
  75 ;             Franz Lisp declaration 'localf' is currently untranslated
  76 (progn)
  77 ;$explose_init
  78 (progn (defvar drapeaux) (defvar modele))
  79 (progn (defvar lvar) (defvar permut))
  80 ;*****************************************************************************
  81 ;                FORME CONTRACTE D'UN POLYNOME SYMETRIQUE
  82 ; appels : contract(pol,lvar) ou tcontract(pol,lvar) t pour test
  83 ;     $lvar = [x1, ..., xp] = ((mlist) x1, ..., xp) au depart
  84 ;     lvar  = (x1, ..., xp)
  85 ;     card = p
  86 ;-----------------------------------------------------------------------------
  87 ;                     AVEC TEST DE SYMETRIE
  88 ; tpartpol ramene des partition sout la forme [part](2). Pour
  89 ; l'ecrivain
  90 ; de polyn\^ome on peut directement utiliser la fonction 2ecrit.
  91 ; Sinon, pour utiliser $distri_ecrit (i.e. ecrit_pol), il faut faire
  92 ; ch1repol avant.
  93 ;-----------------------------------------------------------------------------
  94 (defun $tcontract_init ($pol $lvar)
  95   (if (meval (list '($is) (list '(mequal) $pol 0))) 0
  96       (catch 'rate
  97         (2ecrit (tpartpol $pol $lvar) (cdr $lvar)))))
  98 ;-----------------------------------------------------------------------------
  99 ;                    SANS TEST DE SYMETRIE
 100 ;-----------------------------------------------------------------------------
 101 (defun $contract_init ($pol $lvar)
 102   (if (meval (list '($is) (list '(mequal) $pol 0))) 0
 103       (meval (list
 104                   '($distri_ecrit)
 105                   (meval (list '($partpol) $pol $lvar))
 106                   $lvar))))
 107 ;======================================================================
 108 ;              PASSAGE D'UN POLYNOME SYMETRIQUE MACSYMA
 109 ;                    A SA REPRESENTATION PARTITIONNEE.
 110 ;
 111 ; $partpol_init ramene REP(pol)(1) ............ lexicographique decroissant
 112 ; $partpol      ...... REP(pol)(2) ............ lexicographique decroissant
 113 ; tpartpol si on desire en plus tester la symetrie
 114 ; ----------------------------------------------------------------------------
 115 ;                    AVEC TEST DE SYMETRIE
 116 ; $pol est un polynome macsyma sous la forme expand
 117 ; Si on part d'un polynome non symetrique une erreur est declanchee
 118 ; $p_tpartpol_init ramene REP(pol)(1) dans l'ordre des longueurs decroissant
 119 ; p_tpartpol       ...... REP(pol)(2) ......................................
 120 ;============================================================================
 121 ; lpol est sous la forme((coe . [partition](2))...)  mais n'est pas
 122 ; un polynome partitionne car on n'a encore retire aucun monome.
 123 (defun $p_tpartpol_init ($pol $lvar)
 124   (if (meval (list '($is) (list '(mequal) $pol 0))) (cons '(mlist) nil)
 125       (catch 'rate
 126         (macsy_list (ch1repol (p_tpartpol $pol $lvar))))))
 127
 128 ; appele par pui en drapeau 5 &&&&&&mais si erreur pas de cath !
 129
 130 (defun p_tpartpol ($pol $lvar)
 131    (if (meval (list '($is) (list '(mequal) $pol 0))) (cons '(mlist) nil)
 132        (let ((lpol (sort (prep (lect $pol (cons '(mlist) $lvar)))
 133            '$orlong)))
 134     (tpartpol2 lpol (list-length (cdr $lvar))))))
 135
 136 (defun $tpartpol_init ($pol $lvar)
 137    (if (meval (list '($is) (list '(mequal) $pol 0))) (cons '(mlist) nil)
 138        (catch 'rate
 139               (macsy_list (ch1repol (tpartpol $pol $lvar))))))
 140
 141 ;appele par e_sym (22 ..)  &&&&&&&& pas de catch si erreur
 142 ; pol est sous la forme((coe .I)...) les +petit part devant
 143
 144 (defun tpartpol ($pol $lvar)
 145   (let ((pol (sort (prep (lect $pol $lvar)) '$lexic)))
 146     (tpartpol2 pol (list-length (cdr $lvar)))))
 147
 148 (defun $lexic (p q)
 149   (and (not (equal (cdr p) (cdr q))) ($lex (cdr p) (cdr q))))
 150 ;--------------------------------------------------------------------------
 151 ;                    LE PARTITIONNEUR
 152 ;-----------------------------------------------------------------------------
 153 ; Reunir les partitions identiques, rendre le coefficient de la
 154 ; forme monomiale associee.
 155 (defun tpartpol2 (lpol card)
 156   (tpartpol3 card lpol (card_orbit (cdar lpol) card))
 157   lpol)
 158 ;-----------------------------------------------------------------------
 159 ;              partitionnement avec test de symetrie : tpartpol3
 160 ;                 modification physique sur la forme distribuee
 161 ; Tout les monome ont leur exposants represente par [partition](2)
 162 ; on a perdu l'information : "exposant attache a une variable"
 163 ; dont on n'a pas besoin.
 164 ; Un compteur
 165 ; Le monome de tete de pol = (coe1 . exposants1)
 166 ; le deuxieme monome de tete de pol = (coe2 . exposants2)
 167 ; 1- exposants1 est le (p+1)ieme que
 168 ;      l'on trouve identique a exposants2
 169 ;          (on a reordonne les exposants)
 170 ;     1-1 coe1 distinct de coe2 ==> non symetrique
 171 ;     1-2 coe1 = coe2 ==> (compteur ==> compteur-1) : on a un element
 172 ;                         supplementaire de l'orbite de exposants2
 173 ;                           on l'elimine de pol puisqu'un seul representant
 174 ;                          suffit
 175 ; 2- exposant1 different de exposants2
 176 ;     1-1 compteur > 1 ==> Il manque des elements de l'orbit de exposants2
 177 ;     1-2 compteur = 1 ==> - (coe1 . exposants1) represente toute l'orbite
 178 ;                            contenu dans pol
 179 ;                          - On passe a l'orbite suivante qui
 180 ;                              est (coe2 . exposants2) en mettant
 181 ;                              compteur = card_orbite(exposants2)
 182 ;--------------------------------------------------------------------------
 183 ; pas le cas de la cste
 184 ; les coe sont egaux
 185 ; ex : 3xy + 3yz + 2xz est non symetrique
 186 (defun tpartpol3 (card pol compteur)
 187   (if (null (cdr pol))
 188       (or (eql 1 compteur)
 189           (throw 'rate '|manque des elements de l'orbite|))
 190       (let ((coe1 (caar pol)) (coe2 (caadr pol)) (exp1 (cdar pol))
 191             (exp2 (cdadr pol)))
 192         (if (equal exp1 exp2)
 193             (if (equal coe1 coe2)
 194                 (tpartpol3 card (rplacd pol (cddr pol))
 195                     (1- compteur))
 196                 (throw 'rate
 197                       '|polynome non symetrique par ses coefficients|))
 198             (if (eql 1 compteur)
 199                 (tpartpol3 card (cdr pol) (card_orbit exp2 card))
 200                 (throw 'rate '|manque des monomes|))))))
 201 ;-----------------------------------------------------------------
 202 ;             PARTITIONNEMENT D'UN POLYNOME
 203 ; ENTREES : UN POLYNOME $psym  SOUS FORME DISTRIBUEE OU SINON DONNES PAR UNE
 204 ;           REPRESENTATION MACSYMA ET LA LISTE $lvar DE SES VARIABLES
 205 ; SORTIE  : LE POLYNOME PARTITIONNE AVEC COMME REPRESENTATION :
 206 ;            REP([pol]) = [ppart](1)
 207 ;-----------------------------------------------------------------
 208 (defun $partpol_init ($psym $lvar)
 209  (if (meval (list '($is) (list '(mequal) $psym 0))) (cons '(mlist) nil)
 210      (cons '(mlist)
 211         (mapcan #'(lambda ($exposant)
 212                           (and (apply '>= (cddr $exposant))
 213                                (list $exposant)))
 214                 (if (equal 'mlist (caar $psym))
 215                     (cdr $psym) ; $psym est sous forme distribuee
 216                     (cdr (meval (list '($distri_lect)
 217                                        $psym
 218                                        $lvar))))))))
 219
 220 ;-----------------------------------------------------------------------
 221 ;         FORME CONTRACTEE D'UN POLYNOME ==> REP([pol])(i)
 222 ; $cont2part_init --> [$ppart](1) ($ car liste macsyma)
 223 ; cont2part et partipol --> [ppart](2)
 224 ; appele par $elem_init et $pui_init
 225 ;-----------------------------------------------------------------------
 226 ; dans lect on met $polcontrac sous forme expand
 227 (defun $cont2part_init ($polcontrac $lvar)
 228       (if (meval (list '($is) (list '(mequal) $polcontrac 0)))
 229           (cons '(mlist)  nil)
 230           (macsy_list (cont2part_rep1 $polcontrac $lvar))))
 231
 232 (defun cont2part_rep1 (polcontrac $lvar)
 233   (ordonne_expo (lect polcontrac $lvar)))
 234 ; Rend des partitions de type 2
 235
 236 (defun cont2part (polcontrac lvar)
 237   (ch2repol (cont2part_rep1 polcontrac (cons '(mlist) lvar))))
 238 ; on le met sous forme rat
 239 ; rat = (((rat simp)...) . listebase . 1)
 240
 241 (defun $partipol (polycontracte)
 242   (let ((rat ($rat polycontracte)))
 243     (cons '(mlist) (mapcar 'list2croch (partipol rat)))))
 244
 245 (defun partipol (rat) (partipol1 (cadr rat) (cons nil nil)))
 246
 247 ; on aura lpart = ( (coe a1 m1 a2 m2 ....) ....)  par
 248 ; accrochage de bord
 249
 250 (defun partipol1 (rat lpart)
 251   (partipol2 lpart (cdr rat) nil nil)
 252   (cdr lpart))
 253
 254 ; representation des partitions avec un coe et eventuellement des
 255 ; exposants nuls en tete.
 256 ; rat est la liste elementaire
 257 ;on depile ie on remonte
 258 ; on a enfin une partition solution
 259 ; deplacement en profondeur dans
 260 ; le coefficient qui est une liste de base et auquel on enleve la variable.
 261 ; on met en instance le deplacement en largeur dans rrat
 262 ; on note dans la pile l'etat actuel : la pile la partition
 263 ; en construction et le reste des exposant coe a parcourir
 264 ; en largeur
 265
 266 (defun partipol2 (lpart rat part pile)
 267   (if (null rat) (and pile (apply 'partipol2 (cons lpart pile)))
 268       (if (numberp (cadr rat))
 269           (partipol2
 270               (cdr (rplacd lpart
 271                            (list ($part0 (sort
 272                                             (cons (car rat) part) '<)
 273                                          (cadr rat)))))
 274               (cddr rat) part pile)
 275           (partipol2 lpart (cdadr rat) (cons (car rat) part)
 276               (if (cddr rat) (list (cddr rat) part pile) pile)))))
 277 ;=======================================================================
 278 (defun $part2cont_init ($ppart $lvar)
 279   (if (null (cdr $ppart)) 0
 280        (meval (list '($distri_ecrit) $ppart $lvar))))
 281 ;========================================================================
 282 ;                 RAMENER TOUT LE POLYNOME SYMETRIQUE
 283 ;                   ASSOCIE A UNE FORME CONTRACTEE
 284 ;                   dans k[y1, ... ,yn][x1, ... ,xn]
 285 ; EN SE SERVANT DE permut ECRIT PAR PHILIPPE ESPERET
 286 ; appel : explose(polynome,[x1, ... ,xp])
 287 ;-------------------------------------------------------------------------
 288 ;                     L'APPEL PRINCIPAL
 289 ; lvar = [x1, ... ,xp] on ne demande pas de preciser y1, ..., yn
 290 ; On utilise le lecteur rendant les formes distribuees et
 291 ; l'ecrivain qui se trouvent dans le fichier util.l
 292 ;-------------------------------------------------------------------------
 293
 294 (defun $explose_init ($pc $lvar)
 295          (if (meval (list '($is) (list '(mequal) $pc 0))) 0
 296   (let ((lcoenuplet (lect $pc  $lvar)))
 297     (meval (list '($distri_ecrit)
 298                   (cons '(mlist)
 299                         (mapcan #'(lambda (coenuplet)
 300                                     (let ((coe (car coenuplet)))
 301                                            (mapcar #'(lambda (permu)
 302                                                         (list* '(mlist)
 303                                                                 coe permu))
 304                                                    (permut (cdr coenuplet)))))
 305                                  lcoenuplet))
 306                    $lvar)))))
 307 ;=======================================================================
 308 ;                ORBITE D'UN POLYNOME
 309 ;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 310 ;========================================================================
 311 ;                       UTILITAIRES LOCAUX
 312 ;------------------------------------------------------------------------
 313 ; prep : Fonction ordonnant tout d'abord tout les exposants
 314 ; dans le sens decroissant puis representent ces exposants
 315 ; sous forme [partition](2) (sans les 0)
 316 ; Creons tout d'abords une fonction, ordonne_expo, qui a partir d'un polynome
 317 ; sous forme distribuee range tout les exposants dans l'ordre
 318 ; decroissant
 319
 320 (defun ordonne_expo (pol)
 321   (mapcar #'(lambda (mon)
 322              (cons (car mon)
 323                    (sort (cdr mon) '>)))
 324           pol))
 325
 326 ; Avec changement de repre'sentation
 327 (defun prep (pol) (ch2repol (ordonne_expo pol)))