src/numerical/slatec/dqc25c.lisp

   1 ;;; Compiled by f2cl version:
   2 ;;; ("f2cl1.l,v 46c1f6a93b0d 2012/05/03 04:40:28 toy $"
   3 ;;;  "f2cl2.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
   4 ;;;  "f2cl3.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
   5 ;;;  "f2cl4.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
   6 ;;;  "f2cl5.l,v 46c1f6a93b0d 2012/05/03 04:40:28 toy $"
   7 ;;;  "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $"
   8 ;;;  "macros.l,v fceac530ef0c 2011/11/26 04:02:26 toy $")
   9
  10 ;;; Using Lisp CMU Common Lisp snapshot-2012-04 (20C Unicode)
  11 ;;;
  12 ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  13 ;;;           (:coerce-assigns :as-needed) (:array-type ':array)
  14 ;;;           (:array-slicing t) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :slatec)
  18
  19
  20 (let ((x
  21        (make-array 11
  22                    :element-type 'double-float
  23                    :initial-contents '(0.9914448613738104 0.9659258262890683
  24                                        0.9238795325112867 0.8660254037844386
  25                                        0.7933533402912352 0.7071067811865476
  26                                        0.6087614290087207 0.5
  27                                        0.3826834323650898 0.25881904510252074
  28                                        0.1305261922200516))))
  29   (declare (type (array double-float (11)) x))
  30   (defun dqc25c (f a b c result abserr krul neval)
  31     (declare (type (f2cl-lib:integer4) neval krul)
  32              (type (double-float) abserr result c b a))
  33     (prog ((fval (make-array 25 :element-type 'double-float))
  34            (cheb12 (make-array 13 :element-type 'double-float))
  35            (cheb24 (make-array 25 :element-type 'double-float)) (i 0) (isym 0)
  36            (k 0) (kp 0) (ak22 0.0) (amom0 0.0) (amom1 0.0) (amom2 0.0) (cc 0.0)
  37            (centr 0.0) (hlgth 0.0) (p2 0.0) (p3 0.0) (p4 0.0) (resabs 0.0)
  38            (resasc 0.0) (res12 0.0) (res24 0.0) (u 0.0))
  39       (declare (type (array double-float (25)) fval cheb24)
  40                (type (array double-float (13)) cheb12)
  41                (type (double-float) u res24 res12 resasc resabs p4 p3 p2 hlgth
  42                                     centr cc amom2 amom1 amom0 ak22)
  43                (type (f2cl-lib:integer4) kp k isym i))
  44       (setf cc (/ (- (* 2.0 c) b a) (- b a)))
  45       (if (< (abs cc) 1.1) (go label10))
  46       (setf krul (f2cl-lib:int-sub krul 1))
  47       (multiple-value-bind
  48             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  49              var-11 var-12)
  50           (dqk15w f #'dqwgtc c p2 p3 p4 kp a b result abserr resabs resasc)
  51         (declare (ignore var-0 var-1 var-7 var-8))
  52         (setf c var-2)
  53         (setf p2 var-3)
  54         (setf p3 var-4)
  55         (setf p4 var-5)
  56         (setf kp var-6)
  57         (setf result var-9)
  58         (setf abserr var-10)
  59         (setf resabs var-11)
  60         (setf resasc var-12))
  61       (setf neval 15)
  62       (if (= resasc abserr) (setf krul (f2cl-lib:int-add krul 1)))
  63       (go label50)
  64      label10
  65       (setf hlgth (* 0.5 (- b a)))
  66       (setf centr (* 0.5 (+ b a)))
  67       (setf neval 25)
  68       (setf (f2cl-lib:fref fval (1) ((1 25)))
  69               (* 0.5 (funcall f (+ hlgth centr))))
  70       (setf (f2cl-lib:fref fval (13) ((1 25)))
  71               (multiple-value-bind (ret-val var-0)
  72                   (funcall f centr)
  73                 (declare (ignore))
  74                 (when var-0
  75                   (setf centr var-0))
  76                 ret-val))
  77       (setf (f2cl-lib:fref fval (25) ((1 25)))
  78               (* 0.5 (funcall f (- centr hlgth))))
  79       (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
  80                     ((> i 12) nil)
  81         (tagbody
  82           (setf u
  83                   (* hlgth
  84                      (f2cl-lib:fref x ((f2cl-lib:int-sub i 1)) ((1 11)))))
  85           (setf isym (f2cl-lib:int-sub 26 i))
  86           (setf (f2cl-lib:fref fval (i) ((1 25))) (funcall f (+ u centr)))
  87           (setf (f2cl-lib:fref fval (isym) ((1 25))) (funcall f (- centr u)))
  88          label20))
  89       (dqcheb x fval cheb12 cheb24)
  90       (setf amom0 (f2cl-lib:flog (abs (/ (- 1.0 cc) (+ 1.0 cc)))))
  91       (setf amom1 (+ 2.0 (* cc amom0)))
  92       (setf res12
  93               (+ (* (f2cl-lib:fref cheb12 (1) ((1 13))) amom0)
  94                  (* (f2cl-lib:fref cheb12 (2) ((1 13))) amom1)))
  95       (setf res24
  96               (+ (* (f2cl-lib:fref cheb24 (1) ((1 25))) amom0)
  97                  (* (f2cl-lib:fref cheb24 (2) ((1 25))) amom1)))
  98       (f2cl-lib:fdo (k 3 (f2cl-lib:int-add k 1))
  99                     ((> k 13) nil)
 100         (tagbody
 101           (setf amom2 (- (* 2.0 cc amom1) amom0))
 102           (setf ak22
 103                   (coerce
 104                    (the f2cl-lib:integer4
 105                         (f2cl-lib:int-mul (f2cl-lib:int-sub k 2)
 106                                           (f2cl-lib:int-sub k 2)))
 107                    'double-float))
 108           (if (= (* (the f2cl-lib:integer4 (truncate k 2)) 2) k)
 109               (setf amom2 (+ amom2 (/ -4.0 (- ak22 1.0)))))
 110           (setf res12 (+ res12 (* (f2cl-lib:fref cheb12 (k) ((1 13))) amom2)))
 111           (setf res24 (+ res24 (* (f2cl-lib:fref cheb24 (k) ((1 25))) amom2)))
 112           (setf amom0 amom1)
 113           (setf amom1 amom2)
 114          label30))
 115       (f2cl-lib:fdo (k 14 (f2cl-lib:int-add k 1))
 116                     ((> k 25) nil)
 117         (tagbody
 118           (setf amom2 (- (* 2.0 cc amom1) amom0))
 119           (setf ak22
 120                   (coerce
 121                    (the f2cl-lib:integer4
 122                         (f2cl-lib:int-mul (f2cl-lib:int-sub k 2)
 123                                           (f2cl-lib:int-sub k 2)))
 124                    'double-float))
 125           (if (= (* (the f2cl-lib:integer4 (truncate k 2)) 2) k)
 126               (setf amom2 (+ amom2 (/ -4.0 (- ak22 1.0)))))
 127           (setf res24 (+ res24 (* (f2cl-lib:fref cheb24 (k) ((1 25))) amom2)))
 128           (setf amom0 amom1)
 129           (setf amom1 amom2)
 130          label40))
 131       (setf result res24)
 132       (setf abserr (abs (- res24 res12)))
 133      label50
 134       (go end_label)
 135      end_label
 136       (return (values nil nil nil c result abserr krul neval)))))
 137
 138 (in-package #:cl-user)
 139 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 140 (eval-when (:load-toplevel :compile-toplevel :execute)
 141   (setf (gethash 'fortran-to-lisp::dqc25c
 142                  fortran-to-lisp::*f2cl-function-info*)
 143           (fortran-to-lisp::make-f2cl-finfo
 144            :arg-types '(t (double-float) (double-float) (double-float)
 145                         (double-float) (double-float)
 146                         (fortran-to-lisp::integer4)
 147                         (fortran-to-lisp::integer4))
 148            :return-values '(nil nil nil fortran-to-lisp::c
 149                             fortran-to-lisp::result fortran-to-lisp::abserr
 150                             fortran-to-lisp::krul fortran-to-lisp::neval)
 151            :calls '(fortran-to-lisp::dqcheb fortran-to-lisp::dqk15w))))
 152