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

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