share/sym/direct.lisp

   1 ; Fichier direct.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 ;                CALCULS D'IMAGES DIRECTES, D'ORBITES ...
  21 ;                    DANS LE CAS LE PLUS GENERAL
  22 ;            LA FONCTION resolvante PEUX ETRE AMENEE A UTILISER CE PROGRAMME
  23 ;        LORSQUE LA FONCTION RESOLVANTE N'A AUCUNE PROPRIETE EXPLOITABLE
  24 ; Si CE N'EST SON ARITE EN COMPARAISON DES DEGRES DES POLYNOMES A TRANSFORMER
  25 ;===========================================================================
  26 ;                        INTERFACE
  27
  28 (in-package :maxima)
  29 (macsyma-module directnew)
  30
  31 ;; Fonctions MACSYMA
  32 (mdefprop $direct
  33     ((lambda ()) ((mlist) $list_pol $x $fonction $list_list_var)
  34      ((mprog) (($operation))
  35       (($direct_init) $list_pol $x $fonction $list_list_var)))
  36     mexpr)
  37 (add2lnc '(($direct) $list_pol $x $fonction $list_list_var) $functions)
  38 (mdefprop $orbit
  39     ((lambda ()) ((mlist) $fonction $list_var)
  40      ((mprog) (($operation)) (($orbit_init) $fonction $list_var)))
  41     mexpr)
  42 (add2lnc '(($orbit) $fonction $list_var) $functions)
  43 (mdefprop $multi_orbit
  44     ((lambda ()) ((mlist) $fonction $list_var)
  45      ((mprog) (($operation)) (($multi_orbit_init) $fonction $list_var)))
  46     mexpr)
  47 (add2lnc '(($multi_orbit) $fonction $list_var) $functions)
  48 (mdefprop $pui_direct
  49     ((lambda ()) ((mlist) $multi_orbit $list_list_var $ldegre)
  50      ((mprog) (($operation))
  51       (($pui_direct_init) $multi_orbit $list_list_var $ldegre)))
  52     mexpr)
  53 (add2lnc '(($pui_direct) $multi_orbit $list_list_var $ldegre) $functions)
  54 ; SYM pour Macsyma
  55 ;============================================================================
  56 ;                     IMAGE DIRECTE
  57 ; P1,....,Pp polynomes dans k[X^i] et de degre di respectivement
  58 ; X^i = (x^i_1, x^i_2, ..., x^i_di) pour i= 1 ... p
  59 ; X = (X^1,X^2,...,X^p)
  60 ; P polynome dans k[X]
  61 ;                ON CALCUL P_*(P1,...,Pp)
  62 ;============================================================================
  63 ;                   DECLARATIONS AU COMPILATEUR FRANZLISP
  64 (progn
  65   (defvar $direct)
  66   (defvar $pui)
  67   (defvar $elem)
  68   (defvar sauvedrapeau))
  69 ; $orbit_init
  70 ; $multi_orbit_init
  71 ; $pui_direct_init
  72 ; $direct_init
  73 ;** FTOC. WARNING:
  74 ;             Franz Lisp declaration 'localf' is currently untranslated
  75 (progn)
  76 ;============================================================================
  77 ;                  Orbite d'un polynome dans k[y1,...,yn] sous S_n
  78 ;                      k eventuellemnt anneau de polynomes
  79 ;----------------------------------------------------------------------------
  80 ; On fait permuter ses variables et on l'applique a chacune
  81 ; de ces permutations. Puis on elimine les egaux au fur et
  82 ; a mesure
  83 ;----------------------------------------------------------------------------
  84 (defun $orbit_init ($p $lvar) (cons '(mlist) (orbit $p  $lvar)))
  85 (defun orbit ($p $lvar)
  86   (let ((p_dist (lect $p $lvar))
  87         (list$pol (list $p)))
  88     (orbit2 list$pol p_dist (cdr $lvar) nil)
  89     list$pol))
  90 ; les permutations circulaires ne changeraient rien
  91 (defun orbit2 (list$pol p_dist f_lvar d_lvar)
  92   (and (cdr f_lvar)
  93        (mapc #'(lambda (f_lvar2)
  94                 (let (($pol (ecrit_pol  p_dist (append d_lvar f_lvar2))))
  95                   (or (contient list$pol $pol)
  96                       (flet ((franz.attach (newelt oldlist)
  97                                  "equivalent to Franz Lisp 'attach'."
  98                                  (progn
  99                                    (rplacd oldlist
 100                                     (cons (car oldlist) (cdr oldlist)))
 101                                    (rplaca oldlist newelt))))
 102                         (franz.attach $pol list$pol))))
 103                 (orbit2 list$pol p_dist (cdr f_lvar2)
 104                         (append d_lvar (list (car f_lvar2)))))
 105              (permut_circu (cdr f_lvar) (list (car f_lvar))))
 106        (orbit2 list$pol p_dist (cdr f_lvar)
 107                (append d_lvar (list (car f_lvar))))))
 108 ; on ne ramene pas l'identite
 109 (defun permut_circu (debut fin)
 110   (cond
 111     ((null (cdr debut)) (list (cons (car debut) fin)))
 112     (t (cons (append debut fin)
 113              (permut_circu (cdr debut) (cons (car debut) fin))))))
 114 (defun contient (list$pol $pol)
 115     (catch 'trouve
 116         (progn
 117           (mapc #'(lambda ($pol2)
 118                    (and (meval (list '($is)
 119                                      (list '(mequal) $pol $pol2)))
 120                                (throw 'trouve t)))
 121                 list$pol)
 122           nil)))
 123 ;==========================================================================
 124 ;        CALCUL DE L'ORBITE DU POLYNOME P SOUS S_d1xS_d2x...xS_dp
 125 ;--------------------------------------------------------------------------
 126 (defun $multi_orbit_init ($p $llvar)
 127   (cons '(mlist) (multi_orbit_init $p (cdr $llvar))))
 128 ; sous S_0
 129 (defun multi_orbit_init ($p llvar)
 130   (cond
 131     ((null llvar) (list $p))
 132     (t (multi_orbit (orbit $p (car llvar)) nil (cdr llvar)))))
 133 ; On se deplace en largeur dans l'arbre ie. on fait agir tout S_i avant
 134 ; de passer a S_(i+1).
 135 ; En d'autres termes : On calcul l'orbite du polynome P sous
 136 ; S_1 x ... x S_i et on en deduit son orbite sous S_1 x ... x S_(i+1).
 137 ; Quand on passe a S_(i+1) si un des polynomes generes par l'action de
 138 ; S_(i+1) (sur un polynome q de l'etape S_i ) est egal a un
 139 ; polynome r (distinct de q bien entendu!) genere par
 140 ; l'action de S_i on elimine froidement r. (Pourquoi refaire ce qui vient
 141 ;  d'etre fait ?)
 142 ; au depart i = 1  et  llvar = (X^2 X^3 ... X^p) (cf. probleme general)
 143 ; on a toute l'orbite sous S_1 x ... x S_(i+1).
 144 ; si i+1 =p
 145 ; passe a i+2
 146 ; on ote de lpoli les polynomes communs a orbit
 147 (defun multi_orbit (lpoli lpoli+1 llvar)
 148   (cond
 149     ((null lpoli)
 150      (cond
 151        ((null (cdr llvar)) lpoli+1)
 152        (t (multi_orbit lpoli+1 nil (cdr llvar)))))
 153     (t (let ((orbit (orbit (car lpoli) (car llvar))))
 154          (epure lpoli (cons nil (copy-tree orbit)))
 155          (multi_orbit (cdr lpoli) (nconc orbit lpoli+1) llvar)))))
 156 ; Que fait epure? He bien il enleve physiquement de (cdr l1) tout
 157 ; les polynomes se trouvant eventuellement dans (cdr l2) en les diminuant
 158 ; toutes deux physiquement.
 159
 160 (defun epure (l1 l2)
 161   (and (cdr l1)
 162        (cond
 163             ((catch 'trouve
 164                           (dans l2 (cadr l1)))
 165                              ; car on calcul la difference
 166                              ; on l'a retire de l2 (ne reviendra pas)
 167           (epure (rplacd l1 (cddr l1)) l2)) ; allez! oust!
 168                   ; l2 diminue' physiquement egalement
 169          (t (epure (cdr l1) l2)))))
 170
 171 ; on regarde si l'oppose de $-pol est dans l2
 172 ; si oui on le retire de l2 et on repond : t sinon : nil
 173
 174 (defun dans (l2 $pol)
 175   (and (cdr l2)
 176        (cond
 177          ((meval (list '($is) (list '(mequal) (cadr l2) $pol)))
 178           (rplacd l2 (cddr l2))
 179              ; on en profite pour le retirer de l2
 180             (throw 'trouve t))
 181          (t (dans (cdr l2) $pol)))))
 182
 183 ;=========================================================================
 184 ;        REMARQUE SUR CE QUI PRECEDE
 185 ;=========================================================================
 186 ; On peut se demander : Pourquoi ne pas lire une seule fois
 187 ; le polynome en le mettant sous la forme d'une
 188 ; liste (c m1 m2 ... mp) representant cm1m2...mp ? ou c est
 189 ; un element de k et chaque mi un monome de k[X^i].
 190 ; Ceci n'a pas ete fait pour 3 raisons
 191 ; la premiere etant que la donnee d'entree (le polynome P) est
 192 ; forcement petite (sinon les calculs satureront par la suite)
 193 ; et que le calcul de sa multi_orbite est negligeable devant
 194 ; ce qui l'attend. Alors au vu des difficultes mises en evidence
 195 ; par les deux autres raisons on se dit que ce n'est vraiment
 196 ; pas la peine.
 197 ; la seconde est qu'on est amene a comparer l'egalite des polynomes
 198 ; a chaque etape i de multi_orbit. Et que meme si les monomes
 199 ; de k[X^1,...,X^(i-1)] sont mis en coefficients comment fait-on
 200 ; pour ceux dependant des X^q (q > i)?
 201 ; La troisieme est que le coefficient lie a un monome en X^i est
 202 ; en fait un polynome en les autres groupe de variables et qu'il
 203 ; faudra bien les reunir d'une facon ou d'une autre.
 204 ; Apres maintes considerations j'ai opte pour la version decrite
 205 ; precedemment qui oblige a repasser le lecteur sur le polynome et ses
 206 ; permute's a chaque fois que l'on veut calculer son orbite sous S_di.
 207 ;===========================================================================
 208 ;                CALCUL DES FONCTIONS PUISSANCES
 209 ;          SOUS FORME CONTRACTE SOUS S_d1 x ... x S_dp = S_D
 210 ; SOIT O = {f1,f2, ...,fo} des polynomes en X^1, en X^2, ... et
 211 ; en X^p. On cherche les fonctions puissances P_r(O) (r= 1..o)
 212 ; sur O mais dans une forme contracte sous
 213 ; S_D (O etant bien constitue pour que cela soit possible).
 214 ;(ie. ne prendre qu'un monome par orbite)
 215 ;-----------------------------------------------------------------------------
 216 ;                       EXEMPLE
 217 ; P = a*x + b*y
 218 ;  X^1=(x,y) elementaires : e1, e2, puissances : p1, p2
 219 ;  X^2=(a,b) elementaires : f1, f2, puissances : g1, g2.
 220 ; O = (ax + by , ay + bx)
 221 ; P_1(O) = ax + by + ay + bx = (a + b)(x +y)
 222 ;         forme contracte : ax
 223 ;         P_1(O) = e1*f1 = p1*g1
 224 ; P_2(O) = (ax + by)^2 + (ay + bx)^2
 225 ;        = a^2x^2 + b^2y^2 + a^2y^2 + b^2y^2 + 2(axby + aybx)
 226 ;        = (a^2 + b^y^2)(x^2 + y^2) + 4axby
 227 ;         forme contracte : a^2x^2 + 4axby
 228 ;         P_1(O) = (e1^2 - 2e2)(f1^2 - 2f2) + 4e2f2
 229 ;                = p2g2 + (p1^2 - p2)(g1^2 -g2)
 230 ;-----------------------------------------------------------------------------
 231 ;                    CONTRAINTE
 232 ; SE DEBARASSER SYSTEMATIQUEMENT DE TOUT MONOMES SI ON EN A DEJA
 233 ; UN DANS SA MULTI_ORBITE AFIN D'EVITER AU MIEUX L'EXPLOSION EN ESPACE.
 234 ; CE QUI EXPLIQUE EN PARTIE POURQUOI ON PREFERE LES FONCTIONS PUISSANCES
 235 ; AUX FONCTIONS SYMETRIQUES ELEMENTAIRES SUR O.
 236 ; On ne garde que les monomes representes par des multipartitions.
 237 ; Remarque : il serait plus efficace d'utiliser le logiciel de
 238 ; Jean-Charles Faugere.
 239 ;-----------------------------------------------------------------------------
 240 ; 1_ l'appel et la boucle principale
 241 ; on retire le degre en tete
 242 (defun $pui_direct_init ($or $llvar $ldegre)
 243   (cons '(mlist)
 244         (cdr (pui_direct (cdr $or)
 245                          (mapcar 'cdr (cdr $llvar))
 246                          (cdr $ldegre)))))
 247
 248 (defun pui_direct (or llvar ldegre)
 249   (let* (
 250         (ldegre_arite (mapcar 'list
 251                              ldegre
 252                              (mapcar 'list-length llvar)))
 253         (degre_resol (* (list-length or) ;le degre de P_*(P1,...,Pp)
 254                         (apply '* (mapcar #'(lambda (nb) (apply
 255                                                         'binomial nb))
 256                                            ldegre_arite)))))
 257     (do ((o (and (print degre_resol) (1- degre_resol))
 258             (and (print o) (1- o)))
 259          (listpui (list  (pui_contract 0 or llvar degre_resol ldegre_arite))
 260                   (cons  (pui_contract 0 or llvar o ldegre_arite) listpui)))
 261         ((eql 0 o) (cons degre_resol listpui)))))
 262
 263 ; 2_ Obtention de la rieme fonction puissance
 264 ; dans Or on a des polynomes macsyma
 265 ; dans $pui_contract des polynomes sous formes contractees
 266 ; on ne conserve que les monomes dont les exposants correspondent a des
 267 ; multipartitions
 268 ; Ramene un polynome macsyma
 269 ;-----------------------------------------------------------------------
 270 (defun pui_contract ($pui_cont or llvar r ldegre_arite)
 271   (cond
 272     ((null or) $pui_cont)
 273     (t (pui_contract
 274            ($add_sym (multi_partitions ($exp_sym (car or) r)
 275                                        llvar
 276                                        ldegre_arite)
 277                $pui_cont)
 278            (cdr or) llvar r ldegre_arite))))
 279
 280 ; on jette les momones a exposants non multi_partitionne dans $pol.
 281 ; map applique a toute les sous-listes et rend son deuxieme arguments
 282 ; ie. la premiere liste.
 283
 284 (defun multi_partitions ($pol llvar ldegre_arite)
 285   (do ((rllvar (cdr (reverse llvar)) (cdr rllvar))
 286        (rldegre_arite (cdr (reverse ldegre_arite)) (cdr rldegre_arite))
 287        (pol_multipartitionne
 288            (garde_que_partition_init (cons nil
 289                                            (lect $pol
 290                                                 (cons '(mlist) (car (last
 291                                                                      llvar)))))
 292                                      (car (last ldegre_arite)))))
 293       ((null rllvar)
 294        (if pol_multipartitionne
 295            (multi_ecrit pol_multipartitionne llvar) 0))
 296     (setq pol_multipartitionne
 297           (apply 'nconc
 298                  (mapl #'(lambda (p_m)
 299                           (rplaca p_m
 300                                   (distribu (cdar p_m)
 301                                       (garde_que_partition
 302                                         (cons nil (lect (caar p_m)
 303                                                         (cons '(mlist)
 304                                                             (car rllvar))))
 305                                          (car rldegre_arite)))))
 306                        pol_multipartitionne)))))
 307
 308 ; le coefficient binomial permet de tenir compte de l'arite'
 309 ; de la fonction resolvante.
 310
 311 (defun garde_que_partition_init ($pol_dist degre_arite)
 312   (do ((pol $pol_dist) (degre (car degre_arite)) (arite (cadr degre_arite)))
 313       ((null (cdr pol)) (cdr  $pol_dist))
 314     (let ((exposants (cdadr pol)))
 315          (cond
 316               ((apply #'>= exposants)
 317                (setq pol (cdr pol))
 318                (let ((lg (longueur exposants)))
 319                      (rplaca pol (list ($mult_sym (caar pol)
 320                                               (binomial (- degre lg)
 321                                                         (- arite lg)))
 322                                        exposants))))
 323       (t (rplacd pol (cddr pol)))))))
 324
 325 ; remarque en common lisp : (>= 4 3 2 1 0) ==> true permet de tester
 326 ; si une liste est une partition
 327
 328 (defun garde_que_partition ($pol_dist degre_arite)
 329   (do ((pol $pol_dist)
 330        (degre (car degre_arite))
 331        (arite (cadr degre_arite)))
 332       ((null (cdr pol)) (cdr $pol_dist))
 333       (let ((exposants (cdadr pol)))
 334            (cond
 335                 ((apply '>= exposants)
 336                  (setq pol (cdr pol))
 337                  (let ((lg (longueur exposants)))
 338                       (rplaca (car pol)
 339                               ($mult_sym (caar pol)
 340                                          (binomial (- degre lg)
 341                                                    (- arite lg))))))
 342                (t (rplacd pol (cddr pol)))))))
 343
 344 (defun distribu (multipartition pol_part)
 345   (mapcar #'(lambda (coe_partition)
 346              (cons (car coe_partition)
 347                    (cons (cdr coe_partition) multipartition)))
 348           pol_part))
 349 ;=========================================================================
 350 ;                BOUCLE PINCIPALE DE L'IMAGE DIRECTE
 351 ;=========================================================================
 352 (defun $direct_init ($lpol $x $p $llvar)
 353   (direct (cdr $lpol) $x $p  $llvar))
 354
 355 (defun direct (l$pol $x $p $llvar)
 356   (cond ((equal '$parallele $directnew)
 357          (direct-parallele l$pol $x $p $llvar))
 358         (t
 359            (let* (
 360                   ($multi_base (if (equal '$elementaire $direct)
 361                                    (macsy_list (multi_polynome2ele l$pol $x))
 362                                    (multi_ele2pui (multi_polynome2ele l$pol $x))))
 363                   (pui_* (pui_direct  (multi_orbit_init $p (cdr $llvar))
 364                                       (mapcar 'cdr (cdr $llvar))
 365                                       (mapcar 'cadr (cdr $multi_base) )))
 366                   (degre (car pui_*)))
 367                  (pui2polynome '$y
 368                                (cons degre
 369                                      (mapcar #'(lambda ($pi)
 370                                                    (multi_decomp_init
 371                                                              $pi
 372                                                              $multi_base
 373                                                              $llvar ))
 374                                              (cdr pui_*))))))))
 375
 376 ; Ici on calcule les fonctions puissances des racines de la resolvante
 377 ; au fur et a mesure. Avant nous calculions les fonctions puissances
 378 ; des racines de la resolvante generique sur la base des formes
 379 ; monomiales et ensuite on specialisait. En fait nous n'exploitions
 380 ; pas l'aspect parallele du calcul.
 381
 382 (defun direct-parallele (l$pol $x $p $llvar)
 383 (let* (
 384        (llvar (mapcar 'cdr (cdr $llvar)))
 385          ($multi_base (if (equal '$elementaire $direct)
 386                          (macsy_list (multi_polynome2ele l$pol $x))
 387                          (multi_ele2pui (multi_polynome2ele l$pol $x))))
 388        (multi_orbite (multi_orbit_init $p (cdr $llvar)))
 389        (ldegre (mapcar 'cadr (cdr $multi_base) )) ; degres des polynomes
 390        (ldegre_arite (mapcar 'list
 391                              ldegre
 392                              (mapcar 'list-length llvar)))
 393        (degre_resol (* (list-length multi_orbite) ;le degre de P_*(P1,...,Pp)
 394                         (apply '* (mapcar #'(lambda (nb) (apply
 395                                                         'binomial nb))
 396                                            ldegre_arite)))))
 397            (do ((r (and (print degre_resol) (1- degre_resol))
 398                    (and (print r) (1- r)))
 399             (listpui (list (multi_decomp_init (pui_contract 0
 400                                                     multi_orbite
 401                                                     llvar
 402                                                     degre_resol
 403                                                     ldegre_arite)
 404                                             $multi_base
 405                                             $llvar ) )
 406                     (cons  (multi_decomp_init  (pui_contract 0
 407                                                     multi_orbite
 408                                                     llvar
 409                                                     r
 410                                                     ldegre_arite)
 411                                             $multi_base
 412                                             $llvar)
 413                              listpui)))
 414         ((eql 0 r)
 415          (pui2polynome '$y  (cons degre_resol listpui))))))
 416
 417 ;====================================================================
 418 ;  FONCTIONS SYMETRIQUES ELEMENTAIRES DES RACINES DES POLYNOMES DE
 419 ; l$pol EN LA VARIABLE $y.
 420
 421 (defun multi_polynome2ele (l$pol $y)
 422   (mapcar #'(lambda ($pol) (polynome2ele $pol $y)) l$pol))
 423
 424 ;=========================================================================
 425 ;          DECOMPOSITION D'UN POLYNOME SYMETRIQUE CONTRACTE
 426 ;             EN PLUSIEURS PAQUETS DE VARIABLES
 427 ;           DONT LES FONCTIONS SYMETRIQUES ELEMENTAIRES
 428 ;               RESPECTIVES SONT DONNEES
 429 ;=========================================================================
 430
 431 (defun multi_decomp_init ( $multi_pc $multi_base $llvar)
 432   (cond
 433     ((equal '$elementaires $direct)
 434      (meval (list '($multi_elem)  $multi_base
 435                                   $multi_pc $llvar)))
 436     (t (meval (list '($multi_pui) $multi_base
 437                                    $multi_pc $llvar)))))
 438
 439 ; on a les fonctions symetriques elementaires des racines des differents
 440 ; polynomes. On recupere en fonction d'elles leurs fonctions puissances.
 441
 442 (defun multi_ele2pui (multi_lelem)
 443   (cons '(mlist)
 444        (mapcar #'(lambda (lelem)
 445                      (meval (list '($ele2pui)
 446                                    (car lelem)
 447                                    (cons '(mlist) lelem))))
 448                 multi_lelem)))
 449
 450
 451