Mercurial > hg > octave-nkf

annotate libcruft/slatec-fn/d9lgmc.f @ 12992:0413227e407b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
dMatrix.cc: Fix tests introduced by changeset f217edac2c71. Declare variables Mt and rvt as shared, and increase tolerance for tests which fail on cygwin.
author Marco atzeri <marco.atzeri@gmail.com>
date Tue, 23 Aug 2011 20:58:54 -0400
parents 30c606bec7a8
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
1 *DECK D9LGMC
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
2 DOUBLE PRECISION FUNCTION D9LGMC (X)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
3 C***BEGIN PROLOGUE D9LGMC
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
4 C***SUBSIDIARY
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
5 C***PURPOSE Compute the log Gamma correction factor so that
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
6 C LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
7 C + D9LGMC(X).
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
8 C***LIBRARY SLATEC (FNLIB)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
9 C***CATEGORY C7E
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
10 C***TYPE DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
11 C***KEYWORDS COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
12 C LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
13 C***AUTHOR Fullerton, W., (LANL)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
14 C***DESCRIPTION
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
15 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
16 C Compute the log gamma correction factor for X .GE. 10. so that
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
17 C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
18 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
19 C Series for ALGM on the interval 0. to 1.00000E-02
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
20 C with weighted error 1.28E-31
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
21 C log weighted error 30.89
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
22 C significant figures required 29.81
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
23 C decimal places required 31.48
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
24 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
25 C***REFERENCES (NONE)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
26 C***ROUTINES CALLED D1MACH, DCSEVL, INITDS, XERMSG
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
27 C***REVISION HISTORY (YYMMDD)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
28 C 770601 DATE WRITTEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
29 C 890531 Changed all specific intrinsics to generic. (WRB)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
30 C 890531 REVISION DATE from Version 3.2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
31 C 891214 Prologue converted to Version 4.0 format. (BAB)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
32 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
33 C 900720 Routine changed from user-callable to subsidiary. (WRB)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
34 C***END PROLOGUE D9LGMC
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
35 DOUBLE PRECISION X, ALGMCS(15), XBIG, XMAX, DCSEVL, D1MACH
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
36 LOGICAL FIRST
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
37 SAVE ALGMCS, NALGM, XBIG, XMAX, FIRST
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
38 DATA ALGMCS( 1) / +.1666389480 4518632472 0572965082 2 D+0 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
39 DATA ALGMCS( 2) / -.1384948176 0675638407 3298605913 5 D-4 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
40 DATA ALGMCS( 3) / +.9810825646 9247294261 5717154748 7 D-8 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
41 DATA ALGMCS( 4) / -.1809129475 5724941942 6330626671 9 D-10 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
42 DATA ALGMCS( 5) / +.6221098041 8926052271 2601554341 6 D-13 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
43 DATA ALGMCS( 6) / -.3399615005 4177219443 0333059966 6 D-15 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
44 DATA ALGMCS( 7) / +.2683181998 4826987489 5753884666 6 D-17 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
45 DATA ALGMCS( 8) / -.2868042435 3346432841 4462239999 9 D-19 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
46 DATA ALGMCS( 9) / +.3962837061 0464348036 7930666666 6 D-21 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
47 DATA ALGMCS( 10) / -.6831888753 9857668701 1199999999 9 D-23 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
48 DATA ALGMCS( 11) / +.1429227355 9424981475 7333333333 3 D-24 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
49 DATA ALGMCS( 12) / -.3547598158 1010705471 9999999999 9 D-26 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
50 DATA ALGMCS( 13) / +.1025680058 0104709120 0000000000 0 D-27 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
51 DATA ALGMCS( 14) / -.3401102254 3167487999 9999999999 9 D-29 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
52 DATA ALGMCS( 15) / +.1276642195 6300629333 3333333333 3 D-30 /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
53 DATA FIRST /.TRUE./
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
54 C***FIRST EXECUTABLE STATEMENT D9LGMC
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
55 IF (FIRST) THEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
56 NALGM = INITDS (ALGMCS, 15, REAL(D1MACH(3)) )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
57 XBIG = 1.0D0/SQRT(D1MACH(3))
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
58 XMAX = EXP (MIN(LOG(D1MACH(2)/12.D0), -LOG(12.D0*D1MACH(1))))
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
59 ENDIF
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
60 FIRST = .FALSE.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
61 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
62 IF (X .LT. 10.D0) CALL XERMSG ('SLATEC', 'D9LGMC',
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
63 + 'X MUST BE GE 10', 1, 2)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
64 IF (X.GE.XMAX) GO TO 20
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
65 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
66 D9LGMC = 1.D0/(12.D0*X)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
67 IF (X.LT.XBIG) D9LGMC = DCSEVL (2.0D0*(10.D0/X)**2-1.D0, ALGMCS,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
68 1 NALGM) / X
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
69 RETURN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
70 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
71 20 D9LGMC = 0.D0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
72 CALL XERMSG ('SLATEC', 'D9LGMC', 'X SO BIG D9LGMC UNDERFLOWS', 2,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
73 + 1)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
74 RETURN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
75 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
76 END