src/numerical/slatec/dqc25s.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.9238795325112868 0.8660254037844386
  25                                        0.7933533402912352 0.7071067811865475
  26                                        0.6087614290087205 0.5
  27                                        0.3826834323650898 0.2588190451025208
  28                                        0.1305261922200516))))
  29   (declare (type (array double-float (11)) x))
  30   (defun dqc25s
  31          (f a b bl br alfa beta ri rj rg rh result abserr resasc integr nev)
  32     (declare (type (f2cl-lib:integer4) nev integr)
  33              (type (array double-float (*)) rh rg rj ri)
  34              (type (double-float) resasc abserr result beta alfa br bl b a))
  35     (f2cl-lib:with-multi-array-data
  36         ((ri double-float ri-%data% ri-%offset%)
  37          (rj double-float rj-%data% rj-%offset%)
  38          (rg double-float rg-%data% rg-%offset%)
  39          (rh double-float rh-%data% rh-%offset%))
  40       (prog ((cheb12 (make-array 13 :element-type 'double-float))
  41              (cheb24 (make-array 25 :element-type 'double-float))
  42              (fval (make-array 25 :element-type 'double-float)) (i 0) (isym 0)
  43              (centr 0.0) (dc 0.0) (factor 0.0) (fix 0.0) (hlgth 0.0)
  44              (resabs 0.0) (res12 0.0) (res24 0.0) (u 0.0))
  45         (declare (type (array double-float (25)) fval cheb24)
  46                  (type (array double-float (13)) cheb12)
  47                  (type (double-float) u res24 res12 resabs hlgth fix factor dc
  48                                       centr)
  49                  (type (f2cl-lib:integer4) isym i))
  50         (setf nev 25)
  51         (if (and (= bl a) (or (/= alfa 0.0) (= integr 2) (= integr 4)))
  52             (go label10))
  53         (if (and (= br b) (or (/= beta 0.0) (= integr 3) (= integr 4)))
  54             (go label140))
  55         (multiple-value-bind
  56               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
  57                var-10 var-11 var-12)
  58             (dqk15w f #'dqwgts a b alfa beta integr bl br result abserr resabs
  59              resasc)
  60           (declare (ignore var-0 var-1 var-7 var-8))
  61           (setf a var-2)
  62           (setf b var-3)
  63           (setf alfa var-4)
  64           (setf beta var-5)
  65           (setf integr var-6)
  66           (setf result var-9)
  67           (setf abserr var-10)
  68           (setf resabs var-11)
  69           (setf resasc var-12))
  70         (setf nev 15)
  71         (go label270)
  72        label10
  73         (setf hlgth (* 0.5 (- br bl)))
  74         (setf centr (* 0.5 (+ br bl)))
  75         (setf fix (- b centr))
  76         (setf (f2cl-lib:fref fval (1) ((1 25)))
  77                 (* 0.5 (funcall f (+ hlgth centr)) (expt (- fix hlgth) beta)))
  78         (setf (f2cl-lib:fref fval (13) ((1 25)))
  79                 (*
  80                  (multiple-value-bind (ret-val var-0)
  81                      (funcall f centr)
  82                    (declare (ignore))
  83                    (when var-0
  84                      (setf centr var-0))
  85                    ret-val)
  86                  (expt fix beta)))
  87         (setf (f2cl-lib:fref fval (25) ((1 25)))
  88                 (* 0.5 (funcall f (- centr hlgth)) (expt (+ fix hlgth) beta)))
  89         (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
  90                       ((> i 12) nil)
  91           (tagbody
  92             (setf u
  93                     (* hlgth
  94                        (f2cl-lib:fref x ((f2cl-lib:int-sub i 1)) ((1 11)))))
  95             (setf isym (f2cl-lib:int-sub 26 i))
  96             (setf (f2cl-lib:fref fval (i) ((1 25)))
  97                     (* (funcall f (+ u centr)) (expt (- fix u) beta)))
  98             (setf (f2cl-lib:fref fval (isym) ((1 25)))
  99                     (* (funcall f (- centr u)) (expt (+ fix u) beta)))
 100            label20))
 101         (setf factor (expt hlgth (+ alfa 1.0)))
 102         (setf result 0.0)
 103         (setf abserr 0.0)
 104         (setf res12 0.0)
 105         (setf res24 0.0)
 106         (if (> integr 2) (go label70))
 107         (dqcheb x fval cheb12 cheb24)
 108         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 109                       ((> i 13) nil)
 110           (tagbody
 111             (setf res12
 112                     (+ res12
 113                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 114                           (f2cl-lib:fref ri-%data% (i) ((1 25)) ri-%offset%))))
 115             (setf res24
 116                     (+ res24
 117                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 118                           (f2cl-lib:fref ri-%data% (i) ((1 25)) ri-%offset%))))
 119            label30))
 120         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 121                       ((> i 25) nil)
 122           (tagbody
 123             (setf res24
 124                     (+ res24
 125                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 126                           (f2cl-lib:fref ri-%data% (i) ((1 25)) ri-%offset%))))
 127            label40))
 128         (if (= integr 1) (go label130))
 129         (setf dc (f2cl-lib:flog (- br bl)))
 130         (setf result (* res24 dc))
 131         (setf abserr (abs (* (- res24 res12) dc)))
 132         (setf res12 0.0)
 133         (setf res24 0.0)
 134         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 135                       ((> i 13) nil)
 136           (tagbody
 137             (setf res12
 138                     (+ res12
 139                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 140                           (f2cl-lib:fref rg-%data% (i) ((1 25)) rg-%offset%))))
 141             (setf res24
 142                     (+ res24
 143                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 144                           (f2cl-lib:fref rg-%data% (i) ((1 25)) rg-%offset%))))
 145            label50))
 146         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 147                       ((> i 25) nil)
 148           (tagbody
 149             (setf res24
 150                     (+ res24
 151                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 152                           (f2cl-lib:fref rg-%data% (i) ((1 25)) rg-%offset%))))
 153            label60))
 154         (go label130)
 155        label70
 156         (setf (f2cl-lib:fref fval (1) ((1 25)))
 157                 (* (f2cl-lib:fref fval (1) ((1 25)))
 158                    (f2cl-lib:flog (- fix hlgth))))
 159         (setf (f2cl-lib:fref fval (13) ((1 25)))
 160                 (* (f2cl-lib:fref fval (13) ((1 25))) (f2cl-lib:flog fix)))
 161         (setf (f2cl-lib:fref fval (25) ((1 25)))
 162                 (* (f2cl-lib:fref fval (25) ((1 25)))
 163                    (f2cl-lib:flog (+ fix hlgth))))
 164         (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
 165                       ((> i 12) nil)
 166           (tagbody
 167             (setf u
 168                     (* hlgth
 169                        (f2cl-lib:fref x ((f2cl-lib:int-sub i 1)) ((1 11)))))
 170             (setf isym (f2cl-lib:int-sub 26 i))
 171             (setf (f2cl-lib:fref fval (i) ((1 25)))
 172                     (* (f2cl-lib:fref fval (i) ((1 25)))
 173                        (f2cl-lib:flog (- fix u))))
 174             (setf (f2cl-lib:fref fval (isym) ((1 25)))
 175                     (* (f2cl-lib:fref fval (isym) ((1 25)))
 176                        (f2cl-lib:flog (+ fix u))))
 177            label80))
 178         (dqcheb x fval cheb12 cheb24)
 179         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 180                       ((> i 13) nil)
 181           (tagbody
 182             (setf res12
 183                     (+ res12
 184                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 185                           (f2cl-lib:fref ri-%data% (i) ((1 25)) ri-%offset%))))
 186             (setf res24
 187                     (+ res24
 188                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 189                           (f2cl-lib:fref ri-%data% (i) ((1 25)) ri-%offset%))))
 190            label90))
 191         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 192                       ((> i 25) nil)
 193           (tagbody
 194             (setf res24
 195                     (+ res24
 196                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 197                           (f2cl-lib:fref ri-%data% (i) ((1 25)) ri-%offset%))))
 198            label100))
 199         (if (= integr 3) (go label130))
 200         (setf dc (f2cl-lib:flog (- br bl)))
 201         (setf result (* res24 dc))
 202         (setf abserr (abs (* (- res24 res12) dc)))
 203         (setf res12 0.0)
 204         (setf res24 0.0)
 205         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 206                       ((> i 13) nil)
 207           (tagbody
 208             (setf res12
 209                     (+ res12
 210                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 211                           (f2cl-lib:fref rg-%data% (i) ((1 25)) rg-%offset%))))
 212             (setf res24
 213                     (+ res24
 214                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 215                           (f2cl-lib:fref rg-%data% (i) ((1 25)) rg-%offset%))))
 216            label110))
 217         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 218                       ((> i 25) nil)
 219           (tagbody
 220             (setf res24
 221                     (+ res24
 222                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 223                           (f2cl-lib:fref rg-%data% (i) ((1 25)) rg-%offset%))))
 224            label120))
 225        label130
 226         (setf result (* (+ result res24) factor))
 227         (setf abserr (* (+ abserr (abs (- res24 res12))) factor))
 228         (go label270)
 229        label140
 230         (setf hlgth (* 0.5 (- br bl)))
 231         (setf centr (* 0.5 (+ br bl)))
 232         (setf fix (- centr a))
 233         (setf (f2cl-lib:fref fval (1) ((1 25)))
 234                 (* 0.5 (funcall f (+ hlgth centr)) (expt (+ fix hlgth) alfa)))
 235         (setf (f2cl-lib:fref fval (13) ((1 25)))
 236                 (*
 237                  (multiple-value-bind (ret-val var-0)
 238                      (funcall f centr)
 239                    (declare (ignore))
 240                    (when var-0
 241                      (setf centr var-0))
 242                    ret-val)
 243                  (expt fix alfa)))
 244         (setf (f2cl-lib:fref fval (25) ((1 25)))
 245                 (* 0.5 (funcall f (- centr hlgth)) (expt (- fix hlgth) alfa)))
 246         (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
 247                       ((> i 12) nil)
 248           (tagbody
 249             (setf u
 250                     (* hlgth
 251                        (f2cl-lib:fref x ((f2cl-lib:int-sub i 1)) ((1 11)))))
 252             (setf isym (f2cl-lib:int-sub 26 i))
 253             (setf (f2cl-lib:fref fval (i) ((1 25)))
 254                     (* (funcall f (+ u centr)) (expt (+ fix u) alfa)))
 255             (setf (f2cl-lib:fref fval (isym) ((1 25)))
 256                     (* (funcall f (- centr u)) (expt (- fix u) alfa)))
 257            label150))
 258         (setf factor (expt hlgth (+ beta 1.0)))
 259         (setf result 0.0)
 260         (setf abserr 0.0)
 261         (setf res12 0.0)
 262         (setf res24 0.0)
 263         (if (or (= integr 2) (= integr 4)) (go label200))
 264         (dqcheb x fval cheb12 cheb24)
 265         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 266                       ((> i 13) nil)
 267           (tagbody
 268             (setf res12
 269                     (+ res12
 270                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 271                           (f2cl-lib:fref rj-%data% (i) ((1 25)) rj-%offset%))))
 272             (setf res24
 273                     (+ res24
 274                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 275                           (f2cl-lib:fref rj-%data% (i) ((1 25)) rj-%offset%))))
 276            label160))
 277         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 278                       ((> i 25) nil)
 279           (tagbody
 280             (setf res24
 281                     (+ res24
 282                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 283                           (f2cl-lib:fref rj-%data% (i) ((1 25)) rj-%offset%))))
 284            label170))
 285         (if (= integr 1) (go label260))
 286         (setf dc (f2cl-lib:flog (- br bl)))
 287         (setf result (* res24 dc))
 288         (setf abserr (abs (* (- res24 res12) dc)))
 289         (setf res12 0.0)
 290         (setf res24 0.0)
 291         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 292                       ((> i 13) nil)
 293           (tagbody
 294             (setf res12
 295                     (+ res12
 296                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 297                           (f2cl-lib:fref rh-%data% (i) ((1 25)) rh-%offset%))))
 298             (setf res24
 299                     (+ res24
 300                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 301                           (f2cl-lib:fref rh-%data% (i) ((1 25)) rh-%offset%))))
 302            label180))
 303         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 304                       ((> i 25) nil)
 305           (tagbody
 306             (setf res24
 307                     (+ res24
 308                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 309                           (f2cl-lib:fref rh-%data% (i) ((1 25)) rh-%offset%))))
 310            label190))
 311         (go label260)
 312        label200
 313         (setf (f2cl-lib:fref fval (1) ((1 25)))
 314                 (* (f2cl-lib:fref fval (1) ((1 25)))
 315                    (f2cl-lib:flog (+ fix hlgth))))
 316         (setf (f2cl-lib:fref fval (13) ((1 25)))
 317                 (* (f2cl-lib:fref fval (13) ((1 25))) (f2cl-lib:flog fix)))
 318         (setf (f2cl-lib:fref fval (25) ((1 25)))
 319                 (* (f2cl-lib:fref fval (25) ((1 25)))
 320                    (f2cl-lib:flog (- fix hlgth))))
 321         (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
 322                       ((> i 12) nil)
 323           (tagbody
 324             (setf u
 325                     (* hlgth
 326                        (f2cl-lib:fref x ((f2cl-lib:int-sub i 1)) ((1 11)))))
 327             (setf isym (f2cl-lib:int-sub 26 i))
 328             (setf (f2cl-lib:fref fval (i) ((1 25)))
 329                     (* (f2cl-lib:fref fval (i) ((1 25)))
 330                        (f2cl-lib:flog (+ u fix))))
 331             (setf (f2cl-lib:fref fval (isym) ((1 25)))
 332                     (* (f2cl-lib:fref fval (isym) ((1 25)))
 333                        (f2cl-lib:flog (- fix u))))
 334            label210))
 335         (dqcheb x fval cheb12 cheb24)
 336         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 337                       ((> i 13) nil)
 338           (tagbody
 339             (setf res12
 340                     (+ res12
 341                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 342                           (f2cl-lib:fref rj-%data% (i) ((1 25)) rj-%offset%))))
 343             (setf res24
 344                     (+ res24
 345                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 346                           (f2cl-lib:fref rj-%data% (i) ((1 25)) rj-%offset%))))
 347            label220))
 348         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 349                       ((> i 25) nil)
 350           (tagbody
 351             (setf res24
 352                     (+ res24
 353                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 354                           (f2cl-lib:fref rj-%data% (i) ((1 25)) rj-%offset%))))
 355            label230))
 356         (if (= integr 2) (go label260))
 357         (setf dc (f2cl-lib:flog (- br bl)))
 358         (setf result (* res24 dc))
 359         (setf abserr (abs (* (- res24 res12) dc)))
 360         (setf res12 0.0)
 361         (setf res24 0.0)
 362         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 363                       ((> i 13) nil)
 364           (tagbody
 365             (setf res12
 366                     (+ res12
 367                        (* (f2cl-lib:fref cheb12 (i) ((1 13)))
 368                           (f2cl-lib:fref rh-%data% (i) ((1 25)) rh-%offset%))))
 369             (setf res24
 370                     (+ res24
 371                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 372                           (f2cl-lib:fref rh-%data% (i) ((1 25)) rh-%offset%))))
 373            label240))
 374         (f2cl-lib:fdo (i 14 (f2cl-lib:int-add i 1))
 375                       ((> i 25) nil)
 376           (tagbody
 377             (setf res24
 378                     (+ res24
 379                        (* (f2cl-lib:fref cheb24 (i) ((1 25)))
 380                           (f2cl-lib:fref rh-%data% (i) ((1 25)) rh-%offset%))))
 381            label250))
 382        label260
 383         (setf result (* (+ result res24) factor))
 384         (setf abserr (* (+ abserr (abs (- res24 res12))) factor))
 385        label270
 386         (go end_label)
 387        end_label
 388         (return
 389          (values nil
 390                  a
 391                  b
 392                  nil
 393                  nil
 394                  alfa
 395                  beta
 396                  nil
 397                  nil
 398                  nil
 399                  nil
 400                  result
 401                  abserr
 402                  resasc
 403                  integr
 404                  nev))))))
 405
 406 (in-package #:cl-user)
 407 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 408 (eval-when (:load-toplevel :compile-toplevel :execute)
 409   (setf (gethash 'fortran-to-lisp::dqc25s
 410                  fortran-to-lisp::*f2cl-function-info*)
 411           (fortran-to-lisp::make-f2cl-finfo
 412            :arg-types '(t (double-float) (double-float) (double-float)
 413                         (double-float) (double-float) (double-float)
 414                         (array double-float (*)) (array double-float (*))
 415                         (array double-float (*)) (array double-float (*))
 416                         (double-float) (double-float) (double-float)
 417                         (fortran-to-lisp::integer4)
 418                         (fortran-to-lisp::integer4))
 419            :return-values '(nil fortran-to-lisp::a fortran-to-lisp::b nil nil
 420                             fortran-to-lisp::alfa fortran-to-lisp::beta nil nil
 421                             nil nil fortran-to-lisp::result
 422                             fortran-to-lisp::abserr fortran-to-lisp::resasc
 423                             fortran-to-lisp::integr fortran-to-lisp::nev)
 424            :calls '(fortran-to-lisp::dqcheb fortran-to-lisp::dqk15w))))
 425