solving a $p$-adic integration involving maximum function I am always struggling when it comes to dealing with maximum functions. I am trying to find the solution to this integral
$$\int_{\mathbb{Z}_p^3}||xy,xz,yz||_p^sd\mu(x,y,z),$$
where $||xy,xz,yz||_p^s=\sup\{|xy|_p^s,|xz|_p^s,|yz|_p^s\}$ and $d\mu$ is the normalized Haar measure.
Note: this is the same $p$-adic integration in Igusa local zeta function.
Ex: $\int_{\mathbb{z}_p}|x|_p^sd\mu(x)=\frac{1-p^{-1}}{1-p^{-s-1}}$
 A: I'll solve this by a similar method I use in your other question: problem in $p$-adic integration
$$\int \max(|xy|, |xz|, |yz|)^s d\mu$$
Expand over each first digit,
$$\frac{1}{p^3} \sum_{a=0}^{p-1}\sum_{b=0}^{p-1}\sum_{c=0}^{p-1}\int \max(|(a+px)(b+py)|, |(a+px)(c+pz)|, |(b+py)(c+pz)|)^s d\mu$$
Now we split the sums between $0$ and $1$ to $p-1$, since $0$ as a first digit causes the p-adic absolute value to be smaller, the other digits cause it to be 1.
When there are no 0s or only one 0 as a digit, then the max is simply 1. 
$$\frac{1}{p^3}\left((p-1)^3 + 3*(p-1)^2 + 3*\sum_{a=1}^{p-1}\int \max(|(a+px)py|, |(a+px)pz|, |p^2yz|)^s d\mu + \int\max(|p^2xy|, |p^2xz|, |p^2yz|)^s d\mu\right)$$
If we focus on just on this first integral in the sum, it simplifies substantially because $|a+px|=1$
$$\int \max(|py|, |pz|, |p^2yz|)^s d\mu$$
The 3rd term in the max is always smaller than the first two, and we can factor out a $p^{-s}$
$$p^{-s}\int \max(|y|, |z|)^s d\mu$$
This integral is exactly the one solved in the question linked at the start, so I'll put it here,
$$p^{-s}\frac{1-p^{-2}}{1-p^{-s-2}}$$
Now plugging that back in we have,
$$\frac{1}{p^3}\left((p-1)^3 + 3(p-1)^2 + 3p^{-s}\frac{1-p^{-2}}{1-p^{-s-2}} + p^{-2s}\int\max(|xy|, |xz|, |yz|)^s d\mu\right)$$
At this point we've got the original integral we started with! So all that's left is using the "snake eating it's own tail" trick. Call the integral J to make it cleaner to solve,
$$J = \frac{ (p-1)^3 + 3(p-1)^2}{p^3} + 3p^{-s-3}\frac{1-p^{-2}}{1-p^{-s-2}} + p^{-2s-3}J $$
$$J(1-p^{-2s-3}) = \frac{ (p-1)^3 + 3(p-1)^2}{p^3} + 3p^{-s-3}\frac{1-p^{-2}}{1-p^{-s-2}} $$
$$J= \frac{ (p-1)^3 + 3(p-1)^2}{p^3(1-p^{-2s-3}) } + 3p^{-s-3}\frac{1-p^{-2}}{(1-p^{-s-2})(1-p^{-2s-3}) } $$
It can be simplified more but I'll leave it here.
A: I'll go ahead and add a second answer for the general case since you asked, I think it would muddle the first answer since this stands on its own.
Now I work the general case which I take to be this integral here,
$$J_n(s) = \int \max_{i=1}^n (\prod_{\substack{j=1 \\ j \ne i}}^n |x_j|) ^s d\mu$$
Expand it along the first digits,
$$J_n(s) = \frac{1}{p^n} \sum_{\substack{a_1 \in [0, p-1] \\ \vdots \\ a_n \in [0, p-1]}}\int \max_{i=1}^n (\prod_{\substack{j=1 \\ j \ne i}}^n |a_j+px_j|) ^s d\mu$$
At this point I'm imagining how the p-adic integers break up into $\mathbb{Z}_p = p \mathbb{Z}_p \cup \mathbb{Z}_p^\times$, because this is really what matters for the maximum. So, we can cut the sum apart at the digits, $a=0$ or $a \in [1,p-1]$. In all there are $2^n$ distinct types because of the $n$ variables you're either in $p \mathbb{Z}_p$ or $\mathbb{Z}_p^\times$, and by symmetry there are $\binom{n}{k}$ of each type with $k$ nonzero first-digit variables.
When we look at a specific integral with $k$ in $p \mathbb{Z}_p$ and $n-k$ in $\mathbb{Z}_p^\times$, then we need to simplify the inner maximum function. Everywhere a variable in $\mathbb{Z}_p^\times$ appears, it no longer depends on it, as its p-adic absolute value is 1. I think a specific mini example might help, for instance if we have $k=3$ and $n=4$ then we might have an integral with a max of something like, 
$$\max(|xyz|, |wyz|, |wxz|, |wxy|)^s$$
Since $n-k=1$ I'll say wlog that $|w|=1$, this leaves us with,
$$\max(|xyz|, |yz|, |xz|, |xy|)^s$$
Because these are all p-adic integers, we know that $|x|\le 1$ and so it must be that $|xyz|\le |yz|$ and so can  be safely removed,
$$\max(|yz|, |xz|, |xy|)^s$$
I could go through this more specifically however I think it's much clear to say it this way since it makes it apparent how $k$ variables in $p\mathbb{Z}_p$ leads to the max in $J_k(s)$. So let's return to our problem, and don't forget to then factor out the power of $p$.
$$J_n(s) = \frac{1}{p^n}\left( (p-1)^n+ \sum_{k=1}^n \binom{n}{k}(p-1)^{n-k} p^{-ks+s} J_k(s)\right)$$
For the sake of a cleaner formula, let's choose $J_0(s) = p^{-s}$ by definition since it completes the formula here more nicely, just aesthetic.
$$J_n(s) = \frac{1}{p^n} \sum_{k=0}^n \binom{n}{k}(p-1)^{n-k} p^{-ks+s} J_k(s)$$
At this point we can subtract out the highest term, subtract it to the other side, then divide; now we have a relation for the higher integral strictly in terms of lower ones.
$$J_n(s) = \frac{p^{-n}}{1-p^{-n-ns+s}} \sum_{k=0}^{n-1} \binom{n}{k}(p-1)^{n-k} p^{-ks+s} J_k(s)$$
At this point, you have $J_n(s)$ in terms of strictly lower $J_k(s)$, so you can do row reduction on the system of equations. I don't know if this simplifies to any cute closed form. I realize I'm kind of terse so if you have any questions or want me to elaborate on this answer in any spots just let me know I'm happy to fill in more details for you. You should keep me honest and check to see that this does indeed end up giving you the same answer for this and the previous question.
