share/linearalgebra/mring.lisp

   1 ;; A Maxima ring structure
   2 ;; Copyright (C) 2005, 2007, 2021, 2022 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 (in-package :maxima)
  28
  29 ;; Let's have version numbers 1,2,3,...
  30 (eval-when (:compile-toplevel :load-toplevel :execute)
  31   ($put '$mring 1 '$version))
  32
  33 ;; (1) In maxima-grobner.lisp, there is a structure 'ring.'
  34
  35 ;; (2) Some functions in this structure, for example 'great' might
  36 ;;     not be defined for a ring; when this is the case, a function
  37 ;;     can signal an error.
  38
  39 ;; (3) Floating point addition isn't associative; so a mring needn't
  40 ;;     be a ring.  But a mring is 'close' to being a ring.
  41
  42 ;; Description of the mring fields:
  43
  44 ;; For floating point numbers (both real and complex), the lu_factor code uses
  45 ;; partial pivoting. Thus lu_factor needs to be able to compare the magnitudes.
  46 ;; Thus we include a function `great` in the mring structure. The function 'great`
  47 ;; is always supposed be be called on the magnitude of a ring element.
  48
  49 ;; For nonordered fields (such as generalring) where partial pivoting isn't needed,
  50 ;; we can define
  51 ;;   :great #'(lambda (a b) (declare (ignore a)) (eq t (meqp b 0)))
  52 ;; This inhibits row swapping unless it is absolutely needed to avoid dividing
  53 ;; by zero.
  54
  55 ;; The function coerce-to-lisp-float is only used by lu_factor to estimate the
  56 ;; matrix condition number. This is only important for the floating point case.
  57 ;; Setting :coerce-to-lisp-float to nil inhibits estimating the matrix condition
  58 ;; number. Except for matrices of floating point numbers, it's best to set
  59 ;;coerce-to-lisp-float to nil.
  60 (defstruct mring
  61   name
  62   coerce-to-lisp-float
  63   abs
  64   great
  65   add
  66   div
  67   rdiv
  68   reciprocal
  69   mult
  70   sub
  71   negate
  72   psqrt
  73   add-id
  74   mult-id
  75   fzerop
  76   adjoint
  77   maxima-to-mring
  78   mring-to-maxima)
  79
  80 (eval-when
  81     (:compile-toplevel :load-toplevel :execute)
  82   (defmvar $%mrings `((mlist) $floatfield $complexfield $rationalfield $crering
  83                       $generalring $bigfloatfield $runningerror $noncommutingring ))
  84   (defvar *gf-rings* '(gf-coeff-ring gf-ring ef-ring)) )
  85
  86 (defun $require_ring (ringname pos fun)
  87   (if (or ($member ringname $%mrings) (member ringname *gf-rings*))
  88     (get ringname 'ring)
  89     (merror (intl:gettext "The ~:M argument of the function '~:M' must be the name of a ring") pos fun)))
  90
  91 ;; flonum-epsilon is defined in clmacs
  92 (defparameter *floatfield*
  93   (make-mring
  94    :name '$floatfield
  95    :coerce-to-lisp-float #'cl:identity
  96    :abs #'abs
  97    :great #'>
  98    :add #'+
  99    :div #'/
 100    :rdiv #'/
 101    :reciprocal #'/
 102    :mult #'*
 103    :sub #'-
 104    :negate #'-
 105    :psqrt #'(lambda (s) (if (>= s 0) (cl:sqrt s) nil))
 106    :add-id #'(lambda () 0.0)
 107    :mult-id #'(lambda () 1.0)
 108    :fzerop #'(lambda (s)  (< (abs s) (* 4 +flonum-epsilon+)))
 109    :adjoint #'cl:identity
 110    :mring-to-maxima #'cl:identity
 111    :maxima-to-mring #'(lambda (s)
 112                         (setq s ($float s))
 113                         (if (floatp s) s (merror (intl:gettext "Unable to convert ~:M to a float") s)))))
 114
 115 (setf (get '$floatfield 'ring) *floatfield*)
 116
 117 (defun coerce-expr-to-clcomplex (ex)
 118   (let* (($numer t)
 119          ($float t)
 120          (split (risplit ex))
 121          ;; convert int, rat, and bfloat to float
 122          (re ($float (car split)))
 123          (im ($float (cdr split))))
 124     (declare (special $numer $float))
 125     (if (and (floatp re) (floatp im))
 126         (complex re im)
 127       (merror (intl:gettext "Unable to convert ~:M to a complex float") ex))))
 128
 129 (defparameter *complexfield*
 130   (make-mring
 131    :name '$complexfield
 132    :coerce-to-lisp-float #'cl:identity
 133    :abs #'cl:abs
 134    :great #'>
 135    :add #'+
 136    :div #'/
 137    :rdiv #'/
 138    :reciprocal #'/
 139    :mult #'*
 140    :sub #'-
 141    :negate #'-
 142    :psqrt #'(lambda (s) (if (and (= 0 (imagpart s)) (>= (realpart s) 0)) (cl:sqrt s) nil))
 143    :add-id #'(lambda () 0.0)
 144    :mult-id #'(lambda () 1.0)
 145    :fzerop #'(lambda (s) (< (abs s) (* 4 +flonum-epsilon+)))
 146    :adjoint #'cl:conjugate
 147    :mring-to-maxima #'(lambda (s) (add (realpart s) (mul '$%i (imagpart s)))) ;; was complexify
 148    :maxima-to-mring #'coerce-expr-to-clcomplex))
 149
 150 (setf (get '$complexfield 'ring) *complexfield*)
 151
 152 (defun rational-square-root (z)
 153   (let ((re (realpart z)) (im (imagpart z)) (q))
 154     (setq re (div ($num re) ($denom re)))
 155     (setq im (div ($num im) ($denom im)))
 156     (setq q (risplit (ftake '%sqrt (add re (mul '$%i im)))))
 157     (setq re (car q)
 158           im (cdr q))
 159     (cond ((and ($ratnump re) ($ratnump im))
 160            (complex (/ ($num re) ($denom re))
 161                     (/ ($num im) ($denom im))))
 162           (t nil))))
 163
 164 ;; We allow both real and complex rational numbers.
 165 (defparameter *rationalfield*
 166   (make-mring
 167    :name '$rationalfield
 168    :coerce-to-lisp-float nil
 169    :abs #'cl:abs
 170    :great #'(lambda (a b) (declare (ignore a)) (eql b 0))
 171    :add #'+
 172    :div #'/
 173    :rdiv #'/
 174    :reciprocal #'/
 175    :mult #'*
 176    :sub #'-
 177    :negate #'-
 178    :psqrt #'rational-square-root
 179    :add-id #'(lambda () 0)
 180    :mult-id #'(lambda () 1)
 181    :fzerop #'(lambda (s) (eql s 0))
 182    :adjoint #'cl:identity
 183    :mring-to-maxima #'(lambda (s)
 184       (let ((re (realpart s))
 185             (im (imagpart s)))
 186           (add (div (numerator re) (denominator re))
 187                (mul '$%i (div (numerator im) (denominator im))))))
 188    :maxima-to-mring
 189    #'(lambda (s)
 190         (setq s ($rationalize s))
 191         (let* ((z (trisplit s))
 192                (re (car z))
 193                (im (cdr z)))
 194        (if (and ($ratnump re) ($ratnump im))
 195          (complex (/ ($num re) ($denom re))
 196                   (/ ($num im) ($denom im)))
 197         (merror (intl:gettext "Unable to convert ~:M to a rational number") s))))))
 198 (setf (get '$rationalfield 'ring) *rationalfield*)
 199
 200 (defparameter *crering*
 201   (make-mring
 202    :name '$crering
 203    :coerce-to-lisp-float nil
 204    :abs #'$cabs
 205    :great #'(lambda (a b) (declare (ignore a)) (eq t (meqp b 0)))
 206    :add #'add
 207    :div #'div
 208    :rdiv #'div
 209    :reciprocal #'(lambda (s) (div 1 s))
 210    :mult #'mult
 211    :sub #'sub
 212    :negate #'(lambda (s) (mul -1 s))
 213    :psqrt #'(lambda (s) (ftake 'mexpt s (div 1 2)))
 214    :add-id #'(lambda () 0)
 215    :mult-id #'(lambda () 1)
 216    :fzerop #'(lambda (s) (eq t (meqp s 0)))
 217    :adjoint #'(lambda (s) (ftake '$conjugate s))
 218    :mring-to-maxima #'cl:identity
 219    :maxima-to-mring #'(lambda (s) ($rat s))))
 220
 221 (setf (get '$crering 'ring) *crering*)
 222
 223 (defparameter *generalring*
 224   (make-mring
 225    :name '$generalring
 226    :coerce-to-lisp-float nil
 227    :abs #'$cabs
 228    :great #'(lambda (a b) (declare (ignore a)) (eq t (meqp b 0)))
 229    :add #'(lambda (a b) (add a b))
 230    :div #'(lambda (a b) (div a b))
 231    :rdiv #'(lambda (a b) (div a b))
 232    :reciprocal #'(lambda (s) (div 1 s))
 233    :mult #'(lambda (a b) (mul a b))
 234    :sub #'(lambda (a b) (sub a b))
 235    :negate #'(lambda (a) (mul -1 a))
 236    :psqrt #'(lambda (s) (ftake 'mexpt s (div 1 2)))
 237    :add-id #'(lambda () 0)
 238    :mult-id #'(lambda () 1)
 239    :fzerop #'(lambda (s) (eq t (meqp s 0)))
 240    :adjoint #'(lambda (s) (ftake '$conjugate s))
 241    :mring-to-maxima #'cl:identity
 242    :maxima-to-mring #'cl:identity))
 243
 244 (setf (get '$generalring 'ring) *generalring*)
 245
 246 (defparameter *bigfloatfield*
 247   (make-mring
 248    :name '$bigfloatfield
 249    :coerce-to-lisp-float #'(lambda (s)
 250                              (setq s ($rectform ($float s)))
 251                              (complex ($realpart s) ($imagpart s)))
 252
 253    :abs #'$cabs
 254    :great #'mgrp
 255    :add #'(lambda (a b) ($rectform (add a b)))
 256    :div #'(lambda (a b) ($rectform (div a b)))
 257    :rdiv #'(lambda (a b) ($rectform (div a b)))
 258    :reciprocal #'(lambda (s) (div 1 s))
 259    :mult #'(lambda (a b) ($rectform (mul a b)))
 260    :sub #'(lambda (a b) ($rectform (sub a b)))
 261    :negate #'(lambda (a) (mul -1 a))
 262    :psqrt #'(lambda (s) (if (mlsp s 0) nil (ftake '%sqrt s)))
 263    :add-id #'(lambda () bigfloatzero)
 264    :mult-id #'(lambda () bigfloatone)
 265    :fzerop #'(lambda (s) (like s bigfloatzero))
 266    :adjoint #'cl:identity
 267    :mring-to-maxima #'cl:identity
 268    :maxima-to-mring #'(lambda (s)
 269                         (setq s ($rectform ($bfloat s)))
 270                         (if (complex-number-p s #'bigfloat-or-number-p) s
 271                                 (merror (intl:gettext "Unable to convert matrix entry to a big float"))))))
 272
 273 (setf (get '$bigfloatfield 'ring) *bigfloatfield*)
 274
 275 ;; --- *gf-rings* --- (used by src/numth.lisp) ------------------------------ ;;
 276 ;;                                                                            ;;
 277 (defparameter *gf-coeff-ring*
 278   (make-mring
 279    :name 'gf-coeff-ring
 280    :coerce-to-lisp-float nil
 281    :abs #'gf-cmod
 282    :great #'(lambda (a b) (declare (ignore a)) (= 0 b))
 283    :add #'gf-cplus-b
 284    :div #'(lambda (a b) (gf-ctimes a (gf-cinv b)))
 285    :rdiv #'(lambda (a b) (gf-ctimes a (gf-cinv b)))
 286    :reciprocal #'gf-cinv
 287    :mult #'gf-ctimes
 288    :sub #'(lambda (a b) (gf-cplus-b a (gf-cminus-b b)))
 289    :negate #'gf-cminus-b
 290    :psqrt #'(lambda (a) (let ((rs (zn-nrt a 2 *gf-char*))) (when rs (car rs))))
 291    :add-id #'(lambda () 0)
 292    :mult-id #'(lambda () 1)
 293    :fzerop #'(lambda (s) (= 0 s))
 294    :adjoint nil
 295    :mring-to-maxima #'cl:identity
 296    :maxima-to-mring #'cl:identity ))
 297 ;;
 298 (setf (get 'gf-coeff-ring 'ring) *gf-coeff-ring*)
 299 ;;
 300 (defparameter *gf-ring*
 301   (make-mring
 302    :name 'gf-ring
 303    :coerce-to-lisp-float nil
 304    :abs #'gf-mod
 305    :great #'(lambda (a b) (declare (ignore a)) (null b))
 306    :add #'gf-plus
 307    :div #'(lambda (a b) (gf-times a (gf-inv b *gf-red*) *gf-red*))
 308    :rdiv #'(lambda (a b) (gf-times a (gf-inv b *gf-red*) *gf-red*))
 309    :reciprocal #'(lambda (a) (gf-inv a *gf-red*))
 310    :mult #'(lambda (a b) (gf-times a b *gf-red*))
 311    :sub #'(lambda (a b) (gf-plus a (gf-minus b)))
 312    :negate #'gf-minus
 313    :psqrt #'(lambda (a)
 314               (let ((rs (gf-nrt-exit (gf-nrt a 2 *gf-red* *gf-ord*))))
 315                 (when rs (cadr rs)) ))
 316    :add-id #'(lambda () nil)
 317    :mult-id #'(lambda () '(0 1))
 318    :fzerop #'(lambda (s) (null s))
 319    :adjoint nil
 320    :mring-to-maxima #'gf-x2p
 321    :maxima-to-mring #'gf-p2x ))
 322 ;;
 323 (setf (get 'gf-ring 'ring) *gf-ring*)
 324 ;;
 325 (defparameter *ef-ring*
 326   (make-mring
 327    :name 'ef-ring
 328    :coerce-to-lisp-float nil
 329    :abs #'gf-mod
 330    :great #'(lambda (a b) (declare (ignore a)) (null b))
 331    :add #'gf-plus
 332    :div #'(lambda (a b) (gf-times a (gf-inv b *ef-red*) *ef-red*))
 333    :rdiv #'(lambda (a b) (gf-times a (gf-inv b *ef-red*) *ef-red*))
 334    :reciprocal #'(lambda (a) (gf-inv a *ef-red*))
 335    :mult #'(lambda (a b) (gf-times a b *ef-red*))
 336    :sub #'(lambda (a b) (gf-plus a (gf-minus b)))
 337    :negate #'gf-minus
 338    :psqrt #'(lambda (a)
 339               (let ((rs (gf-nrt-exit (gf-nrt a 2 *ef-red* *ef-ord*))))
 340                 (when rs (cadr rs)) ))
 341    :add-id #'(lambda () nil)
 342    :mult-id #'(lambda () '(0 1))
 343    :fzerop #'(lambda (s) (null s))
 344    :adjoint nil
 345    :mring-to-maxima #'gf-x2p
 346    :maxima-to-mring #'gf-p2x ))
 347
 348 (setf (get 'ef-ring 'ring) *ef-ring*)
 349 ;;                                                                            ;;
 350 ;; -------------------------------------------------------------------------- ;;
 351
 352 (defun fp-abs (a)
 353   (list (abs (first a)) (second a)))
 354
 355 (defun fp+ (a b)
 356   (cond ((= (first a) 0.0) b)
 357         ((= (first b) 0.0) a)
 358         (t
 359          (let ((s (+ (first a) (first b))))
 360            (if (= 0.0 s) (merror (intl:gettext "floating point divide by zero")))
 361            (list s (ceiling (+ 1
 362                                (abs (/ (* (first a) (second a)) s))
 363                                (abs (/ (* (first b) (second b)) s)))))))))
 364
 365 (defun fp- (a b)
 366   (cond ((= (first a) 0.0) (list (- (first b)) (second b)))
 367         ((= (first b) 0.0) a)
 368         (t
 369          (let ((s (- (first a) (first b))))
 370            (if (= 0.0 s) (merror (intl:gettext "floating point divide by zero")))
 371            (list s (ceiling (+ 1
 372                                (abs (/ (* (first a) (second a)) s))
 373                                (abs (/ (* (first b) (second b)) s)))))))))
 374
 375 (defun fp* (a b)
 376   (if (or (= (first a) 0.0) (= (first b) 0.0)) (list 0.0 0)
 377     (list (* (first a) (first b)) (+ 1 (second a) (second b)))))
 378
 379 (defun fp/ (a b)
 380   (if (= (first a) 0) (list 0.0 0)
 381     (list (/ (first a) (first b)) (+ 1 (second a) (second b)))))
 382
 383 (defun $addmatrices(fn &rest m)
 384   (mfuncall '$apply '$matrixmap `((mlist) ,fn ,@m)))
 385
 386 (defparameter *runningerror*
 387   (make-mring
 388    :name '$runningerror
 389    :coerce-to-lisp-float #'(lambda (s) (if (consp s) (first s) s))
 390    :abs #'fp-abs
 391    :great #'(lambda (a b) (> (first a) (first b)))
 392    :add #'fp+
 393    :div #'fp/
 394    :rdiv #'fp/
 395    :reciprocal #'(lambda (s) (fp/ (list 1 0) s))
 396    :mult #'fp*
 397    :sub #'fp-
 398    :negate #'(lambda (s) (list (- (first s)) (second s)))
 399    :psqrt #'(lambda (s) (if (> (first s) 0) (list (cl:sqrt (first s)) (+ 1 (second s))) nil))
 400    :add-id #'(lambda () (list 0 0))
 401    :mult-id #'(lambda () (list 1 0))
 402    :fzerop #'(lambda (s) (like (first s) 0))
 403    :adjoint #'cl:identity
 404    :mring-to-maxima #'(lambda (s) `((mlist) ,@s))
 405    :maxima-to-mring #'(lambda (s) (if ($listp s) (cdr s) (list ($float s) 1)))))
 406
 407 (setf (get '$runningerror 'ring) *runningerror*)
 408
 409 (defparameter *noncommutingring*
 410   (make-mring
 411    :name '$noncommutingring
 412    :coerce-to-lisp-float nil
 413    :abs #'$cabs
 414    :great #'(lambda (a b) (declare (ignore a)) (not (invertible-matrixp b '$noncommutingring)))
 415    :add #'(lambda (a b) (add a b))
 416    :div #'(lambda (a b) (progn
 417                           (let (($matrix_element_mult ".")
 418                                 ($matrix_element_transpose '$transpose))
 419                             (setq b (if ($matrixp b) ($invert_by_lu b '$noncommutingring)
 420                                       (take '(mncexpt) b -1)))
 421                             (take '(mnctimes) a b))))
 422
 423    :rdiv #'(lambda (a b) (progn
 424                            (let (($matrix_element_mult ".")
 425                                  ($matrix_element_transpose '$transpose))
 426                              (setq b (if ($matrixp b) ($invert_by_lu b '$noncommutingring)
 427                                        (take '(mncexpt) b -1)))
 428                              (take  '(mnctimes) b a))))
 429
 430
 431    :reciprocal #'(lambda (s) (progn
 432                                (let (($matrix_element_mult ".")
 433                                      ($matrix_element_transpose '$transpose))
 434                                  (if ($matrixp s) ($invert_by_lu s '$noncommutingring)
 435                                    (take '(mncexpt) s -1)))))
 436
 437    :mult #'(lambda (a b) (progn
 438                            (let (($matrix_element_mult ".")
 439                                  ($matrix_element_transpose '$transpose))
 440                              (take  '(mnctimes) a b))))
 441
 442
 443    :sub #'(lambda (a b) (sub a b))
 444    :negate #'(lambda (a) (mul -1 a))
 445    :add-id #'(lambda () 0)
 446    :psqrt #'(lambda (s) (take '(%sqrt) s))
 447    :mult-id #'(lambda () 1)
 448    :fzerop #'(lambda (s) (not (invertible-matrixp s '$noncommutingring)))
 449    :adjoint #'(lambda (s) ($transpose (take '($conjugate) s)))
 450    :mring-to-maxima #'cl:identity
 451    :maxima-to-mring #'cl:identity))
 452
 453 (setf (get '$noncommutingring 'ring) *noncommutingring*)
 454
 455 (defun ring-eval (e fld)
 456   (let ((fadd (mring-add fld))
 457         (fnegate (mring-negate fld))
 458         (fmult (mring-mult fld))
 459         (fdiv (mring-div fld))
 460         (fabs (mring-abs fld))
 461         (fconvert (mring-maxima-to-mring fld)))
 462
 463     (cond ((or ($numberp e) (symbolp e))
 464            (funcall fconvert (meval e)))
 465
 466           ;; I don't think an empty sum or product is possible here. If it is, append
 467           ;; the appropriate initial-value to reduce. Using the :inital-value isn't
 468           ;; a problem, but (fp* (a b) (1 0)) --> (a (+ b 1)).  A better value is
 469           ;; (fp* (a b) (1 0)) --> (a b).
 470
 471           ((op-equalp e 'mplus)
 472            (reduce fadd (mapcar #'(lambda (s) (ring-eval s fld)) (margs e)) :from-end t))
 473
 474           ((op-equalp e 'mminus)
 475            (funcall fnegate (ring-eval (first (margs e)) fld)))
 476
 477           ((op-equalp e 'mtimes)
 478            (reduce fmult (mapcar #'(lambda (s) (ring-eval s fld)) (margs e)) :from-end t))
 479
 480           ((op-equalp e 'mquotient)
 481            (funcall fdiv (ring-eval (first (margs e)) fld)(ring-eval (second (margs e)) fld)))
 482
 483           ((op-equalp e 'mabs) (funcall fabs (ring-eval (first (margs e)) fld)))
 484
 485           ((and (or (eq (mring-name fld) '$floatfield) (eq (mring-name fld) '$complexfield))
 486                 (consp e) (consp (car e)) (gethash (mop e) *flonum-op*))
 487            (apply (gethash (mop e) *flonum-op*) (mapcar #'(lambda (s) (ring-eval s fld)) (margs e))))
 488
 489           (t (merror (intl:gettext "Unable to evaluate ~:M in the ring '~:M'") e (mring-name fld))))))
 490
 491 (defmspec $ringeval (e)
 492   (let ((fld (get (or (car (member (meval (nth 2 e)) $%mrings)) '$generalring) 'ring)))
 493     (funcall (mring-mring-to-maxima fld) (ring-eval (nth 1 e) fld))))