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
         end if
   *
  -C      write (6,*) phase,m_max
  -      if (mod(abs(phase+m_max),2) .eq. 1) result = -result
  +      if (mod(abs(phase+m_max+ll/2),2) .eq. 1) result = -result
   *
         return
         end
  
positive.code:       implicit none
        integer l1,l2,l3,m1,m2,m3
        integer l_min,l_mid,l_max,m_min,m_mid,m_max,
       &    ll_min,ll_mid,ll_max,l_sum
        integer i,j,k,l,m,n,phase
        double precision result,wa,wb,wc
  *
        result = 0.0d00
        m3 = -(m1+m2)
  C      write (6,'(3I4)') l1,l2,l3
  C      write (6,'(3I4)') m1,m2,m3
        if (abs(m1) .gt. l1) return
        if (abs(m2) .gt. l2) return
        if (abs(m3) .gt. l3) return
  
negative.code: *   Calculate product of 3j symbols which appears in the angular       * 
  *   integral over three spherical tensors. Product is                  *
  *   (l1 l2 l3 / m1 m2 m3) (l1 l2 l3 / 0 0 0).                          *
  *                                                                      *
  *   l1,l2,l3 - l values                                                *
  *   m1,m2 - m values; m3 = -(m1+m2)                                    *
  