Interpolate vertex attributes with $z$ AFTER homogeneous divide Usually, if we want to interpolate some vertex attributes, say $A_1$, $A_2$ and $A_3$, we have to first divide them by their corresponding $z$ in view space, then we interpolate them using barycentric weights and finally we multiply them by the interpolated $z$ at the current pixel. To interpolate this $z$, we have to inverse the 3 $z$ of each vertex, then interpolate them and finally inverse the result to have the interpolated $z$.
In other words, the interpolated attribute $A$ is equal to:
$$
A = z \left(w_0 \times \frac{A_1}{Z_1} + w_1 \times \frac{A_2}{Z_2} + w_2 \times \frac{A_3}{Z_3}\right)
$$
where $w_i$ are the barycentric weights, $Z_i$ are the depths of each vertex and $z$ is the interpolated depth:
$$
z = \frac{1}{w_0 \times \frac{1}{Z_1} + w_1 \times \frac{1}{Z_2} + w_2 \times \frac{1}{Z_3}}
$$
Now my question is: instead of going back and forth between attributes and their reciprocals, can we just use the vertices'$z$ in normalized device coordinates, that is, the $z$s that have been divided by $w$, after the perspective projection?
If I understand correctly, we interpolate the reciprocal of vertex attributes, because it's linear in screen space, but what about the $z$ in NDC space, after the homogeneous divide? If we can manage to use it, we already have a $z$ that we can linearly interpolate, because it has been homogenized?
In the book &quot;Game Engine Architecture&quot;, page 667, they say:

with w-buffering, we cannot linearly interpolate depths directly. Depths must be inverted prior to interpolation and the re-inverted prior to being stored the w-buffer.

In other words, we don't need to invert $z$ back and forth to find the interpolated depth: we just have to use the $z$ in NDC space, that has been divided by view-space $z$. So we can directly interpolate the $z$ in NDC space to have perspective-correct depth, but how can I use this interpolated depth to interpolate other vertex attributes?
The $Z$ of NDC space is related to $1/Z_\text{view}$ but not the same. With a typical projection matrix,  they're related by an affine remapping,
$$
Z_\text{NDC} = a \, \frac{1}{Z_\text{view}} + b
$$
where $a, b$ are some constants related to the near and far planes. So, in general you wouldn't be able to substitute $Z_\text{NDC}$ for $1/Z_\text{view}$ in the interpolation expressions. The $a$ factor isn't a problem, as it would factor out and cancel, but if $b \neq 0$ then it will produce some extra terms and throw off the interpolation.
If I'm hearing you right, it sounds like you're concerned about the idea of needing to prepare every single vertex attribute value (multiplying by $1/Z_\text{view}$) before interpolation and then reverse that afterward. Well, there is an easier way: we can create perspective-corrected barycentric coordinates $w'_0, w'_1, w'_2$, defined by
$$
w'_i \equiv \frac{w_i / Z_i}{w_0/Z_0 + w_1/Z_1 + w_2/Z_2}
$$
Then we can interpolate any attributes by $w'_0 A_0 + w'_1 A_1 + w'_2 A_2$, and it will be perspective-correct interpolation. You can verify this is equivalent to your original interpolation formula in your question. If we have more than three scalar attributes to interpolate (which we usually do), this would be the cheaper way to do it.