integer solution of $\frac{x^3+y^3+z^3-xy(x+y)-yz(y+z)-xz(x+z)-2xyz}{(x+y+z)(x+y-z)(x-y+z)(x-y-z)}=\frac{1}{2016}$ Let $x,y,z$ be positive integers such that $\frac{x^3+y^3+z^3-xy(x+y)-yz(y+z)-xz(x+z)-2xyz}{(x+y+z)(x+y-z)(x-y+z)(x-y-z)}=\frac{1}{2016}$.
How to find all solutions ?
I have no any idea.
Thanks in advance.
 A: Hint:
Since$$\frac{x^{3}+y^{3}+z^{3}-xy(x+y)-yz(y+z)-zx(z+x)-2xyz}{(x+y+z)(x+y-z)(x-y+z)(x-y-z)}=\frac{1}{2}\left(\frac{1}{-x+y+z}+\frac{1}{x-y+z}+\frac{1}{x+y-z}+\frac{1}{x+y+z}\right),$$
 we are trying to solve $$\frac{1}{-x+y+z}+\frac{1}{x-y+z}+\frac{1}{x+y-z}+\frac{1}{x+y+z}=\frac{1}{1008}.$$
 This is equivalent to solving the equation $$\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{a+b+c}=\frac{1}{1008}.$$
A: This method uses Pythagorean triples $x^2+y^2=z^2$ to solve,
$$\frac{1}{-x+y+z}+\frac{1}{x-y+z}+\frac{1}{x+y-z}+\frac{1}{x+y+z}=\frac{1}{N}\tag1$$
It does not give all the solutions, but does give an easy way to tackle general $N$. 
Let $x\leq y\leq z$. The method is based on the fortuitous observation that some $x,y,z$ (of a higher search range) were Pythagorean triples. Substitute, 
$$x,\,y,\,z = a^2-b^2,\;2ab,\;a^2+b^2$$
into $(1)$, and one gets a much simpler equation,
$$(a-b)b = N\tag2$$
Thus, quick solutions can be found using the divisors of $N$. For $N = 1008$, this has thirty, namely,
$$b = \color{brown}{1}, 2, 3, 4, 6, 7, 8, 9, 12, 14, 16, 18, 21, 24, 28, 36, 42, 48, 56, 63, 72, 84, 112, 126, 144, 168, 252, 336, 504, \color{brown}{1008}.$$
For example, using $a,b = 1009,\color{brown}1$ yields,
$$x,\,y,\,z = 1018080,\; 2018,\; 1018082$$
while $a,b = 1009,\color{brown}{1008}$ gives,
$$x,\,y,\,z = 2017,\; 2034144,\; 2034145$$
and so on for thirty solutions to $(1)$ that also obey $x^2+y^2=z^2$.
