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