)], \label{eq:ajl}$$ with $N_s$ being the number of sites and $$\begin{aligned}
& & A(j_1, j_2,..., j_{N_e/2};l_1, l_2,..., l_{N_e/2}) = \nonumber  \\
& &
\left|
   \begin{array}{cccc}
       a(j_1,l_1) & a(j_1,l_2)&\ldots & a(j_1,l_{N_e/2}) \\
       a(j_2,l_1) & a(j_2,l_2)&\ldots & a(j_2,l_{N_e/2}) \\
       \multicolumn{4}{c}{\dotfill}\\
       a(j_{N_e/2},l_1) & a(j_{N_e/2},l_2)&\ldots & a(j_{N_e/2},l_{N_e/2}) \label{eq:det} \\
   \end{array}
\right| ;\ \end{aligned}$$ $j_1,j_2, \ldots$ are arranged in the ascending order; so are $l_1, l_2, \ldots$. Then the ground state energy $$E_g = < H > \equiv <\Psi_s| H |\Psi_s> / <\Psi_s|\Psi_s> \label{eq:average}$$ is obtained using a M.C. procedure \[\]. In order to minimize computation time in the M.C. computation the values of the cofactors of the matrix in eq. (\[eq:det\]) were stored and corrected at each time when the electron distribution was modified. Errors accumulated in the cofactors after corrected many times were avoided by recalculating cofactors as determinants, instead of getting them by correction, after a certain number of M.C. steps. We tested our programs using exact diagonalization results for small systems. We optimized $E_g$ with respect to $g$, $\Delta_{\mib k} $ and $\mu$. In the later stage of work we employed the correlated measurements method \[\] in the process of searching optimal parameter values minimizing $E_g$. It enables to precisely locate 