1
;; Maxima functions for floor, ceiling, and friends
2
;; Copyright (C) 2004, 2005, Barton Willis
5
;; Department of Mathematics,
6
;; University of Nebraska at Kearney
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.
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.
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.
28
(macsyma-module nummod)
30
;; Let's have version numbers 1,2,3,...
32
(eval-when (:compile-toplevel :load-toplevel :execute)
33
(mfuncall '$declare '$integervalued '$feature)
34
($put '$nummod 2 '$version))
36
;; charfun(pred) evaluates to 1 when the predicate 'pred' evaluates
37
;; to true; it evaluates to 0 when 'pred' evaluates to false; otherwise,
38
;; it evaluates to a noun form.
40
(defprop $charfun simp-charfun operators)
42
(defun simp-charfun (e bool z)
44
(setq e (specrepcheck e))
45
(cond ((member 'simp (car e)) e)
47
(let (($prederror nil))
48
(setq e (simplifya ($ratdisrep (nth 1 e)) z))
49
(setq bool (mevalp (mfuncall '$ev e '$nouns)))
52
(t `(($charfun simp) ,e)))))))
54
(defun integer-part-of-sum (e)
55
(let ((i-sum 0) (n-sum 0) (o-sum 0) (n))
58
(cond ((maxima-integerp ei)
59
(setq i-sum (add ei i-sum)))
60
((or (ratnump ei) (floatp ei) ($bfloatp ei))
61
(setq n-sum (add ei n-sum)))
63
(setq o-sum (add ei o-sum)))))
64
(setq n (simplify `(($floor) ,n-sum)))
65
(setq i-sum (add i-sum n))
66
(setq o-sum (add o-sum (sub n-sum n)))
67
(values i-sum o-sum)))
69
(defprop $floor simp-floor operators)
71
(defprop $floor tex-matchfix tex)
72
(defprop $floor (("\\left \\lfloor " ) " \\right \\rfloor") texsym)
74
; These defprops for $entier are copied from orthopoly/orthopoly-init.lisp.
76
(defprop $entier tex-matchfix tex)
77
(defprop $entier (("\\lfloor ") " \\rfloor") texsym)
79
;; When constantp(x) is true, we use bfloat evaluation to try to determine
80
;; the floor. If numerical evaluation of e is ill-conditioned, this function
81
;; can misbehave. I'm somewhat uncomfortable with this, but it is no worse
82
;; than some other stuff. One safeguard -- the evaluation is done with three
83
;; values for fpprec. When the floating point evaluations are
84
;; inconsistent, bailout and return nil. I hope that this function is
85
;; much more likely to return nil than it is to return a bogus value.
87
(defun pretty-good-floor-or-ceiling (x fn &optional (digits 'no-value))
88
(let (($fpprec) ($float2bf t))
89
(cond ((eq digits 'no-value)
91
;; When x doesn't evaluate to a bigfloat, bailout and return nil.
92
;; This happens when, for example, x = asin(2). For now, bfloatp
93
;; evaluates to nil for a complex big float. If this ever changes,
94
;; the freeof(%i, xf) will detect the complex big float and return nil.
96
(meval `((msetq) $fpprec 25))
97
(let (($algebraic t) (xf ($bfloat (setq x ($rectform x)))))
98
(if (and ($bfloatp xf) ($freeof '$%i xf))
99
(pretty-good-floor-or-ceiling x fn
100
(max 25 (ceiling (* 0.4 (nth 2 xf))))) nil)))
102
(let ((f1) (f2) (f3) (xf) (eps) (lb) (ub) (n))
103
(meval `((msetq) $fpprec ,digits))
104
(setq f1 (mfuncall fn ($bfloat x)))
107
(meval `((msetq) $fpprec ,digits))
108
(setq f2 (mfuncall fn ($bfloat x)))
111
(meval `((msetq) $fpprec ,digits))
112
(setq f3 (mfuncall fn (setq xf ($bfloat x))))
114
;; Let's say that the true value of x is in the interval
115
;; [xf - |xf| * eps, xf + |xf| * eps], where eps = 10^(20 - digits).
116
;; Define n to be the number of integers in this interval; we have
118
(setq eps (power ($bfloat 10) (- 20 digits)))
119
(setq lb (sub xf (mult (simplify (list '(mabs) xf)) eps)))
120
(setq ub (add xf (mult (simplify (list '(mabs) xf)) eps)))
121
(setq n (sub (mfuncall '$ceiling ub) (mfuncall '$ceiling lb)))
123
;; Provided f1, f2, and f3 are the same and that n = 0, return f1.
125
(if (and (= f1 f2) (= f2 f3) (= n 0)) f1 nil))))))
127
;; (a) The function fpentier rounds a bigfloat towards zero--we need to
130
;; (b) Mostly for fun, floor(minf) --> und and etc. I suppose floor
131
;; should be undefined for many other arguments---for example
132
;; floor(a < b), floor(true).
134
;; I think floor(complex number) should be undefined too. Some folks
135
;; might like floor(a + %i b) --> floor(a) + %i floor(b). But
136
;; this would violate the integer-valuedness of floor.
138
(defun simp-floor (e e1 z)
140
(setq e (specrepcheck e))
141
(cond ((not (member 'simp (car e)))
142
(setq e (simplifya ($ratdisrep (nth 1 e)) z))
143
(if (ratnump e) (setq e (/ (cadr e) (caddr e))))
144
(cond ((numberp e) (floor e))
147
(setq e1 (fpentier e))
148
(if (and (minusp (cadr e)) (not (zerop1 (sub e1 e))))
151
(($orderlessp e (neg e))
152
(sub* 0 (simplifya `(($ceiling) ,(neg e)) nil)))
154
((maxima-integerp e) e)
156
((and (setq e1 (mget e '$numer)) (floor e1)))
158
((or (member e infinities) (eq e '$und) (eq e '$ind)) '$und)
159
((or (like e '$zerob)) -1)
160
((or (like e '$zeroa)) 0)
162
((and ($constantp e) (pretty-good-floor-or-ceiling e '$floor)))
165
(let ((i-sum) (o-sum))
166
(multiple-value-setq (i-sum o-sum) (integer-part-of-sum e))
167
(setq o-sum (if (like i-sum 0) `(($floor simp) ,o-sum)
168
(simplifya `(($floor) ,o-sum) nil)))
171
;; handle 0 < e < 1 implies floor(e) = 0 and
172
;; -1 < e < 0 implies floor(e) = -1.
174
((and (eq ($compare 0 e) '&<) (eq ($compare e 1) '&<)) 0)
175
((and (eq ($compare -1 e) '&<) (eq ($compare e 0) '&<)) -1)
176
(t `(($floor simp) ,e))))
179
(defprop $ceiling simp-ceiling operators)
181
(defprop $ceiling tex-matchfix tex)
182
(defprop $ceiling (("\\left \\lceil " ) " \\right \\rceil") texsym)
184
(defun simp-ceiling (e e1 z)
186
(setq e ($ratdisrep e))
187
(cond ((not (member 'simp (car e)))
188
(setq e (simplifya ($ratdisrep (nth 1 e)) z))
189
(if (ratnump e) (setq e (/ (cadr e) (caddr e))))
190
(cond ((numberp e) (ceiling e))
192
(sub 0 (simplify `(($floor) ,(sub 0 e)))))
194
(($orderlessp e (neg e))
195
(sub* 0 (simplifya `(($floor) ,(neg e)) nil)))
197
((maxima-integerp e) e)
198
((and (setq e1 (mget e '$numer)) (ceiling e1)))
200
((or (member e infinities) (eq e '$und) (eq e '$ind)) '$und)
201
((or (like e '$zerob)) 0)
202
((or (like e '$zeroa)) 1)
204
((and ($constantp e) (pretty-good-floor-or-ceiling e '$ceiling)))
207
(let ((i-sum) (o-sum))
208
(multiple-value-setq (i-sum o-sum) (integer-part-of-sum e))
209
(setq o-sum (if (like i-sum 0) `(($ceiling simp) ,o-sum)
210
(simplifya `(($ceiling) ,o-sum) nil)))
214
;; handle 0 < e < 1 implies ceiling(e) = 1 and
215
;; -1 < e < 0 implies ceiling(e) = 0.
217
((and (eq ($compare 0 e) '&<) (eq ($compare e 1) '&<)) 1)
218
((and (eq ($compare -1 e) '&<) (eq ($compare e 0) '&<)) 0)
219
(t `(($ceiling simp) ,e))))
222
(defprop $mod simp-nummod operators)
223
(defprop $mod tex-infix tex)
224
(defprop $mod (" \\rm{mod} ") texsym)
225
(defprop $mod 180. tex-rbp)
226
(defprop $mod 180. tex-rbp)
228
;; See "Concrete Mathematics," Section 3.21.
230
(defun simp-nummod (e e1 z)
232
(cond ((not (member 'simp (car e)))
233
(setq e (mapcar #'specrepcheck (margs e)))
234
(setq e (mapcar #'(lambda (s) (simplifya s z)) e))
235
(let ((x (nth 0 e)) (y (nth 1 e)))
236
(cond ((or (like 0 y) (like 0 x)) x)
237
((like 1 y) (sub* x `(($floor) ,x)))
238
((and ($constantp x) ($constantp y) (not (like 0 y)))
239
(sub* x (mul* y `(($floor) ,(div x y)))))
240
((not (like 1 (setq e1 ($gcd x y))))
241
(mul* e1 `(($mod) ,(div* x e1) ,(div* y e1))))
243
(t `(($mod simp) ,x ,y)))))
added
removed
src/nummod.lisp
