share/algebra/grob1.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2
   3 (in-package :maxima)
   4
   5 ;               COMPUTATIONS OF GROEBNER BASES
   6 ;written by D.Lazard (jun. 1988)
   7 ;-----------------------------------------------------------
   8 ;This package contains (with french comments)
   9 ; LIRE : a lisp function which convert a list of multivariate polynomials with
  10 ;       integer coefficients from macsyma CRE form to our internal form; the
  11 ;       whole list has to be in CRE form, not only its members.
  12 ;
  13 ; STAND(CRE_list,ordering) : a macsyma function computing the groebner base of
  14 ;       a list of polynomials acceptable by LIRE ; at the present the result
  15 ;       is left in internal form and bounded to the  (macsyma) variables BASE
  16 ;       and ERREXP1. ORDERING may be ORDLEX, the pure lexicographical ordering
  17 ;       or (better) ORDDEG the degree-reverse-lexicographical ordering; the
  18 ;       ordering on the variables, determined by the ordering in the CRE form,
  19 ;       is reversed, when passing from ORDLEX to ORDDEG.
  20 ;   Sample usage:
  21 ;     (c1)  Stand([x^2-x-u,u^2-1],ordlex);
  22 ;
  23 ;       Some statistics are given by macsyma global variables:
  24 ;               NBSYZ: number of computed syzygies
  25 ;               NBRED: number of reduction steps
  26 ;               NBREDUC: number of reductions of polynomials already in the
  27 ;                       base
  28 ;               NBRED0: number of reductions leading to 0
  29 ;       Moreover, statistics, and the leading monomials of the polynomials in
  30 ;       the base are displayed when a syzygy is computed
  31 ;
  32 ; (ELIMPOL BASE ORDERING VAR): A lisp function which, from the result of
  33 ;       STAND (called with the same ordering), returns
  34 ;               1 if the input of STAND, viewed as a system of equations
  35 ;                       has no solution in an algebraic closed field
  36 ;               "dimension positive" if the set of solutions in an algebraic
  37 ;                       closed field is not finite
  38 ;               or a macsyma  univariate polynomial in the variable VAR, the
  39 ;                       roots of which are the values of the first coordinate
  40 ;                       of the solutions; the multiplicity of a root is the
  41 ;                       total multiplicity of the corresponding solutions.
  42 ;       This function is mainly useful if ORDERING is ORDDEG.
  43 ;
  44 ; ELIMPOL_ERR (ORDERING, VAR): a macsyma function which call ELIMPOL on
  45 ;       ERREXP1.
  46 ;-----------------------------------
  47
  48 ;A FAIRE
  49 ;l'ecrivain pour la sortie de grobner
  50 ;ne trier que les paires non visitees
  51 ;diminuer le nombre de divisions (comme dans le cas lexico ?)
  52
  53 ;DECLARATIONS
  54 (declare-top
  55   (special
  56     $nbsyz $nbred $nbreduc $nbred0 ordexp orddeg ordpairs ordlex
  57     $orddeg $ordlex
  58     ))
  59
  60 ;; Keep the header for the last cre polynomial formed during a LIRE.
  61 (defvar $header nil)
  62
  63 ;(defun $stand ($liste ordre)(grobner (lire $liste ordre) ordre))
  64 (defun $stand ($liste &optional (ordre 'ordlex))
  65   (let* ((lis (lire $liste ordre))
  66          (base  (grobner lis ordre)))
  67     (cons '(mlist)
  68           (loop for w in base collect ($decode_poly w nil)))))
  69
  70
  71
  72
  73 ;LE LECTEUR
  74 ;Il faut traduire une liste de polynomes CRE; on suppose ici les
  75 ;coefficients rationnels
  76 ;Il faut avoir appele' RAT sur la liste de polynomes pour toutes les variables
  77 ;apparaissent dans les en-tetes de tous les polynomes
  78
  79 ;Passer de CRE a polynome en la premiere variable
  80 (defun lire1 (pol var)                  ;en-tete deja enlevee
  81         (cond ((numberp pol)  (list (list pol 0)))
  82               ((equal (car pol) var)
  83                  (let ((res (cons nil nil)))
  84                    (do ((p (cdr pol) (cddr p))
  85                         (r res (cdr r)))
  86                        ((null p) (cdr res))
  87                        (rplacd r (cons (list (cadr p) (car p))
  88                                      nil)))))
  89               (t (list (list pol 0)))))
  90
  91 ;passer de CRE a polynome en n variables
  92 (defun lire2 (pol lvar)
  93         (cond ((null (cdr lvar))
  94                 (lire1 pol (car lvar)))
  95               (t
  96                 (mapcan #'(lambda (u)
  97                            (mapcar #'(lambda (v)
  98                                       (append v (cdr u)))
  99                              (lire2 (car u) (cdr lvar))))
 100                   (lire1 pol (car lvar))))))
 101
 102 ;appliquer lire2 a une liste macsyma de polynomes
 103 (defun lire (mliste ordexp)
 104   (or ($ratp mliste) (setq mliste ($rat mliste)))
 105   (setq $header (car (second mliste)))
 106   (or (eq (car $header) 'mrat) (error "bad rat form"))
 107         (mapcar #'(lambda (u)
 108                    (sort
 109                      (lire2 (cadr u)
 110                             (reverse (cadddr (caadr mliste))))
 111                      #'(lambda (v w) (ordexp (cdr v) (cdr w)))))
 112                 (cdr mliste)))
 113
 114 ;LES ORDRES
 115 ;ordre lexico inverse pour les exposants
 116 (defun orddeg (m1 m2)
 117         (let ((d (- (apply '+ m1) (apply '+ m2))))
 118             (cond ((plusp d) t)
 119                   ((minusp d) nil)
 120                   (t (do ((m1 (cdr m1) (cdr m1))
 121                           (m2 (cdr m2) (cdr m2))
 122                           (a  (- (car m1) (car m2)) (- (car m1) (car m2))))
 123                        ((not (and (zerop a) m1)) (minusp a)))))))
 124 (setq orddeg 'orddeg)
 125 (setq $orddeg 'orddeg)
 126
 127 (defun ordlex (m1 m2)
 128         (do ((m1 (cdr m1) (cdr m1))
 129              (m2 (cdr m2) (cdr m2))
 130              (a (- (car m1) (car m2)) (- (car m1) (car m2))))
 131             ((not (and (zerop a) m1)) (plusp a))))
 132
 133 (setq ordlex 'ordlex)
 134 (setq $ordlex 'ordlex)
 135
 136 (defun ordexp (m1 m2)           ;ancienne version
 137         (let ((d (- (apply '+ m1) (apply '+ m2))))
 138             (cond ((plusp d) t)
 139                   ((minusp d) nil)
 140                   (t (do ((m1 m1 (cdr m1))
 141                           (m2 m2 (cdr m2)))
 142                        ((not (and m1 (= (car m1) (car m2))))
 143                          (and m1 (< (car m1) (car m2)))))))))
 144
 145 (defun ordpairs (p1 p2)         ;nouvelle version moins bonne, mais
 146         (let ((d1 (apply '+ (caar p1)))         ; ? meilleure pour lexico
 147               (d2 (apply '+ (caar p2))))
 148            (or (< d1 d2)
 149                (and (= d1 d2)
 150                     (or (funcall ordexp (caadr p2) (caadr p1))
 151                         (and (equal (caadr p1) (caadr p2))
 152                              (funcall ordexp (cdadr p2) (cdadr p1))))))))
 153
 154 (defun ordpairs (p1 p2)         ;ancienne version
 155         (let ((exp1 (caar p1))
 156               (exp2 (caar p2)))
 157            (or (funcall ordexp exp2 exp1)
 158                (and (equal exp1 exp2)
 159                     (or (funcall ordexp (caadr p2) (caadr p1))
 160                         (and (equal (caadr p1) (caadr p2))
 161                              (funcall ordexp (cdadr p2) (cdadr p1))))))))
 162
 163 (setq ordpairs 'ordpairs)
 164
 165 ;OPERATIONS ELEMENTAIRES PRESERVATIVES (POLYNOMES) COEFFICIENTS GENERAUX
 166
 167 ;definition des multiplications
 168 ;mactimes peut etre change si les coefficients sont restreints
 169 ;par exemple, pour des entiers, on peut prendre "times"
 170
 171 (defun *gpolmon (pol mon)
 172         (mapcar #'(lambda (u)
 173                    (*gterm u mon))
 174             pol))
 175
 176 (defun *gterm (m1 m2)
 177         (cons
 178           (mactimes (car m1) (car m2))
 179           (mapcar '+ (cdr m1) (cdr m2))))
 180
 181
 182 ;(defun +pol (p1 p2)                    ;non destructif
 183 ;       (cond ((null p1) p2)            ;recursion non terminale
 184 ;             ((null p2) p1)            ;non utilise
 185 ;             ((equal (cdar p1) (cdar p2))
 186 ;               (let ((a (macplus (caar p1) (caar p2))))
 187 ;                 (cond ((zerop a) (+pol (cdr p1) (cdr p2)))
 188 ;                       (t (cons
 189 ;                            (cons a (cdar p1))
 190 ;                            (+pol (cdr p1) (cdr p2)))))))
 191 ;             ((funcall ordexp (cdar p1)(cdar p2))
 192 ;                 (cons (car p1) (+pol (cdr p1) p2)))
 193 ;             (t  (cons (car p2) (+pol (cdr p2) p1)))))
 194 ;
 195 ;(defun syz (pair)
 196 ;       (let ((f (caddr pair))
 197 ;             (g (cdddr pair))
 198 ;             (exp (caar pair)))
 199 ;         (-pol (*polmon f (cons (caar g) (mapcar '- exp (cdar f))))
 200 ;               (*polmon (cdr g) (cons (caar f) (mapcar '- exp (cdar g)))))))
 201
 202 ;OPERATIONS ELEMENTAIRES DESTRUCTIVES (POLYNOMES) COEFFICIENTS GENERAUX
 203
 204 (defun *gcofpol (c p)   ;destructif resultat dans p
 205         (mapc #'(lambda(u) (rplaca u (mactimes c (car u))))   p))
 206
 207         ;Cette fonction retourne et lie a p1 la difference (cdr p1) - p2
 208         ;iteratif et destructif
 209         ;le polynome nul est retourne (cons nil nil)
 210 (defun -gpol (p1 p2)
 211         (do ((p p1)
 212              (q p2)
 213              (e1) (e2))
 214             ((null q)
 215               (rplaca p1 (cadr p1))
 216               (rplacd p1 (cddr p1)))
 217             (setq e1 (cdadr p))
 218             (setq e2 (cdar q))
 219             (cond
 220               ((funcall ordexp e1 e2)
 221                 (setq p (cdr p)))
 222               ((equal e1 e2)
 223                 (let ((a (macmoins (caadr p) (caar q))))
 224                   (cond ((zerop a)
 225                           (rplacd p (cddr p)))
 226                         (t
 227                           (setq p (cdr p))
 228                           (rplaca (car p) a))))
 229                 (setq q (cdr q)))
 230               (t
 231                 (setq p (cdr (rplacd p (cons
 232                                           (cons
 233                                             (macmoins 0 (caar q))
 234                                             (cdar q))
 235                                           (cdr p)))))
 236                 (setq q (cdr q))))))
 237
 238 ;On definit maintenant la division d'un polynome par une base unitaire
 239 ;Le polynome est remplace par le resultat
 240
 241 (defun red1 (pol1 p pol2)                               ;destructif
 242         (let ((a (mapcar '- (cdar p) (cdar pol2))))     ;resultat dans pol1
 243           (cond ((minusp (apply 'min a))                ;pol2 est unitaire
 244                   pol1)
 245                 (t
 246                   (-gpol p
 247                         (*gpolmon
 248                             (cdr pol2)
 249                             (cons (caar p) a)))))))
 250
 251 (defun div-terme (p1 p  base)
 252         (do ((a (cdar p) (cdar p))
 253              (b nil a))
 254             ((eq a b) p)
 255             (mapc #'(lambda (u) (red1 p1 p u))       base)))
 256
 257 (defun divise (p1 base)
 258         (do ((p (cdr (div-terme p1 p1 base)) (cdr (div-terme p1 p base))))
 259             ((null p) p1)))
 260
 261 (defun monic (pol)
 262         (let ((c (caar pol)))
 263           (mapc #'(lambda (u)
 264                    (rplaca u (macdiv (car u) c)))
 265             pol)))
 266
 267 ;       ;Comme ce qui precede, mais pseudo division pour eviter
 268 ;       ;les divisions de coefficients
 269 ;(defun red2 (pol1 pol2)                                        ;destructif
 270 ;       (let ((a (mapcar '- (cdar pol1) (cdar pol2)))   ;resultat dans pol1
 271 ;             (c (caar pol1)))
 272 ;         (cond ((minusp (apply 'min a)) nil)           ;sans division
 273 ;               ((-pol (*cofpol (caar pol2) pol1)       ;rend nil si pol1..
 274 ;                      (*polmon    (cdr pol2)           ;..non modifie
 275 ;                                  (cons c  a)))))))
 276 ;
 277 ;(defun div-tete2 (p1 base)
 278 ;       (do ((a (cdar p1) (cdar p1))
 279 ;            (b nil a))
 280 ;           ((eq a b) p1)
 281 ;           (mapc #'(lambda (u) (red2 p1 u))       base)))
 282 ;
 283 ;(defun divise2 (p1 base)
 284 ;       (do ((p (cdr p1) (cdr (div-tete2 p base))))
 285 ;           ((null p) p1)))
 286 ;
 287 ;OPERATIONS ELEMENTAIRES PRESERVATIVES (POLYNOMES) COEFFICIENTS ENTIERS
 288
 289 ;definition des multiplications
 290 ;mactimes peut etre change si les coefficients sont restreints
 291 ; exemple, pour des entiers, on peut prendre "times"
 292
 293 (defun *ipolmon (pol mon)
 294         (mapcar #'(lambda (u) (*iterm u mon)) pol))
 295
 296 (defun *iterm (m1 m2)
 297         (cons
 298           (* (car m1) (car m2))
 299           (mapcar #'+ (cdr m1) (cdr m2))))
 300
 301 (defun isyz (pair)
 302         (let ((f (caddr pair))
 303               (g (cdddr pair))
 304               (exp (caar pair))
 305               (a)(b)(d))
 306           (setq $nbsyz (1+ $nbsyz)
 307                 a (caar g)
 308                 b (caar f)
 309                 d (gcd a b)
 310                 a (quotient a d)
 311                 b (quotient b d))
 312           (-ipol (*ipolmon f (cons a (mapcar '- exp (cdar f))))
 313                 (*ipolmon (cdr g) (cons b (mapcar '- exp (cdar g)))))))
 314
 315 ;OPERATIONS ELEMENTAIRES DESTRUCTIVES (POLYNOMES) COEFFICIENTS ENTIERS
 316
 317 (defun *icofpol (c p)   ;destructif resultat dans p
 318         (mapc #'(lambda(u) (rplaca u (* c (car u))))   p))
 319
 320         ;Cette fonction retourne et lie a p1 la difference (cdr p1) - p2
 321         ;iteratif et destructif
 322         ;le polynome nul est retourne (cons nil nil)
 323 (defun -ipol (p1 p2)
 324         (-ipol2 p1 p2)
 325         (rplaca p1 (cadr p1))
 326         (rplacd p1 (cddr p1)))
 327
 328 (defun -ipol2 (p1 p2)
 329         (do ((p p1)
 330              (q p2)
 331              (e1) (e2))
 332             ((null q))
 333             (setq e1 (cdadr p))
 334             (setq e2 (cdar q))
 335             (cond
 336               ((funcall ordexp e1 e2)
 337                 (setq p (cdr p)))
 338               ((equal e1 e2)
 339                 (let ((a (- (caadr p) (caar q))))
 340                   (cond ((zerop a)
 341                           (rplacd p (cddr p)))
 342                         (t
 343                           (setq p (cdr p))
 344                           (rplaca (car p) a))))
 345                 (setq q (cdr q)))
 346               (t
 347                 (setq p (cdr (rplacd p (cons
 348                                           (cons
 349                                             (- (caar q))
 350                                             (cdar q))
 351                                           (cdr p)))))
 352                 (setq q (cdr q))))))
 353
 354         ;On definit maintenant la division d'un polynome par une base
 355         ;Le polynome est remplace par le resultat
 356         ;Pseudo division pour eviter
 357         ;les divisions de coefficients
 358 (defun ired2 (pol1 pol2)                                ;destructif
 359         (let ((a (mapcar '- (cdar pol1) (cdar pol2)))   ;resultat dans pol1
 360               (c)(b)(d))                                ;sans division
 361           (cond ((minusp (apply 'min a)) nil)           ;rend nil si pol1..
 362                 (t (setq $nbred (1+ $nbred)             ;..non modifie
 363                          b (caar pol2)
 364                          c (caar pol1)
 365                          d (gcd b c)                    ;rend (cons nil nil)
 366                          b (quotient b d)               ;si pol1 devient nul
 367                          c (quotient c d))
 368                    (*icofpol b (cdr pol1))
 369                    (-ipol  pol1
 370                            (*ipolmon  (cdr pol2)
 371                                       (cons c  a)))))))
 372
 373 (defun iredp (p1 p q)
 374         (let ((pp (cdr p))
 375               (a) (b) (c) (d))
 376           (and pp
 377                (cond ((not (minusp (apply 'min
 378                                      (setq a (mapcar '- (cdar pp) (cdar q))))))
 379                        (setq    b (caar q)
 380                                 c (caar pp)
 381                                 d (gcd b c)
 382                                 b (quotient b d)
 383                                 c (quotient c d))
 384                        (*icofpol b p1)
 385                        (-ipol2 pp
 386                           (*ipolmon (cdr q)
 387                                     (cons c a)))
 388                        (rplacd p (cdr pp)))))))
 389
 390 (defun idivp (p1 p base)                ;pour idivise2, il faut multiplier
 391         (do ((a (cdadr p) (cdadr p))    ;tout le polynome dividende
 392              (b nil a))
 393             ((eq a b) (cdr p))
 394             (mapc #'(lambda (u) (iredp p1 p u))       base)))
 395
 396 (defun idivise2 (p1 base)
 397         (do ((p p1 (idivp p1 p base)))
 398             ((null p) (primpart p1))))
 399
 400 ;OPERATIONS SUR LES COEFFICIENTS
 401
 402 (defun mactimes (a b)
 403         (meval (list '(mtimes) a b)))
 404
 405 (defun macplus (a b)
 406         (meval (list '(mplus) a b)))
 407
 408 (defun macmoins (a b)
 409         (meval (list '(mplus) a (list '(mminus) b))))
 410
 411 (defun macdiv (a b)
 412         (meval (list '(mquotient) a b)))
 413
 414 ;(defun *gcd (lnb)
 415 ;       (do ((r (car lnb) (gcd r (car l)))
 416 ;            (l (cdr lnb) (cdr l)))
 417 ;           ((null l) r)))
 418
 419 (defun primpart (p)
 420         (let ((d (do ((p (cdr p) (cdr p))       ;calcul du contenu
 421                       (g (caar p) (gcd g (caar p))))
 422                      ((or (eql (abs g) 1) (null p))
 423                          g) )))
 424           (or (eql (abs d) 1)
 425               (mapc #'(lambda (u)               ;diviser par le contenu
 426                  (rplaca u (quotient (car u) d)))
 427                 p))
 428           p))
 429
 430 ;(defun primpart (p)
 431 ;       (let ((d (*gcd (mapcar 'car
 432 ;                        p))))
 433 ;         (mapc #'(lambda (u)
 434 ;                  (rplaca u (quotient (car u) d)))
 435 ;               p)))
 436
 437 ;LA BASE STANDARD
 438
 439 ;Construction de la base
 440 ;Methode a la Buchberger-Bouzeghoub
 441 ;base: liste de polynomes
 442 ;paires: candidats sygyzies:   ((exp."deja reduite").((exp1.exp2).(f.g)))
 443
 444 ;algorithme:
 445 ;       reduire tous les polynomes
 446 ;       calculer la liste des paires; vivier construit celles qui ne sont
 447 ;pas encore visitees
 448 ;       reduire la premiere paire, inserer le resultat dans la base et
 449 ;rereduire
 450 ;remettre a jour la liste des paires et recommencer
 451
 452 (defun grobner (gener ordexp)
 453           (setq $base (cons nil (sort (copy-tree gener)
 454                                 #'(lambda (u v) (ordexp (cdar v) (cdar u))))))
 455           (let ((paires nil))
 456           (setq $nbsyz 0 $nbred 0 $nbreduc 0 $nbred0 0)
 457           (reduire $base (cdr $base) (cddr $base))
 458           (setq paires
 459             (mapcon #'(lambda (u)
 460                        (let ((exp (cdaar u))
 461                              (f (car u)))
 462                          (mapcar #'(lambda (v) (apparier f v exp (cdar v)))
 463                                  (cdr u))))
 464                     (cdr $base)))
 465           (setq paires (cons
 466                          nil
 467                          (sort paires ordpairs)))
 468           (do ((l (vivier paires))                      ;parcourir les paires
 469                (exp) (exp1) (exp2) (p) (q))
 470               ((null l))
 471               (setq     p (car l)        exp (caar p)
 472                         exp1 (caadr p)   exp2 (cdadr p))
 473               (cond ((cdar p)                   ;paire deja vue
 474                       (setq l (cdr l)))
 475                                 ;on cherche maintenant h dans base tel
 476                                 ;que h divise exp et que les syzygies
 477                                 ;(f h) et (g h) ont ete calculees
 478                     ((do ((l1 (cdr $base) l3)           ;"critere 3"
 479                           (l3 (cddr $base) (cdr l3))
 480                           (h (cadr $base) (car l3))
 481                           (exph (cdaadr $base) (cdaar l3)))
 482                          ((cond ((null l1) t)           ;pas trouve de h
 483                                 ((funcall ordexp exph exp) (not (setq l1 nil)))
 484                                 ((and
 485                                    (not (minusp (apply 'min
 486                                                   (mapcar '- exp exph))))
 487                                    (let ((exph1 (mapcar 'max exph exp1)))
 488                                       (or (funcall ordexp exp exph1)
 489                                           (and (equal exp exph1)
 490                                                (funcall ordexp exp2 exph))))
 491                                    (let ((exph2 (mapcar 'max exph exp2)))
 492                                       (or (funcall ordexp exp exph2)
 493                                           (and (equal exp exph2)
 494                                                (funcall ordexp exp1 exph)))))))
 495                 ;               (t nil)
 496                            l1))
 497                       (setq l (cdr l))
 498                       (rplacd (car p) t))       ;le critere 3 est verifie
 499                     (t
 500                       (rplacd (car p) t)        ;la paire va etre traitee
 501                       (setq p (isyz p)
 502                             q (inserer p $base))
 503                       (cond (q
 504                               (reduire $base q nil)
 505                               (format t "~% nbsyz = ~a    nbred = ~a    nbreduc = ~a     nbred0 = ~a    lbase = ~d~%"
 506                                       $nbsyz $nbred $nbreduc $nbred0 (length (cdr base)))
 507
 508                               (print (escalier (cdr $base)))
 509                               (print "
 510 ")
 511                               (rafraichir  paires p $base)
 512                               (setq l (vivier paires)))
 513                             (t
 514                               (setq l (cdr l))
 515                               (format t "~% nbsyz = ~a    nbred = ~a    nbreduc = ~a     nbred0 = ~a    lbase = ~a~% paires a voir = ~d~%"
 516                                       $nbsyz $nbred $nbreduc $nbred0 (length (cdr base)) (length l))
 517
 518                               )))))     ;paire suivante
 519           (setq $base (cdr $base))
 520           (mapc #'(lambda (u) (idivise2 u $base))
 521                 $base)                          ;reduire completement
 522           (mapc 'monic $base)))
 523
 524 (defun apparier (p q expp expq)
 525           (cons
 526             (cons
 527               (mapcar #'max expp expq)          ;caar
 528               (zerop (apply #'+                 ;cdar ;exposants etrangers
 529                        (mapcar 'min expp expq))))
 530             (cons
 531               (cons expp expq)                  ;caadr  cdadr
 532               (cons p q))))                     ;caddr  cdddr
 533
 534 (defun vivier (paires)
 535         (let ((v (cons nil nil)))
 536           (do ((l (cdr paires) (cdr l))
 537                (r v))
 538               ((null l))
 539               (or  (cdaar l)                    ;paire deja vue
 540                    (setq r (cdr (rplacd r
 541                                   (cons (car l) nil))))))
 542           (sort (cdr v)  ordpairs)))
 543
 544 (defun rafraichir (paires q base)
 545           (let ((nvp (cons nil nil)))
 546             (cond ((car q)                      ;paires p q pour p dans base
 547                     (do ((l (cdr base) (cdr l))
 548                          (p (cadr base) (cadr l))
 549                          (nvpa nvp)             ;liste des nouvelles paires
 550                          (expq (cdar q))
 551                          (bit t))
 552                         ((null l))
 553                         (cond ((eq q p)
 554                                 (setq bit nil))
 555                               (bit
 556                                 (setq nvpa (cdr (rplacd nvpa
 557                                                  (cons (apparier p q
 558                                                                (cdar p) expq)
 559                                                    nil)))))
 560                               (t
 561                                 (setq nvpa (cdr (rplacd nvpa
 562                                                  (cons (apparier q p
 563                                                                expq (cdar p))
 564                                                    nil)))))))))
 565             (do ((l paires)                     ;rafraichir les anciennes
 566                  (pp) (f) (g) (ef) (eg))
 567                 ((null (cdr l)))
 568                 (setq pp (cadr l)
 569                       f (caddr pp) g (cdddr pp)
 570                       ef (cdar f) eg (cdar g))
 571                 (cond ((or (null ef)            ;enlever si un pol est nul
 572                            (null eg))
 573                         (rplacd l (cddr l)))
 574                       ((or (not (equal ef (caadr pp)))  ;un pol a change?
 575                            (not (equal eg (cdadr pp))))
 576                         (setq pp (cond ((funcall ordexp ef eg)
 577                                                 (apparier g f eg ef))
 578                                        (t       (apparier f g ef eg))))
 579                         (rplacd l (cons pp (cddr l)))
 580                         (setq l (cdr l)))
 581                       (t (setq l (cdr l)))))
 582             (nconc paires (cdr nvp))))
 583
 584           ;recherche des reductibles et reduction des elements de base
 585 (defun reduire (base a b)       ;a ---> dernier element reduit
 586           (do ((l0 a)                   ;       (vient d'etre insere)
 587                (l1 (cdr a))     ;b ---> provenance de cet element
 588                (q0 (car a))
 589                (q (cadr a) (car l1))
 590                (bit (eq a (cdr b)) (or bit (eq l1 (cdr b)))))
 591               ((null l1) base)
 592               (cond (bit                ;reduction par tout ce qui precede
 593                       (do ((l (cdr base) (cdr l)))      ;q reductible?
 594                           ((cond
 595                             ((eq l l1) (setq l0 l1))    ;q irreductible
 596                             ((ired2 q (car l))          ;ired2 est un predicat
 597                               (rplacd l0 (cdr l1))      ;enlever q
 598                               (setq $nbreduc (1+ $nbreduc))
 599                               (and (setq l1 (inserer q base))
 600                                    (setq q0 q b l0 l0 l1
 601                                          bit (eq l0 (cdr b))))
 602                               t)))))
 603                     (t                  ;reduction par q0 seulement
 604                       (cond ((ired2 q q0)               ;q reductible
 605                               (rplacd l0 (cdr l1))      ;enlever q
 606                               (setq $nbreduc (1+ $nbreduc))
 607                               (and (setq l1 (inserer q base))
 608                                    (setq q0 q b l0 l0 l1
 609                                          bit (eq l0 (cdr b)))))
 610                             (t (setq l0 l1)))))
 611               (setq l1 (cdr l0))))
 612
 613 (defun inserer (q base)                         ;retourne le lieu d'insertion
 614         (do ((l2 base)                          ;chercher ou`
 615              (l3 (cdr base))                    ;inserer q...
 616              (qq (cadr base) (car l3))
 617              (expq (cdar q)))
 618             ((cond
 619                ((null expq)
 620                  (setq $nbred0 (1+ $nbred0))
 621                  (setq l2 q))   ;si q est nul retourner nil..
 622                ((or                     ;... pour inserer, mais t pour cond
 623                   (null l3)
 624                   (funcall ordexp (cdar qq) expq))
 625                  (rplacd l2 nil)                ;pour accelerer la reduction
 626                  (setq q (primpart q))
 627                  (setq q (idivise2 q (cdr $base)))      ;reduire completement
 628                  (rplacd l2 (cons  q l3))))     ;inserer q
 629               (cdr l2))                                 ;lieu a retourner
 630             (cond ((ired2 q qq)                         ;reduire q
 631                     (setq expq (cdar q))                ;et
 632                     (setq l2 base))                     ;repartir
 633                   (t
 634                     (setq l2 l3)))
 635             (setq l3 (cdr l2))))
 636
 637
 638 ;RESOLUTION DE SYSTEMES
 639 ;On part d'une base standard pour un ordre compatible avec le degre. Elle
 640 ;est representee par une liste de polynomes, chacun d'eux etant une liste de
 641 ;monomes de la forme (coeff e1 ... en).
 642
 643 ;L'escalier est la liste des vecteurs exposants des monomes dominants.
 644 (defun escalier (base)
 645         (mapcar 'cdar base))
 646
 647 ;La avriete est vide(i.e. de dim. -1) si l'escalier est reduit au vecteur nul
 648 ;(monome constant)
 649 (defun dim-1 (esc)
 650         (apply 'and (mapcar 'zerop (car esc))) )
 651
 652 ;La dimension est 0 si il y a un element de l'escalier sur chaque axe.
 653 ;On suppose la base reduite, ce qui implique qu'il n'y a pas 2 elements de
 654 ;l'escalier sur le meme axe
 655 (defun dim0 (esc)
 656          (or (dim-1 esc)
 657              (do ((l esc (cdr l))
 658                   (i 0
 659                      (+ i
 660                         (- 2 (min 2
 661                                (apply '+
 662                                  (mapcar #'(lambda (u) (min u 1))
 663                                          (car l))))))))
 664                  ((null l)(= i (length (car esc)))))))
 665
 666 ;Construction de la liste des monomes sous l'escalier
 667 ;Le predicat sousp indique si un monome est sous l'escalier
 668 (defun sous-esc (esc)
 669         (do ((l ())
 670              (m (mapcar #'(lambda (u)u 0) (car esc)))
 671              (i (cons 0 nil)))
 672             ((null i) l)
 673             (rplaca i (1+ (car i)))
 674             (cond ((sousp m esc)
 675                      (setq l (cons (copy-tree m) l))
 676                      (setq i m))
 677                   (t
 678                     (rplaca i 0)
 679                     (setq i (cdr i))))))
 680
 681 ;Test mon2 ne divise pas mon1 pour 2 vecteurs d'exposants
 682 (defun notdivmon (mon1 mon2)
 683         (minusp (apply 'min (mapcar '- mon1 mon2))))
 684
 685 ;sousp est utilise par sous-esc pour tester si un monome est au sous
 686 ;l'escalier
 687 (defun sousp (mon esc)
 688         (apply 'and
 689                (mapcar #'(lambda (u) (notdivmon mon u)) esc)))
 690
 691 ;ICI, IL FAUT RANGER LE RESULTAT DE LA FONCTION PRECEDENTE POUR L'ORDRE
 692 ;COMPATIBLE CHOISI; CE RESULTAT EST UTILISE SOUS LE NOM DE basli
 693
 694 ;Construction de la matrice dont le polynome caracteristique a pour racine
 695 ;les premieres coordonnees des solutions
 696 ;les lignes et colonnes sont indexees par basli, chaque ligne etant en
 697 ;representation creuse
 698 ;LA FONCTION DIVISE RETOURNE LE RESTE DE LA DIVISION D'UN POLYNOME PAR LA
 699 ;BASE STANDARD
 700 (defun umat (basli base)
 701         (mapcar #'(lambda (u) (divise
 702                                (cons
 703                                  (cons 1 (rplaca (copy-tree u) (1+ (car u))))
 704                                  nil)
 705                                base))
 706                  basli))
 707
 708 ;densification d'une ligne
 709 (defun densif (ligne basli)
 710         (let ((nl (cons '(mlist) nil)))
 711           (do ((b basli (cdr b))
 712                (l ligne)
 713                (c 0 0)
 714                (nl nl (cdr (rplacd nl (cons c nil)))))
 715               ((null b))
 716               (and (equal (car b) (cdar l))
 717                    (setq c (caar l))
 718                    (setq l (cdr l))))
 719           nl))
 720
 721 ;La suite appelle des fonctions macsyma et se trouve dans japsmac.l
 722 ;Partie macsyma de japs.l
 723
 724 ;matrice macsyma; umatr est le resultat de umat
 725 (defun matmac (basli base)
 726         (meval (cons '($matrix) (mapcar #'(lambda (u) (densif u basli))
 727                                          (umat basli base)))))
 728
 729 ;Enfin le polynome eliminant, retourne en macsyma
 730 (defun elimpol (base ordexp var1)
 731         (setq $ratmx t)
 732         (let ((esc (escalier base))
 733               ($ratmx t))
 734           (cond
 735             ((dim-1 esc) 1)
 736             ((dim0 esc)
 737               ($charpoly
 738                    (matmac   (sort (sous-esc esc) ordexp)
 739                              base)
 740                    var1))
 741             (t "dimension positive"))))
 742
 743 (defun $elimpol_err (ordexp var1)
 744         (elimpol $errexp1  ordexp var1))
 745
 746 ;;;; -*- Mode: LISP; Package: Macsyma; Base:10 -*- Saved by dl
 747 ;;;; Macsyma version 309.0
 748 ;(mdefprop $elimpol
 749 ;          ((lambda nil)
 750 ;           ((mlist) $base $var1 $ordre)
 751 ;           ((elimpol) $base $var1 $ordre))
 752 ;          mexpr)
 753 ;(add2lnc '(($elimpol) $base $var1 $ordre) $functions)
 754
 755 (defun decmon (a)
 756         (mapcar
 757             #'(lambda (u v) (list '(mexpt) u v))
 758             lvar a))
 759
 760 (defun decterm (a)
 761         (meval
 762            (cons
 763               '(mtimes)
 764               (cons
 765                  (car a)
 766                  (decmon (cdr a))))))
 767
 768 (defun decpol (a)
 769         (meval
 770            (cons '(mplus)
 771               (mapcar 'decterm a))))
 772
 773 (defun decode (lpol lvar)
 774         (meval
 775            (cons '(mlist)
 776               (mapcar 'decpol lpol))))
 777
 778 (defun $decode_poly(poly sample &aux header)
 779   (cond ((and sample (consp sample)) (setq header (car sample)))
 780         (t  (setq header $header)))
 781   (let ((monoms (loop for v in (fourth header)
 782                        collect (list v 1 1))))
 783     (loop for v in poly
 784            with ans = 0
 785            do (setq ans
 786                     (pplus
 787                       ans
 788                       (loop for deg in (cdr v)
 789                              for mon in monoms
 790                              with term = (car v)
 791                              when (not (eql 0 deg))
 792                              do (setq term (ptimes term
 793                                                 (pexpt mon deg)))
 794                              finally (return term)
 795                              )))
 796                     finally (return
 797                               (cons header
 798                                     (cons ans 1))))))
 799
 800 (defun show-lazard (term)
 801   (cond ((or (atom term) (atom (car term)))(error "bad term"))
 802         ((consp (caar term)) (loop for v in term do (show-lazard v)))
 803         (t (displa ($decode_poly term nil)))))