← Back to branch summary

~ubuntu-branches/ubuntu/karmic/maxima/karmic

« back to all changes in this revision

Viewing changes to src/csimp2.lisp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Sauthier
  • Date: 2009-07-13 15:38:41 UTC
  • mfrom: (3.1.3 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090713153841-gtux06oun30kuuo7
Tags: 5.17.1-1ubuntu1
https://launchpad.net/bugs/296643
https://launchpad.net/bugs/242243
* Merge from debian unstable, remaining changes (LP: #296643, LP: #242243):
   - debian/maxima-doc.doc-base.{tips, plotting}:
    + Use .shtml instead of .html to fix lintian errors.
   - debian/maxima-emacs.emacsen-install:
    + Install symlinks for source files rather than copying them.  This
      makes find-function work.
    + Install symlink for *.lisp so that we don't need to add
      /usr/share/emacs/site-lisp/maxima to load-path.
  - debian/maxima-emacs.emacsen-startup:
    + Remove use of /usr/share/emacs/site-lisp/maxima, since this
      causes load-path shadows and is not needed anymore.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/csimp2.lisp
14
14
 
15
15
(load-macsyma-macros rzmac)
16
16
 
17
 
(declare-top (special var %p%i varlist plogabs half%pi nn* dn*))
 
17
(declare-top (special var %p%i varlist plogabs half%pi nn* dn* $factlim))
 
18
 
 
19
(defmvar $gammalim 10000
 
20
  "Controls simplification of gamma for rational number arguments.")
 
21
 
 
22
(defvar $gamma_expand nil
 
23
  "Expand gamma(z+n) for n an integer when T.") 
18
24
 
19
25
(defmfun simpplog (x vestigial z)
20
26
  (declare (ignore vestigial))
176
182
(defmfun simpgamma (x vestigial z)
177
183
  (declare (ignore vestigial))
178
184
  (oneargcheck x)
179
 
  (let* ((j (simpcheck (cadr x) z))
180
 
         (jr ($realpart j))
181
 
         (ji ($imagpart j)))
182
 
    (cond ((floatp j) (gammafloat j))
183
 
          ((and (numberp jr) 
184
 
                (numberp ji)
185
 
                (or $numer (floatp jr) (floatp ji)))
186
 
           (complexify (gamma-lanczos (complex (float jr)
187
 
                                               (float ji)))))
188
 
          ((or (not (mnump j))
189
 
               (ratgreaterp (simpabs (list '(%abs) j) 1 t) $gammalim))
190
 
           (eqtest (list '(%gamma) j) x))
 
185
  (let ((j (simpcheck (cadr x) z)))
 
186
    (cond ((and (floatp j)
 
187
                (or (zerop j)
 
188
                    (and (< j 0)
 
189
                         (zerop (nth-value 1 (truncate j))))))
 
190
           (merror "gamma(~:M) is undefined." j))
 
191
          ((floatp j) (gammafloat j))
 
192
          ((and ($bfloatp j)
 
193
                (or (zerop1 j)
 
194
                    (and (eq ($sign j) '$neg)
 
195
                         (zerop1 (sub j ($truncate j))))))
 
196
           (merror "gamma(~:M) is undefined." j))
 
197
          (($bfloatp j) 
 
198
           ;; Adding 4 digits in the call to bffac. For $fpprec up to about 256
 
199
           ;; and an argument up to about 500.0 the accuracy of the result is
 
200
           ;; better than 10^(-$fpprec).
 
201
           (mfuncall '$bffac (m+ j -1) (+ $fpprec 4)))
 
202
          ((and (complex-number-p j 'float-or-rational-p)
 
203
                (or $numer (floatp ($realpart j)) (floatp ($imagpart j))))
 
204
           (complexify (gamma-lanczos (complex ($float ($realpart j))
 
205
                                               ($float ($imagpart j))))))
 
206
          ((and (complex-number-p j 'bigfloat-or-number-p)
 
207
                (or $numer ($bfloatp ($realpart j)) 
 
208
                           ($bfloatp ($imagpart j))))
 
209
           ;; Adding 4 digits in the call to cbffac. See comment above.
 
210
           (mfuncall '$cbffac 
 
211
                     (add -1 ($bfloat ($realpart j)) 
 
212
                             (mul '$%i ($bfloat ($imagpart j))))
 
213
                     (+ $fpprec 4)))
 
214
          ((eq j '$inf) '$inf) ; Simplify to $inf to be more consistent.
 
215
          ((and $gamma_expand
 
216
                (mplusp j) 
 
217
                (integerp (cadr j)))
 
218
           ;; Expand gamma(z+n) for n an integer.
 
219
           (let ((n (cadr j))
 
220
                 (z (simplify (cons '(mplus) (cddr j)))))
 
221
             (cond 
 
222
               ((> n 0)
 
223
                (mul (simplify (list '($pochhammer) z n))
 
224
                     (simplify (list '(%gamma) z))))
 
225
               ((< n 0)
 
226
                (setq n (- n))
 
227
                (div (mul (power -1 n) (simplify (list '(%gamma) z)))
 
228
                     ;; We factor to get the order (z-1)*(z-2)*...
 
229
                     ;; and not (1-z)*(2-z)*... 
 
230
                     ($factor
 
231
                       (simplify (list '($pochhammer) (sub 1 z) n))))))))
191
232
          ((integerp j)
192
 
           (cond ((> j 0) (simpfact (list '(mfactorial) (1- j)) 1 nil))
 
233
           (cond ((> j 0)
 
234
                  (cond ((<= j $factlim)
 
235
                         ;; Positive integer less than $factlim. Evaluate.
 
236
                         (simplify (list '(mfactorial) (1- j))))
 
237
                         ;; Positive integer greater $factlim. Noun form.
 
238
                        (t (eqtest (list '(%gamma) j) x))))
 
239
                 ;; Negative integer. Throw a Maxima error.
193
240
                 (errorsw (throw 'errorsw t))
194
241
                 (t (merror "gamma(~:M) is undefined" j))))
195
242
          ($numer (gammafloat (fpcofrat j)))
196
243
          ((alike1 j '((rat) 1 2))
197
244
           (list '(mexpt simp) '$%pi j))
198
 
          ((or (ratgreaterp j 1) (ratgreaterp 0 j)) (gammared j))
 
245
          ((and (mnump j)
 
246
                (ratgreaterp $gammalim (simplify (list '(mabs) j)))
 
247
                (or (ratgreaterp j 1) (ratgreaterp 0 j)))
 
248
           ;; Expand for rational numbers less than $gammalim.
 
249
           (gammared j))
199
250
          (t (eqtest (list '(%gamma) j) x)))))
200
251
 
201
252
(defun gamma (y) ;;; numerical evaluation for 0 < y < 1
202
253
  (prog (sum coefs)
203
 
     (setq coefs (list 0.035868343d0 -0.193527817d0 0.48219939d0
204
 
                       -0.75670407d0 0.91820685d0 -0.89705693d0
205
 
                       0.98820588d0 -0.57719165d0))
 
254
     (setq coefs (list 0.035868343 -0.193527817 0.48219939
 
255
                       -0.75670407 0.91820685 -0.89705693
 
256
                       0.98820588 -0.57719165))
206
257
     (unless (atom y) (setq y (fpcofrat y)))
207
258
     (setq sum (car coefs) coefs (cdr coefs))
208
259
     loop (setq sum (+ (* sum y) (car coefs)))
246
297
;;
247
298
;; to handle the case of Re(z) <= 0.
248
299
;;
249
 
;; See http://winnie.fit.edu/~gabdo/gamma.m for some matlab code to
250
 
;; compute this and http://winnie.fit.edu/~gabdo/gamma.txt for a nice
 
300
;; See http://my.fit.edu/~gabdo/gamma.m for some matlab code to
 
301
;; compute this and http://my.fit.edu/~gabdo/gamma.txt for a nice
251
302
;; discussion of Lanczos method and an improvement of Lanczos method.
252
303
;;
253
304
;;
257
308
;; accuracy.
258
309
 
259
310
(defun gamma-lanczos (z)
260
 
  (declare (type (complex double-float) z)
 
311
  (declare (type (complex flonum) z)
261
312
           (optimize (safety 3)))
262
 
  (let ((c (make-array 15 :element-type 'double-float
 
313
  (let ((c (make-array 15 :element-type 'flonum
263
314
                       :initial-contents
264
 
                       '(0.99999999999999709182d0
265
 
                         57.156235665862923517d0
266
 
                         -59.597960355475491248d0
267
 
                         14.136097974741747174d0
268
 
                         -0.49191381609762019978d0
269
 
                         .33994649984811888699d-4
270
 
                         .46523628927048575665d-4
271
 
                         -.98374475304879564677d-4
272
 
                         .15808870322491248884d-3
273
 
                         -.21026444172410488319d-3
274
 
                         .21743961811521264320d-3
275
 
                         -.16431810653676389022d-3
276
 
                         .84418223983852743293d-4
277
 
                         -.26190838401581408670d-4
278
 
                         .36899182659531622704d-5))))
279
 
    (declare (type (simple-array double-float (15)) c))
 
315
                       '(0.99999999999999709182
 
316
                         57.156235665862923517
 
317
                         -59.597960355475491248
 
318
                         14.136097974741747174
 
319
                         -0.49191381609762019978
 
320
                         .33994649984811888699e-4
 
321
                         .46523628927048575665e-4
 
322
                         -.98374475304879564677e-4
 
323
                         .15808870322491248884e-3
 
324
                         -.21026444172410488319e-3
 
325
                         .21743961811521264320e-3
 
326
                         -.16431810653676389022e-3
 
327
                         .84418223983852743293e-4
 
328
                         -.26190838401581408670e-4
 
329
                         .36899182659531622704e-5))))
 
330
    (declare (type (simple-array flonum (15)) c))
280
331
    (if (minusp (realpart z))
281
332
        ;; Use the reflection formula
282
333
        ;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
284
335
        ;; Gamma(z) = pi/z/sin(pi*z)/Gamma(-z)
285
336
        ;;
286
337
        ;; If z is a negative integer, Gamma(z) is infinity.  Should
287
 
        ;; we test for this?  Throw an error?  What
288
 
        (/ (float pi 1d0)
289
 
           (* (- z) (sin (* (float pi 1d0) z))
 
338
        ;; we test for this?  Throw an error?
 
339
        ;; The test must be done by the calling routine.
 
340
        (/ (float pi)
 
341
           (* (- z) (sin (* (float pi) z))
290
342
              (gamma-lanczos (- z))))
291
343
        (let* ((z (- z 1))
292
344
               (zh (+ z 1/2))
293
345
               (zgh (+ zh 607/128))
294
 
               (zp (expt zgh (/ zh 2)))
295
346
               (ss 
296
 
                (do ((sum 0d0)
 
347
                (do ((sum 0.0)
297
348
                     (pp (1- (length c)) (1- pp)))
298
349
                    ((< pp 1)
299
350
                     sum)
300
351
                  (incf sum (/ (aref c pp) (+ z pp))))))
301
 
          (* (sqrt (float (* 2 pi) 1d0))
302
 
             (+ ss (aref c 0))
303
 
             (* zp (exp (- zgh)) zp))))))
 
352
          (let ((result 
 
353
                 ;; We check for an overflow. The last positive value in 
 
354
                 ;; double-float precicsion for which Maxima can calculate 
 
355
                 ;; gamma is ~171.6243 (CLISP 2.46 and GCL 2.6.8)
 
356
                 (ignore-errors
 
357
                   (let ((zp (expt zgh (/ zh 2))))
 
358
                     (* (sqrt (float (* 2 pi)))
 
359
                        (+ ss (aref c 0))
 
360
                        (* (/ zp (exp zgh)) zp))))))
 
361
            (cond ((null result)
 
362
                   ;; No result. Overflow.
 
363
                   (merror "Overflow in `gamma-lanczos'."))
 
364
                  ((or (float-nan-p (realpart result))
 
365
                       (float-inf-p (realpart result)))
 
366
                   ;; Result, but beyond extreme values. Overflow.
 
367
                   (merror "Overflow in `gamma-lanczos'."))
 
368
                  (t result)))))))
304
369
 
305
370
(defun gammafloat (a)
306
 
  (realpart (gamma-lanczos (complex a 0d0))))
 
371
  (let ((a (float a)))
 
372
    (cond ((minusp a)
 
373
           (/ (float (- pi))
 
374
              (* a (sin (* (float pi) a)))
 
375
              (gammafloat (- a))))
 
376
          ((< a 10)
 
377
           (slatec::dgamma a))
 
378
          (t
 
379
           (let ((result
 
380
                  (let ((c (* (sqrt (* 2 (float pi)))
 
381
                              (exp (slatec::d9lgmc a)))))
 
382
                    (let ((v (expt a (- (* .5d0 a) 0.25d0))))
 
383
                      (* v
 
384
                         (/ v (exp a))
 
385
                         c)))))
 
386
             (if (or (float-nan-p result)
 
387
                     (float-inf-p result))
 
388
                 (merror "Overflow in `gammafloat'")
 
389
                 result))))))
 
390
 
307
391
 
308
392
(declare-top (special $numer $trigsign))
309
393
 
310
 
(defmfun simperf (x vestigial z &aux y)
311
 
  (declare (ignore vestigial))
312
 
  (oneargcheck x)
313
 
  (setq y (simpcheck (cadr x) z))
314
 
  (cond ((zerop1 y) y)
315
 
        ((or (floatp y) (and $numer (integerp y))) (erf (float y)))
316
 
        ((eq y '$inf) 1)
317
 
        ((eq y '$minf) -1)
318
 
        ((and $trigsign (mminusp* y)) (neg (list '(%erf simp) (neg y))))
319
 
        (t (eqtest (list '(%erf) y) x))))
320
 
 
321
 
 
322
 
(defmfun erf (y)
323
 
  (slatec:derf (float y 1d0)))
324
 
 
325
 
(defmfun erfc (y)
326
 
  (slatec:derfc (float y 1d0)))
 
394
;;; The Error functions Erf, Erfc, Erfi and the function Generalized Erf have 
 
395
;;; been implemented as simplifying functions. The code is moved to the file
 
396
;;; gamma.lisp and commented out at this place. 10/2008 DK.
 
397
 
 
398
;(defmfun simperf (x vestigial z &aux y)
 
399
;  (declare (ignore vestigial))
 
400
;  (oneargcheck x)
 
401
;  (setq y (simpcheck (cadr x) z))
 
402
;  (cond ((zerop1 y) y)
 
403
;       ((or (floatp y) (and $numer (integerp y))) (erf (float y)))
 
404
;       ((eq y '$inf) 1)
 
405
;       ((eq y '$minf) -1)
 
406
;       ;;((and $trigsign (mminusp* y)) (neg (list '(%erf simp) (neg y))))
 
407
;       ((and $trigsign (great (neg y) y)) (neg (take '(%erf) (neg y))))
 
408
;       (t (eqtest (list '(%erf) y) x))))
 
409
 
 
410
 
 
411
;(defmfun erf (y)
 
412
;  (slatec:derf (float y)))
 
413
 
 
414
;(defmfun erfc (y)
 
415
;  (slatec:derfc (float y)))
327
416
 
328
417
(defmfun $zeromatrix (m n) ($ematrix m n 0 1 1))
329
418