If I understood your question right, you would like to obtain a Linear Constant Coefficient Difference Equation (LCCDE) representation of a given system from its given System Transfer Function, $H(z)$, when it's expressed in a pole-zero product form such as:
$$
H(z) = g \frac {(z - d_0)(z - d_1)} {(z - c_0)(z - c_1)}
\tag{1}
$$
instead of a, more directly apparent, power series form such as:
$$
H(z) = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 + a_1 z^{-1} + a_2 z^{-2}}
\tag{2}
$$
As known, the latter form can be directly converted to its corresponding LCCDE representation from the coefficients ${b_k}$ and ${a_k}$ as follows:
$$y[n] + a_1 y[n-1] + a_2 y[n-2] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2]$$
Where the pattern of conversion speaks for itself, without a need of further description.
Now, we'd like to know whether one can obtain a similar conversion from the pole-zero form into the LCCDE without expanding the product?
Well yes and no! Depending on how you consider the process of computing the coefficents. If your intention is on the fact that you do not want to algebricaly perform the paranthesis expansion operation, but instead use a machanism that computes each coefficient separately by a set of arithmetic operations; then yes you can devise such a method to, for example, compute $a_3$ or $b_5$ without computing any other terms. And No, eventually you must always compute the coefficient $a_k$ that multiplies $z^{-k}$ so as to find the coefficient that multiplies $y[n-k]$, hence you need to always expand the brackets and compute multipliers of $z^{-k}$ in one or the other way.
Then I want to present here a simple method to find those coefficients :
given a polynomial in the product of zeros form: $$p(x) = (x-a)(x-b)(x-c)$$ we can perform the following operation to get its coefficients in the expanded form $$p(x) = d_0 x^3 + d_1 x^2 + d_2 x + d_3$$
$$d[k] = (1 -a)\star(1 -b)\star(1-c)$$ where $k =0,1,2,3$ in this particular case.
As clearly seen this is a discrete convolution operation, where each multiplier bracket is turned into a respective convolution operand. For which a simple MATLAB line to compute $d[k]$ would be:
d = conv([1 -a],conv([1 -b],[1 -c]));
where a
, b
, and c
would be replaced by their numeric values or if you're competent in the symbolic math capability, you could instead obtain a pure symbolic result for the coeffients $d_k$, which would therefore give you a way to compute LCCDE coefficients from pole-zero products of $H(z)$ in a rather roundabout way.
Finally, the computations for a power of $z^{-1}$ products is basicly the same if you describe the product as $$P(z) = (1 - az^{-1})(1-bz^{-1})(1-cz^{-1})$$ and,
$$P(z) = d_0 + d_1 z^{-1} + d_2 z^{-2} + d_3 z^{-3} $$
where again $$d[k] = (1 -a)\star(1 -b)\star(1-c)$$