1;; A Maxima ring structure
2;; Copyright (C) 2005, 2007, Barton Willis
3
4;; Barton Willis
5;; Department of Mathematics
6;; University of Nebraska at Kearney
7;; Kearney NE 68847
8;; willisb@unk.edu
9
10;; This source code is licensed under the terms of the Lisp Lesser
11;; GNU Public License (LLGPL). The LLGPL consists of a preamble, published
12;; by Franz Inc. (http://opensource.franz.com/preamble.html), and the GNU
13;; Library General Public License (LGPL), version 2, or (at your option)
14;; any later version.  When the preamble conflicts with the LGPL,
15;; the preamble takes precedence.
16
17;; This library is distributed in the hope that it will be useful,
18;; but WITHOUT ANY WARRANTY; without even the implied warranty of
19;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20;; Library General Public License for details.
21
22;; You should have received a copy of the GNU Library General Public
23;; License along with this library; if not, write to the
24;; Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
25;; Boston, MA  02110-1301, USA.
26
27;; Let's have version numbers 1,2,3,...
28
29(eval-when (:compile-toplevel :load-toplevel :execute)
30  ($put '$mring 1 '$version))
31
32;; (1) In maxima-grobner.lisp, there is a structure 'ring.'
33
34;; (2) Some functions in this structure, for example 'great' might
35;;     not be defined for a ring; when this is the case, a function
36;;     can signal an error.
37
38;; (3) Floating point addition isn't associative; so a mring needn't
39;;     be a ring.  But a mring is 'close' to being a ring.
40
41;; Description of the mring fields:
42
43(defstruct mring
44  name
45  coerce-to-lisp-float
46  abs
47  great
48  add
49  div
50  rdiv
51  reciprocal
52  mult
53  sub
54  negate
55  psqrt
56  add-id
57  mult-id
58  fzerop
59  adjoint
60  maxima-to-mring
61  mring-to-maxima)
62
63(eval-when
64#-gcl (:compile-toplevel :load-toplevel :execute)
65#+gcl (compile load eval)
66  (defmvar $%mrings `((mlist) $floatfield $complexfield $rationalfield $crering
67                       $generalring $bigfloatfield $runningerror $noncommutingring ))
68  (defvar *gf-rings* '(gf-coeff-ring gf-ring ef-ring)) )
69
70(defun $require_ring (ringname pos fun)
71  (if (or ($member ringname $%mrings) (member ringname *gf-rings*))
72    (get ringname 'ring)
73    (merror (intl:gettext "The ~:M argument of the function '~:M' must be the name of a ring") pos fun)))
74
75
76(defparameter *floatfield*
77  (make-mring
78   :name '$floatfield
79   :coerce-to-lisp-float #'cl:identity
80   :abs #'abs
81   :great #'>
82   :add #'+
83   :div #'/
84   :rdiv #'/
85   :reciprocal #'/
86   :mult #'*
87   :sub #'-
88   :negate #'-
89   :psqrt #'(lambda (s) (if (>= s 0) (cl:sqrt s) nil))
90   :add-id #'(lambda () 0.0)
91   :mult-id #'(lambda () 1.0)
92   :fzerop #'(lambda (s)  (< (abs s) (* 4 flonum-epsilon)))
93   :adjoint #'cl:identity
94   :mring-to-maxima #'cl:identity
95   :maxima-to-mring #'(lambda (s)
96			(setq s ($float s))
97			(if (floatp s) s (merror "Unable to convert ~:M to a long float" s)))))
98
99(setf (get '$floatfield 'ring) *floatfield*)
100
101(defparameter *complexfield*
102  (make-mring
103   :name '$complexfield
104   :coerce-to-lisp-float #'cl:identity
105   :abs #'abs
106   :great #'>
107   :add #'+
108   :div #'/
109   :rdiv #'/
110   :reciprocal #'/
111   :mult #'*
112   :sub #'-
113   :negate #'-
114   :psqrt #'(lambda (s) (if (and (= 0 (imagpart s)) (>= (realpart s) 0)) (cl:sqrt s) nil))
115   :add-id #'(lambda () 0.0)
116   :mult-id #'(lambda () 1.0)
117   :fzerop #'(lambda (s) (< (abs s) (* 4 flonum-epsilon)))
118   :adjoint #'cl:conjugate
119   :mring-to-maxima #'(lambda (s) (add (realpart s) (mult '$%i (imagpart s)))) ;; was complexify
120   :maxima-to-mring #'(lambda (s)
121			(progn
122			  (setq s (coerce-expr-to-clcomplex ($rectform (meval s))))
123			  (if (complexp s)
124			      s
125			      (merror "Unable to convert ~:M to a complex long float" s))))))
126
127(defun coerce-expr-to-clcomplex (s)
128  (complex (funcall (coerce-float-fun ($realpart s))) (funcall (coerce-float-fun ($imagpart s)))))
129
130(setf (get '$complexfield 'ring) *complexfield*)
131
132(defparameter *rationalfield*
133  (make-mring
134   :name '$rationalfield
135   :coerce-to-lisp-float #'(lambda (s) ($float s))
136   :abs #'abs
137   :great #'>
138   :add #'+
139   :div #'/
140   :rdiv #'/
141   :reciprocal #'/
142   :mult #'*
143   :sub #'-
144   :negate #'-
145   :psqrt #'(lambda (s) (let ((x))
146			  (cond ((>= s 0)
147				 (setq x (isqrt (numerator s)))
148				 (setq x (/ x (isqrt (denominator s))))
149				 (if (= s (* x x)) x nil))
150				(t nil))))
151   :add-id #'(lambda () 0)
152   :mult-id #'(lambda () 1)
153   :fzerop #'(lambda (s) (= s 0))
154   :adjoint #'cl:identity
155   :mring-to-maxima #'(lambda (s) (simplify `((rat) ,(numerator s) ,(denominator s))))
156   :maxima-to-mring
157   #'(lambda (s)
158       (if (or (floatp s) ($bfloatp s)) (setq s ($rationalize s)))
159       (if ($ratnump s) (if (integerp s) s (/ ($num s) ($denom s)))
160	 (merror "Unable to convert ~:M to a rational number" s)))))
161
162(setf (get '$rationalfield 'ring) *rationalfield*)
163
164(defparameter *crering*
165  (make-mring
166   :name '$crering
167   :coerce-to-lisp-float nil
168   :abs #'(lambda (s) (simplify (mfuncall '$cabs s)))
169   :great #'(lambda (a b) (declare (ignore a)) (eq t (meqp b 0)))
170   :add #'add
171   :div #'div
172   :rdiv #'div
173   :reciprocal #'(lambda (s) (div 1 s))
174   :mult #'mult
175   :sub #'sub
176   :negate #'(lambda (s) (mult -1 s))
177   :psqrt #'(lambda (s) (if (member (csign ($ratdisrep s)) `($pos $pz $zero)) (take '(%sqrt) s) nil))
178   :add-id #'(lambda () 0)
179   :mult-id #'(lambda () 1)
180   :fzerop #'(lambda (s) (eq t (meqp s 0)))
181   :adjoint #'(lambda (s) (take '($conjugate) s))
182   :mring-to-maxima #'(lambda (s) s)
183   :maxima-to-mring #'(lambda (s) ($rat s))))
184
185(setf (get '$crering 'ring) *crering*)
186
187(defparameter *generalring*
188  (make-mring
189   :name '$generalring
190   :coerce-to-lisp-float nil
191   :abs #'(lambda (s) (simplify (mfuncall '$cabs s)))
192   :great #'(lambda (a b) (declare (ignore a)) (eq t (meqp b 0)))
193   :add #'(lambda (a b) (add a b))
194   :div #'(lambda (a b) (div a b))
195   :rdiv #'(lambda (a b) (div a b))
196   :reciprocal #'(lambda (s) (div 1 s))
197   :mult #'(lambda (a b) (mult a b))
198   :sub #'(lambda (a b) (sub a b))
199   :negate #'(lambda (a) (mult -1 a))
200   :psqrt #'(lambda (s) (if (member (csign s) `($pos $pz $zero)) (take '(%sqrt) s) nil))
201   :add-id #'(lambda () 0)
202   :mult-id #'(lambda () 1)
203   :fzerop #'(lambda (s) (eq t (meqp (sratsimp s) 0)))
204   :adjoint #'(lambda (s) (take '($conjugate) s))
205   :mring-to-maxima #'(lambda (s) s)
206   :maxima-to-mring #'(lambda (s) s)))
207
208(setf (get '$generalring 'ring) *generalring*)
209
210(defparameter *bigfloatfield*
211  (make-mring
212   :name '$bigfloatfield
213   :coerce-to-lisp-float #'(lambda (s)
214			     (setq s ($rectform ($float s)))
215			     (complex ($realpart s) ($imagpart s)))
216
217   :abs #'(lambda (s) (simplify (mfuncall '$cabs s)))
218   :great #'mgrp
219   :add #'(lambda (a b) ($rectform (add a b)))
220   :div #'(lambda (a b) ($rectform (div a b)))
221   :rdiv #'(lambda (a b) ($rectform (div a b)))
222   :reciprocal #'(lambda (s) (div 1 s))
223   :mult #'(lambda (a b) ($rectform (mult a b)))
224   :sub #'(lambda (a b) ($rectform (sub a b)))
225   :negate #'(lambda (a) (mult -1 a))
226   :psqrt #'(lambda (s) (if (mlsp s 0) nil (take '(%sqrt) s)))
227   :add-id #'(lambda () bigfloatzero)
228   :mult-id #'(lambda () bigfloatone)
229   :fzerop #'(lambda (s) (like s bigfloatzero))
230   :adjoint #'cl:identity
231   :mring-to-maxima #'(lambda (s) s)
232   :maxima-to-mring #'(lambda (s)
233			(setq s ($rectform ($bfloat s)))
234			(if (or (eq s '$%i) (complex-number-p s 'bigfloat-or-number-p)) s
235				(merror "Unable to convert matrix entry to a big float")))))
236
237(setf (get '$bigfloatfield 'ring) *bigfloatfield*)
238
239;; --- *gf-rings* --- (used by src/numth.lisp) ------------------------------ ;;
240;;                                                                            ;;
241(defparameter *gf-coeff-ring*
242  (make-mring
243   :name 'gf-coeff-ring
244   :coerce-to-lisp-float nil
245   :abs #'gf-cmod
246   :great #'(lambda (a b) (declare (ignore a)) (null b))
247   :add #'gf-cplus-b
248   :div #'(lambda (a b) (gf-ctimes a (gf-cinv b)))
249   :rdiv #'(lambda (a b) (gf-ctimes a (gf-cinv b)))
250   :reciprocal #'gf-cinv
251   :mult #'gf-ctimes
252   :sub #'(lambda (a b) (gf-cplus-b a (gf-cminus-b b)))
253   :negate #'gf-cminus-b
254   :psqrt #'(lambda (a) (let ((rs (zn-nrt a 2 *gf-char*))) (when rs (car rs))))
255   :add-id #'(lambda () 0)
256   :mult-id #'(lambda () 1)
257   :fzerop #'(lambda (s) (= 0 s))
258   :adjoint nil
259   :mring-to-maxima #'cl:identity
260   :maxima-to-mring #'cl:identity ))
261;;
262(setf (get 'gf-coeff-ring 'ring) *gf-coeff-ring*)
263;;
264(defparameter *gf-ring*
265  (make-mring
266   :name 'gf-ring
267   :coerce-to-lisp-float nil
268   :abs #'gf-mod
269   :great #'(lambda (a b) (declare (ignore a)) (null b))
270   :add #'gf-plus
271   :div #'(lambda (a b) (gf-times a (gf-inv b *gf-red*) *gf-red*))
272   :rdiv #'(lambda (a b) (gf-times a (gf-inv b *gf-red*) *gf-red*))
273   :reciprocal #'(lambda (a) (gf-inv a *gf-red*))
274   :mult #'(lambda (a b) (gf-times a b *gf-red*))
275   :sub #'(lambda (a b) (gf-plus a (gf-minus b)))
276   :negate #'gf-minus
277   :psqrt #'(lambda (a)
278              (let ((rs (gf-nrt-exit (gf-nrt a 2 *gf-red* *gf-ord*))))
279                (when rs (cadr rs)) ))
280   :add-id #'(lambda () nil)
281   :mult-id #'(lambda () '(0 1))
282   :fzerop #'(lambda (s) (null s))
283   :adjoint nil
284   :mring-to-maxima #'gf-x2p
285   :maxima-to-mring #'gf-p2x ))
286;;
287(setf (get 'gf-ring 'ring) *gf-ring*)
288;;
289(defparameter *ef-ring*
290  (make-mring
291   :name 'ef-ring
292   :coerce-to-lisp-float nil
293   :abs #'gf-mod
294   :great #'(lambda (a b) (declare (ignore a)) (null b))
295   :add #'gf-plus
296   :div #'(lambda (a b) (gf-times a (gf-inv b *ef-red*) *ef-red*))
297   :rdiv #'(lambda (a b) (gf-times a (gf-inv b *ef-red*) *ef-red*))
298   :reciprocal #'(lambda (a) (gf-inv a *ef-red*))
299   :mult #'(lambda (a b) (gf-times a b *ef-red*))
300   :sub #'(lambda (a b) (gf-plus a (gf-minus b)))
301   :negate #'gf-minus
302   :psqrt #'(lambda (a)
303              (let ((rs (gf-nrt-exit (gf-nrt a 2 *ef-red* *ef-ord*))))
304                (when rs (cadr rs)) ))
305   :add-id #'(lambda () nil)
306   :mult-id #'(lambda () '(0 1))
307   :fzerop #'(lambda (s) (null s))
308   :adjoint nil
309   :mring-to-maxima #'gf-x2p
310   :maxima-to-mring #'gf-p2x ))
311
312(setf (get 'ef-ring 'ring) *ef-ring*)
313;;                                                                            ;;
314;; -------------------------------------------------------------------------- ;;
315
316(defun fp-abs (a)
317  (list (abs (first a)) (second a)))
318
319(defun fp+ (a b)
320  (cond ((= (first a) 0.0) b)
321	((= (first b) 0.0) a)
322	(t
323	 (let ((s (+ (first a) (first b))))
324	   (if (= 0.0 s) (merror "floating point divide by zero"))
325	   (list s (ceiling (+ 1
326			       (abs (/ (* (first a) (second a)) s))
327			       (abs (/ (* (first b) (second b)) s)))))))))
328
329(defun fp- (a b)
330  (cond ((= (first a) 0.0) (list (- (first b)) (second b)))
331	((= (first b) 0.0) a)
332	(t
333	 (let ((s (- (first a) (first b))))
334	   (if (= 0.0 s) (merror "floating point divide by zero"))
335	   (list s (ceiling (+ 1
336			       (abs (/ (* (first a) (second a)) s))
337			       (abs (/ (* (first b) (second b)) s)))))))))
338
339(defun fp* (a b)
340  (if (or (= (first a) 0.0) (= (first b) 0.0)) (list 0.0 0)
341    (list (* (first a) (first b)) (+ 1 (second a) (second b)))))
342
343(defun fp/ (a b)
344  (if (= (first a) 0) (list 0.0 0)
345    (list (/ (first a) (first b)) (+ 1 (second a) (second b)))))
346
347(defun $addmatrices(fn &rest m)
348  (mfuncall '$apply '$matrixmap `((mlist) ,fn ,@m)))
349
350(defparameter *runningerror*
351  (make-mring
352   :name '$runningerror
353   :coerce-to-lisp-float #'(lambda (s) (if (consp s) (first s) s))
354   :abs #'fp-abs
355   :great #'(lambda (a b) (> (first a) (first b)))
356   :add #'fp+
357   :div #'fp/
358   :rdiv #'fp/
359   :reciprocal #'(lambda (s) (fp/ (list 1 0) s))
360   :mult #'fp*
361   :sub #'fp-
362   :negate #'(lambda (s) (list (- (first s)) (second s)))
363   :psqrt #'(lambda (s) (if (> (first s) 0) (list (cl:sqrt (first s)) (+ 1 (second s))) nil))
364   :add-id #'(lambda () (list 0 0))
365   :mult-id #'(lambda () (list 1 0))
366   :fzerop #'(lambda (s) (like (first s) 0))
367   :adjoint #'cl:identity
368   :mring-to-maxima #'(lambda (s) `((mlist) ,@s))
369   :maxima-to-mring #'(lambda (s) (if ($listp s) (cdr s) (list ($float s) 1)))))
370
371(setf (get '$runningerror 'ring) *runningerror*)
372
373(defparameter *noncommutingring*
374  (make-mring
375   :name '$noncommutingring
376   :coerce-to-lisp-float nil
377   :abs #'(lambda (s) (simplify (mfuncall '$cabs s)))
378   :great #'(lambda (a b) (declare (ignore a)) (eq t (meqp b 0)))
379   :add #'(lambda (a b) (add a b))
380   :div #'(lambda (a b) (progn
381			  (let (($matrix_element_mult ".")
382				($matrix_element_transpose '$transpose))
383			    (setq b (if ($matrixp b) ($invert_by_lu b '$noncommutingring)
384				      (take '(mncexpt) b -1)))
385			    (take '(mnctimes) a b))))
386
387   :rdiv #'(lambda (a b) (progn
388			   (let (($matrix_element_mult ".")
389				 ($matrix_element_transpose '$transpose))
390			     (setq b (if ($matrixp b) ($invert_by_lu b '$noncommutingring)
391				       (take '(mncexpt) b -1)))
392			     (take  '(mnctimes) b a))))
393
394
395   :reciprocal #'(lambda (s) (progn
396			       (let (($matrix_element_mult ".")
397				     ($matrix_element_transpose '$transpose))
398				 (if ($matrixp s) ($invert_by_lu s '$noncommutingring)
399				   (take '(mncexpt) s -1)))))
400
401   :mult #'(lambda (a b) (progn
402			   (let (($matrix_element_mult ".")
403				 ($matrix_element_transpose '$transpose))
404			     (take  '(mnctimes) a b))))
405
406
407   :sub #'(lambda (a b) (sub a b))
408   :negate #'(lambda (a) (mult -1 a))
409   :add-id #'(lambda () 0)
410   :psqrt #'(lambda (s) (take '(%sqrt) s))
411   :mult-id #'(lambda () 1)
412   :fzerop #'(lambda (s) (eq t (meqp s 0)))
413   :adjoint #'(lambda (s) ($transpose (take '($conjugate) s)))
414   :mring-to-maxima #'cl:identity
415   :maxima-to-mring #'cl:identity))
416
417(setf (get '$noncommutingring 'ring) *noncommutingring*)
418
419(defun ring-eval (e fld)
420  (let ((fadd (mring-add fld))
421	(fnegate (mring-negate fld))
422	(fmult (mring-mult fld))
423	(fdiv (mring-div fld))
424	(fabs (mring-abs fld))
425	(fconvert (mring-maxima-to-mring fld)))
426
427    (cond ((or ($numberp e) (symbolp e))
428	   (funcall fconvert (meval e)))
429
430	  ;; I don't think an empty sum or product is possible here. If it is, append
431	  ;; the appropriate initial-value to reduce. Using the :inital-value isn't
432	  ;; a problem, but (fp* (a b) (1 0)) --> (a (+ b 1)).  A better value is
433	  ;; (fp* (a b) (1 0)) --> (a b).
434
435	  ((op-equalp e 'mplus)
436	   (reduce fadd (mapcar #'(lambda (s) (ring-eval s fld)) (margs e)) :from-end t))
437
438	  ((op-equalp e 'mminus)
439	   (funcall fnegate (ring-eval (first (margs e)) fld)))
440
441	  ((op-equalp e 'mtimes)
442	   (reduce fmult (mapcar #'(lambda (s) (ring-eval s fld)) (margs e)) :from-end t))
443
444	  ((op-equalp e 'mquotient)
445	   (funcall fdiv (ring-eval (first (margs e)) fld)(ring-eval (second (margs e)) fld)))
446
447	  ((op-equalp e 'mabs) (funcall fabs (ring-eval (first (margs e)) fld)))
448
449	  ((and (or (eq (mring-name fld) '$floatfield) (eq (mring-name fld) '$complexfield))
450		(consp e) (consp (car e)) (gethash (mop e) *flonum-op*))
451	   (apply (gethash (mop e) *flonum-op*) (mapcar #'(lambda (s) (ring-eval s fld)) (margs e))))
452
453	  (t (merror "Unable to evaluate ~:M in the ring '~:M'" e (mring-name fld))))))
454
455(defmspec $ringeval (e)
456  (let ((fld (get (or (car (member (nth 2 e) $%mrings)) '$generalring) 'ring)))
457    (funcall (mring-mring-to-maxima fld) (ring-eval (nth 1 e) fld))))
458