src/numerical/slatec/dai.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 ':simple-array)
  14 ;;;           (:array-slicing nil) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :slatec)
  18
  19
  20 (let ((naif 0)
  21       (naig 0)
  22       (x3sml 0.0)
  23       (xmax 0.0)
  24       (aifcs
  25        (make-array 13
  26                    :element-type 'double-float
  27                    :initial-contents '(-0.03797135849667 0.05919188853726364
  28                                        9.862928057727998e-4
  29                                        6.848843819076567e-6
  30                                        2.5942025962194713e-8
  31                                        6.176612774081375e-11
  32                                        1.0092454172466118e-13
  33                                        1.2014792511179938e-16
  34                                        1.0882945588716992e-19
  35                                        7.751377219668488e-23
  36                                        4.4548112037175636e-26
  37                                        2.1092845231692343e-29
  38                                        8.370173591074134e-33)))
  39       (aigcs
  40        (make-array 13
  41                    :element-type 'double-float
  42                    :initial-contents '(0.018152365581161272
  43                                        0.021572563166010757
  44                                        2.567835698748325e-4
  45                                        1.4265214119792405e-6
  46                                        4.572114920018043e-9
  47                                        9.52517084356471e-12
  48                                        1.3925634605771398e-14
  49                                        1.5070999142762378e-17
  50                                        1.2559148312567778e-20
  51                                        8.306307377082133e-24
  52                                        4.465753849371857e-27
  53                                        1.9900855034518868e-30
  54                                        7.4702885256533335e-34)))
  55       (first$ nil))
  56   (declare (type (f2cl-lib:integer4) naif naig)
  57            (type (double-float) x3sml xmax)
  58            (type (simple-array double-float (13)) aifcs aigcs)
  59            (type f2cl-lib:logical first$))
  60   (setq first$ f2cl-lib:%true%)
  61   (defun dai (x)
  62     (declare (type (double-float) x))
  63     (prog ((theta 0.0) (xm 0.0) (z 0.0) (xmaxt 0.0) (dai 0.0))
  64       (declare (type (double-float) dai xmaxt z xm theta))
  65       (cond
  66         (first$
  67          (setf naif
  68                  (initds aifcs 13
  69                   (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3)))))
  70          (setf naig
  71                  (initds aigcs 13
  72                   (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3)))))
  73          (setf x3sml (expt (f2cl-lib:d1mach 3) 0.3334))
  74          (setf xmaxt (expt (* -1.5 (f2cl-lib:flog (f2cl-lib:d1mach 1))) 0.6667))
  75          (setf xmax
  76                  (-
  77                   (+ xmaxt
  78                      (/ (* (- xmaxt) (f2cl-lib:flog xmaxt))
  79                         (+ (* 4.0 (f2cl-lib:fsqrt xmaxt)) 1.0)))
  80                   0.01))))
  81       (setf first$ f2cl-lib:%false%)
  82       (if (>= x -1.0) (go label20))
  83       (multiple-value-bind (var-0 var-1 var-2)
  84           (d9aimp x xm theta)
  85         (declare (ignore var-0))
  86         (setf xm var-1)
  87         (setf theta var-2))
  88       (setf dai (* xm (cos theta)))
  89       (go end_label)
  90      label20
  91       (if (> x 1.0) (go label30))
  92       (setf z 0.0)
  93       (if (> (abs x) x3sml) (setf z (expt x 3)))
  94       (setf dai
  95               (+ 0.375
  96                  (- (dcsevl z aifcs naif)
  97                     (* x (+ 0.25 (dcsevl z aigcs naig))))))
  98       (go end_label)
  99      label30
 100       (if (> x xmax) (go label40))
 101       (setf dai (* (daie x) (exp (/ (* -2.0 x (f2cl-lib:fsqrt x)) 3.0))))
 102       (go end_label)
 103      label40
 104       (setf dai 0.0)
 105       (xermsg "SLATEC" "DAI" "X SO BIG AI UNDERFLOWS" 1 1)
 106       (go end_label)
 107      end_label
 108       (return (values dai nil)))))
 109
 110 (in-package #:cl-user)
 111 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 112 (eval-when (:load-toplevel :compile-toplevel :execute)
 113   (setf (gethash 'fortran-to-lisp::dai fortran-to-lisp::*f2cl-function-info*)
 114           (fortran-to-lisp::make-f2cl-finfo :arg-types '((double-float))
 115                                             :return-values '(nil)
 116                                             :calls '(fortran-to-lisp::xermsg
 117                                                      fortran-to-lisp::daie
 118                                                      fortran-to-lisp::dcsevl
 119                                                      fortran-to-lisp::d9aimp
 120                                                      fortran-to-lisp::initds
 121                                                      fortran-to-lisp::d1mach))))
 122