Blame view

fvn_fnlib/zlngam.f 2.65 KB
38581db0c   daniau   git-svn-id: https...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
        complex(8) function zlngam (zin)
        implicit none
  c august 1980 edition.  w. fullerton c3, los alamos scientific lab.
  c eventually clngam should make use of c8lgmc for all z except for
  c z in the vicinity of 1 and 2.
        complex(8) zin, z, corr, zlnrel, z9lgmc
        real(8) d1mach,pi,sq2pil,bound,dxrel,rmax,cabsz
        real(8) x,y,argsum,zarg
        integer irold,ir,i,n
        external z9lgmc, zarg, zlnrel,d1mach
        data pi / 3.1415926535 8979324e0 /
        data sq2pil / 0.9189385332 0467274e0 /
  c
        data bound, dxrel, rmax / 3*0.0 /
  c
        if (bound.ne.0.) go to 10
        n = -0.30*log(d1mach(3))
  c bound = n*(0.1*eps)**(-1/(2*n-1))/(pi*exp(1))
        bound = 0.1171*dble(n)*(0.1*d1mach(3))**(-1./(2.*dble(n)-1.))
        dxrel = sqrt (d1mach(4))
        rmax = d1mach(2)/log(d1mach(2))
  c
   10   z = zin
        x = real(zin)
        y = aimag(zin)
  c
        corr = (0.0, 0.0)
        cabsz = abs(z)
        if (cabsz.gt.rmax) call seteru (
       1  44hzlngam   abs(z) so large result may overflow, 44, 2, 2)
        if (x.ge.0.0 .and. cabsz.gt.bound) go to 50
        if (x.lt.0.0 .and. abs(y).gt.bound) go to 50
  c
        if (cabsz.lt.bound) go to 20
  c
  c use the reflection formula for real(z) negative, cabs(z) large, and
  c abs(aimag(y)) small.
  c
        call entsrc (irold, 1)
        if (y.gt.0.0) z = conjg (z)
        corr = exp (-dcmplx(0.0,2.0*pi)*z)
        if (real(corr).eq.1.0 .and. aimag(corr).eq.0.0) call seteru (
       1  31hzlngam  z is a negative integer, 31, 3, 2)
  c
        zlngam = sq2pil + 1.0 - dcmplx(0.0,pi)*(z-0.5) - zlnrel(-corr)
       1  + (z-0.5)*log(1.0-z) - z - z9lgmc(1.0-z)
        if (y.gt.0.0) zlngam = conjg (zlngam)
  c
        call erroff
        call entsrc (ir, irold)
        if (abs(y).gt.dxrel) return
  c  X in sequence field on next line is to preserve trailing blank
        if (0.5* abs((1.-corr)*zlngam/z).lt.dxrel) call seteru (68hzlngam    X
       1 answer lt half precision because z too near negative integer, 68,
       2 1, 1)
        return
  c
  c use the recursion relation for cabs(z) small.
  c
   20   if (x.ge.(-0.5) .or. abs(y).gt.dxrel) go to 30
        if ( abs((z-aint(x-0.5))/x).lt.dxrel) call seteru ( 68hzlngam  ans
       1wer lt half precision because z too near negative integer, 68,1,1)
  c
   30   n = sqrt (bound**2 - y**2) - x + 1.0
        argsum = 0.0
        corr = (1.0, 0.0)
        do 40 i=1,n
          argsum = argsum + zarg(z)
          corr = z*corr
          z = 1.0 + z
   40   continue
  c
        if (real(corr).eq.0.0 .and. aimag(corr).eq.0.0) call seteru (
       1  31hzlngam  z is a negative integer, 31, 3, 2)
        corr = -dcmplx (log(abs(corr)), argsum)
  c
  c use stirling-s approximation for large z.
  c
   50   zlngam = sq2pil + (z-0.5)*log(z) - z + corr + z9lgmc(z)
        return
  c
        end