src/nummod.lisp

   1 ;; Maxima functions for floor, ceiling, and friends
   2 ;; Copyright (C) 2004, 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 (in-package :maxima)
  28
  29 (macsyma-module nummod)
  30
  31 ;; Let's have version numbers 1,2,3,...
  32
  33 (eval-when (:load-toplevel :execute)
  34   ($put '$nummod 3 '$version)
  35   ;; Let's remove built-in symbols from list for user-defined properties.
  36   (setq $props (remove '$nummod $props)))
  37
  38 (defmacro opcons (op &rest args)
  39   `(simplify (list (list ,op) ,@args)))
  40
  41 ;; charfun(pred) evaluates to 1 when the predicate 'pred' evaluates
  42 ;; to true; it evaluates to 0 when 'pred' evaluates to false; otherwise,
  43 ;; it evaluates to a noun form.
  44
  45 (defprop $charfun simp-charfun operators)
  46 (defprop %charfun simp-charfun operators)
  47
  48 (defun simp-charfun (e yy z)
  49   (declare (ignore yy))
  50   (oneargcheck e)
  51   (setq e (take '($is) (simplifya (specrepcheck (second e)) z)))
  52   (let* (($prederror nil)
  53          (bool (mevalp e)))
  54     (cond ((eq t bool) 1)
  55           ((eq nil bool) 0)
  56           ((op-equalp e '$is) `(($charfun simp) ,(second e)))
  57           (t `(($charfun simp) ,e)))))
  58
  59 (defun integer-part-of-sum (e)
  60   (let ((i-sum 0) (n-sum 0) (o-sum 0) (n))
  61     (setq e (margs e))
  62     (dolist (ei e)
  63       (cond ((maxima-integerp ei)
  64              (setq i-sum (add ei i-sum)))
  65             ((or (ratnump ei) (floatp ei) ($bfloatp ei))
  66              (setq n-sum (add ei n-sum)))
  67             (t
  68              (setq o-sum (add ei o-sum)))))
  69     (setq n (opcons '$floor n-sum))
  70     (setq i-sum (add i-sum n))
  71     (setq o-sum (add o-sum (sub n-sum n)))
  72     (values i-sum o-sum)))
  73
  74 (defprop $floor simp-floor operators)
  75
  76 (defprop $floor tex-matchfix tex)
  77 (defprop $floor (("\\left \\lfloor " ) " \\right \\rfloor") texsym)
  78
  79 ;; $floor distributes over lists, matrices, and equations
  80 (setf (get '$floor 'distribute_over) '(mlist $matrix mequal))
  81
  82 ; These defprops for $entier are copied from orthopoly/orthopoly-init.lisp.
  83
  84 (defprop $entier tex-matchfix tex)
  85 (defprop $entier (("\\lfloor ") " \\rfloor") texsym)
  86
  87 ;; $entier and $fix distributes over lists, matrices, and equations
  88 (setf (get '$entier 'distribute_over) '(mlist $matrix mequal))
  89 (setf (get '$fix 'distribute_over) '(mlist $matrix mequal))
  90
  91 ;; For an example, see pretty-good-floor-or-ceiling. Code courtesy of Stavros Macrakis.
  92
  93 (defmacro bind-fpprec (val &body exprs)
  94   `(let (fpprec bigfloatzero bigfloatone bfhalf bfmhalf)
  95      (let (($fpprec (fpprec1 nil ,val)))
  96        ,@exprs)))
  97
  98 ;; Return true if the expression can be formed using rational numbers, logs, mplus, mexpt, or mtimes.
  99
 100 (defun use-radcan-p (e)
 101   (or ($ratnump e) (and (op-equalp e '%log 'mexpt 'mplus 'mtimes) (every 'use-radcan-p (cdr e)))))
 102
 103 ;; When constantp(x) is true, we use bfloat evaluation to try to determine
 104 ;; the ceiling or floor. If numerical evaluation of e is ill-conditioned, this function
 105 ;; can misbehave.  I'm somewhat uncomfortable with this, but it is no worse
 106 ;; than some other stuff. One safeguard -- the evaluation is done with three
 107 ;; values for fpprec.  When the floating point evaluations are
 108 ;; inconsistent, bailout and return nil.  I hope that this function is
 109 ;; much more likely to return nil than it is to return a bogus value.
 110
 111 (defun pretty-good-floor-or-ceiling (x fn &optional digits)
 112   (let (($float2bf t) ($algebraic t) (f1) (f2) (f3) (eps) (lb) (ub) (n))
 113
 114     (setq digits (if (and (integerp digits) (> 0 digits)) digits 25))
 115     (catch 'done
 116
 117       ;; To handle ceiling(%i^%i), we need to apply rectform. If bfloat
 118       ;; is improved, it might be possible to remove this call to rectform.
 119
 120       (setq x ($rectform x))
 121       (if (not ($freeof '$%i '$minf '$inf '$und '$infinity x)) (throw 'done nil))
 122
 123       ;; When x doesn't evaluate to a bigfloat, bailout and return nil.
 124       ;; This happens when, for example, x = asin(2). For now, bfloatp
 125       ;; evaluates to nil for a complex big float. If this ever changes,
 126       ;; this code might need to be repaired.
 127
 128       (bind-fpprec digits
 129         (setq f1 ($bfloat x))
 130         (if (not ($bfloatp f1)) (throw 'done nil)))
 131
 132       (incf digits 20)
 133       (setq f2 (bind-fpprec digits ($bfloat x)))
 134       (if (not ($bfloatp f2)) (throw 'done nil))
 135
 136       (incf digits 20)
 137       (bind-fpprec digits
 138         (setq f3 ($bfloat x))
 139         (if (not ($bfloatp f3)) (throw 'done nil))
 140
 141         ;; Let's say that the true value of x is in the interval
 142         ;; [f3 - |f3| * eps, f3 + |f3| * eps], where eps = 10^(20 - digits).
 143         ;; Define n to be the number of integers in this interval; we have
 144
 145         (setq eps (power ($bfloat 10) (- 20 digits)))
 146         (setq lb (sub f3 (mult (take '(mabs) f3) eps)))
 147         (setq ub (add f3 (mult (take '(mabs) f3) eps)))
 148
 149         ;; If n > 1, we're going to give up. This is true iff ceiling(ub) -
 150         ;; ceiling(lb) > 1. However, we don't want to blindly try to calculate
 151         ;; that: if ub and lb are enormous, we probably won't have enough
 152         ;; precision in the bigfloats to calculate either ceiling. Start by
 153         ;; taking the difference: if it's at least 2, then we know that n >= 2.
 154         (when (fpgreaterp (cdr (sub ub lb)) (cdr ($bfloat 2)))
 155           (throw 'done nil))
 156
 157         (setq n (sub (take '($ceiling) ub)
 158                      (take '($ceiling) lb))))
 159
 160       (setq f1 (take (list fn) f1))
 161       (setq f2 (take (list fn) f2))
 162       (setq f3 (take (list fn) f3))
 163
 164       ;; Provided f1 = f2 = f3 and n = 0, return f1.
 165       ;; If n = 1 and (use-radcan-p e) and ($radcan e) is a $ratnump, return
 166       ;; floor / ceiling of radcan(x).
 167
 168       (cond ((and (= f1 f2 f3) (= n 0)) f1)
 169             ((and (=  f1 f2 f3) (= n 1) (use-radcan-p x))
 170              (setq x ($radcan x))
 171              (if ($ratnump x) (take (list fn) x) nil))
 172             (t nil)))))
 173
 174
 175 ;; (a) The function fpentier rounds a bigfloat towards zero--we need to
 176 ;;     check for that.
 177
 178 ;; (b) Since limit(f(x))=(m)inf implies limit(floor(f(x)))=(m)inf,
 179 ;;     floor(inf/minf) = inf/minf.  Since minf<limit(f(x))<inf implies
 180 ;;     minf<limit(floor(f(x)))<inf, floor(ind)=ind
 181
 182 ;;     I think floor(complex number) should be undefined too.  Some folks
 183 ;;     might like floor(a + %i b) --> floor(a) + %i floor(b). But
 184 ;;     this would violate the integer-valuedness of floor.
 185 ;;     So floor(infinity) remains unsimplified
 186
 187 (defun simp-floor (e e1 z)
 188   (oneargcheck e)
 189   (setq e (simplifya (specrepcheck (nth 1 e)) z))
 190
 191   (cond ((numberp e) (floor e))
 192
 193         ((ratnump e) (floor (cadr e) (caddr e)))
 194
 195         (($bfloatp e)
 196          (setq e1 (fpentier e))
 197          (if (and (minusp (cadr e)) (not (zerop1 (sub e1 e))))
 198              (1- e1)
 199              e1))
 200
 201         ((maxima-integerp e) e)
 202
 203         (($orderlessp e (neg e))
 204          (sub 0 (opcons '$ceiling (neg e))))
 205
 206         ((and (setq e1 (mget e '$numer)) (floor e1)))
 207
 208         ((or (member e infinities) (eq e '$und) (eq e '$ind)) '$und)
 209         ;;((member e '($inf $minf $ind $und)) e) ; Proposed code
 210         ;; leave floor(infinity) as is
 211         ((or (eq e '$zerob)) -1)
 212         ((or (eq e '$zeroa)) 0)
 213
 214         ((and ($constantp e) (pretty-good-floor-or-ceiling e '$floor)))
 215
 216         ((mplusp e)
 217          (let ((i-sum) (o-sum))
 218            (multiple-value-setq (i-sum o-sum) (integer-part-of-sum e))
 219            (setq o-sum (if (equal i-sum 0) `(($floor simp) ,o-sum)
 220                          (opcons '$floor o-sum)))
 221            (add i-sum o-sum)))
 222
 223         ;; handle 0 <= e < 1 implies floor(e) = 0 and
 224         ;; -1 <= e < 0 implies floor(e) = -1.
 225
 226         ((and (member ($compare 0 e) '("<" "<=") :test #'equal)
 227               (equal ($compare e 1) "<"))
 228          0)
 229         ((and (member ($compare -1 e) '("<" "<=") :test #'equal)
 230               (equal ($compare e 0) "<"))
 231          -1)
 232         (t `(($floor simp) ,e))))
 233
 234 (defun floor-integral (x)
 235   (let ((f (take '($floor) x)))
 236     (mul f (div 1 2) (add (mul 2 x) -1 (neg f)))))
 237
 238 (defun ceiling-integral (x)
 239   (let ((f (take '($ceiling) x)))
 240     (mul f (div 1 2) (add 1 (mul 2 x) (neg f)))))
 241
 242 (putprop '$floor `((x) ,#'floor-integral) 'integral)
 243 (putprop '$ceiling `((x) ,#'ceiling-integral) 'integral)
 244
 245 (defprop $ceiling simp-ceiling operators)
 246
 247 (defprop $floor simplim%floor simplim%function)
 248
 249 (defun simplim%floor (expr var val)
 250   (let* ((arg (cadr expr))
 251          (b (behavior arg var val))
 252          (arglimab (limit arg var val 'think)) ; with $zeroa $zerob
 253          (arglim (ridofab arglimab)))
 254     (cond ((and (= b -1)
 255                 (maxima-integerp arglim))
 256            (m- arglim 1))
 257           ((and (= b 1)
 258                 (maxima-integerp arglim))
 259            arglim)
 260           ((and ($constantp arglim)
 261                 (not (maxima-integerp arglim)))
 262            (simplify (list '($floor) arglim)))
 263           (t
 264            (throw 'limit nil)))))
 265
 266 (defprop $ceiling tex-matchfix tex)
 267 (defprop $ceiling (("\\left \\lceil " ) " \\right \\rceil") texsym)
 268
 269 ;; $ceiling distributes over lists, matrices, and equations.
 270 (setf (get '$ceiling 'distribute_over) '(mlist $matrix mequal))
 271
 272 (defun simp-ceiling (e e1 z)
 273   (oneargcheck e)
 274   (setq e (simplifya (specrepcheck (nth 1 e)) z))
 275   (cond ((numberp e) (ceiling e))
 276
 277         ((ratnump e) (ceiling (cadr e) (caddr e)))
 278
 279         (($bfloatp e)
 280          (sub 0 (opcons '$floor (sub 0 e))))
 281
 282         ((maxima-integerp e) e)
 283
 284         (($orderlessp e (neg e))
 285          (sub* 0 (opcons '$floor (neg e))))
 286
 287         ((and (setq e1 (mget e '$numer)) (ceiling e1)))
 288
 289         ((or (member e infinities) (eq e '$und) (eq e '$ind)) '$und)
 290         ((or (eq e '$zerob)) 0)
 291         ((or (eq e '$zeroa)) 1)
 292
 293         ((and ($constantp e) (pretty-good-floor-or-ceiling e '$ceiling)))
 294
 295         ((mplusp e)
 296          (let ((i-sum) (o-sum))
 297            (multiple-value-setq (i-sum o-sum) (integer-part-of-sum e))
 298            (setq o-sum (if (equal i-sum 0) `(($ceiling simp) ,o-sum)
 299                          (opcons '$ceiling o-sum)))
 300            (add i-sum o-sum)))
 301
 302         ;; handle 0 < e <= 1 implies ceiling(e) = 1 and
 303         ;; -1 < e <= 0 implies ceiling(e) = 0.
 304
 305         ((and (equal ($compare 0 e) "<")
 306               (member ($compare e 1) '("<" "<=") :test #'equal))
 307          1)
 308         ((and (equal ($compare -1 e) "<")
 309               (member ($compare e 0) '("<" "<=") :test #'equal))
 310          0)
 311         (t `(($ceiling simp) ,e))))
 312
 313 (defprop $ceiling simplim%ceiling simplim%function)
 314
 315 (defun simplim%ceiling (expr var val)
 316   (let* ((arg (cadr expr))
 317          (b (behavior arg var val))
 318          (arglimab (limit arg var val 'think)) ; with $zeroa $zerob
 319          (arglim (ridofab arglimab)))
 320     (cond ((and (= b -1)
 321                 (maxima-integerp arglim))
 322            arglim)
 323           ((and (= b 1)
 324                 (maxima-integerp arglim))
 325            (m+ arglim 1))
 326           ((and ($constantp arglim)
 327                 (not (maxima-integerp arglim)))
 328            (simplify (list '($ceiling) arglim)))
 329           (t
 330            (throw 'limit nil)))))
 331
 332
 333 (defprop $mod simp-nummod operators)
 334 (defprop $mod tex-infix tex)
 335 (defprop $mod (" \\rm{mod} ") texsym)
 336 (defprop $mod 180. tex-lbp)
 337 (defprop $mod 180. tex-rbp)
 338
 339 ;; $mod distributes over lists, matrices, and equations
 340 (setf (get '$mod 'distribute_over) '(mlist $matrix mequal))
 341
 342 ;; See "Concrete Mathematics," Section 3.21.
 343
 344 (defun simp-nummod (e e1 z)
 345   (twoargcheck e)
 346   (let ((x (simplifya (specrepcheck (cadr e)) z))
 347         (y (simplifya (specrepcheck (caddr e)) z)))
 348     (cond ((or (equal 0 y) (equal 0 x)) x)
 349           ((equal 1 y) (sub x (opcons '$floor x)))
 350           ((and ($constantp x) ($constantp y))
 351            (sub x (mul y (opcons '$floor (div x y)))))
 352           ((not (equal 1 (setq e1 ($gcd x y))))
 353            (mul e1 (opcons '$mod (div x e1) (div y e1))))
 354           (t `(($mod simp) ,x ,y)))))
 355
 356 ;; The function 'round' rounds a number to the nearest integer. For a tie, round to the
 357 ;; nearest even integer.
 358
 359 (defprop %round simp-round operators)
 360 (setf (get '%round 'integer-valued) t)
 361 (setf (get '%round 'reflection-rule) #'odd-function-reflect)
 362 (setf (get '$round 'alias) '%round)
 363 (setf (get '$round 'verb) '%round)
 364 (setf (get '%round 'noun) '$round)
 365
 366 ;; round distributes over lists, matrices, and equations.
 367 (setf (get '%round 'distribute_over) '(mlist $matrix mequal))
 368
 369 (defun simp-round (e yy z)
 370   (oneargcheck e)
 371   (setq yy (caar e)) ;; find a use for the otherwise unused YY.
 372   (setq e (simplifya (specrepcheck (second e)) z))
 373   (cond (($featurep e '$integer) e) ;; takes care of round(round(x)) --> round(x).
 374         ((member e '($inf $minf $und $ind) :test #'eq) e)
 375         ((eq e '$zerob) 0)
 376         ((eq e '$zeroa) 0)
 377         (t
 378          (let* ((lb (take '($floor) e))
 379                 (ub (take '($ceiling) e))
 380                 (sgn (csign (sub (sub ub e) (sub e lb)))))
 381            (cond ((eq sgn t) `((,yy simp) ,e))
 382                  ((eq sgn '$neg) ub)
 383                  ((eq sgn '$pos) lb)
 384                  ((alike lb ub) lb) ;; For floats that are integers, this can happen. Maybe featurep should catch this.
 385                  ((and (eq sgn '$zero) ($featurep lb '$even)) lb)
 386                  ((and (eq sgn '$zero) ($featurep ub '$even)) ub)
 387                  ((apply-reflection-simp yy e t))
 388                  (t `((,yy simp) ,e)))))))
 389
 390 (defprop %round simplim%round simplim%function)
 391
 392 (defun simplim%round (expr var val)
 393   (let* ((arg (cadr expr))
 394          (b (behavior arg var val))
 395          (arglimab (limit arg var val 'think)) ; with $zeroa $zerob
 396          (arglim (ridofab arglimab)))
 397     (cond ((and (= b -1)
 398                 (maxima-integerp (m+ 1//2 arglim)))
 399            (m- arglim 1//2))
 400           ((and (= b 1)
 401                 (maxima-integerp (m+ 1//2 arglim)))
 402            (m+ arglim 1//2))
 403           ((and (mnump arglim)
 404                 (not (maxima-integerp (m+ 1//2 arglim))))
 405            (simplify (list '(%round) arglim)))
 406           (t
 407            (throw 'limit nil)))))
 408
 409 ;; Round a number towards zero.
 410
 411 (defprop %truncate simp-truncate operators)
 412 (setf (get '%truncate 'integer-valued) t)
 413 (setf (get '%truncate 'reflection-rule) #'odd-function-reflect)
 414 (setf (get '$truncate 'alias) '%truncate)
 415 (setf (get '$truncate 'verb) '%truncate)
 416 (setf (get '%truncate 'noun) '$truncate)
 417
 418 ;; truncate distributes over lists, matrices, and equations.
 419 (setf (get '%truncate 'distribute_over) '(mlist $matrix mequal))
 420
 421 (defun simp-truncate (e yy z)
 422   (oneargcheck e)
 423   (setq yy (caar e)) ;; find a use for the otherwise unused YY.
 424   (setq e (simplifya (specrepcheck (second e)) z))
 425   (cond (($featurep e '$integer) e) ;; takes care of truncate(truncate(x)) --> truncate(x).
 426         ((member e '($inf $minf $und $ind) :test #'eq) e)
 427         ((eq e '$zerob) 0)
 428         ((eq e '$zeroa) 0)
 429         (t
 430          (let ((sgn (csign e)))
 431            (cond ((member sgn '($neg $nz) :test #'eq) (take '($ceiling) e))
 432                  ((member sgn '($zero $pz $pos) :test #'eq) (take '($floor) e))
 433                  ((apply-reflection-simp yy e t))
 434                  (t `((,yy simp) ,e)))))))
 435
 436 ;;; integration for signum, unit_step, and mod.
 437
 438 ;; integrate(signum(x),x) = |x|
 439 (defun signum-integral (x)
 440                  (take '(mabs) x))
 441
 442 ;; integrate(unit_step(x),x) = (x + |x|)/2
 443 (defun unit-step-integral (x)
 444                 (div (add x (take '(mabs) x)) 2))
 445
 446 ;; We have mod(x,a) = x-a*floor(x/a). This gives
 447 ;;   integrate(mod(x,a),x) = (a^2 floor(x/a)^2 + (a^2 - 2 a x) floor(x/a) + x^2)/2
 448 ;; In terms of mod(x,a), an antiderivative is
 449 ;;   integrate(mod(x,a),x) = (mod(x,a)^2-a*mod(x,a)+a*x)/2
 450 ;; Before this function is called, Maxima checks if a explicitly depends on x. So
 451 ;; this function doesn't need to do that check.
 452 (defun mod-integral (x a)
 453                 (let ((q (take '($mod) x a)))
 454                    (div (add (mul q q) (mul -1 a q) (mul a x)) 2)))
 455
 456 (putprop '%signum (list (list 'x) #'signum-integral) 'integral)
 457 (putprop '$unit_step (list (list 'x) #'unit-step-integral) 'integral)
 458
 459 ;; integrate(mod(x,a),a) doesn't have representation in terms of functions
 460 ;; known to Maxima, I think. (Barton Willis, 2020).
 461 (putprop '$mod (list (list 'x 'y) #'mod-integral nil) 'integral)