*
* $Id: dgamma.F,v 1.1.1.1 1996/02/15 17:48:17 mclareni Exp $
*
* $Log: dgamma.F,v $
* Revision 1.1.1.1  1996/02/15 17:48:17  mclareni
* Kernlib
*
*
#include "kernnum/pilot.h"
      DOUBLE PRECISION FUNCTION DGAMMA(X)
      LOGICAL MFLAG,RFLAG
      REAL SX
      DOUBLE PRECISION X,U,F,ZERO,ONE,THREE,FOUR,PI
      DOUBLE PRECISION C(0:24),H,ALFA,B0,B1,B2
      DATA ZERO /0.0D0/, ONE /1.0D0/, THREE /3.0D0/, FOUR /4.0D0/
#if defined(CERNLIB_NUMHIPRE)
      DATA NC /24/
      DATA PI    /3.14159 26535 89793 23846 26433 83D0/
      DATA C( 0) /3.65738 77250 83382 43849 88068 39D0/
      DATA C( 1) /1.95754 34566 61268 26928 33742 26D0/
      DATA C( 2) / .33829 71138 26160 38915 58510 73D0/
      DATA C( 3) / .04208 95127 65575 49198 51083 97D0/
      DATA C( 4) / .00428 76504 82129 08770 04289 08D0/
      DATA C( 5) / .00036 52121 69294 61767 02198 22D0/
      DATA C( 6) / .00002 74006 42226 42200 27170 66D0/
      DATA C( 7) / .00000 18124 02333 65124 44603 05D0/
      DATA C( 8) / .00000 01096 57758 65997 06993 06D0/
      DATA C( 9) / .00000 00059 87184 04552 00046 95D0/
      DATA C(10) / .00000 00003 07690 80535 24777 71D0/
      DATA C(11) / .00000 00000 14317 93029 61915 76D0/
      DATA C(12) / .00000 00000 00651 08773 34803 70D0/
      DATA C(13) / .00000 00000 00025 95849 89822 28D0/
      DATA C(14) / .00000 00000 00001 10789 38922 59D0/
      DATA C(15) / .00000 00000 00000 03547 43620 17D0/
      DATA C(16) / .00000 00000 00000 00168 86075 04D0/
      DATA C(17) / .00000 00000 00000 00002 73543 58D0/
      DATA C(18) / .00000 00000 00000 00000 30297 74D0/
      DATA C(19) /-.00000 00000 00000 00000 00571 22D0/
      DATA C(20) / .00000 00000 00000 00000 00090 77D0/
      DATA C(21) /-.00000 00000 00000 00000 00005 05D0/
      DATA C(22) / .00000 00000 00000 00000 00000 41D0/
      DATA C(23) /-.00000 00000 00000 00000 00000 03D0/
      DATA C(24) / .00000 00000 00000 00000 00000 01D0/
#endif
#if defined(CERNLIB_NUMLOPRE)
      DATA NC /15/
      DATA PI    /3.14159 26535 89793 24D0/
      DATA C( 0) /3.65738 77250 83382 44D0/
      DATA C( 1) /1.95754 34566 61268 27D0/
      DATA C( 2) / .33829 71138 26160 39D0/
      DATA C( 3) / .04208 95127 65575 49D0/
      DATA C( 4) / .00428 76504 82129 09D0/
      DATA C( 5) / .00036 52121 69294 62D0/
      DATA C( 6) / .00002 74006 42226 42D0/
      DATA C( 7) / .00000 18124 02333 65D0/
      DATA C( 8) / .00000 01096 57758 66D0/
      DATA C( 9) / .00000 00059 87184 05D0/
      DATA C(10) / .00000 00003 07690 81D0/
      DATA C(11) / .00000 00000 14317 93D0/
      DATA C(12) / .00000 00000 00651 09D0/
      DATA C(13) / .00000 00000 00025 96D0/
      DATA C(14) / .00000 00000 00001 11D0/
      DATA C(15) / .00000 00000 00000 04D0/
#endif
      U=X
      IF(X .LE. ZERO) THEN
       IF(X .EQ. INT(X)) THEN
        CALL KERMTR('C305.1',LGFILE,MFLAG,RFLAG)
        IF(MFLAG) THEN
         SX=X
         IF(LGFILE .EQ. 0) THEN
          WRITE(*,100) SX
         ELSE
          WRITE(LGFILE,100) SX
         END IF
        END IF
        IF(.NOT.RFLAG) CALL ABEND
        DGAMMA=ZERO
        RETURN
       ELSE
        U=ONE-U
       END IF
      END IF
      F=ONE
      IF(U .LT. THREE) THEN
       DO 1 I = 1,INT(FOUR-U)
       F=F/U
    1  U=U+ONE
      ELSE
       DO 2 I = 1,INT(U-THREE)
       U=U-ONE
    2  F=F*U
      END IF
      U=U-THREE
      H=U+U-ONE
      ALFA=H+H
      B1=ZERO
      B2=ZERO
      DO 3 I = NC,0,-1
      B0=C(I)+ALFA*B1-B2
      B2=B1
    3 B1=B0
      U=F*(B0-H*B2)
      IF(X .LT. ZERO) U=PI/(SIN(PI*X)*U)
      DGAMMA=U
      RETURN
  100 FORMAT(1X,'DGAMMA ... ARGUMENT IS NON-POSITIVE INTEGER = ',E15.1)
      END