src/numerical/slatec/dasyik.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 nil) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :slatec)
  18
  19
  20 (let ((con
  21        (make-array 2
  22                    :element-type 'double-float
  23                    :initial-contents '(0.3989422804014327 1.2533141373155003)))
  24       (c
  25        (make-array 65
  26                    :element-type 'double-float
  27                    :initial-contents '(-0.208333333333333 0.125
  28                                        0.334201388888889 -0.401041666666667
  29                                        0.0703125 -1.02581259645062
  30                                        1.84646267361111 -0.8912109375
  31                                        0.0732421875 4.66958442342625
  32                                        -11.207002616223 8.78912353515625
  33                                        -2.3640869140625 0.112152099609375
  34                                        -28.2120725582002 84.6362176746007
  35                                        -91.81824154324 42.5349987453885
  36                                        -7.36879435947963 0.227108001708984
  37                                        212.570130039217 -765.252468141182
  38                                        1059.990452528 -699.579627376133
  39                                        218.190511744212 -26.4914304869516
  40                                        0.572501420974731 -1919.45766231841
  41                                        8061.72218173731 -13586.5500064341
  42                                        11655.3933368645 -5305.6469786134
  43                                        1200.90291321635 -108.090919788395
  44                                        1.72772750258446 20204.2913309661
  45                                        -96980.5983886375 192547.001232532
  46                                        -203400.177280416 122200.464983017
  47                                        -41192.6549688976 7109.51430248936
  48                                        -493.915304773088 6.07404200127348
  49                                        -242919.187900551 1311763.61466298
  50                                        -2998015.91853811 3763271.2976564
  51                                        -2813563.22658653 1268365.27332162
  52                                        -331645.172484564 45218.7689813627
  53                                        -2499.83048181121 24.3805296995561
  54                                        3284469.85307204 -1.97068191184322e7
  55                                        5.09526024926646e7 -7.41051482115327e7
  56                                        6.6344512274729e7 -3.75671766607634e7
  57                                        1.32887671664218e7 -2785618.12808645
  58                                        308186.404612662 -13886.089753717
  59                                        110.017140269247))))
  60   (declare (type (array double-float (2)) con)
  61            (type (array double-float (65)) c))
  62   (defun dasyik (x fnu kode flgik ra arg in y)
  63     (declare (type (array double-float (*)) y)
  64              (type (f2cl-lib:integer4) in kode)
  65              (type (double-float) arg ra flgik fnu x))
  66     (f2cl-lib:with-multi-array-data
  67         ((y double-float y-%data% y-%offset%))
  68       (prog ((ak 0.0) (ap 0.0) (coef 0.0) (etx 0.0) (fn 0.0) (gln 0.0) (s1 0.0)
  69              (s2 0.0) (t$ 0.0) (tol 0.0) (t2 0.0) (z 0.0) (j 0) (jn 0) (k 0)
  70              (kk 0) (l 0))
  71         (declare (type (f2cl-lib:integer4) l kk k jn j)
  72                  (type (double-float) z t2 tol t$ s2 s1 gln fn etx coef ap ak))
  73         (setf tol (f2cl-lib:d1mach 3))
  74         (setf tol (max tol 1.0e-15))
  75         (setf fn fnu)
  76         (setf z (/ (- 3.0 flgik) 2.0))
  77         (setf kk (f2cl-lib:int z))
  78         (f2cl-lib:fdo (jn 1 (f2cl-lib:int-add jn 1))
  79                       ((> jn in) nil)
  80           (tagbody
  81             (if (= jn 1) (go label10))
  82             (setf fn (- fn flgik))
  83             (setf z (/ x fn))
  84             (setf ra (f2cl-lib:fsqrt (+ 1.0 (* z z))))
  85             (setf gln (f2cl-lib:flog (/ (+ 1.0 ra) z)))
  86             (setf etx
  87                     (coerce (the f2cl-lib:integer4 (f2cl-lib:int-sub kode 1))
  88                             'double-float))
  89             (setf t$ (+ (* ra (- 1.0 etx)) (/ etx (+ z ra))))
  90             (setf arg (* fn (- t$ gln) flgik))
  91            label10
  92             (setf coef (exp arg))
  93             (setf t$ (/ 1.0 ra))
  94             (setf t2 (* t$ t$))
  95             (setf t$ (/ t$ fn))
  96             (setf t$ (f2cl-lib:sign t$ flgik))
  97             (setf s2 1.0)
  98             (setf ap 1.0)
  99             (setf l 0)
 100             (f2cl-lib:fdo (k 2 (f2cl-lib:int-add k 1))
 101                           ((> k 11) nil)
 102               (tagbody
 103                 (setf l (f2cl-lib:int-add l 1))
 104                 (setf s1 (f2cl-lib:fref c (l) ((1 65))))
 105                 (f2cl-lib:fdo (j 2 (f2cl-lib:int-add j 1))
 106                               ((> j k) nil)
 107                   (tagbody
 108                     (setf l (f2cl-lib:int-add l 1))
 109                     (setf s1 (+ (* s1 t2) (f2cl-lib:fref c (l) ((1 65)))))
 110                    label20))
 111                 (setf ap (* ap t$))
 112                 (setf ak (* ap s1))
 113                 (setf s2 (+ s2 ak))
 114                 (if (< (max (abs ak) (abs ap)) tol) (go label40))
 115                label30))
 116            label40
 117             (setf t$ (abs t$))
 118             (setf (f2cl-lib:fref y-%data% (jn) ((1 *)) y-%offset%)
 119                     (* s2
 120                        coef
 121                        (f2cl-lib:fsqrt t$)
 122                        (f2cl-lib:fref con (kk) ((1 2)))))
 123            label50))
 124         (go end_label)
 125        end_label
 126         (return (values nil nil nil nil ra arg nil nil))))))
 127
 128 (in-package #:cl-user)
 129 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 130 (eval-when (:load-toplevel :compile-toplevel :execute)
 131   (setf (gethash 'fortran-to-lisp::dasyik
 132                  fortran-to-lisp::*f2cl-function-info*)
 133           (fortran-to-lisp::make-f2cl-finfo
 134            :arg-types '((double-float) (double-float)
 135                         (fortran-to-lisp::integer4) (double-float)
 136                         (double-float) (double-float)
 137                         (fortran-to-lisp::integer4) (array double-float (*)))
 138            :return-values '(nil nil nil nil fortran-to-lisp::ra
 139                             fortran-to-lisp::arg nil nil)
 140            :calls '(fortran-to-lisp::d1mach))))
 141