archive/share/trash/qq.lisp

   1 ;;;;;;;;;;;;;;;;;;; -*- Mode: Lisp; Package: Macsyma -*- ;;;;;;;;;;;;;;;;;;;
   2 ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
   3 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   4
   5 (in-package :maxima)
   6
   7 (macsyma-module qq)
   8
   9 ;; (load-macsyma-macros Numerm)
  10
  11 ;;; 10:55 pm Sept 25, 1981 - LPH & GJC
  12 ;;; added $numer and $float set to T so that quanc8(x,x,0,1/2); works. this
  13 ;;; is done inside the prog for $quanc8 so that they are put back when done.
  14 ;;; 3:30 pm Feb 12, 1982 - JPG, LPH, & GJC
  15 ;;; removed the $numer:$float:true and replaced with mto-float which is
  16 ;;; defined in maxsrc;numer > .
  17 ;;; 3:45 am April 11, 1982 - LPH
  18 ;;; fixed an error in estimating TOLERR due to neglect of sub-interval size
  19 ;;; when abserr check is made.
  20
  21 (defmvar $quanc8_flag 0.0
  22   "integer part-how many regions failed to converge. fractional part-
  23   how much of the interval was left when it failed")
  24 (defmvar $quanc8_errest 0.0
  25   "an estimated error based on the difference between one- and two-panel
  26   calculation of an interval, summed over the whole region")
  27 (defmvar $quanc8_abserr 1.0e-8
  28   "absolute error tolerance for reasonable single-precision answers")
  29 (defmvar $quanc8_relerr 1.0e-4
  30   "relative error condition for reasonable run-times")
  31
  32 ;; (DEFUN ($QUANC8 TRANSLATED-MMACRO) (F A B &OPTIONAL (C NIL C-P))
  33 ;;   (IF (NOT C-P)
  34 ;;       `((CALL-QUANC8) ,F ,A ,B)
  35 ;;       `((CALL-QUANC8) ((LAMBDA) ((MLIST) ,A) ,F) ,B ,C)))
  36
  37 (defmspec $quanc8 (form)
  38   (if (cdddr (setq form (cdr form)))
  39       (apply #'call-quanc8
  40              (meval `((lambda) ((mlist) ,(cadr form)) ,(car form)))
  41              (mapcar #'meval (cddr form)))
  42       (apply #'call-quanc8 (mapcar #'meval form))))
  43
  44 (def%tr $quanc8 (form)
  45   `($float
  46     ,@(cdr (tr-lisp-function-call
  47             (if (cdddr (setq form (cdr form)))
  48                 `((call-quanc8)
  49                   ((lambda) ((mlist) ,(cadr form)) ,(car form)) ,@(cddr form))
  50                 `((call-quanc8) ,@form)) nil))))
  51
  52 (defvar quanc8-free-list ()
  53   "For efficient calls to quanc8 keep arrays we need on a free list.")
  54
  55 (defvar quanc8-|^]| ())
  56 (defun quanc8-|^]| () (setq quanc8-|^]| t))
  57
  58 (defun call-quanc8 (fun a b)
  59   (bind-tramp1$
  60    fun fun
  61    (let ((vals (if quanc8-free-list
  62                    (pop quanc8-free-list)
  63                    (list (*array nil 'flonum 17.)
  64                          (*array nil 'flonum 17.)
  65                          (*array nil 'flonum 9. 31.)
  66                          (*array nil 'flonum 9. 31.)
  67                          (*array nil 'flonum 32.))))
  68          (user-timesofar (cons #'quanc8-|^]| user-timesofar)))
  69      (prog1
  70       (quanc8 fun
  71               (mto-float a)
  72               (mto-float b)
  73               (car vals)
  74               (car (cdr vals))
  75               (car (cddr vals))
  76               (car (cdddr vals))
  77               (car (cddddr vals)))
  78       (push vals quanc8-free-list)))))
  79
  80 (defun quanc8 (fun a b x-arr f-arr xsave-arr fsave-arr qright-arr)
  81   (declare (type (simple-array cl:float)
  82                  x-arr f-arr xsave-arr fsave-arr qright-arr))
  83   ;; local macros for typing convenience.
  84   (macrolet ((x (j) `(arraycall flonum x-arr ,j))
  85              (f (j) `(arraycall flonum f-arr ,j))
  86              (xsave (j k) `(arraycall flonum xsave-arr ,j ,k))
  87              (fsave (j k) `(arraycall flonum fsave-arr ,j ,k))
  88              (qright (j) `(arraycall flonum qright-arr ,j)))
  89     ;; Rudimentary (non-ansi GCL compatible) error handling.
  90     (let (errset)
  91       (or (car (errset
  92     (prog ((levmin 1.) (levmax 30.) (levout 6.) (nomax 5000.) (nofin 0)
  93          (w0 (//$ 3956.0 14175.0)) (w1 (//$ 23552.0 14175.0))
  94          (w2 (//$ -3712.0 14175.0)) (w3 (//$ 41984.0 14175.0))
  95          (w4 (//$ -18160.0 14175.0))
  96          (result 0.0) (cor11 0.0) (area 0.0)
  97          (nofun 0.) (lev 0.) (nim 1.) (qprev 0.0)
  98          (stone (//$ (-$ b a) 16.0))
  99          (i 0.)
 100          (step 0.0) (qleft 0.0) (qnow 0.0) (qdiff 0.0)
 101          (esterr 0.0) (tolerr 0.0) (temp 0.0)
 102          ($numer t) ($float t))
 103
 104         (declare (cl:float w0 w1 w2 w3 w4 result cor11 area qprev stone step
 105                          qleft qnow qdiff esterr tolerr temp)
 106                  (fixnum i levmin levmax levout nomax nofin nofun lev nim)
 107                  (boolean $numer $float))
 108
 109         (setq
 110          $quanc8_flag 0.0
 111          $quanc8_errest 0.0
 112          nofin (- nomax (* 8 (+ levmax (* -1. levout) (expt 2. (1+ levout))))))
 113
 114         (cond ((= a b)
 115                (return 0.0)))
 116
 117         (setf (x 0.)
 118               a)
 119
 120         (setf (x 16.)
 121               b)
 122
 123         (setf (x 8.)
 124               (*$ 0.5 (+$ (x 0.)
 125                           (x 16.))))
 126
 127         (setf (x 4.)
 128               (*$ 0.5 (+$ (x 0.)
 129                           (x 8.))))
 130
 131         (setf (x 12.)
 132               (*$ 0.5 (+$ (x 16.)
 133                           (x 8.))))
 134
 135         (setf (x 2.)
 136               (*$ 0.5 (+$ (x 0.)
 137                           (x 4.))))
 138
 139         (setf (x 6.)
 140               (*$ 0.5 (+$ (x 4.)
 141                           (x 8.))))
 142
 143         (setf (x 10.)
 144               (*$ 0.5 (+$ (x 12.)
 145                           (x 8.))))
 146
 147         (setf (x 14.)
 148               (*$ 0.5 (+$ (x 12.)
 149                           (x 16.))))
 150
 151
 152   do-25
 153         (when quanc8-|^]|
 154               (mtell "QUANC8 calculating at X= ~S" (x i))
 155               (setq quanc8-|^]| nil))
 156
 157         (setf (f i)
 158               (fcall$ fun (x i)))
 159
 160         (setq i (+ 2. i))
 161
 162         (if (> i 16.)
 163             (go do-25-done))
 164
 165         (go do-25)
 166
 167 do-25-done
 168         (setq nofun 9.)
 169
 170  tag-30
 171         (setq i 1.)
 172
 173   do-30
 174         (setf (x i)
 175               (*$ 0.5 (+$ (x (1- i))
 176                           (x (1+ i)))))
 177         (when quanc8-|^]|
 178               (mtell "QUANC8 calculating at X= ~S" (x i))
 179               (setq quanc8-|^]| nil))
 180
 181         (setf (f i)
 182               (fcall$ fun (x i)))
 183
 184         (setq i (+ 2. i))
 185
 186         (if (> i 15.)
 187             (go do-30-done))
 188
 189         (go do-30)
 190
 191 do-30-done
 192         (setq nofun (+ 8. nofun))
 193
 194         (setq step (//$ (-$ (x 16.)
 195                             (x 0.))
 196                         16.0))
 197
 198         (setq qleft (*$ step (+$ (*$ w0 (+$ (f 0.)
 199                                             (f 8.)))
 200                                  (*$ w1 (+$ (f 1.)
 201                                             (f 7.)))
 202                                  (*$ w2 (+$ (f 2.)
 203                                             (f 6.)))
 204                                  (*$ w3 (+$ (f 3.)
 205                                             (f 5.)))
 206                                  (*$ w4 (f 4.)))))
 207
 208         (setf (qright (1+ lev))
 209               (*$ step (+$ (*$ w0 (+$ (f 8.)
 210                                       (f 16.)))
 211                            (*$ w1 (+$ (f 9.)
 212                                       (f 15.)))
 213                            (*$ w2 (+$ (f 10.)
 214                                       (f 14.)))
 215                            (*$ w3 (+$ (f 11.)
 216                                       (f 13.)))
 217                            (*$ w4 (f 12.)))))
 218
 219         (setq qnow (+$ qleft (qright (1+ lev))))
 220
 221         (setq qdiff (-$ qnow qprev))
 222
 223         (setq area (+$ area qdiff))
 224
 225         (setq esterr (//$ (abs qdiff)
 226                           1023.0))
 227
 228         (setq tolerr (*$
 229                       (//$ step stone)
 230                       (max $quanc8_abserr
 231                            (*$ (abs area)
 232                                $quanc8_relerr))))
 233
 234         (if (< lev levmin)
 235             (go tag-50))
 236
 237         (if (or (> lev levmax)
 238                 (= lev levmax))
 239             (go tag-62))
 240
 241         (if (> nofun nofin)
 242             (go tag-60))
 243
 244         (if (< esterr tolerr)
 245             (go tag-70))
 246
 247  tag-50
 248         (setq nim (* 2. nim)
 249               lev (1+ lev))
 250
 251         (setq i 1.)
 252
 253   do-52
 254         (setf (fsave i lev)
 255               (f (+ 8. i)))
 256
 257         (setf (xsave i lev)
 258               (x (+ 8. i)))
 259
 260         (setq i (1+ i))
 261
 262         (if (> i 8.)
 263             (go do-52-done))
 264
 265         (go do-52)
 266
 267 do-52-done
 268         (setq qprev qleft)
 269
 270         (setq i 1.)
 271
 272   do-55
 273         (setf (f (+ 18. (* -2. i)))
 274               (f (+ 9. (* -1. i))))
 275
 276         (setf (x (+ 18. (* -2. i)))
 277               (x (+ 9. (* -1. i))))
 278
 279         (setq i (1+ i))
 280
 281         (if (> i 8.)
 282             (go tag-30))
 283
 284         (go do-55)
 285
 286  tag-60
 287         (setq nofin (* 2. nofin))
 288
 289         (setq levmax levout)
 290
 291         (setq $quanc8_flag (+$ $quanc8_flag (//$ (-$ b (x 0.))
 292                                  (-$ b a))))
 293
 294         (go tag-70)
 295
 296  tag-62
 297         (setq $quanc8_flag (+$ $quanc8_flag 1.0))
 298
 299  tag-70
 300         (setq result (+$ result qnow)
 301               $quanc8_errest (+$ $quanc8_errest esterr)
 302               cor11 (+$ cor11 (//$ qdiff 1023.0)))
 303
 304  tag-72
 305         (if (= nim (* 2. (// nim 2.)))
 306             (go tag-75))
 307
 308         (setq nim (// nim 2.))
 309
 310         (setq lev (1- lev))
 311
 312         (go tag-72)
 313
 314  tag-75
 315         (setq nim (1+ nim))
 316
 317         (if (or (< lev 0.)
 318                 (= lev 0.))
 319             (go tag-80))
 320
 321         (setq qprev (qright lev))
 322
 323         (setq i 1.)
 324
 325         (setf (x 0.)
 326               (x 16.))
 327
 328         (setf (f 0.)
 329               (f 16.))
 330
 331   do-78
 332         (setf (f (* 2. i))
 333               (fsave i lev))
 334
 335         (setf (x (* 2. i))
 336               (xsave i lev))
 337
 338         (setq i (1+ i))
 339
 340         (if (> i 8.)
 341             (go tag-30))
 342
 343         (go do-78)
 344
 345  tag-80
 346         (setq result (+$ result cor11))
 347
 348         (if (= $quanc8_errest 0.0)
 349             (return result))
 350
 351  tag-82
 352     (setq temp (+$ $quanc8_errest (abs result)))
 353
 354         (if (not (= temp (abs result)))
 355             (return result))
 356
 357         (setq $quanc8_errest (*$ 2.0 $quanc8_errest))
 358
 359         (go tag-82))))
 360
 361           ;; For whatever reason and with a suitable
 362           ;; discrete meaning of convergence...
 363           (merror "QUANC8 failed to converge.")))))