src/numerical/slatec/dqc25f.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 dqc25f
  31          (f a b omega integr nrmom maxp1 ksave result abserr neval resabs
  32           resasc momcom chebmo)
  33     (declare (type (array double-float (*)) chebmo)
  34              (type (f2cl-lib:integer4) momcom neval ksave maxp1 nrmom integr)
  35              (type (double-float) resasc resabs abserr result omega b a))
  36     (f2cl-lib:with-multi-array-data
  37         ((chebmo double-float chebmo-%data% chebmo-%offset%))
  38       (prog ((cheb12 (make-array 13 :element-type 'double-float))
  39              (cheb24 (make-array 25 :element-type 'double-float))
  40              (d (make-array 25 :element-type 'double-float))
  41              (d1 (make-array 25 :element-type 'double-float))
  42              (d2 (make-array 25 :element-type 'double-float))
  43              (fval (make-array 25 :element-type 'double-float))
  44              (v (make-array 28 :element-type 'double-float)) (i 0) (iers 0)
  45              (isym 0) (j 0) (k 0) (m 0) (noequ 0) (noeq1 0) (ac 0.0) (an 0.0)
  46              (an2 0.0) (as 0.0) (asap 0.0) (ass 0.0) (centr 0.0) (conc 0.0)
  47              (cons$ 0.0) (cospar 0.0) (estc 0.0) (ests 0.0) (hlgth 0.0)
  48              (oflow 0.0) (parint 0.0) (par2 0.0) (par22 0.0) (p2 0.0) (p3 0.0)
  49              (p4 0.0) (resc12 0.0) (resc24 0.0) (ress12 0.0) (ress24 0.0)
  50              (sinpar 0.0))
  51         (declare (type (array double-float (28)) v)
  52                  (type (array double-float (25)) fval d2 d1 d cheb24)
  53                  (type (array double-float (13)) cheb12)
  54                  (type (double-float) sinpar ress24 ress12 resc24 resc12 p4 p3
  55                                       p2 par22 par2 parint oflow hlgth ests
  56                                       estc cospar cons$ conc centr ass asap as
  57                                       an2 an ac)
  58                  (type (f2cl-lib:integer4) noeq1 noequ m k j isym iers i))
  59         (setf oflow (f2cl-lib:d1mach 2))
  60         (setf centr (* 0.5 (+ b a)))
  61         (setf hlgth (* 0.5 (- b a)))
  62         (setf parint (* omega hlgth))
  63         (if (> (abs parint) 2.0) (go label10))
  64         (multiple-value-bind
  65               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
  66                var-10 var-11 var-12)
  67             (dqk15w f #'dqwgtf omega p2 p3 p4 integr a b result abserr resabs
  68              resasc)
  69           (declare (ignore var-0 var-1 var-7 var-8))
  70           (setf omega var-2)
  71           (setf p2 var-3)
  72           (setf p3 var-4)
  73           (setf p4 var-5)
  74           (setf integr var-6)
  75           (setf result var-9)
  76           (setf abserr var-10)
  77           (setf resabs var-11)
  78           (setf resasc var-12))
  79         (setf neval 15)
  80         (go label170)
  81        label10
  82         (setf conc (* hlgth (cos (* centr omega))))
  83         (setf cons$ (* hlgth (sin (* centr omega))))
  84         (setf resasc oflow)
  85         (setf neval 25)
  86         (if (or (< nrmom momcom) (= ksave 1)) (go label120))
  87         (setf m (f2cl-lib:int-add momcom 1))
  88         (setf par2 (* parint parint))
  89         (setf par22 (+ par2 2.0))
  90         (setf sinpar (sin parint))
  91         (setf cospar (cos parint))
  92         (setf (f2cl-lib:fref v (1) ((1 28))) (/ (* 2.0 sinpar) parint))
  93         (setf (f2cl-lib:fref v (2) ((1 28)))
  94                 (/
  95                  (+ (* 8.0 cospar) (/ (* (- (+ par2 par2) 8.0) sinpar) parint))
  96                  par2))
  97         (setf (f2cl-lib:fref v (3) ((1 28)))
  98                 (/
  99                  (+ (* 32.0 (- par2 12.0) cospar)
 100                     (/ (* 2.0 (+ (* (- par2 80.0) par2) 192.0) sinpar) parint))
 101                  (* par2 par2)))
 102         (setf ac (* 8.0 cospar))
 103         (setf as (* 24.0 parint sinpar))
 104         (if (> (abs parint) 24.0) (go label30))
 105         (setf noequ 25)
 106         (setf noeq1 (f2cl-lib:int-sub noequ 1))
 107         (setf an 6.0)
 108         (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 109                       ((> k noeq1) nil)
 110           (tagbody
 111             (setf an2 (* an an))
 112             (setf (f2cl-lib:fref d (k) ((1 25)))
 113                     (* -2.0 (- an2 4.0) (- par22 an2 an2)))
 114             (setf (f2cl-lib:fref d2 (k) ((1 25)))
 115                     (* (- an 1.0) (- an 2.0) par2))
 116             (setf (f2cl-lib:fref d1 ((f2cl-lib:int-add k 1)) ((1 25)))
 117                     (* (+ an 3.0) (+ an 4.0) par2))
 118             (setf (f2cl-lib:fref v ((f2cl-lib:int-add k 3)) ((1 28)))
 119                     (- as (* (- an2 4.0) ac)))
 120             (setf an (+ an 2.0))
 121            label20))
 122         (setf an2 (* an an))
 123         (setf (f2cl-lib:fref d (noequ) ((1 25)))
 124                 (* -2.0 (- an2 4.0) (- par22 an2 an2)))
 125         (setf (f2cl-lib:fref v ((f2cl-lib:int-add noequ 3)) ((1 28)))
 126                 (- as (* (- an2 4.0) ac)))
 127         (setf (f2cl-lib:fref v (4) ((1 28)))
 128                 (+ (f2cl-lib:fref v (4) ((1 28)))
 129                    (* -56.0 par2 (f2cl-lib:fref v (3) ((1 28))))))
 130         (setf ass (* parint sinpar))
 131         (setf asap
 132                 (/
 133                  (-
 134                   (/
 135                    (+
 136                     (-
 137                      (/
 138                       (+
 139                        (-
 140                         (/
 141                          (- (* (- (* 210.0 par2) 1.0) cospar)
 142                             (* (- (* 105.0 par2) 63.0) ass))
 143                          an2)
 144                         (* (+ 1.0 (* -1 15.0 par2)) cospar))
 145                        (* 15.0 ass))
 146                       an2)
 147                      cospar)
 148                     (* 3.0 ass))
 149                    an2)
 150                   cospar)
 151                  an2))
 152         (setf (f2cl-lib:fref v ((f2cl-lib:int-add noequ 3)) ((1 28)))
 153                 (+ (f2cl-lib:fref v ((f2cl-lib:int-add noequ 3)) ((1 28)))
 154                    (* -2.0 asap par2 (- an 1.0) (- an 2.0))))
 155         (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 156             (dgtsl noequ d1 d d2
 157              (f2cl-lib:array-slice v double-float (4) ((1 28))) iers)
 158           (declare (ignore var-0 var-1 var-2 var-3 var-4))
 159           (setf iers var-5))
 160         (go label50)
 161        label30
 162         (setf an 4.0)
 163         (f2cl-lib:fdo (i 4 (f2cl-lib:int-add i 1))
 164                       ((> i 13) nil)
 165           (tagbody
 166             (setf an2 (* an an))
 167             (setf (f2cl-lib:fref v (i) ((1 28)))
 168                     (/
 169                      (+
 170                       (* (- an2 4.0)
 171                          (-
 172                           (* 2.0
 173                              (- par22 an2 an2)
 174                              (f2cl-lib:fref v
 175                                             ((f2cl-lib:int-sub i 1))
 176                                             ((1 28))))
 177                           ac))
 178                       as
 179                       (* (- par2)
 180                          (+ an 1.0)
 181                          (+ an 2.0)
 182                          (f2cl-lib:fref v ((f2cl-lib:int-sub i 2)) ((1 28)))))
 183                      (* par2 (- an 1.0) (- an 2.0))))
 184             (setf an (+ an 2.0))
 185            label40))
 186        label50
 187         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 188                       ((> j 13) nil)
 189           (tagbody
 190             (setf (f2cl-lib:fref chebmo-%data%
 191                                  (m
 192                                   (f2cl-lib:int-sub (f2cl-lib:int-mul 2 j) 1))
 193                                  ((1 maxp1) (1 25))
 194                                  chebmo-%offset%)
 195                     (f2cl-lib:fref v (j) ((1 28))))
 196            label60))
 197         (setf (f2cl-lib:fref v (1) ((1 28)))
 198                 (/ (* 2.0 (- sinpar (* parint cospar))) par2))
 199         (setf (f2cl-lib:fref v (2) ((1 28)))
 200                 (+ (/ (* (+ 18.0 (/ -48.0 par2)) sinpar) par2)
 201                    (/ (* (- (/ 48.0 par2) 2.0) cospar) parint)))
 202         (setf ac (* -24.0 parint cospar))
 203         (setf as (* -8.0 sinpar))
 204         (if (> (abs parint) 24.0) (go label80))
 205         (setf an 5.0)
 206         (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 207                       ((> k noeq1) nil)
 208           (tagbody
 209             (setf an2 (* an an))
 210             (setf (f2cl-lib:fref d (k) ((1 25)))
 211                     (* -2.0 (- an2 4.0) (- par22 an2 an2)))
 212             (setf (f2cl-lib:fref d2 (k) ((1 25)))
 213                     (* (- an 1.0) (- an 2.0) par2))
 214             (setf (f2cl-lib:fref d1 ((f2cl-lib:int-add k 1)) ((1 25)))
 215                     (* (+ an 3.0) (+ an 4.0) par2))
 216             (setf (f2cl-lib:fref v ((f2cl-lib:int-add k 2)) ((1 28)))
 217                     (+ ac (* (- an2 4.0) as)))
 218             (setf an (+ an 2.0))
 219            label70))
 220         (setf an2 (* an an))
 221         (setf (f2cl-lib:fref d (noequ) ((1 25)))
 222                 (* -2.0 (- an2 4.0) (- par22 an2 an2)))
 223         (setf (f2cl-lib:fref v ((f2cl-lib:int-add noequ 2)) ((1 28)))
 224                 (+ ac (* (- an2 4.0) as)))
 225         (setf (f2cl-lib:fref v (3) ((1 28)))
 226                 (+ (f2cl-lib:fref v (3) ((1 28)))
 227                    (* -42.0 par2 (f2cl-lib:fref v (2) ((1 28))))))
 228         (setf ass (* parint cospar))
 229         (setf asap
 230                 (/
 231                  (-
 232                   (/
 233                    (-
 234                     (/
 235                      (-
 236                       (+
 237                        (/
 238                         (+ (* (- (* 105.0 par2) 63.0) ass)
 239                            (* (- (* 210.0 par2) 1.0) sinpar))
 240                         an2)
 241                        (* (- (* 15.0 par2) 1.0) sinpar))
 242                       (* 15.0 ass))
 243                      an2)
 244                     (* 3.0 ass)
 245                     sinpar)
 246                    an2)
 247                   sinpar)
 248                  an2))
 249         (setf (f2cl-lib:fref v ((f2cl-lib:int-add noequ 2)) ((1 28)))
 250                 (+ (f2cl-lib:fref v ((f2cl-lib:int-add noequ 2)) ((1 28)))
 251                    (* -2.0 asap par2 (- an 1.0) (- an 2.0))))
 252         (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 253             (dgtsl noequ d1 d d2
 254              (f2cl-lib:array-slice v double-float (3) ((1 28))) iers)
 255           (declare (ignore var-0 var-1 var-2 var-3 var-4))
 256           (setf iers var-5))
 257         (go label100)
 258        label80
 259         (setf an 3.0)
 260         (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 261                       ((> i 12) nil)
 262           (tagbody
 263             (setf an2 (* an an))
 264             (setf (f2cl-lib:fref v (i) ((1 28)))
 265                     (/
 266                      (+
 267                       (* (- an2 4.0)
 268                          (+
 269                           (* 2.0
 270                              (- par22 an2 an2)
 271                              (f2cl-lib:fref v
 272                                             ((f2cl-lib:int-sub i 1))
 273                                             ((1 28))))
 274                           as))
 275                       ac
 276                       (* (- par2)
 277                          (+ an 1.0)
 278                          (+ an 2.0)
 279                          (f2cl-lib:fref v ((f2cl-lib:int-sub i 2)) ((1 28)))))
 280                      (* par2 (- an 1.0) (- an 2.0))))
 281             (setf an (+ an 2.0))
 282            label90))
 283        label100
 284         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 285                       ((> j 12) nil)
 286           (tagbody
 287             (setf (f2cl-lib:fref chebmo-%data%
 288                                  (m (f2cl-lib:int-mul 2 j))
 289                                  ((1 maxp1) (1 25))
 290                                  chebmo-%offset%)
 291                     (f2cl-lib:fref v (j) ((1 28))))
 292            label110))
 293        label120
 294         (if (< nrmom momcom) (setf m (f2cl-lib:int-add nrmom 1)))
 295         (if (and (< momcom (f2cl-lib:int-sub maxp1 1)) (>= nrmom momcom))
 296             (setf momcom (f2cl-lib:int-add momcom 1)))
 297         (setf (f2cl-lib:fref fval (1) ((1 25)))
 298                 (* 0.5 (funcall f (+ centr hlgth))))
 299         (setf (f2cl-lib:fref fval (13) ((1 25)))
 300                 (multiple-value-bind (ret-val var-0)
 301                     (funcall f centr)
 302                   (declare (ignore))
 303                   (when var-0
 304                     (setf centr var-0))
 305                   ret-val))
 306         (setf (f2cl-lib:fref fval (25) ((1 25)))
 307                 (* 0.5 (funcall f (- centr hlgth))))
 308         (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
 309                       ((> i 12) nil)
 310           (tagbody
 311             (setf isym (f2cl-lib:int-sub 26 i))
 312             (setf (f2cl-lib:fref fval (i) ((1 25)))
 313                     (funcall f
 314                              (+
 315                               (* hlgth
 316                                  (f2cl-lib:fref x
 317                                                 ((f2cl-lib:int-sub i 1))
 318                                                 ((1 11))))
 319                               centr)))
 320             (setf (f2cl-lib:fref fval (isym) ((1 25)))
 321                     (funcall f
 322                              (- centr
 323                                 (* hlgth
 324                                    (f2cl-lib:fref x
 325                                                   ((f2cl-lib:int-sub i 1))
 326                                                   ((1 11)))))))
 327            label130))
 328         (dqcheb x fval cheb12 cheb24)
 329         (setf resc12
 330                 (* (f2cl-lib:fref cheb12 (13) ((1 13)))
 331                    (f2cl-lib:fref chebmo-%data%
 332                                   (m 13)
 333                                   ((1 maxp1) (1 25))
 334                                   chebmo-%offset%)))
 335         (setf ress12 0.0)
 336         (setf k 11)
 337         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 338                       ((> j 6) nil)
 339           (tagbody
 340             (setf resc12
 341                     (+ resc12
 342                        (* (f2cl-lib:fref cheb12 (k) ((1 13)))
 343                           (f2cl-lib:fref chebmo-%data%
 344                                          (m k)
 345                                          ((1 maxp1) (1 25))
 346                                          chebmo-%offset%))))
 347             (setf ress12
 348                     (+ ress12
 349                        (*
 350                         (f2cl-lib:fref cheb12
 351                                        ((f2cl-lib:int-add k 1))
 352                                        ((1 13)))
 353                         (f2cl-lib:fref chebmo-%data%
 354                                        (m (f2cl-lib:int-add k 1))
 355                                        ((1 maxp1) (1 25))
 356                                        chebmo-%offset%))))
 357             (setf k (f2cl-lib:int-sub k 2))
 358            label140))
 359         (setf resc24
 360                 (* (f2cl-lib:fref cheb24 (25) ((1 25)))
 361                    (f2cl-lib:fref chebmo-%data%
 362                                   (m 25)
 363                                   ((1 maxp1) (1 25))
 364                                   chebmo-%offset%)))
 365         (setf ress24 0.0)
 366         (setf resabs (abs (f2cl-lib:fref cheb24 (25) ((1 25)))))
 367         (setf k 23)
 368         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 369                       ((> j 12) nil)
 370           (tagbody
 371             (setf resc24
 372                     (+ resc24
 373                        (* (f2cl-lib:fref cheb24 (k) ((1 25)))
 374                           (f2cl-lib:fref chebmo-%data%
 375                                          (m k)
 376                                          ((1 maxp1) (1 25))
 377                                          chebmo-%offset%))))
 378             (setf ress24
 379                     (+ ress24
 380                        (*
 381                         (f2cl-lib:fref cheb24
 382                                        ((f2cl-lib:int-add k 1))
 383                                        ((1 25)))
 384                         (f2cl-lib:fref chebmo-%data%
 385                                        (m (f2cl-lib:int-add k 1))
 386                                        ((1 maxp1) (1 25))
 387                                        chebmo-%offset%))))
 388             (setf resabs
 389                     (+ resabs
 390                        (abs (f2cl-lib:fref cheb24 (k) ((1 25))))
 391                        (abs
 392                         (f2cl-lib:fref cheb24
 393                                        ((f2cl-lib:int-add k 1))
 394                                        ((1 25))))))
 395             (setf k (f2cl-lib:int-sub k 2))
 396            label150))
 397         (setf estc (abs (- resc24 resc12)))
 398         (setf ests (abs (- ress24 ress12)))
 399         (setf resabs (* resabs (abs hlgth)))
 400         (if (= integr 2) (go label160))
 401         (setf result (- (* conc resc24) (* cons$ ress24)))
 402         (setf abserr (+ (abs (* conc estc)) (abs (* cons$ ests))))
 403         (go label170)
 404        label160
 405         (setf result (+ (* conc ress24) (* cons$ resc24)))
 406         (setf abserr (+ (abs (* conc ests)) (abs (* cons$ estc))))
 407        label170
 408         (go end_label)
 409        end_label
 410         (return
 411          (values nil
 412                  nil
 413                  nil
 414                  omega
 415                  integr
 416                  nil
 417                  nil
 418                  nil
 419                  result
 420                  abserr
 421                  neval
 422                  resabs
 423                  resasc
 424                  momcom
 425                  nil))))))
 426
 427 (in-package #:cl-user)
 428 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 429 (eval-when (:load-toplevel :compile-toplevel :execute)
 430   (setf (gethash 'fortran-to-lisp::dqc25f
 431                  fortran-to-lisp::*f2cl-function-info*)
 432           (fortran-to-lisp::make-f2cl-finfo
 433            :arg-types '(t (double-float) (double-float) (double-float)
 434                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 435                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 436                         (double-float) (double-float)
 437                         (fortran-to-lisp::integer4) (double-float)
 438                         (double-float) (fortran-to-lisp::integer4)
 439                         (array double-float (*)))
 440            :return-values '(nil nil nil fortran-to-lisp::omega
 441                             fortran-to-lisp::integr nil nil nil
 442                             fortran-to-lisp::result fortran-to-lisp::abserr
 443                             fortran-to-lisp::neval fortran-to-lisp::resabs
 444                             fortran-to-lisp::resasc fortran-to-lisp::momcom
 445                             nil)
 446            :calls '(fortran-to-lisp::dqcheb fortran-to-lisp::dgtsl
 447                     fortran-to-lisp::dqk15w fortran-to-lisp::d1mach))))
 448