anchor.diff:  
   contains
   
  -  subroutine compute_multiplicities(counts, nc_)
  +  subroutine initialize_dirichlet(counts, nc_)
  +    ! set alphabet_size, n_data
       ! set n_multi, multi_z, multi
       integer(int32), intent(in) :: counts(:)
       integer(int32), intent(in) :: nc_
  +    alphabet_size = nc_
  +    n_data = sum(counts)
  +    call compute_multiplicities(counts)
  +  end subroutine initialize_dirichlet
  +
  +  subroutine compute_multiplicities(counts)
  +    ! set n_multi, multi_z, multi
  +    integer(int32), intent(in) :: counts(:)
       integer(int32)              :: nbins
       integer(int32)              :: i_,k_,ni_
       integer(int32)              :: err
       integer(int32)              :: nmax
       integer(int32), allocatable :: multi0(:)
  -
  -    alphabet_size = nc_
  -    n_data = sum(counts)
   
       ! compute multiplicities 
       ! nmax is the largest number of samples in a bin
  
positive.code:      return
    end if
  
    call compute_multiplicities(counts, nc)
  
    call compute_integration_range()
  
negative.code:      return
    end if
  
    call initialize_dirichlet(counts, nc)
    call compute_multiplicities(counts, nc)
  
    estimate = hdir(alpha)
  