1;; Fichier util.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(in-package :maxima) 20(macsyma-module util macros) 21;--------------------------------------------------------------------------- 22; DECLARATION DES MACROS 23; pour le type 2 des polynomes partitionnes avec en tete de chaque 24; terme partitionne sa longueur 25;--------------------------------------------------------------------------- 26 27;---------------------------------------------------------------------------- 28; LES UTILITAIRES 29;---------------------------------------------------------------------------- 30; On a des coefficients dans k[y1, ...,yn] 31; p=(t1 t2 ... tn) t1 > t2 > ... 32; t = (longueur coe . partition) 33 34(progn (defvar d) (defvar lvar) (defvar permut)) 35 36(progn) 37(progn) 38 39 40(defun $estpartition (l) 41 (apply '>= (cdr l))) 42 43;-------------------------------------------------------------------------- 44; MERGE PHYSIQUE AVEC SOMME SUR LES GRANDS ENTIERS 45; ON UTILISE LE PREDICAT SUR LES TERMES ET NON SUR LES MONOMES(comme merge) 46;--------------------------------------------------------------------------- 47(defun somme (l1 l2 pr) 48 (cond 49 ((null l1) l2) 50 ((null l2) l1) 51 (t (let ((t1 (termi l1)) (t2 (termi l2))) 52 (cond 53 ((equal (tmon t1) (tmon t2)) 54 (chcoeterm t1 ($add_sym (tcoe t1) (tcoe t2))) 55 (cond 56 ((and (numberp (tcoe t1)) 57 (zerop (tcoe t1))) 58 (somme (cdr l1) (cdr l2) pr)) 59 (t 60 (somme2 (cdr l2) l1 pr) l1))) 61 ((funcall pr t1 t2) (somme2 l2 l1 pr) l1) 62 (t (somme2 l1 l2 pr) l2)))))) 63 64(defun somme2 (l2 l1 pr) 65 (do ((l1 l1) (l2 l2) (ll2 nil) (t1 (termi (cdr l1)) (termi (cdr l1))) 66 (t2 (termi l2))) 67 ((or (null l2) (and (null (cdr l1)) (nconc l1 l2)))) 68 (cond 69 ((equal (tmon t1) (tmon t2)) 70 (chcoeterm t1 ($add_sym (tcoe t1) (tcoe t2))) (setq l2 (cdr l2)) 71 (setq t2 (termi l2)) 72 (if (and (numberp (tcoe t1)) (zerop (tcoe t1))) 73 (setq l1 (rplacd l1 (cddr l1))) 74 (setq l1 (cdr l1)))) 75 ((funcall pr t1 t2) (setq l1 (cdr l1))) 76 (t (setq ll2 (cdr l1)) (setq l1 (cdr (rplacd l1 l2))) 77 (setq l2 ll2) (setq t2 (termi l2)))))) 78 79;====================================================================== 80; CREATION D'UNE LISTE LISP DE nb VARIABLES GENERIQUES : 81; ($x1 ... $x(nb)) 82; (lvar 2 '(r)) = ; ($X1 $X2 R) 83;====================================================================== 84 85(defun lvar (nb lvar) 86 (cond 87 ((eql 0 nb) lvar) 88 (t (lvar (1- nb) 89 (cons (flet ((franz.concat (&rest args) 90 "equivalent to Franz Lisp 'concat'." 91 (values (intern 92 (format nil "~{~A~}" args))))) 93 (franz.concat '$x nb)) 94 lvar))))) 95 96;(lvar_lettre 2 '(r) 'x) 97; (X1 X2 R) 98 99(defun lvar_lettre (nb lvar lettre) 100 (cond 101 ((eql 0 nb) lvar) 102 (t (lvar_lettre (1- nb) 103 (cons (flet ((franz.concat (&rest args) 104 "equivalent to Franz Lisp 'concat'." 105 (values (intern (format nil "~{~A~}" args))))) 106 (franz.concat lettre nb)) 107 lvar) 108 lettre)))) 109 110;=========================================================================== 111; Calcul du degre d'un polynome symetrique 112; avec REP([pol]) = [lppart](2) 113 114(defun $degrep (pol) 115 (setq d 0) 116 (mapc #'(lambda (di) 117 (and (< d di) 118 (setq d di))) 119 (mapcar #'(lambda (mon) ($degre (cddr mon))) pol )) 120 d) 121 122; Calcul du degre d'une forme monomiale avec REP([forme mon])=[partition](2) 123; mon = (lgI coeI . I) 124 125(defun $degre (mon) 126 (if (or (constantp mon) (null mon)) 0 127 (+ (* (car mon) (cadr mon)) 128 ($degre (cddr mon))))) 129 130;--------------------------------------------------------------------------- 131; TESTE SI ON A AFFAIRE A UNE CONSTANTE APRES LE LECTEUR 132; termpart = REP([somme orbitale]) 133 134; avec [somme orbitale] = (coe.[partition]) 135 136(defun constante (termpart) 137 (or (null (cdr termpart)) 138 (eval (cons 'and 139 (mapcar #'(lambda (exposant) (eql 0 exposant)) 140 (cdr termpart)))))) 141 142; avec [somme orbitale] = (longueur coe.[partition]) 143; il suffit de tester si la longueur est nulle 144 145(defun lconstante (ltermpart) (eql 0 (car ltermpart))) 146 147; Calcul des longueurs de chaque partition contenue dans la liste listparts 148; sous forme [partition](2) 149 150(defun lgparts (ppart) 151 (mapcar #'(lambda (mon) (cons ($calculvar (cdr mon)) mon)) ppart)) 152 153; Calcul de la longueur d'une partition I. 154; Pour [partition](1) 155 156(defun longueur (i) 157 (if (or (null i) (eql 0 (car i))) 0 158 (1+ (longueur (cdr i))))) 159 160; Pour [partition](2), pouvant se terminer par des 0. 161 162(defun $calculvar (i) 163 (if (or (null i) (eql 0 (car i))) 0 164 (+ (cadr i) ($calculvar (cddr i))))) 165 166;************************************************************************** 167; PREDICATS 168;------------------------------------------------------------------------- 169;term est un [partition](2) c'est a dire sous la forme : 170; dans la forme (a1 m1 a2 m2...) ou mi est la multiplicite' de ai (> a(i+1)) 171;------------------------------------------------------------------------- 172;(2 1 ...) < (2 2 ...) < (3 1 2 1 ...) < (3 1 2 2 ...) 173 174(defun $lex (term1 term2) 175 (cond 176 ((null term1) t) 177 ((null term2) nil) 178 (t (let ((pui1 (car term1)) (nb1 (cadr term1)) (rest1 (cddr term1)) 179 (pui2 (car term2)) (nb2 (cadr term2)) 180 (rest2 (cddr term2))) 181 (cond 182 ((or (< pui1 pui2) 183 (and (eql pui1 pui2) 184 (< nb1 nb2))) 185 t) 186 ((or (< pui2 pui1) 187 (< nb2 nb1)) 188 nil) 189 (t ($lex rest1 rest2))))))) 190 191; q inferieur a p pour l'ordre des longueurs ou p et q sont des 192; sommes orbitales represente'es par des [terme partionne](2) avec 193; la longueur en plus en tete 194 195(defun orlongsup (p q) 196 (cond 197 ((equal (cddr p) (cddr q)) nil) 198 ((> (car p) (car q))) 199 ((eql (car p) (car q)) ($lex (cddr q) (cddr p))) 200 (t nil))) 201;---------------------------------------------------- 202; le vrai ordre des longueurs q inferieur a p pour cet ordre : ; p { q 203 204;(orlongsup '(2 a 2 1 3 1) '(1 a 4 1)) 205;T 206 207;>(orlong '(2 a 2 1 3 1) '(1 a 4 1)) 208;NIL 209;---------------------------------------------------- 210 211(defun orlong (p q) 212 (cond 213 ((equal (cddr p) (cddr q)) nil) 214 ((< (car p) (car q))) 215 ((eql (car p) (car q)) ($lex (cddr p) (cddr q))) 216 (t nil))) 217 218(defun $orlong_cst (p q) 219 (cond 220 ((lconstante p)) 221 ((lconstante q) nil) 222 (t (orlong p q)))) 223 224; p > q 225 226(defun $e_lexinv_cst (mon1 mon2) 227 (cond 228 ((constante mon1)) 229 ((constante mon2) nil) 230 (t ($e_lexinv mon1 mon2)))) 231 232; p=(lg(I) coeI .I) 233 234(defun $e_lexinv (p q) 235 (and (not (equal (cddr p) (cddr q))) ($lex (cddr q) (cddr p)))) 236 237; p < q 238; les constantes sont les + petites 239 240(defun $e_lex_cst (p q) 241 (cond 242 ((constante p)) 243 ((constante q) nil) 244 (t ($e_lex p q)))) 245 246(defun $e_lex (p q) 247 (and (not (equal (cddr p) (cddr q))) ($lex (cddr p) (cddr q)))) 248 249; teste sur deux monomes en representation distribuee (i1 i2 ...) 250 251(defun lex_term (term1 term2) 252 (lex_mon (cdr term1) (cdr term2))) 253 254(defun lex_mon (m1 m2) 255 (and (not (equal m1 m2)) 256 (catch 'trouve 257 (mapc #'(lambda (e1 e2) 258 259 (or (eql e1 e2) 260 (cond 261 ((> e1 e2) 262 (throw 'trouve t)) 263 (t (throw 'trouve nil))))) 264 m1 m2)))) 265 266;*************************************************************************** 267; INTERFACE 268 269 270; Le lecteur utilise la fonction $distri_lect qui appelle distri_lecteur 271 272(defun lect ($pol $lvar) 273 (mapcar 'cdr 274 (cdr (meval (list '($distri_lect) $pol $lvar))))) 275 276;-------------------------------------------------------------------------- 277; [ppart](i) lisp ==> [ppart](i) macsyma 278 279(defun macsy_list (llist) 280 (cons '(mlist) (mapcar #'(lambda (list) (cons '(mlist) list)) llist))) 281 282; sa recipropque : 283 284(defun mac2lisp (list) (mapcar 'cdr (cdr list))) 285;-------------------------------------------------------------------------- 286; [ppart](i) == > polynome macsyma 287 288; Si REP([pol]) =[ppart](1) 289; Pour une liste de polynomes. 290; Mais attention! Si on veut l'utiliser sous Macsyma, il faut 291; rajouter (MLIST SIMP) en car de la liste resultat. 292;-------------------------------------------------------------------------- 293 294(defun ecrit_listpol (listpol lvar) 295 (mapcar #'(lambda (pol) (ecrit_pol pol lvar)) listpol)) 296 297;-------------------------------------------------------------------------- 298; Pour un polynome de plusieurs groupes de variables 299; en representation distribuee : 300; (c m1 m2 ... mp) ou mi est un monome en X^(i) . 301; Par exemple : mi=(3 4 1) represente U^3*V^4*W si X^(i)=(U,V,W). 302; llvar = (X^(1), ..., X^(p)) est une liste de listes de variables 303;-------------------------------------------------------------------------- 304 305(defun multi_ecrit (pol llvar) 306 (cond 307 ((null (cdr pol)) 308 (multi_ecrit_mon (caar pol) (cdr (car pol)) llvar)) 309 (t ($fadd_sym 310 (mapcar #'(lambda (terme) 311 (multi_ecrit_mon (car terme) (cdr terme) 312 llvar)) 313 pol))))) 314 315(defun multi_ecrit_mon (coe llexposants llvar) 316 (cond 317 ((null llvar) coe) 318 (t (mapc #'(lambda (lexposants lvar) 319 (setq coe (ecrit_mon lexposants lvar coe))) 320 llexposants llvar) 321 coe))) 322;-------------------------------------------------------------------------- 323; Pour un polynome a un groupe de variables dont le representaion 324; partitionnee est de type 1, on considere que le polynome 325; est sous forme distribue'e et on se sert de l'ecrivain de polynomes 326; destine' a ce cas afin d'obtenir un polyn\^ome maxima. 327; la fonction au niveau maxima est $distri_ecrit. Mais c'est ennnuyeux 328; de mettre de mlist pour les retirer ensuite, alors j'appelle 329; directement 330; la fonction interne : $ecrivain_sym 331;-------------------------------------------------------------------------- 332 333(defun ecrit_pol (ppart lvar) 334 (meval (list '($distri_ecrit) 335 (macsy_list ppart) (cons '(mlist) lvar)))) 336 337;-------------------------------------------------------------------------- 338; Si REP([pol]) = [ppart](2) on se ramene au cas precedent 339;-------------------------------------------------------------------------- 340 341(defun 2ecrit (ppart lvar) (ecrit_pol (ch1repol ppart) lvar)) 342 343;************************************************************************** 344; CHANGEMENTS DES REPRESENTATIONS DE PARTITIONS 345;----------------------------------------------------------------------- 346; Fonction passant de [partition](1) a [partition](2) 347 348(defun ch2rep (partition1) 349 (and partition1 (not (eql 0 (car partition1))) 350 (ch2rep2 (cdr partition1) (list 1 (car partition1))))) 351 352(defun ch2rep2 (partition1 partition2) 353 (if (or (null partition1) (eql 0 (car partition1))) 354 (nreverse partition2) 355 (if (eql (car partition1) (cadr partition2)) 356 (ch2rep2 (cdr partition1) 357 (rplaca partition2 358 (1+ (car partition2)))) 359 (ch2rep2 (cdr partition1) 360 (cons 1 (cons (car partition1) partition2)))))) 361 362; Passer d'un polynome partitionne avec [partition](1) a [partition](2) 363 364(defun ch2repol (ppart) 365 (mapcar #'(lambda (tpart) (cons (car tpart) (ch2rep (cdr tpart)))) 366 ppart)) 367 368; PASSAGE DE [partition](2) a [partition](1) 369 370(defun ch1rep (partition2) 371 (and partition2 372 (ch1rep2 (cddr partition2) 373 (make-list (cadr partition2) 374 :initial-element (car partition2))))) 375 376(defun ch1rep2 (partition2 partition1) 377 (cond 378 ((null partition2) (nreverse partition1)) 379 (t (ch1rep2 (cddr partition2) 380 (nconc (make-list (cadr partition2) 381 :initial-element (car partition2)) 382 partition1))))) 383 384; Maintenant passer d'un polynome partitionne avec [partition](2) 385; a un polynome partitionne avec [partition](1) 386 387(defun ch1repol (ppart) 388 (mapcar #'(lambda (tpart) (cons (car tpart) (ch1rep (cdr tpart)))) 389 ppart)) 390 391;------- Seconde methode 392; Passage de la premiere represensentation des partitions a la seconde 393; on ramene les cst sans partition associee 394 395; listM =(i1 i2 i3 ...ip) avec 0< i1 <= i2 <= ... <=ip 396 397(defun part (listm) 398 ($part2 (cdr listm) (cons (car listm) (cons 1 nil)))) 399 400(defun $part2 (listm lpartm) 401 (if (null listm) lpartm 402 ($part2 (cdr listm) 403 (if (eql (car listm) (car lpartm)) 404 (progn 405 (rplaca (cdr lpartm) 406 (1+ (cadr lpartm))) 407 lpartm) 408 (cons (car listm) (cons 1 lpartm)))))) 409 410(defun $cherchepui (mmon) (if (atom mmon) 1 (car (last mmon)))) 411 412;========================================================================= 413; CALCUL DU CARDINAL DE L' ORBITE D'UN MONOME 414; DONT LA LISTE DES EXPOSANTS EST DONNE PAR UNE 415; PARTITION DONT LA REPRESENTATION EST [partition](2)=(a1 m1 a2 m2...) 416; qui est n!/(m0!m1!...) ou n=somme_{i=0} mi 417; Ou du coefficient multinomial associe a une partition type [partition](1) 418; qui est |I|!/(i1!i2!....) 419;--------------------------------------------------------------------------- 420(defun $card_orbit ($partition $card) 421 (card_orbit (cdr $partition) $card)) 422 423(defun card_orbit (partition card) 424 (nbperm0 card 425 (- card ($calculvar partition)); le nombre a la puissance 0 426 partition)) 427 428(defun $multinomial (poids $partition) 429 (multinomial (nbperm2 1 poids 1 (cadr $partition)) (cddr $partition))) 430 431(defun multinomial (prec_multinomial partition) 432 (cond 433 ((or (null partition) (= 0 (car partition))) (car prec_multinomial)) 434 (t (multinomial 435 (nbperm2 (car prec_multinomial) (cadr prec_multinomial) 1 436 (car partition)) 437 (cdr partition))))) 438 439(defun nbperm0 (card m0 part) (nbperm (nbperm2 1 card 1 m0) part)) 440 441(defun nbperm (lpermn part) 442 (if (null part) (car lpermn) 443 (nbperm (nbperm2 (car lpermn) (cadr lpermn) 1 (cadr part)) 444 (cddr part)))) 445 446(defun nbperm2 (perm n i mi) 447 (if (< mi i) 448 (list perm n) 449 (nbperm2 (/ (mult perm n) i) 450 (1- n) 451 (1+ i) 452 mi))) 453 454;------------------------------------------------------------------------ 455; Calcul du cardinal du stabilisateur d'une liste ordonnee 456; de paires ou de nombres dans l'ordre lexicographique croissant. 457 458(defmfun $card_stab ($part $egal) ($card_stab_init $part $egal)) 459 460(mdefprop $card_stab 461 ((lambda ()) ((mlist) $part $egal) 462 ((mprog) (($operation)) (($card_stab_init) $part $egal))) 463 mexpr) 464(add2lnc '(($card_stab) $part $egal) $functions) 465;------------------------------------------------------------------------ 466 467(defun $card_stab_init ($part $egal) 468 (card_stab (cdr $part) 469 (find-symbol (string $egal)))) 470 471(defun card_stab (s egal) 472 (let ((lmultip (sort (multiplicites s egal) '<))) 473 (prod_factor lmultip))) 474 475(defun multiplicites (s egal) 476 (multiplicites2 (cdr s) (car s) 1 nil egal)) 477 478(defun multiplicites2 (s ai mi lmultip egal) 479 (cond 480 ((null s) (cons mi lmultip)) 481 ((funcall egal (car s) ai) 482 (multiplicites2 (cdr s) ai 483 (1+ mi) 484 lmultip egal)) 485 (t (multiplicites2 (cdr s) (car s) 1 (cons mi lmultip) egal)))) 486 487; l = (m1 m2 ... mp) croissante , on veut m1!m2!...mp! 488 489(defun prod_factor (l) 490 (apply '* (list_factor l (list (factorielle (car l)))))) 491; l = (mi m(i+1) ... mp) et lfactor = (mi!, ..., m2!, m1!) 492 493(defun list_factor (l lfactor) 494 (cond 495 ((null (cdr l)) lfactor) 496 (t (list_factor (cdr l) 497 (cons (fact_recur (car l) (cadr l) (car lfactor)) lfactor))))) 498 499(defun fact_recur (m1 m2 factm1) 500 (cond 501 ((eql m1 m2) factm1) 502 (t (* (finfact (1+ m1) 503 m2 504 (1+ m1)) 505 factm1)))) 506 507(defun finfact (i m2 finfactm2) 508 (cond 509 ((eql i m2) finfactm2) 510 (t (finfact (1+ i) 511 m2 512 (* (1+ i) 513 finfactm2))))) 514 515(defun factorielle (n) 516 (cond 517 ((eql 0 n) 1) 518 (t (* n (factorielle (1- n)))))) 519 520;________________________________________________________________________ 521 522; OBTENIR TOUTES LES PERMUTATIONS D'UN NUPLET D'ENTIER 523; (Philippe Esperet) remarque : 524; VOIR FONCTIONS permutations et permutations_lex de MAXIMA 525;--------------------------------------------------------------------------- 526 527(defun $lpermut (nuplet) 528 (cons '(mlist) 529 (mapcar #'(lambda (permu) (cons '(mlist) permu)) 530 (permut (cdr nuplet))))) 531 532;====================================================================== 533; EXPRESSION D'UN POLYNOME 534; DONT ON CONNAIT LES FONCTIONS SYMETRIQUES 535; ELEMENTAIRES DES RACINES 536; $fct_elem =[cardinal, e_1,e_2,...,e_cardinal,...] 537;====================================================================== 538 539(defun $ele2polynome ($fct_elem $z) 540 (ele2polynome (cdr $fct_elem) $z)) 541 542(defun ele2polynome (l_degre_$elem $z) 543 (genpoly2 544 (1- (car l_degre_$elem)) 545 -1 ($exp_sym $z (car l_degre_$elem)) (cdr l_degre_$elem) $z)) 546 547(defun genpoly2 (exp sign $pol l_$elem $z) 548 (cond 549 ((null l_$elem) $pol $pol) 550 (t (genpoly2 551 (1- exp) 552 (* -1 sign) 553 ($add_sym $pol 554 ($mult_sym ($mult_sym sign (car l_$elem)) 555 ($exp_sym $z exp))) 556 (cdr l_$elem) $z)))) 557;========================================================================= 558; OBTENIR UN POLYNOME A PARTIR DES FONCTIONS PUISSANCES 559; DE SES RACINES 560 561; fct_pui = (card p1 p2 ...) 562;========================================================================= 563 564(defun $pui2polynome ($var $fct_pui) 565 (pui2polynome $var (cdr $fct_pui))) 566 567(defun pui2polynome (variable fct_pui) 568 (let (($pui2ele '$girard)) 569 (ele2polynome (cdr (meval (list '($pui2ele) (car fct_pui) 570 (cons '(mlist) fct_pui)))) 571 variable))) 572 573;========================================================================= 574; CALCUL DES FONCTIONS SYMETRIQUES ELEMENTAIRES 575; DES RACINES D'UN POLYNOME 576; entrees : $p un polynome en la variable $var 577; sortie : [d , e1, ...,ed] ou d est le degre du polynome 578;========================================================================== 579(defun $polynome2ele ($p $var) 580 (cons '(mlist) (polynome2ele $p $var))) 581 582(defun polynome2ele ($p $var) 583 (let* ((alt -1) 584 (n (meval (list '($HIPOW) $p $var))) 585 (an (meval (list '($COEFF) $p $var n)))) 586 (do ((alt alt (* -1 alt )) 587 (i 1 (1+ i)) 588 (elem nil (cons ($divi_sym ($mult_sym alt 589 (meval (list '($COEFF) 590 $P $var (- n i)))) 591 an) 592 elem))) 593 ((= (1+ n) i) (cons n (nreverse elem)))))) 594 595 596; Obtenir tout les coefficients, meme les nuls, 597; (cn c(n-1)...c0) ou ci coefficient de x**i. 598(defun lcoe2 (precedexp p lcoe) 599 (if (null p) 600 (or (eql 0 precedexp) 601 (rplacd lcoe (make-list precedexp :initial-element 0))) 602 (let ((exp (car p)) (coe (cadr p))) 603 (if (eql precedexp 604 (1+ exp)) 605 (lcoe2 exp (cddr p) (cdr (rplacd lcoe (list coe)))) 606 (lcoe2 exp (cddr p) 607 (last (rplacd lcoe 608 (append 609 (make-list (- precedexp 610 (1+ exp)) 611 :initial-element 0) 612 (list coe))))))))) 613 614;====================================================================== 615 616(defun binomial (n p) 617 (meval (list '(%binomial) n p))) 618 619;====================================================================== 620 621; la fonction maxote est commune a : treillis.lsp , resolvante.lsp, kak.lsp 622; voir dans util.lsp 623 624; ici difference avec le common-lisp de macsyma : 625; / a la place de /! pour la division 626 627(defun maxote (a b) 628 (and (plusp b) 629 (if (eql 1 b) 0 630 (if (eql 0 (rem a b)) 631 (- a (div a b)) 632 (- a (1+ (/ a b))))))) 633