share/numeric/brmbrg.lisp

   1 ;;; -*- Lisp -*-
   2
   3 (defvar $brombergit 11)
   4 (defvar $brombergmin 0)
   5 (defvar $brombergtol '((bigfloat simp 56.) 59029581035870565. -13.)) ;1.b-4
   6 (defvar $brombergabs '((bigfloat simp 56.) 0 0)) ; 0.0b0
   7
   8 (defun fpscale (x m)
   9   (if (equal (car x) 0)
  10       x
  11       (list (car x) (+ (cadr x) m))))
  12
  13 (defun bfmeval3 (x1)
  14   (cond (($bfloatp (setq x1 ($bfloat (meval x1)))) (cdr x1))
  15         (t (displa x1) (merror "bromberg: encountered a non-bigfloat."))))
  16
  17 (defun bqeval3 (f x)
  18   (let*
  19     ((bf-x (bcons x))
  20      (result (funcall f bf-x)))
  21     (if (atom result)
  22       (merror "bromberg: integrand evaluates to a nonnumeric value at ~M" bf-x)
  23       (cdr result))))
  24
  25 (defun $bromberg (&rest l1)
  26   (or (= (length l1) 4) (= (length l1) 3)
  27       (merror "bromberg: wrong number of arguments."))
  28   (let ((fun) (a) (b) (x)
  29         ($bfloat t)
  30         ($float2bf t)
  31         (lim (cdr ($bfloat $brombergtol)))
  32         (limabs (cdr ($bfloat $brombergabs)))
  33         (tt (make-array $brombergit))
  34         (rr (make-array $brombergit))
  35         (zero (intofp 0))
  36         (one (intofp 1))
  37         (three (intofp 3)))
  38     (declare (special $bfloat $float2bf))
  39     (cond ((= (length l1) 4)
  40            (when ($constantp (cadr l1))
  41              (merror
  42                "bromberg: variable of integration cannot be a constant; found ~M"
  43                (cadr l1)))
  44            (setq fun (coerce-bfloat-fun (car l1) `((mlist) ,(cadr l1)))
  45                  l1 (cdr l1)))
  46           (t
  47            (setq fun (coerce-bfloat-fun (car l1)))))
  48     (setq a (bfmeval3 (cadr l1))
  49           b (bfmeval3 (caddr l1))
  50           x (fpdifference b a))
  51     (setf (aref tt 0)
  52           (fpscale (fptimes* x (fpplus (bqeval3 fun b)
  53                                        (bqeval3 fun a)))
  54                    -1))
  55     (setf (aref rr 0) (fptimes* x (bqeval3 fun (fpscale (fpplus b a) -1))))
  56     (do ((l 1 (1+ l))
  57          (m 4 (* m 2))
  58          (y) (z) (cerr))
  59         ((= l $brombergit) (merror "bromberg: failed to converge."))
  60       (setq y (intofp m) z (fpquotient x y))
  61       (setf (aref tt l)
  62             (fpscale (fpplus (aref tt (1- l)) (aref rr (1- l))) -1))
  63       (setf (aref rr l) zero)
  64       (do ((i 1 (+ i 2)))
  65           ((> i m))
  66         (setf (aref rr l)
  67               (fpplus (bqeval3 fun (fpplus (fptimes* z (intofp i)) a))
  68                       (aref rr l))))
  69       (setf (aref rr l) (fpscale (fptimes* z (aref rr l)) 1))
  70       (setq y zero)
  71       (do ((k l (1- k)))
  72           ((zerop k))
  73         (setq y (fpplus (fpscale y 2) three))
  74         (setf (aref tt (1- k))
  75               (fpplus (fpquotient
  76                        (fpdifference (aref tt k) (aref tt (1- k)))
  77                        y)
  78                       (aref tt k)))
  79         (setf (aref rr (1- k))
  80               (fpplus (fpquotient
  81                        (fpdifference (aref rr k) (aref rr (1- k)))
  82                        y)
  83                       (aref rr k))))
  84       (setq y (fpscale (fpplus (aref tt 0) (aref rr 0)) -1))
  85       (when (and
  86              (or (not
  87                   (fplessp
  88                    limabs
  89                    (setq cerr
  90                          (fpabs (fpdifference (aref tt 0)
  91                                               (aref rr 0))))))
  92                  (not (fplessp lim
  93                                (fpquotient
  94                                 cerr
  95                                 (cond ((equal y '(0 0)) one)
  96                                       (t (fpabs y)))))))
  97              (> l $brombergmin))
  98         (return (bcons y))))))