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