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 |
|
---|
23 | ;;----------------------------------------------------------------
|
---|
24 | ;; This package implements BASIC OPERATIONS ON MONOMIALS
|
---|
25 | ;;----------------------------------------------------------------
|
---|
26 | ;; DATA STRUCTURES: Monomials are represented as lists:
|
---|
27 | ;;
|
---|
28 | ;; monom: (n1 n2 ... nk) where ni are non-negative integers
|
---|
29 | ;;
|
---|
30 | ;; However, lists may be implemented as other sequence types,
|
---|
31 | ;; so the flexibility to change the representation should be
|
---|
32 | ;; maintained in the code to use general operations on sequences
|
---|
33 | ;; whenever possible. The optimization for the actual representation
|
---|
34 | ;; should be left to declarations and the compiler.
|
---|
35 | ;;----------------------------------------------------------------
|
---|
36 | ;; EXAMPLES: Suppose that variables are x and y. Then
|
---|
37 | ;;
|
---|
38 | ;; Monom x*y^2 ---> (1 2)
|
---|
39 | ;;
|
---|
40 | ;;----------------------------------------------------------------
|
---|
41 |
|
---|
42 | (deftype exponent ()
|
---|
43 | "Type of exponent in a monomial."
|
---|
44 | 'fixnum)
|
---|
45 |
|
---|
46 | (deftype monom (&optional dim)
|
---|
47 | "Type of monomial."
|
---|
48 | `(simple-array exponent (,dim)))
|
---|
49 |
|
---|
50 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
51 | ;;
|
---|
52 | ;; Construction of monomials
|
---|
53 | ;;
|
---|
54 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
55 |
|
---|
56 | (defmacro make-monom (dim &key (initial-contents nil initial-contents-supplied-p)
|
---|
57 | (initial-element 0 initial-element-supplied-p))
|
---|
58 | "Make a monomial with DIM variables. Additional argument
|
---|
59 | INITIAL-CONTENTS specifies the list of powers of the consecutive
|
---|
60 | variables. The alternative additional argument INITIAL-ELEMENT
|
---|
61 | specifies the common power for all variables."
|
---|
62 | ;;(declare (fixnum dim))
|
---|
63 | `(make-array ,dim
|
---|
64 | :element-type 'exponent
|
---|
65 | ,@(when initial-contents-supplied-p `(:initial-contents ,initial-contents))
|
---|
66 | ,@(when initial-element-supplied-p `(:initial-element ,initial-element))))
|
---|
67 |
|
---|
68 | |
---|
69 |
|
---|
70 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
71 | ;;
|
---|
72 | ;; Operations on monomials
|
---|
73 | ;;
|
---|
74 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
75 |
|
---|
76 | (defmacro monom-elt (m index)
|
---|
77 | "Return the power in the monomial M of variable number INDEX."
|
---|
78 | `(elt ,m ,index))
|
---|
79 |
|
---|
80 | (defun monom-dimension (m)
|
---|
81 | "Return the number of variables in the monomial M."
|
---|
82 | (length m))
|
---|
83 |
|
---|
84 | (defun monom-total-degree (m &optional (start 0) (end (length m)))
|
---|
85 | "Return the todal degree of a monomoal M. Optinally, a range
|
---|
86 | of variables may be specified with arguments START and END."
|
---|
87 | (declare (type monom m) (fixnum start end))
|
---|
88 | (reduce #'+ m :start start :end end))
|
---|
89 |
|
---|
90 | (defun monom-sugar (m &aux (start 0) (end (length m)))
|
---|
91 | "Return the sugar of a monomial M. Optinally, a range
|
---|
92 | of variables may be specified with arguments START and END."
|
---|
93 | (declare (type monom m) (fixnum start end))
|
---|
94 | (monom-total-degree m start end))
|
---|
95 |
|
---|
96 | (defun monom-div (m1 m2 &aux (result (copy-seq m1)))
|
---|
97 | "Divide monomial M1 by monomial M2."
|
---|
98 | (declare (type monom m1 m2 result))
|
---|
99 | (map-into result #'- m1 m2))
|
---|
100 |
|
---|
101 | (defun monom-mul (m1 m2 &aux (result (copy-seq m1)))
|
---|
102 | "Multiply monomial M1 by monomial M2."
|
---|
103 | (declare (type monom m1 m2 result))
|
---|
104 | (map-into result #'+ m1 m2))
|
---|
105 |
|
---|
106 | (defun monom-divides-p (m1 m2)
|
---|
107 | "Returns T if monomial M1 divides monomial M2, NIL otherwise."
|
---|
108 | (declare (type monom m1 m2))
|
---|
109 | (every #'<= m1 m2))
|
---|
110 |
|
---|
111 | (defun monom-divides-monom-lcm-p (m1 m2 m3)
|
---|
112 | "Returns T if monomial M1 divides MONOM-LCM(M2,M3), NIL otherwise."
|
---|
113 | (declare (type monom m1 m2 m3))
|
---|
114 | (every #'(lambda (x y z) (declare (type exponent x y z)) (<= x (max y z))) m1 m2 m3))
|
---|
115 |
|
---|
116 | (defun monom-lcm-divides-monom-lcm-p (m1 m2 m3 m4)
|
---|
117 | "Returns T if monomial MONOM-LCM(M1,M2) divides MONOM-LCM(M3,M4), NIL otherwise."
|
---|
118 | (declare (type monom m1 m2 m3 m4))
|
---|
119 | (every #'(lambda (x y z w) (declare (type exponent x y z w)) (<= (max x y) (max z w))) m1 m2 m3 m4))
|
---|
120 |
|
---|
121 | (defun monom-lcm-equal-monom-lcm-p (m1 m2 m3 m4)
|
---|
122 | "Returns T if monomial MONOM-LCM(M1,M2) equals MONOM-LCM(M3,M4), NIL otherwise."
|
---|
123 | (declare (type monom m1 m2 m3 m4))
|
---|
124 | (every #'(lambda (x y z w) (declare (type exponent x y z w)) (= (max x y) (max z w))) m1 m2 m3 m4))
|
---|
125 |
|
---|
126 | (defun monom-divisible-by-p (m1 m2)
|
---|
127 | "Returns T if monomial M1 is divisible by monomial M2, NIL otherwise."
|
---|
128 | (declare (type monom m1 m2))
|
---|
129 | (every #'>= m1 m2))
|
---|
130 |
|
---|
131 | (defun monom-rel-prime-p (m1 m2)
|
---|
132 | "Returns T if two monomials M1 and M2 are relatively prime (disjoint)."
|
---|
133 | (declare (type monom m1 m2))
|
---|
134 | (every #'(lambda (x y) (declare (type exponent x y)) (zerop (min x y))) m1 m2))
|
---|
135 |
|
---|
136 | (defun monom-equal-p (m1 m2)
|
---|
137 | "Returns T if two monomials M1 and M2 are equal."
|
---|
138 | (declare (type monom m1 m2))
|
---|
139 | (every #'= m1 m2))
|
---|
140 |
|
---|
141 | (defun monom-lcm (m1 m2 &aux (result (copy-seq m1)))
|
---|
142 | "Returns least common multiple of monomials M1 and M2."
|
---|
143 | (declare (type monom m1 m2))
|
---|
144 | (map-into result #'max m1 m2))
|
---|
145 |
|
---|
146 | (defun monom-gcd (m1 m2 &aux (result (copy-seq m1)))
|
---|
147 | "Returns greatest common divisor of monomials M1 and M2."
|
---|
148 | (declare (type monom m1 m2))
|
---|
149 | (map-into result #'min m1 m2))
|
---|
150 |
|
---|
151 | (defun monom-depends-p (m k)
|
---|
152 | "Return T if the monomial M depends on variable number K."
|
---|
153 | (declare (type monom m) (fixnum k))
|
---|
154 | (plusp (elt m k)))
|
---|
155 |
|
---|
156 | (defmacro monom-map (fun m &rest ml &aux (result `(copy-seq ,m)))
|
---|
157 | `(map-into ,result ,fun ,m ,@ml))
|
---|
158 |
|
---|
159 | (defmacro monom-append (m1 m2)
|
---|
160 | `(concatenate 'monom ,m1 ,m2))
|
---|
161 |
|
---|
162 | (defmacro monom-contract (k m)
|
---|
163 | `(subseq ,m ,k))
|
---|
164 |
|
---|
165 | (defun monom-exponents (m)
|
---|
166 | (declare (type monom m))
|
---|
167 | (coerce m 'list))
|
---|