anchor.diff:        double precision xq(100),wq(100)
   
         u=(a+b)/(2.d0*dsqrt(gam))
  -      factor=dexp(u**2-arg)/dsqrt(gam)
  +      factor=dexp(u*u-arg)/dsqrt(gam)
   
   xq(1)= 5.38748089001123    
   xq(2)= 4.60368244955074   
  
positive.code: 
          !Compute the s_q+1_0
  !        s_q_0=s_q_0*(2.d0*q+nlm+1)*b**2/((2.d0*(m+q)+3)*4.d0*(q+1)*gam)
          s_q_0=s_q_0*(2.d0*q+nlm+1)*b*b/((8.d0*(m+q)+12.d0)*(q+1)*gam)
          
          if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi*.5d0)
          ! Increment q
  
negative.code: double precision g,dble_fact,expo
  double precision, parameter :: sq_pi_ov_2=dsqrt(dacos(-1.d0)*0.5d0)
  expo=0.5d0*dfloat(n+1)
  crochet=dble_fact(n-1)/(g+g)**expo
  if(mod(n,2).eq.0)crochet=crochet*sq_pi_ov_2
  end
  
  