share/linearalgebra/mring.lisp

   1 ;; A Maxima ring structure
   2 ;; Copyright (C) 2005, 2007, 2021, 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 (intl:gettext "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 (intl:gettext "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 (intl:gettext "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 (intl:gettext "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)) (= 0 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 (intl:gettext "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 (intl:gettext "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)) (not (invertible-matrixp b '$noncommutingring)))
 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) (not (invertible-matrixp s '$noncommutingring)))
 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 (intl:gettext "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))))