$b, c$ are $n$-dimensional column vectors. $A$ is a $n \times n$ matrix. $D(z)$ is a $n \times n$ diagonal matrix, whose terms are powers of z. $d$ is a constant. The goal is to find the zeros and poles of the formula $c^T[D(z^{-1}) - A]^{-1}b + d$.
The equation for zeros is given as $\det[A - \frac1dbc^T - D(z^{-1})] = 0$. The equation for poles is given as $\det[A - D(z^{-1})] = 0$.
Note that the formula is one-dimensional. So if any valid one-dimensional equation is given, it has a chance of being ok (thus justifying the determinant operation).
I have not been able to produce these equations. It seems to me that the authors are trying to "invert" $b$ and $c$. Then, to make the dimensions match, they take a determinant operation.
However, I can plug in numbers to their equations. We'll focus on the equation for zeros. Here's an instance in which it works: let $b = c = (1, 1)$, $A = \left(\begin{matrix}1 & 0 \\0 & 2 \end{matrix}\right)$, $D(z) = \left(\begin{matrix}z^{-1} & 0 \\0 & z^{-2} \end{matrix}\right)$, $d = 1$. Their formula gives $\det\left(\begin{matrix}-z & -1 \\-1 & 1 - z^{2} \end{matrix}\right) = 0$, whose solutions are here. These results agree with the correct solution, which is $\frac1{z-1} + \frac1{z^2 - 2} + 1 = 0$, and whose solutions are here.
Some numbers don't work. Let $b= (0, 0)$, then their equation for zeros gives strangely placed zeros. Same if $c=(0, 0)$. Also, $d = 0$ makes it fail.
For a more complex example of failure, let $b = c = (1, 1)$, $A = I$, $D(z) = \left(\begin{matrix}z^{-1} & 0 \\0 & z^{-1} \end{matrix}\right)$. Their formula gives $z = \pm 1$, but only $z = -1$ is correct.
So why does this formula sometimes work and sometimes not? How was it derived?
These equations come from equations 4 and 5 in this paper, which are duplicated from equations 16 and 19 from this paywalled paper.
The examples you gave are all cases where there is pole-zero-cancellation, so those zeros are there but there is also a pole at the exact same location so cancel each other out when simplifying the formula.
I assume that the equation for poles is clear, since the inverse of $D(z^{-1}) - A$ is used in the formula. The equation for the zeros can be derived by starting with the state space equation
\begin{align} D(z^{-1})\,x &= A\,x + b\,u, \tag{1a} \\ y &= c^\top x + d\,u. \tag{1b} \end{align}
Namely solving $(\text{1a})$ for $x$ and substituting that into $(\text{1b})$ yields
$$ y = \left(c^\top [D(z^{-1}) - A]^{-1} b + d\right) u, \tag{2} $$
which contains the formula you started with. Now a zero occurs when $y=0$ for a non-zero $x,u$, which can be written as
$$ \begin{bmatrix} 0 \\ y \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} = \underbrace{ \begin{bmatrix} A - D(z^{-1}) & b \\ c^\top & d \end{bmatrix} }_{M(z)} \begin{bmatrix} x \\ u \end{bmatrix}. \tag{3} $$
There is a non-zero $x,u$ which solves this if $M(z)$ is singular and would have a nullspace, so when $\det(M(z))=0$. By now using the Schur complement this can also be written as
$$ \det(M(z)) = \det(d) \det(A - D(z^{-1}) - b\,d^{-1}\,c^\top) = d \det(A - \tfrac{1}{d}b\,c^\top - D(z^{-1})). \tag{4} $$
You can only use this Schur complement when $d\neq0$. If $d=0$ you can still use $\det(M(z))=0$, but this might also add solutions for $z$ at infinity.