Where comes the +1 from in this formula? I'm working on two papers ([1] equation 8, [2] equation 2.3) and I can't figure out why there is an identic formula I can't explain on both.
$p(z) \in \mathbb{C}[z]$ is a monic polynomial with simple roots, degree $n$ and complex coefficients.
I fix n distinct points $z_1, z_2, \dots, z_n$ (not necessarily the roots of $p$). The Lagrange's basis for these points are the polynomials
$$L_i(z) := \prod_{j=1, j \neq i}^n \frac{z - z_j}{z_i - z_j} \qquad i=1, \dots, n$$
so I can write $q$, an interpolating polynomial for $p$ (this is the only linear combination of the $L_i$'s that takes the same values as $p$ on the points $z_i$):
$$q(z) = \sum_{i=1}^n p(z_i) L_i(z) = \sum_{i=1}^n p(z_i) \prod_{j=1, j \neq i}^n \frac{z - z_j}{z_i - z_j}$$
(As a side question: in both the papers it is stated that this is exactly $p$, namely that $q=p$. Why does it have to be?)
By the way, here comes my real question: defining the Weierstrass' corrections
$$W_i(z) := \frac{p(z)}{ \prod_{j=1, j \neq i}^n (z - z_j)} \quad , \quad W_i := W_i(z_i) = \frac{p(z_i)}{ \prod_{j=1, j \neq i}^n (z_i - z_j)}$$
I calculated (in many different ways, the following is the simplest I found)
$$\sum_{i=1}^n p(z_i) \prod_{j=1, j \neq i}^n \frac{z - z_j}{z_i - z_j} = \sum_{i=1}^n p(z_i) \frac{ \prod_{j=1, j \neq i}^n (z - z_j) } { \prod_{j=1, j \neq i}^n (z_i - z_j) } = \sum_{i=1}^n \frac{p(z_i)}{z-z_i} \cdot \frac{ \prod_{j=1}^n (z - z_j) } { \prod_{j=1, j \neq i}^n (z_i - z_j) } = \\
= \left( \prod_{j=1}^n (z - z_j) \right) \cdot \left( \sum_{i=1}^n \frac{1}{z - z_i} \cdot \frac{p(z_i)}{\prod_{j=1, j \neq i}^n (z_i - z_j)} \right) = \left( \prod_{j=1}^n (z - z_j) \right) \cdot \left( \sum_{i=1}^n \frac{W_i}{z - z_i} \right)$$
so that isolating the $i$th term of the summation I get
$$\frac{W_i}{z - z_i} \cdot \prod_{j=1}^n (z - z_j) + \left( \sum_{j=1, j \neq i}^n \frac{W_j}{z - z_j} \right) \cdot \prod_{j=1}^n (z - z_j) = W_i \cdot \prod_{j=1, j \neq i}^n (z - z_j) + \left( \sum_{j=1, j \neq i}^n \frac{W_j}{z - z_j} \right) \cdot \prod_{j=1}^n (z - z_j)$$
while the authors in both papers get
$$W_i \prod_{j \neq i} (z - z_j) + \prod_{j=1}^n (z-z_j) \left( \sum_{j \neq i} \frac{W_j}{z - z_j} + 1 \right)$$
After one entire day of headache I can't still figure out where that $+1$ comes from. Any help will be really, really appreciated. Thanks in advance.
References
[1]On a simultaneous method of Newton–Weierstrass’ type for finding
all zeros of a polynomial - M.S. Petkovic ́, Ð. Herceg,  I. Petkovic ́- Applied Mathematics and Computation 215 (2009) 2456–2463 (doi:10.1016/j.amc.2009.08.048)
[2]POINT ESTIMATION OF CUBICALLY CONVERGENT ROOT FINDING METHOD OF WEIERSTRASS’ TYPE - Lidija Z. Ranˇ ci ́c - FACTA UNIVERSITATIS (NIS) Ser. Math. Inform. Vol. 28, No 4 (2013), 417–428 - http://facta.junis.ni.ac.rs/mai/mai2804/fumi2804_07.pdf
 A: I finally found how the formula with the $+1$ is correct. I've done some mistakes in my reasoning above.
I moved backward, starting from the final form where for each $i = 1, \dots, n$
$$W_i \prod_{j \neq i} (z - z_j) + \prod_{j=1}^n (z-z_j) \left( \sum_{j \neq i} \frac{W_j}{z - z_j} + 1 \right)$$
and showing it is , in fact, the same as $p(z)$.
We already know (see question above) that this formula equals
$$\sum_{i=1}^n p(z_i) L_i(z) + \prod_{i=1}^n (z - z_i) =: \tilde{p}(z)$$
which first summand belongs to a subspace $V \subseteq \mathbb{C} [z]_{\leq n-1}$ (namely $V$ has the $L_i$s as a basis) and this is precisely why this summand alone cannot equal $p$: it has degree $n-1$ while $\deg{p} = n$.
The degree $n$ part of $p$ comes from the right summand of $\tilde{p}$. The key fact here is that $p$ is monic, so it can be represented as an element in a coset of $V$, whose representative is the right summand above[1].
To see that $p = \tilde{p}$ as polynomials observe they are both monic of degree $n$ and they take the same value on the $n$ distinct points $z_1, \dots, z_n$.
[1] (using standard Lagrange interpolation on $n$ distinct points one approximates a function using a degree $n-1$ polynomial, in this way the functions we can approximate exactly are polynomials of degree $n-1$ at most; here the function we interpolate is a polynomial forced to have $1$ as coefficient for $z^n$ so we can hope to interpolate it exactly traslating the space spanned by our interpolating polynomials without adding another interpolation point nor raising the degree of the Lagrange's polynomials).
