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@ 995

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

* empty log message *

File size: 12.8 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 :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 (list (cons (make-list n :initial-element 0) number-or-poly))
78 number-or-poly))
79
80(defun $poly+ (p q n order ring)
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+ (convert-number p n) (convert-number q n) order ring))
85
86(defun $poly- (p q n order ring)
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- (convert-number p n) (convert-number q n) order ring))
91
92(defun $minus-poly (p n ring)
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 (minus-poly (convert-number p n) ring))
97
98(defun $poly* (p q n order ring)
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* (convert-number p n) (convert-number q n) order ring))
103
104(defun $poly/ (p q ring)
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 (common-lisp:/ q) p ring)))
110
111(defun $poly-expt (p l n order ring)
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 (convert-number p n) l order ring))
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 (expr vars order ring
254 &aux
255 (n (length vars))
256 (basis (monom-basis (length vars))))
257 "Evaluate an expression EXPR as polynomial by substituting operators
258+ - * expt with corresponding polynomial operators and variables VARS
259with monomials (1 0 ... 0), (0 1 ... 0) etc. We use special versions
260of binary operators $poly+, $poly-, $minus-poly, $poly* and $poly-expt
261which work like the corresponding functions in the POLY package, but
262accept scalars as arguments as well."
263 (cond
264 ((numberp expr)
265 (cond
266 ((zerop expr) NIL)
267 (t (list (cons (make-list n :initial-element 0) expr)))))
268 ((symbolp expr)
269 (nth (position expr vars) basis))
270 ((consp expr)
271 (case (car expr)
272 (expt
273 (if (= (length expr) 3)
274 ($poly-expt (poly-eval-1 (cadr expr) vars order ring)
275 (caddr expr)
276 n order ring)
277 (error "Too many arguments to EXPT")))
278 (/
279 (if (and (= (length expr) 3)
280 (numberp (caddr expr)))
281 ($poly/ (cadr expr) (caddr expr) ring)
282 (error "The second argument to / must be a number")))
283 (otherwise
284 (let ((r (mapcar
285 #'(lambda (e) (poly-eval-1 e vars order ring))
286 (cdr expr))))
287 (ecase (car expr)
288 (+ (reduce #'(lambda (p q) ($poly+ p q n order ring)) r))
289 (-
290 (if (endp (cdr r))
291 ($minus-poly (car r) n ring)
292 ($poly- (car r)
293 (reduce #'(lambda (p q) ($poly+ p q n order ring)) (cdr r))
294 n
295 order ring)))
296 (*
297 (reduce #'(lambda (p q) ($poly* p q n order ring)) r))
298 )))))))
299
300
301(defun poly-eval (expr vars &optional (order #'lex>) (ring *coefficient-ring*))
302 "Evaluate an expression EXPR, which should be a polynomial
303expression or a list of polynomial expressions (a list of expressions
304marked by prepending keyword :[ to it) given in lisp prefix notation,
305in variables VARS, which should be a list of symbols. The result of
306the evaluation is a polynomial or a list of polynomials (marked by
307prepending symbol '[) in the internal alist form. This evaluator is
308used by the PARSE package to convert input from strings directly to
309internal form."
310 (cond
311 ((numberp expr)
312 (unless (zerop expr)
313 (list (cons (make-list (length vars) :initial-element 0) expr))))
314 ((or (symbolp expr) (not (eq (car expr) :[)))
315 (poly-eval-1 expr vars order ring))
316 (t (cons '[ (mapcar #'(lambda (p) (poly-eval-1 p vars order ring)) (rest expr))))))
317
318
319;; Return the standard basis of the monomials in n variables
320(defun monom-basis (n &aux
321 (basis
322 (copy-tree
323 (make-list n :initial-element
324 (list (cons
325 (make-list
326 n
327 :initial-element 0)
328 1))))))
329 "Generate a list of monomials ((1 0 ... 0) (0 1 0 ... 0) ... (0 0 ... 1)
330which correspond to linear monomials X1, X2, ... XN."
331 (dotimes (i n basis)
332 (setf (elt (caar (elt basis i)) i) 1)))
Note: See TracBrowser for help on using the repository browser.