I give a geometric derivation, we want to count all integral points lying on surface $|x|+|y|+|z|=P$. Actually in 3D space, in every octant, the shape of surface is triangle shape.
For example, if $P=4$, then the shape in the first octant would be
$$
\begin{array}[ccccccccc]
\ & & & &Q(0,0,4)& & & &\\
& & &D(0,1,3))& &D(1,0,3)& & &\\
& &D(0,2,2)& &S(1,1,2)& &D(2,0,2)& &\\
&D(0,3,1)& &S(1,2,1)& &S(2,1,1)& &D(3,0,1)&\\
Q(0,4,0)& &D(1,3,0)& &D(2,2,0)& &D(3,1,0)& &Q(4,0,0)\\
\end{array}
$$
where $Q$ denotes the points that should be shared by 4 octants, $D$ denotes the points that should be shared by 2 octants and $S$ denotes the points belongs only to this octants.
So for the total octants, the number of points with S is
$$n_S=8\cdot\frac{(P-1)(P-2)}{2}=4P^2-12P+8$$
the number of points with D is
$$n_D=8\cdot\frac{3(P-1)}{2}=12P-12$$
the number of points with Q is
$$n_Q=8\cdot\frac{3}{4}=6$$
So the total number would be
$$n=n_S+n_D+n_Q=4P^2+2$$