src/numerical/slatec/dgamma.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 ((ngam 0)
  21       (xmin 0.0)
  22       (xmax 0.0)
  23       (dxrel 0.0)
  24       (gamcs
  25        (make-array 42
  26                    :element-type 'double-float
  27                    :initial-contents '(0.00857119559098933 0.004415381324841007
  28                                        0.05685043681599363 -0.00421983539641856
  29                                        0.0013268081812124603
  30                                        -1.8930245297988805e-4
  31                                        3.606925327441245e-5
  32                                        -6.056761904460864e-6
  33                                        1.0558295463022833e-6
  34                                        -1.811967365542384e-7
  35                                        3.117724964715322e-8
  36                                        -5.354219639019687e-9
  37                                        9.193275519859589e-10
  38                                        -1.5779412802883398e-10
  39                                        2.7079806229349544e-11
  40                                        -4.64681865382573e-12
  41                                        7.97335019200742e-13
  42                                        -1.368078209830916e-13
  43                                        2.3473194865638007e-14
  44                                        -4.027432614949067e-15
  45                                        6.910051747372101e-16
  46                                        -1.185584500221993e-16
  47                                        2.034148542496374e-17
  48                                        -3.490054341717406e-18
  49                                        5.987993856485306e-19
  50                                        -1.027378057872228e-19
  51                                        1.7627028160605298e-20
  52                                        -3.024320653735306e-21
  53                                        5.188914660218398e-22
  54                                        -8.902770842456576e-23
  55                                        1.5274740684933426e-23
  56                                        -2.620731256187363e-24
  57                                        4.496464047830539e-25
  58                                        -7.714712731336878e-26
  59                                        1.323635453126044e-26
  60                                        -2.2709994129429287e-27
  61                                        3.8964189980039913e-28
  62                                        -6.685198115125953e-29
  63                                        1.1469986631400244e-29
  64                                        -1.9679385863451348e-30
  65                                        3.376448816585338e-31
  66                                        -5.793070335782136e-32)))
  67       (pi$ 3.141592653589793)
  68       (sq2pil 0.9189385332046728)
  69       (first$ nil))
  70   (declare (type (f2cl-lib:integer4) ngam)
  71            (type (double-float) xmin xmax dxrel pi$ sq2pil)
  72            (type (simple-array double-float (42)) gamcs)
  73            (type f2cl-lib:logical first$))
  74   (setq first$ f2cl-lib:%true%)
  75   (defun dgamma (x)
  76     (declare (type (double-float) x))
  77     (prog ((sinpiy 0.0) (y 0.0) (dgamma 0.0) (i 0) (n 0))
  78       (declare (type (f2cl-lib:integer4) n i)
  79                (type (double-float) dgamma y sinpiy))
  80       (cond
  81         (first$
  82          (setf ngam
  83                  (initds gamcs 42
  84                   (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3)))))
  85          (multiple-value-bind (var-0 var-1)
  86              (dgamlm xmin xmax)
  87            (declare (ignore))
  88            (setf xmin var-0)
  89            (setf xmax var-1))
  90          (setf dxrel (f2cl-lib:fsqrt (f2cl-lib:d1mach 4)))))
  91       (setf first$ f2cl-lib:%false%)
  92       (setf y (abs x))
  93       (if (> y 10.0) (go label50))
  94       (setf n (f2cl-lib:int x))
  95       (if (< x 0.0) (setf n (f2cl-lib:int-sub n 1)))
  96       (setf y (- x n))
  97       (setf n (f2cl-lib:int-sub n 1))
  98       (setf dgamma (+ 0.9375 (dcsevl (- (* 2.0 y) 1.0) gamcs ngam)))
  99       (if (= n 0) (go end_label))
 100       (if (> n 0) (go label30))
 101       (setf n (f2cl-lib:int-sub n))
 102       (if (= x 0.0) (xermsg "SLATEC" "DGAMMA" "X IS 0" 4 2))
 103       (if (and (< x 0.0f0) (= (- (+ x n) 2) 0.0))
 104           (xermsg "SLATEC" "DGAMMA" "X IS A NEGATIVE INTEGER" 4 2))
 105       (if
 106        (and (< x -0.5) (< (abs (/ (- x (f2cl-lib:aint (- x 0.5))) x)) dxrel))
 107        (xermsg "SLATEC" "DGAMMA"
 108         "ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER" 1 1))
 109       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 110                     ((> i n) nil)
 111         (tagbody (setf dgamma (/ dgamma (- (+ x i) 1))) label20))
 112       (go end_label)
 113      label30
 114       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 115                     ((> i n) nil)
 116         (tagbody (setf dgamma (* (+ y i) dgamma)) label40))
 117       (go end_label)
 118      label50
 119       (if (> x xmax) (xermsg "SLATEC" "DGAMMA" "X SO BIG GAMMA OVERFLOWS" 3 2))
 120       (setf dgamma 0.0)
 121       (if (< x xmin)
 122           (xermsg "SLATEC" "DGAMMA" "X SO SMALL GAMMA UNDERFLOWS" 2 1))
 123       (if (< x xmin) (go end_label))
 124       (setf dgamma
 125               (exp
 126                (+ (- (* (- y 0.5) (f2cl-lib:flog y)) y) sq2pil (d9lgmc y))))
 127       (if (> x 0.0) (go end_label))
 128       (if (< (abs (/ (- x (f2cl-lib:aint (- x 0.5))) x)) dxrel)
 129           (xermsg "SLATEC" "DGAMMA"
 130            "ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER" 1 1))
 131       (setf sinpiy (sin (* pi$ y)))
 132       (if (= sinpiy 0.0)
 133           (xermsg "SLATEC" "DGAMMA" "X IS A NEGATIVE INTEGER" 4 2))
 134       (setf dgamma (/ (- pi$) (* y sinpiy dgamma)))
 135       (go end_label)
 136      end_label
 137       (return (values dgamma nil)))))
 138
 139 (in-package #:cl-user)
 140 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 141 (eval-when (:load-toplevel :compile-toplevel :execute)
 142   (setf (gethash 'fortran-to-lisp::dgamma
 143                  fortran-to-lisp::*f2cl-function-info*)
 144           (fortran-to-lisp::make-f2cl-finfo :arg-types '((double-float))
 145                                             :return-values '(nil)
 146                                             :calls '(fortran-to-lisp::d9lgmc
 147                                                      fortran-to-lisp::xermsg
 148                                                      fortran-to-lisp::dcsevl
 149                                                      fortran-to-lisp::dgamlm
 150                                                      fortran-to-lisp::initds
 151                                                      fortran-to-lisp::d1mach))))
 152