function poch1 (a, x) c august 1980 edition. w. fullerton, c3, los alamos scientific lab. c c evaluate a generalization of pochhammer-s symbol for special c situations that require especially accurate values when x is small in c poch1(a,x) = (poch(a,x)-1)/x c = (gamma(a+x)/gamma(a) - 1.0)/x . c this specification is particularly suited for stably computing c expressions such as c (gamma(a+x)/gamma(a) - gamma(b+x)/gamma(b))/x c = poch1(a,x) - poch1(b,x) c note that poch1(a,0.0) = psi(a) c c when abs(x) is so small that substantial cancellation will occur if c the straightforward formula is used, we use an expansion due c to fields and discussed by y. l. luke, the special functions and their c approximations, vol. 1, academic press, 1969, page 34. c c the ratio poch(a,x) = gamma(a+x)/gamma(a) is written by luke as c (a+(x-1)/2)**x * polynomial in (a+(x-1)/2)**(-2) . c in order to maintain significance in poch1, we write for positive a c (a+(x-1)/2)**x = exp(x*alog(a+(x-1)/2)) = exp(q) c = 1.0 + q*exprel(q) . c likewise the polynomial is written c poly = 1.0 + x*poly1(a,x) . c thus, c poch1(a,x) = (poch(a,x) - 1) / x c = exprel(q)*(q/x + q*poly1(a,x)) + poly1(a,x) c dimension bern(9), gbern(10) external cot, exprel, poch, psi, r1mach c c bern(i) is the 2*i bernoulli number divided by factorial(2*i). data bern( 1) / .8333333333 3333333e-01 / data bern( 2) / -.1388888888 8888889e-02 / data bern( 3) / .3306878306 8783069e-04 / data bern( 4) / -.8267195767 1957672e-06 / data bern( 5) / .2087675698 7868099e-07 / data bern( 6) / -.5284190138 6874932e-09 / data bern( 7) / .1338253653 0684679e-10 / data bern( 8) / -.3389680296 3225829e-12 / data bern( 9) / .8586062056 2778446e-14 / c data pi / 3.1415926535 8979324 e0 / data sqtbig, alneps / 2*0.0 / c if (sqtbig.ne.0.0) go to 10 sqtbig = 1.0/sqrt(24.0*r1mach(1)) alneps = alog(r1mach(3)) c 10 if (x.eq.0.0) poch1 = psi(a) if (x.eq.0.0) return c absx = abs(x) absa = abs(a) if (absx.gt.0.1*absa) go to 70 if (absx*alog(amax1(absa,2.0)).gt.0.1) go to 70 c bp = a if (a.lt.(-0.5)) bp = 1.0 - a - x incr = 0 if (bp.lt.10.0) incr = 11.0 - bp b = bp + float(incr) c var = b + 0.5*(x-1.0) alnvar = alog(var) q = x*alnvar c poly1 = 0.0 if (var.ge.sqtbig) go to 40 var2 = (1.0/var)**2 c rho = 0.5*(x+1.0) gbern(1) = 1.0 gbern(2) = -rho/12.0 term = var2 poly1 = gbern(2)*term c nterms = -0.5*alneps/alnvar + 1.0 if (nterms.gt.9) call seteru ( 1 49hpoch1 nterms is too big, maybe r1mach(3) is bad, 49, 1, 2) if (nterms.lt.2) go to 40 c do 30 k=2,nterms gbk = 0.0 do 20 j=1,k ndx = k - j + 1 gbk = gbk + bern(ndx)*gbern(j) 20 continue gbern(k+1) = -rho*gbk/float(k) c term = term * (float(2*k-2)-x)*(float(2*k-1)-x)*var2 poly1 = poly1 + gbern(k+1)*term 30 continue c 40 poly1 = (x-1.0)*poly1 poch1 = exprel(q)*(alnvar + q*poly1) + poly1 c if (incr.eq.0) go to 60 c c we have poch1(b,x). but bp is small, so we use backwards recursion c to obtain poch1(bp,x). c do 50 ii=1,incr i = incr - ii binv = 1.0/(bp+float(i)) poch1 = (poch1-binv)/(1.0+x*binv) 50 continue c 60 if (bp.eq.a) return c c we have poch1(bp,x), but a is lt -0.5. we therefore use a reflection c formula to obtain poch1(a,x). c sinpxx = sin(pi*x)/x sinpx2 = sin(0.5*pi*x) trig = sinpxx*cot(pi*b) - 2.0*sinpx2*(sinpx2/x) c poch1 = trig + (1.0 + x*trig) * poch1 return c 70 call entsrc (irold, 1) poch1 = poch (a, x) if (poch1.eq.0.0) call erroff call entsrc (irold2, irold) c poch1 = (poch1 - 1.0) / x return c end