function psi (x)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension psics(23), apsics(16)
      external aint, alog, cot, csevl, inits, r1mach, sqrt
c
c series for psi        on the interval  0.          to  1.00000d+00
c                                        with weighted error   2.03e-17
c                                         log weighted error  16.69
c                               significant figures required  16.39
c                                    decimal places required  17.37
c
      data psi cs( 1) /   -.0380570808 35217922e0 /
      data psi cs( 2) /    .4914153930 2938713e0 /
      data psi cs( 3) /   -.0568157478 21244730e0 /
      data psi cs( 4) /    .0083578212 25914313e0 /
      data psi cs( 5) /   -.0013332328 57994342e0 /
      data psi cs( 6) /    .0002203132 87069308e0 /
      data psi cs( 7) /   -.0000370402 38178456e0 /
      data psi cs( 8) /    .0000062837 93654854e0 /
      data psi cs( 9) /   -.0000010712 63908506e0 /
      data psi cs(10) /    .0000001831 28394654e0 /
      data psi cs(11) /   -.0000000313 53509361e0 /
      data psi cs(12) /    .0000000053 72808776e0 /
      data psi cs(13) /   -.0000000009 21168141e0 /
      data psi cs(14) /    .0000000001 57981265e0 /
      data psi cs(15) /   -.0000000000 27098646e0 /
      data psi cs(16) /    .0000000000 04648722e0 /
      data psi cs(17) /   -.0000000000 00797527e0 /
      data psi cs(18) /    .0000000000 00136827e0 /
      data psi cs(19) /   -.0000000000 00023475e0 /
      data psi cs(20) /    .0000000000 00004027e0 /
      data psi cs(21) /   -.0000000000 00000691e0 /
      data psi cs(22) /    .0000000000 00000118e0 /
      data psi cs(23) /   -.0000000000 00000020e0 /
c
c series for apsi       on the interval  0.          to  2.50000d-01
c                                        with weighted error   5.54e-17
c                                         log weighted error  16.26
c                               significant figures required  14.42
c                                    decimal places required  16.86
c
      data apsics( 1) /   -.0204749044 678185e0 /
      data apsics( 2) /   -.0101801271 534859e0 /
      data apsics( 3) /    .0000559718 725387e0 /
      data apsics( 4) /   -.0000012917 176570e0 /
      data apsics( 5) /    .0000000572 858606e0 /
      data apsics( 6) /   -.0000000038 213539e0 /
      data apsics( 7) /    .0000000003 397434e0 /
      data apsics( 8) /   -.0000000000 374838e0 /
      data apsics( 9) /    .0000000000 048990e0 /
      data apsics(10) /   -.0000000000 007344e0 /
      data apsics(11) /    .0000000000 001233e0 /
      data apsics(12) /   -.0000000000 000228e0 /
      data apsics(13) /    .0000000000 000045e0 /
      data apsics(14) /   -.0000000000 000009e0 /
      data apsics(15) /    .0000000000 000002e0 /
      data apsics(16) /   -.0000000000 000000e0 /
c
      data pi     / 3.1415926535 8979324e0/
      data ntpsi, ntapsi, xbig, dxrel /0, 0, 0., 0./
c
      if (ntpsi.ne.0) go to 10
      ntpsi = inits (psics, 23, 0.1*r1mach(3))
      ntapsi = inits (apsics, 16, 0.1*r1mach(3))
c
      xbig = 1.0/sqrt(r1mach(3))
      dxrel = sqrt (r1mach(4))
c
 10   y = abs(x)
      if (y.ge.2.0) go to 30
c
c psi(x) for -2. .le. x .le. 2.
c
      n = x
      if (x.lt.0.) n = n - 1
      y = x - float(n)
      n = n - 1
      psi = csevl (2.*y-1., psics, ntpsi)
      if (n.eq.0) return
c
      n = -n
      if (x.eq.0.) call seteru (14hpsi     x is 0, 14, 2, 2)
      if (x.lt.0. .and. x+float(n-2).eq.0.) call seteru (
     1  31hpsi     x is a negative integer, 31, 3, 2)
      if (x.lt.(-0.5) .and. abs((x-aint(x-0.5))/x).lt.dxrel) call seteru
     1  (68hpsi     answer lt half precision because x too near negative
     1 integer, 68, 1, 1)
c
      do 20 i=1,n
        psi = psi - 1.0/(x+float(i-1))
 20   continue
      return
c
c psi(x) for abs(x) .ge. 2.
c
 30   aux = 0.
      if (y.lt.xbig) aux = csevl (8./y**2-1., apsics, ntapsi)
      if (x.lt.0.) psi = alog(abs(x)) - 0.5/x + aux - pi*cot(pi*x)
      if (x.gt.0.) psi = alog(x) - 0.5/x + aux
      return
c
      end