anchor.diff:  C       Routine called: GAMMA2 for computing gamma function
   C       ======================================================
   C
  +C       DLMF 13.2.42 with prefactors rewritten according to
  +C       DLMF 5.5.3, M(a, b, x) with DLMF 13.2.2
  +C
           IMPLICIT DOUBLE PRECISION (A-H,O-Z)
           ID=-100
           PI=3.141592653589793D0
  
positive.code:         NM=0
          IF (IL1) NM=ABS(A)
          IF (IL2) NM=ABS(AA)
          IF (IL1.OR.IL2) THEN
             HU=1.0D0
             R=1.0D0
  
negative.code: 20      CONTINUE
  25      CALL GAMMA2(A,GA)
          HU1=HU1/GA
  C       DLMF 13.4.4 with substitution t=C/(1-u)
  C       integration u from 0 to 1, i.e. t from C=12/X to infinity
          DO 40 M=2,10,2
             HU2=0.0D0
             G=0.5D0/M
  