Relation between Diophantine equations involving multiple sums with the modulo function and greatest common divisor function.

Question

This one of my questions that I have as I am trying to learn about Diophantine equations, and about what I found through an OEIS search here: https://oeis.org/A360323

To generalize:

Let the function $f$ be the function:

$$f(x,y,z,...)=a_1 x + a_2 x^2+a_3 x^3+ a_4 x^4 +...a_k x^k + b_1 y + b_2 y^2+b_3 y^3+ b_4 y^4 +...b_k y^k+ c_1 z + c_2 z^2+c_3 z^3+ c_4 z^4 +...c_k z^k,...$$

where the coefficients:

$$a_1,a_2,a_3,a_4,...,a_k$$ $$b_1,b_2,b_3,b_4,...,b_k$$ $$b_1,b_2,b_3,b_4,...,b_k$$ are random integers in some chosen interval.

Then count $\text{#}$ the numbers of solutions to:

$$f(x,y,z,...) \bmod p_n = 0$$, where the variables $x,y,z,...$ ranges between: $$x=1..p_n$$ $$y=1..p_n$$ $$z=1..p_n$$ $$\vdots$$

where $p_n$ is the $n$-th prime.

Subtract $$(p_n)^{w}-\text{#}$$

where

$w$ = the number of variables $x,y,z,...$ used in $f(x,y,z,...)$

and compare the following $\left[f \bmod p_n=0\right]$ and $\left[\gcd(f,p_n)=1\right]$ sums, where $\left[ \text{ } \right]$ is the Iverson bracket.

$$(p_n)^w-\sum _{ö=1}^{p_n} \cdots \left(\sum _{z=1}^{p_n} \left(\sum _{x=1}^{p_n} \left(\sum _{y=1}^{p_n} \left[f \bmod p_n=0\right]\right)\right)\right)$$ $$= \sum _{ö=1}^{p_n} \cdots \left(\sum _{z=1}^{p_n} \left(\sum _{y=1}^{p_n} \left(\sum _{x=1}^{p_n} \left[\gcd(f,p_n)=1\right]\right)\right)\right) \tag{$\ast$}$$

Question:

Can you prove this conjectured relation $(\ast)$? It seems a bit tautological since the modulo function and greatest common divisor function are related. But the $(p_n)^w$ has the variable $w$ in it.

Hopefully nothing is wrong in the latex. Here is however a Mathematica program that demonstrates what I am trying to do:

"Mathematica start"
Clear[f, nn, x, y, z, a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4]
nn = 10;(*nn=10 is the number of terms in the sequences*)imax = 4;

"w = number of variables x,y,z,..."
w = 3

"a1,a2,a3,a4"
a1 = RandomInteger[{-imax, imax}];
a2 = RandomInteger[{-imax, imax}];
a3 = RandomInteger[{-imax, imax}];
a4 = RandomInteger[{-imax, imax}];

"b1,b2,b3,b4"
b1 = RandomInteger[{-imax, imax}];
b2 = RandomInteger[{-imax, imax}];
b3 = RandomInteger[{-imax, imax}];
b4 = RandomInteger[{-imax, imax}];

"c1,c2,c3,c4"
c1 = RandomInteger[{-imax, imax}];
c2 = RandomInteger[{-imax, imax}];
c3 = RandomInteger[{-imax, imax}];
c4 = RandomInteger[{-imax, imax}];

f = a1*x + a2*x^2 + a3*x^3 + a4*x^4 + b1*y + b2*y^2 + b3*y^3 + 
  b4*y^4 + c1*z + c2*z^2 + c3*z^3 + c4*z^4
Table[(Prime[n]^w - 
   Sum[Sum[Sum[If[Mod[f, Prime[n]] == 0, 1, 0], {y, 1, Prime[n]}], {x,
       1, Prime[n]}], {z, 1, Prime[n]}]), {n, 1, nn}]
Table[Sum[
  Sum[Sum[If[GCD[f, Prime[n]] == 1, 1, 0], {x, 1, Prime[n]}], {y, 1, 
    Prime[n]}], {z, 1, Prime[n]}], {n, 1, nn}]
%% - %
"The difference between the left and the right hand side appear to \
zero."
"Mathematica end"