Let $(x_1,y_1,x_2,y_2)$ satisfy the divisibility condition, then clearly $p \not| \gcd(x_1,x_2) $ (we call such pairs good) .
We have $ 2(p^n - p^{n-1})p^n= 2p^{2n-1}(p-1)$ such good pairs $(x_1,x_2)$. We say $(c,d)$ is associated to $(a,b)$ if $p^n | ac + bd +1$. Letting $(a,b)$ be a good pair, we count its associated pairs. We may assume by symmetry that $a$ is not divisible by $p$, then letting $u$ be the inverse of $a$ module $p^n$. Then if $(c,d)$ is an associated pair, we have $ c \equiv -u(bd+1) \pmod{p^n}$ and hence $c$ is uniquely determined by $d$ as each module class intersect $\{0,1, \cdots , p^n-1\}$ at exactly one element. So there is exactly $p^n$ associated pairs for every good pair, whence there is $2p^{3n-1}(p-1)$ quadruples $(a_1,a_2,a_3,a_4)$ satisfying the \[ p^n\mid (a_1a_2+a_3a_4+1). \]