anchor.diff:          wa = m
           wb = n
           result = result*sqrt(wa/wb)
  -C        write (6,*) wa,wb,result
   *
  -        if (l_max .ne. l_min+l_mid) then
  +*     Remaining factor
  +*
  +        ll = l_min+l_mid-l_max
  +        if (ll .gt. 0) then
             i = l_min+m_min
             j = l_min-m_min
  -          k = l_mid+m_mid
  -          l = l_mid-m_mid
  -          wc = i*j*k*l-(i*(i-1)*l*(l-1)+j*(j-1)*k*(k-1))/2
  -          result = result*wc
  -C          write (6,*) wc,result
  +          m = l_mid+m_mid
  +          n = l_mid-m_mid
  +          ll_min = max(0,ll-i,ll-n)
  +          ll_max = min(ll,j,m)
  +          k = 1
  +          do l = 1,ll_min
  +            k = -k*j*m/l
  +            j = j-1
  +            m = m-1
  +          end do
  +          do l = 1,ll-ll_min
  +            k = k*i*n/l
  +            i = i-1
  +            n = n-1
  +          end do
  +          l_sum = k
  +          do l = ll_min+1,ll_max
  +            i = i+1
  +            n = n+1
  +            k = -k*j*m*(ll-l+1)/(i*n*l)
  +            j = j-1
  +            m = m-1
  +            l_sum = l_sum+k
  +          end do
  +          wa = l_sum
  +          result = result*wa
           end if
  +	phase = phase+ll/2
         end if
   *
  -C      write (6,*) phase,m_max
         if (mod(abs(phase+m_max),2) .eq. 1) result = -result
   *
         return
  
positive.code:         do l = 3,l_sum-ll_mid-1,2
            m = m*l
          end do
          wb = m
          result = wb/wa
  C        write (6,*) wb,wa,result
  *
          m = 1
          n = 1
  
negative.code: C$Id: ecp_3j_prod.f,v 1.5 1997-05-22 23:52:50 mg201 Exp $
  ************************************************************************
  *                                                                      *
        subroutine ecp_3j_prod (l1,l2,l3,m1,m2,result)
  