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