anchor.diff:        if(N_states > 1) stop "mrcc_stoch : N_states == 1"
         do i=1, int(cps_N(cur_cp))
           !if(.not. actually_computed(i)) stop "not computed"
  -        call get_comb_val(comb(i), mrcc_detail, cp_first_tooth(cur_cp), val)
  +        !call get_comb_val(comb(i), mrcc_detail, cp_first_tooth(cur_cp), val)
  +        call get_comb_val(comb(i), mrcc_detail, cur_cp, val)
           !val = mrcc_detail(1, i) * mrcc_weight_inv(i) * comb_step
           su += val ! cps(i, cur_cp) * val
           su2 += val**2 ! cps(i, cur_cp) * val**2
  
positive.code:   logical, allocatable :: actually_computed(:)
    integer :: total_computed
    
    allocate(delta_det(N_states, N_det_non_ref, 0:comb_teeth, 2))
    allocate(cp(N_states, N_det_non_ref, N_cp, 2), mrcc_detail(N_states, N_det_generators))
    allocate(delta_loc(N_states, N_det_non_ref, 2))
    mrcc_detail = 0d0
  
negative.code: 
    curs = 1d0 - stato
    val = 0d0
    first = cp_first_tooth(cur_cp) 
  
    do j = comb_teeth, first, -1
      !DIR$ FORCEINLINE
      k = mrcc_find(curs, mrcc_cweight,size(mrcc_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
      !if(k < first_det_of_teeth(first)) exit
      if(k == first_det_of_teeth(first)) then
       val += detail(1, k) * mrcc_weight_inv(k) * comb_step * fractage(first)
      else
       val += detail(1, k) * mrcc_weight_inv(k) * comb_step
      end if
  
      curs -= comb_step
    end do
  end subroutine
  