src/numerical/slatec/quadpack.lisp

   1 ;;; Maxima interface to quadpack integration
   2
   3 (in-package :maxima)
   4
   5 ;; This should really be necessary, but the order of compilation in
   6 ;; maxima.system is messed up so that this file is compiled before
   7 ;; plot.lisp.  Rearranging the order of the modules to compile
   8 ;; plot.lisp before quadpack.lisp causes other problems.
   9 (declaim (special *plot-realpart*))
  10
  11 (defvar *debug-quadpack*
  12   nil
  13   "Set to non-NIL to enable printing of the error object when the
  14   Slatec routines throw an error.")
  15
  16 #-(or gcl ecl)
  17 (defmacro get-integrand (fun var)
  18   `(compile nil (coerce-float-fun ,fun `((mlist) ,,var))))
  19
  20 #+(or gcl ecl)
  21 (defmacro get-integrand (fun var)
  22   `(coerce-float-fun ,fun `((mlist) ,,var)))
  23
  24 (defun float-or-lose (val)
  25   (let ((v ($float val)))
  26     (if (numberp v)
  27         v
  28         (error (intl:gettext "~A is not a real floating-point number") val))))
  29
  30 ;; Convert the maxima function FUN of the one variable VAR to a
  31 ;; compiled lisp function returning a float that is suitable for use
  32 ;; as the integrand for the quadpack routines.  If the integrand
  33 ;; cannot be evaluated at the specified point, signal an error.
  34 (defmacro float-integrand-or-lose (name fun var)
  35   (let ((f (gensym "COMPILED-FUN-"))
  36         (x (gensym "VAR-"))
  37         (result (gensym "RESULT-")))
  38     `(let ((,f (get-integrand ,fun ,var)))
  39        (lambda (,x)
  40          (let ((,result (funcall ,f ,x)))
  41            (unless (cl:floatp ,result)
  42              (merror (intl:gettext "~M: Numerical evaluation of ~M at ~M is not a float or is not defined") ,name ,fun ,x))
  43            (float ,result))))))
  44
  45 (defmfun $quad_qag (fun var a b key &key
  46                         (epsrel 1e-8)
  47                         (limit 200)
  48                         (epsabs 0.0))
  49   (quad_argument_check %%pretty-fname fun var a b)
  50   (unless (and (integerp key) (>= key 1) (<= key 6))
  51     (merror (intl:gettext "~M: Key must be 1, 2, 3, 4, 5, or 6; found: ~M")
  52             %%pretty-fname key))
  53   (let* ((lenw (* 4 limit))
  54          (work (make-array lenw :element-type 'flonum))
  55          (iwork (make-array limit :element-type 'f2cl-lib:integer4))
  56          (*plot-realpart* nil))
  57     (handler-case
  58         (multiple-value-bind (junk z-a z-b z-epsabs z-epsrel z-key result abserr neval ier
  59                                    z-limit z-lenw last)
  60             (slatec:dqag (float-integrand-or-lose '$quad_qag fun var)
  61                          (float-or-lose a)
  62                          (float-or-lose b)
  63                          (float-or-lose epsabs)
  64                          (float-or-lose epsrel)
  65                          key
  66                          0.0 0.0 0 0
  67                          limit lenw 0 iwork work)
  68           (declare (ignore junk z-a z-b z-epsabs z-epsrel z-key z-limit z-lenw last))
  69           (list '(mlist) result abserr neval ier))
  70       (error (e)
  71         (when *debug-quadpack*
  72           (format t "~S" e))
  73         `(($quad_qag) ,fun ,var ,a ,b ,key
  74           ((mequal) $epsrel ,epsrel)
  75           ((mequal) $epsabs ,epsabs)
  76           ((mequal) $limit ,limit))))))
  77
  78 (defmfun $quad_qags (fun var a b &key
  79                   (epsrel 1e-8)
  80                   (limit 200)
  81                   (epsabs 0.0))
  82   (quad_argument_check %%pretty-fname fun var a b)
  83   (let* ((lenw (* 4 limit))
  84          (work (make-array lenw :element-type 'flonum))
  85          (iwork (make-array limit :element-type 'f2cl-lib:integer4))
  86          (*plot-realpart* nil))
  87     (handler-case
  88         (multiple-value-bind (junk z-a z-b z-epsabs z-epsrel result abserr neval ier
  89                               z-limit z-lenw last)
  90             (slatec:dqags (float-integrand-or-lose '$quad_qags fun var)
  91                           (float-or-lose a)
  92                           (float-or-lose b)
  93                           (float-or-lose epsabs)
  94                           (float-or-lose epsrel)
  95                           0.0 0.0 0 0
  96                           limit lenw 0 iwork work)
  97           (declare (ignore junk z-a z-b z-epsabs z-epsrel z-limit z-lenw last))
  98           (list '(mlist) result abserr neval ier))
  99       (error (e)
 100         (when *debug-quadpack*
 101           (format t "~S" e))
 102         `(($quad_qags) ,fun ,var ,a ,b
 103           ((mequal) $epsrel ,epsrel)
 104           ((mequal) $epsabs ,epsabs)
 105           ((mequal) $limit ,limit))))))
 106
 107 (defmfun $quad_qagi (fun var a b &key
 108                       (epsrel 1e-8)
 109                       (limit 200)
 110                       (epsabs 0.0))
 111   (quad_argument_check %%pretty-fname fun var a b)
 112   ;; Massage the limits a and b into what Quadpack QAGI wants.
 113   (flet ((fixup (low high)
 114            (let (bnd inf)
 115              ;; Cases to handle: (minf, x), (x, inf), (minf, inf).
 116              ;; Everything else is an error.
 117              (cond ((eq ($limit low) '$minf)
 118                     (cond ((eq ($limit high) '$inf)
 119                            (setf bnd 0)
 120                            (setf inf 2))
 121                           (t
 122                            (setq bnd ($float high))
 123                            (setq inf ($limit low)))))
 124                    ((eq ($limit high) '$inf)
 125                     (setq bnd ($float low))
 126                     (setq inf ($limit high)))
 127                    (t
 128                     (merror "~M: Unexpected limits of integration: ~M, ~M~%"
 129                             %%pretty-fname low high)))
 130              (values bnd inf))))
 131
 132     (multiple-value-bind (bound inf-type)
 133         (fixup a b)
 134       (let* ((lenw (* 4 limit))
 135              (work (make-array lenw :element-type 'flonum))
 136              (iwork (make-array limit :element-type 'f2cl-lib:integer4))
 137              (infinity (ecase inf-type
 138                          ((1 $inf)
 139                           ;; Interval is [bound, infinity]
 140                           1)
 141                          ((-1 $minf)
 142                           ;; Interval is [-infinity, bound]
 143                           -1)
 144                          ((2 $both)
 145                           ;; Interval is [-infinity, infinity]
 146                           2)))
 147              (*plot-realpart* nil))
 148         (handler-case
 149             (multiple-value-bind (junk z-bound z-inf z-epsabs z-epsrel
 150                                        result abserr neval ier
 151                                        z-limit z-lenw last)
 152                 (slatec:dqagi (float-integrand-or-lose '$quad_qagi fun var)
 153                               (float-or-lose bound)
 154                               infinity
 155                               (float-or-lose epsabs)
 156                               (float-or-lose epsrel)
 157                               0.0 0.0 0 0
 158                               limit lenw 0 iwork work)
 159               (declare (ignore junk z-bound z-inf z-epsabs z-epsrel
 160                                z-limit z-lenw last))
 161               (list '(mlist) result abserr neval ier))
 162           (error (e)
 163             (when *debug-quadpack*
 164               (format t "~S" e))
 165             `(($quad_qagi) ,fun ,var ,a ,b
 166               ((mequal) $epsrel ,epsrel)
 167               ((mequal) $epsabs ,epsabs)
 168               ((mequal) $limit ,limit))))))))
 169
 170 (defmfun $quad_qawc (fun var c a b &key
 171                   (epsrel 1e-8)
 172                   (limit 200)
 173                   (epsabs 0.0))
 174   (quad_argument_check %%pretty-fname fun var a b)
 175   (let* ((lenw (* 4 limit))
 176          (work (make-array lenw :element-type 'flonum))
 177          (iwork (make-array limit :element-type 'f2cl-lib:integer4))
 178          (*plot-realpart* nil))
 179     (handler-case
 180         (multiple-value-bind (junk z-a z-b z-c z-epsabs z-epsrel result abserr neval ier
 181                                    z-limit z-lenw last)
 182             (slatec:dqawc (float-integrand-or-lose '$quad_qawc fun var)
 183                           (float-or-lose a)
 184                           (float-or-lose b)
 185                           (float-or-lose c)
 186                           (float-or-lose epsabs)
 187                           (float-or-lose epsrel)
 188                           0.0 0.0 0 0
 189                           limit lenw 0 iwork work)
 190           (declare (ignore junk z-a z-b z-c z-epsabs z-epsrel z-limit z-lenw last))
 191           (list '(mlist) result abserr neval ier))
 192       (error (e)
 193         (when *debug-quadpack*
 194           (format t "~S" e))
 195         `(($quad_qawc) ,fun ,var ,c ,a ,b
 196           ((mequal) $epsrel ,epsrel)
 197           ((mequal) $epsabs ,epsabs)
 198           ((mequal) $limit ,limit))))))
 199
 200 (defmfun $quad_qawf (fun var a omega trig &key
 201                   (epsabs 1e-10)
 202                   (limit 200)
 203                   (maxp1 100)
 204                   (limlst 10))
 205   (let* ((leniw limit)
 206          (lenw (+ (* 2 leniw) (* 25 maxp1)))
 207          (work (make-array lenw :element-type 'flonum))
 208          (iwork (make-array leniw :element-type 'f2cl-lib:integer4))
 209          (integr (case trig
 210                    ((1 %cos $cos) 1)
 211                    ((2 %sin $sin) 2)
 212                    (otherwise
 213                     (merror "~M: the name of the trig function should be sin or cos, not ~M."
 214                             %%pretty-fname
 215                             trig))))
 216          (*plot-realpart* nil))
 217     (handler-case
 218         (multiple-value-bind (junk z-a z-omega z-integr
 219                                    epsabs result abserr neval ier
 220                                    z-limlst z-lst
 221                                    z-leniw z-maxp1 z-lenw)
 222             (slatec:dqawf (float-integrand-or-lose '$quad_qawf fun var)
 223                           (float-or-lose a)
 224                           (float-or-lose omega)
 225                           integr
 226                           (float-or-lose epsabs)
 227                           0.0 0.0 0 0
 228                           limlst 0 leniw maxp1 lenw iwork work)
 229           (declare (ignore junk z-a z-omega z-integr epsabs z-limlst z-lst
 230                            z-leniw z-maxp1 z-lenw))
 231           (list '(mlist) result abserr neval ier))
 232       (error (e)
 233         (when *debug-quadpack*
 234           (format t "~S" e))
 235         `(($quad_qawf) ,fun ,var ,a ,omega ,trig
 236           ((mequal) $epsabs ,epsabs)
 237           ((mequal) $limit ,limit)
 238           ((mequal) $maxp1 ,maxp1)
 239           ((mequal) $limlst ,limlst))))))
 240
 241 (defmfun $quad_qawo (fun var a b omega trig_name &key
 242                   (epsrel 1e-8)
 243                   (limit 200)
 244                   (maxp1 100)
 245                   (epsabs 0.0))
 246   (quad_argument_check %%pretty-fname fun var a b)
 247   (let* ((leniw limit)
 248          (lenw (+ (* 2 leniw) (* 25 maxp1)))
 249          (work (make-array lenw :element-type 'flonum))
 250          (iwork (make-array leniw :element-type 'f2cl-lib:integer4))
 251          (integr (case trig_name
 252                    ((%cos $cos) 1)
 253                    ((%sin $sin) 2)
 254                    (otherwise
 255                     (merror "~M: the name of the trig function should be sin or cos, not ~M."
 256                             %%pretty-fname
 257                             trig_name))))
 258          (*plot-realpart* nil))
 259     (handler-case
 260         (multiple-value-bind (junk z-a z-b z-omega z-integr z-epsabs z-epsrel
 261                                    result abserr neval ier
 262                                    z-leniw z-maxp1 z-lenw z-lst)
 263             (slatec:dqawo (float-integrand-or-lose '$quad_qawo fun var)
 264                           (float-or-lose a)
 265                           (float-or-lose b)
 266                           (float-or-lose omega)
 267                           integr
 268                           (float-or-lose epsabs)
 269                           (float-or-lose epsrel)
 270                           0.0 0.0 0 0
 271                           leniw maxp1 lenw 0 iwork work)
 272           (declare (ignore junk z-a z-b z-omega z-integr z-epsabs z-epsrel
 273                            z-lst z-leniw z-maxp1 z-lenw))
 274           (list '(mlist) result abserr neval ier))
 275       (error (e)
 276         (when *debug-quadpack*
 277           (format t "~S" e))
 278         `(($quad_qawo) ,fun ,var ,a ,b ,omega ,trig_name
 279           ((mequal) $epsrel ,epsrel)
 280           ((mequal) $epsabs ,epsabs)
 281           ((mequal) $limit ,limit)
 282           ((mequal) $maxp1 ,maxp1))))))
 283
 284 (defmfun $quad_qaws (fun var a b alfa beta wfun &key
 285                   (epsrel 1e-8)
 286                   (limit 200)
 287                   (epsabs 0.0))
 288   (quad_argument_check %%pretty-fname fun var a b)
 289   (let* ((lenw (* 4 limit))
 290          (work (make-array lenw :element-type 'flonum))
 291          (iwork (make-array limit :element-type 'f2cl-lib:integer4))
 292          (*plot-realpart* nil))
 293     (handler-case
 294         (multiple-value-bind (junk z-a z-b z-alfa z-beta z-int z-epsabs z-epsrel
 295                                    result abserr neval ier
 296                                    z-limit z-lenw last)
 297             (slatec:dqaws (float-integrand-or-lose '$quad_qaws fun var)
 298                           (float-or-lose a)
 299                           (float-or-lose b)
 300                           (float-or-lose alfa)
 301                           (float-or-lose beta)
 302                           wfun
 303                           (float-or-lose epsabs)
 304                           (float-or-lose epsrel)
 305                           0.0 0.0 0 0
 306                           limit lenw 0 iwork work)
 307           (declare (ignore junk z-a z-b z-alfa z-beta z-int z-epsabs z-epsrel
 308                            z-limit z-lenw last))
 309           (list '(mlist) result abserr neval ier))
 310       (error (e)
 311         (when *debug-quadpack*
 312           (format t "~S" e))
 313         `(($quad_qaws) ,fun ,var ,a ,b ,alfa ,beta ,wfun
 314           ((mequal) $epsrel ,epsrel)
 315           ((mequal) $epsabs ,epsabs)
 316           ((mequal) $limit ,limit))))))
 317
 318 (defmfun $quad_qagp (fun var a b points
 319                   &key (epsrel 1e-8) (epsabs 0.0) (limit 200))
 320   (quad_argument_check %%pretty-fname fun var a b)
 321   (let* ((npts2 (+ 2 (length (cdr points))))
 322          (p (make-array npts2 :element-type 'flonum))
 323          (leniw (max limit (- (* 3 npts2) 2)))
 324          (lenw (- (* 2 leniw) npts2))
 325          (work (make-array lenw :element-type 'flonum))
 326          (iwork (make-array limit :element-type 'f2cl-lib:integer4))
 327          (*plot-realpart* nil))
 328     ;; Check the list of singular points.  Each point must be a real
 329     ;; in the open interval (a, b).  If the point is valid, save it to
 330     ;; the array P.  Otherwise accumulate a list of the invalid points
 331     ;; so we can print a nice error message with a list of the invalid
 332     ;; points.
 333     (let* ((a (float-or-lose a))
 334            (b (float-or-lose b))
 335            (invalid-points
 336              (loop for xp in (cdr points)
 337                    for x = (float-or-lose xp)
 338                    for k from 0
 339                    if (< a x b)
 340                      do (setf (aref p k) x)
 341                    else
 342                      collect x)))
 343       ;; If there are invalid points, throw an error.
 344       (when invalid-points
 345         (merror
 346          (intl:gettext "~M: singular points must be in the open interval (~M, ~M):  ~M")
 347          %%pretty-fname a b (make-mlist-l invalid-points)))
 348
 349       (handler-case
 350           (multiple-value-bind (junk z-a z-b z-npts z-points z-epsabs z-epsrel
 351                                 result abserr neval ier
 352                                 z-leniw z-lenw last)
 353               (slatec:dqagp (float-integrand-or-lose '$quad_qagp fun var)
 354                             (float-or-lose a)
 355                             (float-or-lose b)
 356                             npts2
 357                             p
 358                             (float-or-lose epsabs)
 359                             (float-or-lose epsrel)
 360                             0.0 0.0 0 0
 361                             leniw lenw 0 iwork work)
 362             (declare (ignore junk z-a z-b z-npts z-points z-epsabs z-epsrel
 363                              z-leniw z-lenw last))
 364             (make-mlist result abserr neval ier))
 365         (error (e)
 366           (when *debug-quadpack*
 367             (format t "~S" e))
 368           `(($quad_qagp) ,fun ,var ,a ,b ,points
 369             ((mequal) $epsrel ,epsrel)
 370             ((mequal) $epsabs ,epsabs)
 371             ((mequal) $limit ,limit)))))))
 372
 373 ;; error checking similar to that done by $defint
 374 (defun quad_argument_check (name exp var ll ul)
 375   (setq exp (ratdisrep exp))
 376   (setq var (ratdisrep var))
 377   (setq ll (ratdisrep ll))
 378   (setq ul (ratdisrep ul))
 379   (cond (($constantp var)
 380          (merror "~M: Variable of integration not a variable: ~M" name var)))
 381   (cond ((not (or ($subvarp var) (atom var)))
 382          (merror "~M: Improper variable of integration: ~M" name var))
 383         ((or (among var ul)
 384              (among var ll))
 385          (merror "~M: Limit contains variable of integration: ~M" name var))))
 386
 387 (defun quad-control (parameter &optional new-value)
 388   (values
 389    (slatec:j4save (case parameter
 390                     ($current_error 1)
 391                     ($control 2)
 392                     ($max_message 4)
 393                     (otherwise
 394                      (merror "Parameter should be current_error, control, or max_message")))
 395                   (or new-value 0)
 396                   (if new-value t nil))))
 397
 398 (defmfun $quad_control (parameter &rest new-value)
 399   (quad-control parameter (if new-value (car new-value))))
 400
 401 ;; Tests
 402 ;;
 403 ;; These tests were taken from the QUADPACK book.
 404 ;;
 405 \f
 406 ;; Test 1
 407 ;; integrate(x^alpha*log(1/x),x,0,1)
 408 ;; => (1+alpha)^(-2)
 409 ;;
 410 ;; alpha = 0.9(0.1)0(0.2)2.6
 411 ;;
 412 ;; QAG with key 1, 3, 6
 413 ;;
 414 ;; For key = 1, 3, 6: fails for alpha = -0.9 (ier = 3)
 415 ;;
 416 ;; quad_qag(x^2*log(1/x),x,0,1,3)
 417 ;;
 418 ;;  for alpha : -0.9 thru 0 step 0.1 do print(alpha, float((1+alpha)^(-2)), quad_qag(x^alpha*log(1/x),x,0,1,3));
 419 ;;  for alpha : 0.0 thru 2.6 step 0.2 do print(alpha, float((1+alpha)^(-2)), quad_qag(x^alpha*log(1/x),x,0,1,3));
 420 \f
 421 ;; Test 2
 422 ;; integrate(4^(-alpha)/((x - %pi/4)^2 + 16^(-alpha)), x, 0, 1)
 423 ;; => atan((4-%pi)*4^(alpha-1)) + atan(%pi*4^(alpha-1))
 424 ;;
 425 ;; alpha = 0(1)20
 426 ;; QAG with key = 1, 3, 6
 427 ;;
 428 ;; Fails for key = 1: alpha >= 18 (ier = 2)
 429 ;; Fails for key = 3, 6: alpha >= 19 (ier = 2)
 430 ;;
 431 ;; for alpha : 0.0 thru 20 step 1 do print(alpha, float(atan((4-%pi)*4^(alpha-1)) + atan(%pi*4^(alpha-1))), quad_qag(4^(-alpha)/((x - %pi/4)^2 + 16^(-alpha)),x,0,1,3));
 432 \f
 433 ;; Test 3
 434 ;; integrate(cos(2^alpha*sin(x)), x, 0, %pi)
 435 ;; => %pi * J0(2^alpha)
 436 ;;
 437 ;; alpha = 0(1)10
 438 ;;
 439 ;; QAG with Key 1, 3, 6
 440 ;;
 441 ;; for alpha : 0.0 thru 10 step 1 do print(alpha, float(%pi * bessel_j(0,2^alpha)), quad_qag(cos(2^alpha*sin(x)),x,0,float(%pi),3));
 442 \f
 443 ;; Test 4 (same integral as 1)
 444 ;; integrate(x^alpha*log(1/x),x,0,1)
 445 ;; => (1+alpha)^(-2)
 446 ;;
 447 ;; DQNG, DQAGS, DQAG (key = 1)
 448 ;;
 449 ;; Failures:
 450 ;; DQNG: alpha <= 1.0 (ier = 1)
 451 ;; DQAG: alpha = -0.9 (ier = 3)
 452 ;;
 453 ;;  for alpha : -0.9 thru 0 step 0.1 do print(alpha, float((1+alpha)^(-2)), quad_qags(x^alpha*log(1/x),x,0,1));
 454 ;;  for alpha : 0.0 thru 2.6 step 0.2 do print(alpha, float((1+alpha)^(-2)), quad_qags(x^alpha*log(1/x),x,0,1));
 455 ;;
 456 ;;  for alpha : -0.9 thru 0 step 0.1 do print(alpha, float((1+alpha)^(-2)), quad_qag(x^alpha*log(1/x),x,0,1, 1));
 457 ;;  for alpha : 0.0 thru 2.6 step 0.2 do print(alpha, float((1+alpha)^(-2)), quad_qag(x^alpha*log(1/x),x,0,1, 1));
 458
 459 \f
 460 ;; Test 5
 461 ;; Same integral as 2
 462 ;;
 463 ;; DQNG, DQAGS, DQAG (key = 1)
 464 ;;
 465 ;; Failures:
 466 ;; DQNG:  alpha >= 2 (ier = 1)
 467 ;; DQAGS: alpha >= 10 (ier = 5)
 468 ;; DQAG:  alpha >= 18 (ier = 2)
 469 ;;
 470 ;; for alpha : 0.0 thru 20 step 1 do print(alpha, float(atan((4-%pi)*4^(alpha-1)) + atan(%pi*4^(alpha-1))), quad_qag(4^(-alpha)/((x - %pi/4)^2 + 16^(-alpha)),x,0,1,1));
 471 ;; for alpha : 0.0 thru 20 step 1 do print(alpha, float(atan((4-%pi)*4^(alpha-1)) + atan(%pi*4^(alpha-1))), quad_qags(4^(-alpha)/((x - %pi/4)^2 + 16^(-alpha)),x,0,1));
 472
 473 \f
 474 ;; Test 6
 475 ;; Same integral as test 3
 476 ;;
 477 ;; DQNG, DQAGS, DQAG (key = 6)
 478 ;;
 479 ;; for alpha : 0.0 thru 10 step 1 do print(alpha, float(%pi * bessel_j(0,2^alpha)), quad_qags(cos(2^alpha*sin(x)),x,0,float(%pi),3));
 480 ;; Failures:
 481 ;; DQNG: alpha >= 7 (ier = 1)
 482 \f
 483 ;; Test 7
 484 ;; integrate(|x - 1/3|^alpha, x, 0, 1)
 485 ;; => ((2/3)^(alpha + 1) + (1/3)^(alpha + 1))/(alpha + 1)
 486 ;;
 487 ;; alpha = -0.8(0.1)2.1
 488 ;; DQAGS, DQAGP (point of singularity supplied)
 489 ;;
 490 ;;
 491 ;; for alpha : -8 thru 21 do print(float(alpha/10), float(((2/3)^(alpha/10 + 1) + (1/3)^(alpha/10 + 1))/(alpha/10 + 1)), quad_qags(abs(x-1/3)^(alpha/10), x, 0, 1));
 492 ;; No failures.
 493 \f
 494 ;; Test 8
 495 ;; integrate(abs(x - pi/4)^alpha, x, 0, 1)
 496 ;; => ((1-pi/4)^(alpha+1) + (pi/4)^(alpha + 1))/(alpha + 1)
 497 ;;
 498 ;; alpha = -0.8(0.1)2.1
 499 ;;
 500 ;; DQAGS, DQAGP
 501 ;;
 502 ;; for alpha : -8 thru 21 do print(float(alpha/10), float(((1-%pi/4)^(alpha/10+1) + (%pi/4)^(alpha/10 + 1))/(alpha/10 + 1)), quad_qags(abs(x-%pi/4)^(alpha/10), x, 0, 1));
 503 ;; Failures:
 504 ;; DQAGS: alpha <= -0.5 (ier = 3)
 505 \f
 506 ;; Test 9
 507 ;;
 508 ;; integrate((1-x*x)^(-1/2)/(x+1+2^(-alpha)),x, -1, 1)
 509 ;; => %pi*((1+2^(-alpha))^2-1)^(-1/2)
 510 ;;
 511 ;; alpha = 1(1)20
 512 ;;
 513 ;; for alpha : 1 thru 20 do print(alpha, float(%pi*((1+2^(-alpha))^2-1)^(-1/2)),quad_qaws(1/(x+1+2^(-alpha)), x, -1, 1, -0.5, -0.5, 1));
 514 \f
 515 ;; Test 10
 516 ;; integrate((sin(x))^(alpha - 1), x, 0, %pi/2) =
 517 ;; integrate(x^(alpha - 1)*(sin(x)/x)^(alpha-1), x, 0, %pi/2)
 518 ;; => 2^(alpha - 2)*(Gamma(alpha/2))^2/Gamma(alpha)
 519 ;;
 520 ;; alpha = 0.1(0.1)2
 521 ;;
 522 ;; DQAGS, DQAWS
 523 ;; Failures:
 524 ;; None.
 525 ;; for alpha : 1 thru 20 do print(alpha/10.0, float(2^(alpha/10 - 2)*(gamma(alpha/20))^2/gamma(alpha/10)),quad_qags((sin(x))^(alpha/10 - 1), x, 0, %pi/2));
 526 ;;
 527 ;; for alpha : 1 thru 20 do print(alpha/10.0, float(2^(alpha/10 - 2)*(gamma(alpha/20))^2/gamma(alpha/10)),quad_qaws(if equal(x,0.0) then 1 else (sin(x)/x)^(alpha/10 - 1), x, 0, %pi/2, alpha/10-1, 0, 1));
 528 \f
 529 ;; Test 11
 530 ;; integrate((log(1/x))^(alpha-1), x, 0, 1) =
 531 ;; integrate((1-x)^(alpha - 1)*(log(1/x)/(1-x))^(alpha-1), x, 0, 1)
 532 ;; => Gamma(alpha)
 533 ;;
 534 ;; alpha = 0.1(0.1)2
 535 ;;
 536 ;; DQAGS, DQAWS
 537 ;; for alpha : 1 thru 20 do print(alpha/10.0, float(gamma(alpha/10)),quad_qags(log(1/x)^(alpha/10-1),x,0,1));
 538 ;;
 539 ;; for alpha : 1 thru 20 do print(alpha/10.0, float(gamma(alpha/10)),quad_qaws(if equal(x,1) then 1 else (log(1/x)/(1-x))^(alpha/10-1),x,0,1,0,alpha/10-1,1));
 540
 541 \f
 542 ;; Test 12
 543 ;; integrate(exp(20*(x-1))*sin(2^alpha*x), x, 0, 1)
 544 ;; => (20*sin(2^alpha) - 2^alpha*cos(2^alpha) + 2^alpha*exp(-20))/(400 + 4^alpha)
 545 ;;
 546 ;; alpha = 0(1)9
 547 ;;
 548 ;; DQAG (key = 6), DQAWO
 549 ;; Failures:
 550 ;; None
 551 ;;
 552 ;; for alpha : 0 thru 9 do print(alpha, float((20*sin(2^alpha) - 2^alpha*cos(2^alpha) + 2^alpha*exp(-20))/(400 + 4^alpha)), quad_qag(exp(20*(x-1))*sin(2^alpha*x), x, 0, 1, 6));
 553 ;;
 554 ;; for alpha : 0 thru 9 do print(alpha, float((20*sin(2^alpha) - 2^alpha*cos(2^alpha) + 2^alpha*exp(-20))/(400 + 4^alpha)), quad_qawo(exp(20*(x-1)), x, 0, 1, 2^alpha, 'sin));
 555 ;;
 556 \f
 557 ;; Test 13
 558 ;; integrate((x*(1-x))^(-1/2)*cos(2^alpha*x), x, 0, 1)
 559 ;; => cos(2^(alpha-1))*J0(2^(alpha - 1))
 560 ;;
 561 ;; alpha = 0(1)8
 562 ;;
 563 ;; DQAGS, DQAWO, DQAWS
 564 ;;
 565 ;; Failures:
 566 ;; DQAGS: alpha = 4 (ier = 5)
 567 ;;
 568 ;; for alpha : 0 thru 8 do print(alpha, float(cos(2^(alpha-1))*bessel_j(0,2^(alpha - 1))), quad_qags((x*(1-x))^(-1/2)*cos(2^alpha*x), x, 0, 1));
 569 ;;
 570 ;; This doesn't work:
 571 ;; for alpha : 0 thru 8 do print(alpha, float(cos(2^(alpha-1))*bessel_j(0,2^(alpha - 1))), quad_qawo(if equal(x,0) or equal(x,1) then 0 else (x*(1-x))^(-1/2), x, 0, 1, 2^alpha, 'cos));
 572 ;;
 573 ;; for alpha : 0 thru 8 do print(alpha, float(cos(2^(alpha-1))*bessel_j(0,2^(alpha - 1))), quad_qaws(cos(2^alpha*x), x, 0, 1, -1/2, -1/2, 1));
 574 \f
 575 ;; Test 14
 576 ;;
 577 ;; integrate(x^(-1/2)*exp(-2^(-alpha)*x) * cos(x), x, 0, inf);
 578 ;; => sqrt(%pi)*(1-4^(-alpha))^(-1/4)*cos(atan(2^alpha)/2)
 579 ;;
 580 ;; quad_qawf(x^(-1/2)*exp(-2^(-2)*x), x, 0, 1, cos)
 581 ;; quad_qawf(x^(-1/2)*exp(-2^(-2)*x), x, 1e-8, 1, cos)
 582 ;; quad_qawo(x^(-1/2)*exp(-2^(-2)*x), x, 1e-8, 20*2^2, 1, cos)
 583 ;;
 584 \f
 585 ;; Test 15
 586 ;;
 587 ;; integrate(x^2*exp(-2^(-alpha)*x), x, 0, inf)
 588 ;; => 2^(3*alpha + 1)
 589 ;;
 590 ;; alpha = 0(1)5
 591 ;;
 592 ;; quad_qagi(x^2*exp(-2^(-alpha)*x), x, 0, inf)
 593 ;;
 594 ;; for alpha : 0 thru 5 do print(alpha, float(2^(3*alpha+1)), quad_qagi(x^2*exp(-2^(-alpha)*x), x, 0, inf));
 595 ;;
 596 \f
 597 ;; Test 16
 598 ;; integrate(x^(alpha - 1)/(1+10*x)^2, 0, inf)
 599 ;; => 10^(-alpha)*(1-alpha)*pi/sin(pi*alpha)
 600 ;; if alpha /= 1.  Otherwise result = 1/10 when alpha = 1.
 601 ;;
 602 ;; alpha = 0.1(0.1)1.9
 603 ;;
 604 ;; DQAGI
 605 ;;
 606 ;; Failures: None
 607 ;;
 608 ;; for alpha : 1 thru 19 do print(alpha/10.0, float(10^(-alpha/10)*(1-alpha/10)*%pi/sin(%pi*alpha/10)), quad_qagi(x^(alpha/10 - 1)/(1+10*x)^2, x, 0, inf));
 609 \f
 610 ;; Test 17
 611 ;;
 612 ;; Cauchy principal value of
 613 ;;
 614 ;; integrate(2^(-alpha)*(((x-1)^2 + 4^(-alpha))*(x-2))^(-1), x, 0, 5);
 615 ;; => (2^(-alpha)*ln(3/2) - 2^(-alpha-1)*ln((16 + 4^(-alpha))/(1+4^(-alpha))) - atan(2^(alpha +  2)) - atan(2^alpha))/(1 + 4^(-alpha))
 616 ;;
 617 ;; alpha = 0(1)10
 618 ;;
 619 ;; quad_qawc(2^(-alpha)*((x-1)^2 + 4^(-alpha))^(-1), 2, 0, 5)
 620 ;;
 621 ;; for alpha : 0 thru 10 do print(alpha, float((2^(-alpha)*log(3/2) - 2^(-alpha-1)*log((16 + 4^(-alpha))/(1+4^(-alpha))) - atan(2^(alpha +  2)) - atan(2^alpha))/(1 + 4^(-alpha))), quad_qawc(2^(-alpha)*((x-1)^2 + 4^(-alpha))^(-1), x, 2, 0, 5));