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/.junk/poly-eval.lisp

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

* empty log message *

  • Property svn:mime-type set to application/x-elc
File size: 8.3 KB
Line 
1;; Translate output from parse to a pure list form
2;; assuming variables are VARS
3
4(defun alist-form (plist vars)
5 "Translates an expression PLIST, which should be a list of polynomials
6in variables VARS, to an alist representation of a polynomial.
7It returns the alist. See also PARSE-TO-ALIST."
8 (cond
9 ((endp plist) nil)
10 ((eql (first plist) '[)
11 (cons '[ (mapcar #'(lambda (x) (alist-form x vars)) (rest plist))))
12 (t
13 (assert (eql (car plist) '+))
14 (alist-form-1 (rest plist) vars))))
15
16(defun alist-form-1 (p vars
17 &aux (ht (make-hash-table
18 :test #'equal :size 16))
19 stack)
20 (dolist (term p)
21 (assert (eql (car term) '*))
22 (incf (gethash (powers (cddr term) vars) ht 0) (second term)))
23 (maphash #'(lambda (key value) (unless (zerop value)
24 (push (cons key value) stack))) ht)
25 stack)
26
27(defun powers (monom vars
28 &aux (tab (pairlis vars (make-list (length vars)
29 :initial-element 0))))
30 (dolist (e monom (mapcar #'(lambda (v) (cdr (assoc v tab))) vars))
31 (assert (equal (first e) '^))
32 (assert (integerp (third e)))
33 (assert (= (length e) 3))
34 (let ((x (assoc (second e) tab)))
35 (if (null x) (error "Variable ~a not in the list of variables."
36 (second e))
37 (incf (cdr x) (third e))))))
38
39
40
41(defun convert-number (number-or-poly n)
42 "Returns NUMBER-OR-POLY, if it is a polynomial. If it is a number,
43it converts it to the constant monomial in N variables. If the result
44is a number then convert it to a polynomial in N variables."
45 (if (numberp number-or-poly)
46 (make-poly-from-termlist (list (make-term :monom (make-monom :dimension n) :coeff number-or-poly)))
47 number-or-poly))
48
49(defun $poly+ (ring-and-order p q n)
50 "Add two polynomials P and Q, where each polynomial is either a
51numeric constant or a polynomial in internal representation. If the
52result is a number then convert it to a polynomial in N variables."
53 (poly-add ring-and-order (convert-number p n) (convert-number q n)))
54
55(defun $poly- (ring-and-order p q n)
56 "Subtract two polynomials P and Q, where each polynomial is either a
57numeric constant or a polynomial in internal representation. If the
58result is a number then convert it to a polynomial in N variables."
59 (poly-sub ring-and-order (convert-number p n) (convert-number q n)))
60
61(defun $minus-poly (ring p n)
62 "Negation of P is a polynomial is either a numeric constant or a
63polynomial in internal representation. If the result is a number then
64convert it to a polynomial in N variables."
65 (poly-uminus ring (convert-number p n)))
66
67(defun $poly* (ring-and-order p q n)
68 "Multiply two polynomials P and Q, where each polynomial is either a
69numeric constant or a polynomial in internal representation. If the
70result is a number then convert it to a polynomial in N variables."
71 (poly-mul ring-and-order (convert-number p n) (convert-number q n)))
72
73(defun $poly/ (ring p q)
74 "Divide a polynomials P which is either a numeric constant or a
75polynomial in internal representation, by a number Q."
76 (if (numberp p)
77 (common-lisp:/ p q)
78 (scalar-times-poly ring (common-lisp:/ q) p)))
79
80(defun $poly-expt (ring-and-order p l n)
81 "Raise polynomial P, which is a polynomial in internal
82representation or a numeric constant, to power L. If P is a number,
83convert the result to a polynomial in N variables."
84 (poly-expt ring-and-order (convert-number p n) l))
85
86
87(defun variable-basis (ring n &aux (basis (make-list n)))
88 "Generate a list of polynomials X[i], i=0,1,...,N-1."
89 (dotimes (i n basis)
90 (setf (elt basis i) (make-poly-variable ring n i))))
91
92
93(defun poly-eval-1 (expr vars &optional (ring *ring-of-integers*) (order #'lex>)
94 &aux
95 (ring-and-order (make-ring-and-order :ring ring :order order))
96 (n (length vars))
97 (basis (variable-basis ring (length vars))))
98 "Evaluate an expression EXPR as polynomial by substituting operators
99+ - * expt with corresponding polynomial operators and variables VARS
100with the corresponding polynomials in internal form. We use special
101versions of binary operators $poly+, $poly-, $minus-poly, $poly* and
102$poly-expt which work like the corresponding functions in the POLY
103package, but accept scalars as arguments as well. The result is a
104polynomial in internal form. This operation is somewhat similar to
105the function EXPAND in CAS."
106 (cond
107 ((numberp expr)
108 (cond
109 ((zerop expr) NIL)
110 (t (make-poly-from-termlist (list (make-term :monom (make-monom :dimension n) :coeff expr))))))
111 ((symbolp expr)
112 (nth (position expr vars) basis))
113 ((consp expr)
114 (case (car expr)
115 (expt
116 (if (= (length expr) 3)
117 ($poly-expt ring-and-order
118 (poly-eval-1 (cadr expr) vars ring order)
119 (caddr expr)
120 n)
121 (error "Too many arguments to EXPT")))
122 (/
123 (if (and (= (length expr) 3)
124 (numberp (caddr expr)))
125 ($poly/ ring (cadr expr) (caddr expr))
126 (error "The second argument to / must be a number")))
127 (otherwise
128 (let ((r (mapcar
129 #'(lambda (e) (poly-eval-1 e vars ring order))
130 (cdr expr))))
131 (ecase (car expr)
132 (+ (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) r))
133 (-
134 (if (endp (cdr r))
135 ($minus-poly ring (car r) n)
136 ($poly- ring-and-order
137 (car r)
138 (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) (cdr r))
139 n)))
140 (*
141 (reduce #'(lambda (p q) ($poly* ring-and-order p q n)) r))
142 )))))))
143
144
145
146(defun poly-eval (expr vars &optional (order #'lex>) (ring *ring-of-integers*))
147 "Evaluate an expression EXPR, which should be a polynomial
148expression or a list of polynomial expressions (a list of expressions
149marked by prepending keyword :[ to it) given in Lisp prefix notation,
150in variables VARS, which should be a list of symbols. The result of
151the evaluation is a polynomial or a list of polynomials (marked by
152prepending symbol '[) in the internal alist form. This evaluator is
153used by the PARSE package to convert input from strings directly to
154internal form."
155 (cond
156 ((numberp expr)
157 (unless (zerop expr)
158 (make-poly-from-termlist
159 (list (make-term :monom (make-monom :dimension (length vars)) :coeff expr)))))
160 ((or (symbolp expr) (not (eq (car expr) :[)))
161 (poly-eval-1 expr vars ring order))
162 (t (cons '[ (mapcar #'(lambda (p) (poly-eval-1 p vars ring order)) (rest expr))))))
163
164
165(defun parse-string-to-alist (str vars)
166 "Parse string STR and return a polynomial as a sorted association
167list of pairs (MONOM . COEFFICIENT). For example:
168(parse-string-to-alist \"[x^2-y^2+(-4/3)*u^2*w^3-5,y]\" '(x y u w))
169 ([ (((0 0 2 3) . -4/3) ((0 2 0 0) . -1) ((2 0 0 0) . 1)
170 ((0 0 0 0) . -5))
171 (((0 1 0 0) . 1)))
172The functions PARSE-TO-SORTED-ALIST and PARSE-STRING-TO-SORTED-ALIST
173sort terms by the predicate defined in the ORDER package."
174 (with-input-from-string (stream str)
175 (parse-to-alist vars stream)))
176
177
178(defun parse-to-sorted-alist (vars &optional (order #'lex>) (stream t))
179 "Parses streasm STREAM and returns a polynomial represented as
180a sorted alist. For example:
181(WITH-INPUT-FROM-STRING (S \"X^2-Y^2+(-4/3)*U^2*W^3-5\")
182 (PARSE-TO-SORTED-ALIST '(X Y U W) S))
183returns
184(((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 2 3) . -4/3) ((0 0 0 0) . -5))
185and
186(WITH-INPUT-FROM-STRING (S \"X^2-Y^2+(-4/3)*U^2*W^3-5\")
187 (PARSE-TO-SORTED-ALIST '(X Y U W) T #'GRLEX>) S)
188returns
189(((0 0 2 3) . -4/3) ((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 0 0) . -5))"
190 (sort-poly (parse-to-alist vars stream) order))
191
192(defun parse-string-to-sorted-alist (str vars &optional (order #'lex>))
193 "Parse a string to a sorted alist form, the internal representation
194of polynomials used by our system."
195 (with-input-from-string (stream str)
196 (parse-to-sorted-alist vars order stream)))
197
198(defun sort-poly-1 (p order)
199 "Sort the terms of a single polynomial P using an admissible monomial order ORDER.
200Returns the sorted polynomial. Destructively modifies P."
201 (sort p order :key #'first))
202
203;; Sort a polynomial or polynomial list
204(defun sort-poly (poly-or-poly-list &optional (order #'lex>))
205 "Sort POLY-OR-POLY-LIST, which could be either a single polynomial
206or a list of polynomials in internal alist representation, using
207admissible monomial order ORDER. Each polynomial is sorted using
208SORT-POLY-1."
209 (cond
210 ((eql poly-or-poly-list :syntax-error) nil)
211 ((null poly-or-poly-list) nil)
212 ((eql (car poly-or-poly-list) '[)
213 (cons '[ (mapcar #'(lambda (p) (sort-poly-1 p order))
214 (rest poly-or-poly-list))))
215 (t (sort-poly-1 poly-or-poly-list order))))
216
Note: See TracBrowser for help on using the repository browser.