1 | ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
|
---|
2 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
3 | ;;;
|
---|
4 | ;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
|
---|
5 | ;;;
|
---|
6 | ;;; This program is free software; you can redistribute it and/or modify
|
---|
7 | ;;; it under the terms of the GNU General Public License as published by
|
---|
8 | ;;; the Free Software Foundation; either version 2 of the License, or
|
---|
9 | ;;; (at your option) any later version.
|
---|
10 | ;;;
|
---|
11 | ;;; This program is distributed in the hope that it will be useful,
|
---|
12 | ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
13 | ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
14 | ;;; GNU General Public License for more details.
|
---|
15 | ;;;
|
---|
16 | ;;; You should have received a copy of the GNU General Public License
|
---|
17 | ;;; along with this program; if not, write to the Free Software
|
---|
18 | ;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
---|
19 | ;;;
|
---|
20 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
21 |
|
---|
22 | (in-package :ngrobner)
|
---|
23 |
|
---|
24 | ;;----------------------------------------------------------------
|
---|
25 | ;; This package implements BASIC OPERATIONS ON MONOMIALS
|
---|
26 | ;;----------------------------------------------------------------
|
---|
27 | ;; DATA STRUCTURES: Monomials are represented as lists:
|
---|
28 | ;;
|
---|
29 | ;; monom: (n1 n2 ... nk) where ni are non-negative integers
|
---|
30 | ;;
|
---|
31 | ;; However, lists may be implemented as other sequence types,
|
---|
32 | ;; so the flexibility to change the representation should be
|
---|
33 | ;; maintained in the code to use general operations on sequences
|
---|
34 | ;; whenever possible. The optimization for the actual representation
|
---|
35 | ;; should be left to declarations and the compiler.
|
---|
36 | ;;----------------------------------------------------------------
|
---|
37 | ;; EXAMPLES: Suppose that variables are x and y. Then
|
---|
38 | ;;
|
---|
39 | ;; Monom x*y^2 ---> #(1 2)
|
---|
40 | ;;
|
---|
41 | ;;----------------------------------------------------------------
|
---|
42 |
|
---|
43 | (deftype exponent ()
|
---|
44 | "Type of exponent in a monomial."
|
---|
45 | 'fixnum)
|
---|
46 |
|
---|
47 | (deftype monom (&optional dim)
|
---|
48 | "Type of monomial."
|
---|
49 | `(simple-array exponent (,dim)))
|
---|
50 |
|
---|
51 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
52 | ;;
|
---|
53 | ;; Construction of monomials
|
---|
54 | ;;
|
---|
55 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
56 |
|
---|
57 | (defmacro make-monom (dim &key (initial-contents nil initial-contents-supplied-p)
|
---|
58 | (initial-element 0 initial-element-supplied-p))
|
---|
59 | "Make a monomial with DIM variables. Additional argument
|
---|
60 | INITIAL-CONTENTS specifies the list of powers of the consecutive
|
---|
61 | variables. The alternative additional argument INITIAL-ELEMENT
|
---|
62 | specifies the common power for all variables."
|
---|
63 | ;;(declare (fixnum dim))
|
---|
64 | `(make-array ,dim
|
---|
65 | :element-type 'exponent
|
---|
66 | ,@(when initial-contents-supplied-p `(:initial-contents ,initial-contents))
|
---|
67 | ,@(when initial-element-supplied-p `(:initial-element ,initial-element))))
|
---|
68 |
|
---|
69 | |
---|
70 |
|
---|
71 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
72 | ;;
|
---|
73 | ;; Operations on monomials
|
---|
74 | ;;
|
---|
75 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
76 |
|
---|
77 | (defmacro monom-elt (m index)
|
---|
78 | "Return the power in the monomial M of variable number INDEX."
|
---|
79 | `(elt ,m ,index))
|
---|
80 |
|
---|
81 | (defun monom-dimension (m)
|
---|
82 | "Return the number of variables in the monomial M."
|
---|
83 | (length m))
|
---|
84 |
|
---|
85 | (defun monom-total-degree (m &optional (start 0) (end (length m)))
|
---|
86 | "Return the todal degree of a monomoal M. Optinally, a range
|
---|
87 | of variables may be specified with arguments START and END."
|
---|
88 | (declare (type monom m) (fixnum start end))
|
---|
89 | (reduce #'+ m :start start :end end))
|
---|
90 |
|
---|
91 | (defun monom-sugar (m &aux (start 0) (end (length m)))
|
---|
92 | "Return the sugar of a monomial M. Optinally, a range
|
---|
93 | of variables may be specified with arguments START and END."
|
---|
94 | (declare (type monom m) (fixnum start end))
|
---|
95 | (monom-total-degree m start end))
|
---|
96 |
|
---|
97 | (defun monom-div (m1 m2 &aux (result (copy-seq m1)))
|
---|
98 | "Divide monomial M1 by monomial M2."
|
---|
99 | (declare (type monom m1 m2 result))
|
---|
100 | (map-into result #'- m1 m2))
|
---|
101 |
|
---|
102 | (defun monom-mul (m1 m2 &aux (result (copy-seq m1)))
|
---|
103 | "Multiply monomial M1 by monomial M2."
|
---|
104 | (declare (type monom m1 m2 result))
|
---|
105 | (map-into result #'+ m1 m2))
|
---|
106 |
|
---|
107 | (defun monom-divides-p (m1 m2)
|
---|
108 | "Returns T if monomial M1 divides monomial M2, NIL otherwise."
|
---|
109 | (declare (type monom m1 m2))
|
---|
110 | (every #'<= m1 m2))
|
---|
111 |
|
---|
112 | (defun monom-divides-monom-lcm-p (m1 m2 m3)
|
---|
113 | "Returns T if monomial M1 divides MONOM-LCM(M2,M3), NIL otherwise."
|
---|
114 | (declare (type monom m1 m2 m3))
|
---|
115 | (every #'(lambda (x y z) (declare (type exponent x y z)) (<= x (max y z))) m1 m2 m3))
|
---|
116 |
|
---|
117 | (defun monom-lcm-divides-monom-lcm-p (m1 m2 m3 m4)
|
---|
118 | "Returns T if monomial MONOM-LCM(M1,M2) divides MONOM-LCM(M3,M4), NIL otherwise."
|
---|
119 | (declare (type monom m1 m2 m3 m4))
|
---|
120 | (every #'(lambda (x y z w) (declare (type exponent x y z w)) (<= (max x y) (max z w))) m1 m2 m3 m4))
|
---|
121 |
|
---|
122 | (defun monom-lcm-equal-monom-lcm-p (m1 m2 m3 m4)
|
---|
123 | "Returns T if monomial MONOM-LCM(M1,M2) equals MONOM-LCM(M3,M4), NIL otherwise."
|
---|
124 | (declare (type monom m1 m2 m3 m4))
|
---|
125 | (every #'(lambda (x y z w) (declare (type exponent x y z w)) (= (max x y) (max z w))) m1 m2 m3 m4))
|
---|
126 |
|
---|
127 | (defun monom-divisible-by-p (m1 m2)
|
---|
128 | "Returns T if monomial M1 is divisible by monomial M2, NIL otherwise."
|
---|
129 | (declare (type monom m1 m2))
|
---|
130 | (every #'>= m1 m2))
|
---|
131 |
|
---|
132 | (defun monom-rel-prime-p (m1 m2)
|
---|
133 | "Returns T if two monomials M1 and M2 are relatively prime (disjoint)."
|
---|
134 | (declare (type monom m1 m2))
|
---|
135 | (every #'(lambda (x y) (declare (type exponent x y)) (zerop (min x y))) m1 m2))
|
---|
136 |
|
---|
137 | (defun monom-equal-p (m1 m2)
|
---|
138 | "Returns T if two monomials M1 and M2 are equal."
|
---|
139 | (declare (type monom m1 m2))
|
---|
140 | (every #'= m1 m2))
|
---|
141 |
|
---|
142 | (defun monom-lcm (m1 m2 &aux (result (copy-seq m1)))
|
---|
143 | "Returns least common multiple of monomials M1 and M2."
|
---|
144 | (declare (type monom m1 m2))
|
---|
145 | (map-into result #'max m1 m2))
|
---|
146 |
|
---|
147 | (defun monom-gcd (m1 m2 &aux (result (copy-seq m1)))
|
---|
148 | "Returns greatest common divisor of monomials M1 and M2."
|
---|
149 | (declare (type monom m1 m2))
|
---|
150 | (map-into result #'min m1 m2))
|
---|
151 |
|
---|
152 | (defun monom-depends-p (m k)
|
---|
153 | "Return T if the monomial M depends on variable number K."
|
---|
154 | (declare (type monom m) (fixnum k))
|
---|
155 | (plusp (elt m k)))
|
---|
156 |
|
---|
157 | (defmacro monom-map (fun m &rest ml &aux (result `(copy-seq ,m)))
|
---|
158 | `(map-into ,result ,fun ,m ,@ml))
|
---|
159 |
|
---|
160 | (defmacro monom-append (m1 m2)
|
---|
161 | `(concatenate 'monom ,m1 ,m2))
|
---|
162 |
|
---|
163 | (defmacro monom-contract (k m)
|
---|
164 | `(subseq ,m ,k))
|
---|
165 |
|
---|
166 | (defun monom-exponents (m)
|
---|
167 | (declare (type monom m))
|
---|
168 | (coerce m 'list))
|
---|