Commutation of bosonic operators on finite Hilbert space

I understand the common commutation relations for creationd and annihilation operators, given by: $$[b_i, b_j^\dagger]|n_1 n_2\dots n_N\rangle = b_ib_j^\dagger|n_1 n_2\dots n_N\rangle-b_j^\dagger b_i|n_1 n_2\dots n_N\rangle\\ = \sqrt{n_i(n_j+1)}|n_1\dots n_i-1\ n_j+1\dots n_N\rangle\\ -\sqrt{n_i(n_j+1)}|n_1\dots n_i-1\ n_j+1\dots n_N\rangle=0$$ My question is about equation (2.57) on page 32. They restrict the Fock space to contain at most $$N_P$$ bosons per site, and then they affirm that the following relations hold: $$\left[\bar{b}_{i}, \bar{b}_{j}\right]=0,\left[\bar{b}_{i}, \bar{b}_{j}^{\dagger}\right]=\delta_{i j}\left[1-\frac{N_{P}+1}{N_{P} !}\left(\bar{b}_{i}^{\dagger}\right)^{N_{P}}\left(\bar{b}_{i}\right)^{N_{P}}\right]$$ Where does that come from? I started considering the general Fock state $$|n_1 n_2\dots n_N\rangle=\prod_{k=1}^{N}\frac{1}{\sqrt{n_k!}}{b_k^\dagger}^{n_k}|0\rangle$$ and trying to compute the commutaor expression but obtained nothing useful. Can anyone help me?
You can start with the exact matrix expression for the ladder operators, which we use for numerics, $$b^\dagger = \sum_{i=1}^N \sqrt{i} \mid i \rangle \langle i -1 \mid \quad b = \sum_{i=1}^N \sqrt{i} \mid i - 1 \rangle \langle i \mid,$$ where you will find that these two commute to $$1$$ for any state other than the final state $$\mid N \rangle$$, for which they yield $$[b,b^\dagger]\mid N\rangle = -N \mid N \rangle$$. The main point is that the second term is a kronecker delta function $$\delta_{\hat{n},N} = \frac{(b^\dagger)^N b^N}{N!}$$ because the ladder operators will destroy any state other than the highest allowed state, which is $$\mid N \rangle$$. In order to prove these things analytically it might help to know the following identity $$(b^\dagger)^m b^m \mid n \rangle = \frac{n!}{(n-m)!} \mid n \rangle,$$ and to remember that $$(b^\dagger)^{N+1} = b^{N+1} = 0$$.