close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

source: branches/f4grobner/parse.lisp@ 1013

Last change on this file since 1013 was 1013, checked in by Marek Rychlik, 9 years ago

* empty log message *

File size: 13.1 KB
Line 
1;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
2
3;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4;;;
5;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
6;;;
7;;; This program is free software; you can redistribute it and/or modify
8;;; it under the terms of the GNU General Public License as published by
9;;; the Free Software Foundation; either version 2 of the License, or
10;;; (at your option) any later version.
11;;;
12;;; This program is distributed in the hope that it will be useful,
13;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
14;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15;;; GNU General Public License for more details.
16;;;
17;;; You should have received a copy of the GNU General Public License
18;;; along with this program; if not, write to the Free Software
19;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20;;;
21;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
22
23;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
24;;
25;; Parser of infix notation. This package enables input
26;; of polynomials in human-readable notation outside of Maxima,
27;; which is very useful for debugging.
28;;
29;; NOTE: This package is adapted from CGBLisp.
30;;
31;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
32
33(defpackage "PARSE"
34 (:use :cl :order :monom :term :polynomial :ring)
35 (:export "PARSE PARSE-TO-ALIST" "PARSE-STRING-TO-ALIST"
36 "PARSE-TO-SORTED-ALIST" "PARSE-STRING-TO-SORTED-ALIST" "^" "["
37 "POLY-EVAL"
38 ))
39
40(in-package "PARSE")
41
42(proclaim '(optimize (speed 0) (space 0) (safety 3) (debug 3)))
43
44;; The function PARSE yields the representations as above. The two functions
45;; PARSE-TO-ALIST and PARSE-STRING-TO-ALIST parse polynomials to the alist
46;; representations. For example
47;;
48;; >(parse)x^2-y^2+(-4/3)*u^2*w^3-5 --->
49;; (+ (* 1 (^ X 2)) (* -1 (^ Y 2)) (* -4/3 (^ U 2) (^ W 3)) (* -5))
50;;
51;; >(parse-to-alist '(x y u w))x^2-y^2+(-4/3)*u^2*w^3-5 --->
52;; (((0 0 2 3) . -4/3) ((0 2 0 0) . -1) ((2 0 0 0) . 1) ((0 0 0 0) . -5))
53;;
54;; >(parse-string-to-alist "x^2-y^2+(-4/3)*u^2*w^3-5" '(x y u w)) --->
55;; (((0 0 2 3) . -4/3) ((0 2 0 0) . -1) ((2 0 0 0) . 1) ((0 0 0 0) . -5))
56;;
57;; >(parse-string-to-alist "[x^2-y^2+(-4/3)*u^2*w^3-5,y]" '(x y u w))
58;; ([ (((0 0 2 3) . -4/3) ((0 2 0 0) . -1) ((2 0 0 0) . 1)
59;; ((0 0 0 0) . -5))
60;; (((0 1 0 0) . 1)))
61;; The functions PARSE-TO-SORTED-ALIST and PARSE-STRING-TO-SORTED-ALIST
62;; in addition sort terms by the predicate defined in the ORDER package
63;; For instance:
64;; >(parse-to-sorted-alist '(x y u w))x^2-y^2+(-4/3)*u^2*w^3-5
65;; (((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 2 3) . -4/3) ((0 0 0 0) . -5))
66;; >(parse-to-sorted-alist '(x y u w) t #'grlex>)x^2-y^2+(-4/3)*u^2*w^3-5
67;; (((0 0 2 3) . -4/3) ((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 0 0) . -5))
68
69;;(eval-when (compile)
70;; (proclaim '(optimize safety)))
71
72(defun convert-number (number-or-poly n)
73 "Returns NUMBER-OR-POLY, if it is a polynomial. If it is a number,
74it converts it to the constant monomial in N variables. If the result
75is a number then convert it to a polynomial in N variables."
76 (if (numberp number-or-poly)
77 (make-poly-from-termlist (list (make-term (make-monom :dimension n) number-or-poly)))
78 number-or-poly))
79
80(defun $poly+ (ring-and-order p q n)
81 "Add two polynomials P and Q, where each polynomial is either a
82numeric constant or a polynomial in internal representation. If the
83result is a number then convert it to a polynomial in N variables."
84 (poly-add ring-and-order (convert-number p n) (convert-number q n)))
85
86(defun $poly- (ring-and-order p q n)
87 "Subtract two polynomials P and Q, where each polynomial is either a
88numeric constant or a polynomial in internal representation. If the
89result is a number then convert it to a polynomial in N variables."
90 (poly-sub ring-and-order (convert-number p n) (convert-number q n)))
91
92(defun $minus-poly (ring p n)
93 "Negation of P is a polynomial is either a numeric constant or a
94polynomial in internal representation. If the result is a number then
95convert it to a polynomial in N variables."
96 (poly-uminus ring (convert-number p n)))
97
98(defun $poly* (ring-and-order p q n)
99 "Multiply two polynomials P and Q, where each polynomial is either a
100numeric constant or a polynomial in internal representation. If the
101result is a number then convert it to a polynomial in N variables."
102 (poly-mul ring-and-order (convert-number p n) (convert-number q n)))
103
104(defun $poly/ (ring p q)
105 "Divide a polynomials P which is either a numeric constant or a
106polynomial in internal representation, by a number Q."
107 (if (numberp p)
108 (common-lisp:/ p q)
109 (scalar-times-poly ring (common-lisp:/ q) p)))
110
111(defun $poly-expt (ring-and-order p l n)
112 "Raise polynomial P, which is a polynomial in internal
113representation or a numeric constant, to power L. If P is a number,
114convert the result to a polynomial in N variables."
115 (poly-expt ring-and-order (convert-number p n) l))
116
117(defun parse (&optional stream)
118 "Parser of infis expressions with integer/rational coefficients
119The parser will recognize two kinds of polynomial expressions:
120
121- polynomials in fully expanded forms with coefficients
122 written in front of symbolic expressions; constants can be optionally
123 enclosed in (); for example, the infix form
124 X^2-Y^2+(-4/3)*U^2*W^3-5
125 parses to
126 (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))
127
128- lists of polynomials; for example
129 [X-Y, X^2+3*Z]
130 parses to
131 (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
132 where the first symbol [ marks a list of polynomials.
133
134-other infix expressions, for example
135 [(X-Y)*(X+Y)/Z,(X+1)^2]
136parses to:
137 (:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
138Currently this function is implemented using M. Kantrowitz's INFIX package."
139 (read-from-string
140 (concatenate 'string
141 "#I("
142 (with-output-to-string (s)
143 (loop
144 (multiple-value-bind (line eof)
145 (read-line stream t)
146 (format s "~A" line)
147 (when eof (return)))))
148 ")")))
149
150;; Translate output from parse to a pure list form
151;; assuming variables are VARS
152
153(defun alist-form (plist vars)
154 "Translates an expression PLIST, which should be a list of polynomials
155in variables VARS, to an alist representation of a polynomial.
156It returns the alist. See also PARSE-TO-ALIST."
157 (cond
158 ((endp plist) nil)
159 ((eql (first plist) '[)
160 (cons '[ (mapcar #'(lambda (x) (alist-form x vars)) (rest plist))))
161 (t
162 (assert (eql (car plist) '+))
163 (alist-form-1 (rest plist) vars))))
164
165(defun alist-form-1 (p vars
166 &aux (ht (make-hash-table
167 :test #'equal :size 16))
168 stack)
169 (dolist (term p)
170 (assert (eql (car term) '*))
171 (incf (gethash (powers (cddr term) vars) ht 0) (second term)))
172 (maphash #'(lambda (key value) (unless (zerop value)
173 (push (cons key value) stack))) ht)
174 stack)
175
176(defun powers (monom vars
177 &aux (tab (pairlis vars (make-list (length vars)
178 :initial-element 0))))
179 (dolist (e monom (mapcar #'(lambda (v) (cdr (assoc v tab))) vars))
180 (assert (equal (first e) '^))
181 (assert (integerp (third e)))
182 (assert (= (length e) 3))
183 (let ((x (assoc (second e) tab)))
184 (if (null x) (error "Variable ~a not in the list of variables."
185 (second e))
186 (incf (cdr x) (third e))))))
187
188
189;; New implementation based on the INFIX package of Mark Kantorowitz
190(defun parse-to-alist (vars &optional stream)
191 "Parse an expression already in prefix form to an association list form
192according to the internal CGBlisp polynomial syntax: a polynomial is an
193alist of pairs (MONOM . COEFFICIENT). For example:
194 (WITH-INPUT-FROM-STRING (S \"X^2-Y^2+(-4/3)*U^2*W^3-5\")
195 (PARSE-TO-ALIST '(X Y U W) S))
196evaluates to
197(((0 0 2 3) . -4/3) ((0 2 0 0) . -1) ((2 0 0 0) . 1) ((0 0 0 0) . -5))"
198 (poly-eval (parse stream) vars))
199
200
201(defun parse-string-to-alist (str vars)
202 "Parse string STR and return a polynomial as a sorted association
203list of pairs (MONOM . COEFFICIENT). For example:
204(parse-string-to-alist \"[x^2-y^2+(-4/3)*u^2*w^3-5,y]\" '(x y u w))
205 ([ (((0 0 2 3) . -4/3) ((0 2 0 0) . -1) ((2 0 0 0) . 1)
206 ((0 0 0 0) . -5))
207 (((0 1 0 0) . 1)))
208The functions PARSE-TO-SORTED-ALIST and PARSE-STRING-TO-SORTED-ALIST
209sort terms by the predicate defined in the ORDER package."
210 (with-input-from-string (stream str)
211 (parse-to-alist vars stream)))
212
213
214(defun parse-to-sorted-alist (vars &optional (order #'lex>) (stream t))
215 "Parses streasm STREAM and returns a polynomial represented as
216a sorted alist. For example:
217(WITH-INPUT-FROM-STRING (S \"X^2-Y^2+(-4/3)*U^2*W^3-5\")
218 (PARSE-TO-SORTED-ALIST '(X Y U W) S))
219returns
220(((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 2 3) . -4/3) ((0 0 0 0) . -5))
221and
222(WITH-INPUT-FROM-STRING (S \"X^2-Y^2+(-4/3)*U^2*W^3-5\")
223 (PARSE-TO-SORTED-ALIST '(X Y U W) T #'GRLEX>) S)
224returns
225(((0 0 2 3) . -4/3) ((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 0 0) . -5))"
226 (sort-poly (parse-to-alist vars stream) order))
227
228(defun parse-string-to-sorted-alist (str vars &optional (order #'lex>))
229 "Parse a string to a sorted alist form, the internal representation
230of polynomials used by our system."
231 (with-input-from-string (stream str)
232 (parse-to-sorted-alist vars order stream)))
233
234(defun sort-poly-1 (p order)
235 "Sort the terms of a single polynomial P using an admissible monomial order ORDER.
236Returns the sorted polynomial. Destructively modifies P."
237 (sort p order :key #'first))
238
239;; Sort a polynomial or polynomial list
240(defun sort-poly (poly-or-poly-list &optional (order #'lex>))
241 "Sort POLY-OR-POLY-LIST, which could be either a single polynomial
242or a list of polynomials in internal alist representation, using
243admissible monomial order ORDER. Each polynomial is sorted using
244SORT-POLY-1."
245 (cond
246 ((eql poly-or-poly-list :syntax-error) nil)
247 ((null poly-or-poly-list) nil)
248 ((eql (car poly-or-poly-list) '[)
249 (cons '[ (mapcar #'(lambda (p) (sort-poly-1 p order))
250 (rest poly-or-poly-list))))
251 (t (sort-poly-1 poly-or-poly-list order))))
252
253(defun poly-eval-1 (ring-and-order expr vars
254 &aux
255 (ring (ro-ring ring-and-order))
256 (n (length vars))
257 (basis (monom-basis (length vars))))
258 "Evaluate an expression EXPR as polynomial by substituting operators
259+ - * expt with corresponding polynomial operators and variables VARS
260with monomials (1 0 ... 0), (0 1 ... 0) etc. We use special versions
261of binary operators $poly+, $poly-, $minus-poly, $poly* and $poly-expt
262which work like the corresponding functions in the POLY package, but
263accept scalars as arguments as well."
264 (cond
265 ((numberp expr)
266 (cond
267 ((zerop expr) NIL)
268 (t (make-poly-from-termlist (list (make-monom :dimension n) expr)))))
269 ((symbolp expr)
270 (make-monom :initial-exponents (nth (position expr vars) basis)) )
271 ((consp expr)
272 (case (car expr)
273 (expt
274 (if (= (length expr) 3)
275 ($poly-expt (poly-eval-1 ring-and-order (cadr expr) vars)
276 (caddr expr)
277 n)
278 (error "Too many arguments to EXPT")))
279 (/
280 (if (and (= (length expr) 3)
281 (numberp (caddr expr)))
282 ($poly/ ring (cadr expr) (caddr expr))
283 (error "The second argument to / must be a number")))
284 (otherwise
285 (let ((r (mapcar
286 #'(lambda (e) (poly-eval-1 ring-and-order e vars))
287 (cdr expr))))
288 (ecase (car expr)
289 (+ (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) r))
290 (-
291 (if (endp (cdr r))
292 ($minus-poly ring (car r) n)
293 ($poly- ring-and-order
294 (car r)
295 (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) (cdr r))
296 n)))
297 (*
298 (reduce #'(lambda (p q) ($poly* ring-and-order p q n)) r))
299 )))))))
300
301
302(defun poly-eval (expr vars
303 &optional (order #'lex>) (ring *coefficient-ring*)
304 &aux
305 (ring-and-order (make-ring-and-order ring order)))
306 "Evaluate an expression EXPR, which should be a polynomial
307expression or a list of polynomial expressions (a list of expressions
308marked by prepending keyword :[ to it) given in lisp prefix notation,
309in variables VARS, which should be a list of symbols. The result of
310the evaluation is a polynomial or a list of polynomials (marked by
311prepending symbol '[) in the internal alist form. This evaluator is
312used by the PARSE package to convert input from strings directly to
313internal form."
314 (cond
315 ((numberp expr)
316 (unless (zerop expr)
317 (make-poly-from-termlist
318 (list (make-term (make-monom :dimension (length vars)) expr)))))
319 ((or (symbolp expr) (not (eq (car expr) :[)))
320 (poly-eval-1 ring-and-order expr vars))
321 (t (cons '[ (mapcar #'(lambda (p) (poly-eval-1 ring-and-order p vars)) (rest expr))))))
322
323
324;; Return the standard basis of the monomials in n variables
325(defun monom-basis (n &aux
326 (basis
327 (copy-tree
328 (make-list n :initial-element
329 (list (cons
330 (make-list
331 n
332 :initial-element 0)
333 1))))))
334 "Generate a list of monomials ((1 0 ... 0) (0 1 0 ... 0) ... (0 0 ... 1)
335which correspond to linear monomials X1, X2, ... XN."
336 (dotimes (i n basis)
337 (setf (elt (caar (elt basis i)) i) 1)))
Note: See TracBrowser for help on using the repository browser.