GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
d9lgmc.f
Go to the documentation of this file.
1 *DECK D9LGMC
2  DOUBLE PRECISION FUNCTION d9lgmc (X)
3 C***BEGIN PROLOGUE D9LGMC
4 C***SUBSIDIARY
5 C***PURPOSE Compute the log Gamma correction factor so that
6 C LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X
7 C + D9LGMC(X).
8 C***LIBRARY SLATEC (FNLIB)
9 C***CATEGORY C7E
10 C***TYPE DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
11 C***KEYWORDS COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
12 C LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
13 C***AUTHOR Fullerton, W., (LANL)
14 C***DESCRIPTION
15 C
16 C Compute the log gamma correction factor for X .GE. 10. so that
17 C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X)
18 C
19 C Series for ALGM on the interval 0. to 1.00000E-02
20 C with weighted error 1.28E-31
21 C log weighted error 30.89
22 C significant figures required 29.81
23 C decimal places required 31.48
24 C
25 C***REFERENCES (NONE)
26 C***ROUTINES CALLED D1MACH, DCSEVL, INITDS, XERMSG
27 C***REVISION HISTORY (YYMMDD)
28 C 770601 DATE WRITTEN
29 C 890531 Changed all specific intrinsics to generic. (WRB)
30 C 890531 REVISION DATE from Version 3.2
31 C 891214 Prologue converted to Version 4.0 format. (BAB)
32 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
33 C 900720 Routine changed from user-callable to subsidiary. (WRB)
34 C***END PROLOGUE D9LGMC
35  DOUBLE PRECISION x, algmcs(15), xbig, xmax, dcsevl, d1mach
36  LOGICAL first
37  SAVE algmcs, nalgm, xbig, xmax, first
38  DATA algmcs( 1) / +.1666389480 4518632472 0572965082 2 d+0 /
39  DATA algmcs( 2) / -.1384948176 0675638407 3298605913 5 d-4 /
40  DATA algmcs( 3) / +.9810825646 9247294261 5717154748 7 d-8 /
41  DATA algmcs( 4) / -.1809129475 5724941942 6330626671 9 d-10 /
42  DATA algmcs( 5) / +.6221098041 8926052271 2601554341 6 d-13 /
43  DATA algmcs( 6) / -.3399615005 4177219443 0333059966 6 d-15 /
44  DATA algmcs( 7) / +.2683181998 4826987489 5753884666 6 d-17 /
45  DATA algmcs( 8) / -.2868042435 3346432841 4462239999 9 d-19 /
46  DATA algmcs( 9) / +.3962837061 0464348036 7930666666 6 d-21 /
47  DATA algmcs( 10) / -.6831888753 9857668701 1199999999 9 d-23 /
48  DATA algmcs( 11) / +.1429227355 9424981475 7333333333 3 d-24 /
49  DATA algmcs( 12) / -.3547598158 1010705471 9999999999 9 d-26 /
50  DATA algmcs( 13) / +.1025680058 0104709120 0000000000 0 d-27 /
51  DATA algmcs( 14) / -.3401102254 3167487999 9999999999 9 d-29 /
52  DATA algmcs( 15) / +.1276642195 6300629333 3333333333 3 d-30 /
53  DATA first /.true./
54 C***FIRST EXECUTABLE STATEMENT D9LGMC
55  IF (first) THEN
56  nalgm = initds(algmcs, 15, REAL(D1MACH(3)) )
57  xbig = 1.0d0/sqrt(d1mach(3))
58  xmax = exp(min(log(d1mach(2)/12.d0), -log(12.d0*d1mach(1))))
59  ENDIF
60  first = .false.
61 C
62  IF (x .LT. 10.d0) CALL xermsg('SLATEC', 'D9LGMC',
63  + 'X MUST BE GE 10', 1, 2)
64  IF (x.GE.xmax) go to 20
65 C
66  d9lgmc = 1.d0/(12.d0*x)
67  IF (x.LT.xbig) d9lgmc = dcsevl(2.0d0*(10.d0/x)**2-1.d0, algmcs,
68  1 nalgm) / x
69  RETURN
70 C
71  20 d9lgmc = 0.d0
72  CALL xermsg('SLATEC', 'D9LGMC', 'X SO BIG D9LGMC UNDERFLOWS', 2,
73  + 1)
74  RETURN
75 C
76  END