The number $\pi$ and summation by $SL(2,\mathbb Z)$ Let $f(a,b,c,d)=\sqrt{a^2+b^2}+\sqrt{c^2+d^2}-\sqrt{(a+c)^2+(b+d)^2}$. (it is the defect in the triangle inequality) 
Then, we discovered by heuristic arguments and then verified by computer that
$$\sum f(a,b,c,d)^n = 2-\pi/2$$
where the sum runs over all $a,b,c,d\in\mathbb Z$ such that $a\geq 1,b,c\geq 0, ad-bc=1$ and $n=2$.
It seems that when $n=1$, we obtain $2$ in the right hand side. We have failed to guess the result for $n>2$.
So the question is: can you prove the result for $n=2$? (We can, but in a rather unnatural way. We will write this later and now we want to tease the community, probably somebody can find a beautiful proof.) 
Added: we have two answers for the above question, so the rest is:
Question: Can you guess the result for $n>2$? I tried http://mrob.com/pub/ries/ but nothing interesting was revealed.
PS. In case it can help someone, below are these sums of powers ($n=1,2,3,4,5$), calculated by computer:
$1.9955289122768913 = 2$
$0.4292036731309361 = 2 - \pi/2$
$0.21349025954227965 = $ ?
$0.11983665032283052 = $ ?
$0.06933955916793563 = $ ?
Added: this and something more can be found in https://arxiv.org/abs/1701.07584 and https://arxiv.org/abs/1711.02089
 A: Here goes $n=2$ a la Fedor Petrov.
Notice that his argument is based on the smart identities $x_z+y_z=|z|$, $z_x=|x|+y_x$ and $z_y=|y|+x_y$ where $x,y$ are 2 vectors in $\mathbb R^2$, $z=x+y$ and $x_y$ is the (signed in general, but in our setting everything is non-negative) length of the projection of $x$ to $y$. Fedor's idea can be made into a one-liner: consider $\Phi_1(x,y)=|x|+|y|-x_y-y_x$, check the identity $\Phi_1(x,y)-\Phi_1(z,x)-\Phi_1(z,y)=|x|+|y|-|z|$ and telescope.
For $n=2$ the right function to consider is 
$$
\Phi_2(x,y)=2|x||y|-|x|y_x-|y|x_y
$$ 
I wrote it in the symmetric form, but actually $|x|y_x=|y|x_y=|x||y|\cos\theta$ where $\theta$ is the angle between $x$ and $y$. The difficulty is that in this case the decay "at infinity" is not fast enough to send the sum of faraway terms to $0$: the typical term there is $2|x||y|(1-\cos\theta)=(2|x||y|\sin\frac{\theta }2)2\sin\frac{\theta}2\approx A\theta$ where $A$ is the (invariant) area of the parallelogram spanned by $x$ and $y$. The sum of these expressions at infinity adds an extra term $A$ times the angle between $x$ and $y$ with $-$ sign (the precision is now good enough to ignore the rest; like Fedor, I'll leave the "routine convergence checks" to the reader). So, for two perpendicular unit vectors, we get $2\cdot1\cdot1-0-0-1\cdot\frac \pi 2$ as requested.
I wonder if we can continue with that a bit. The first task would be to find $\Phi_3$ but, even if we are lucky and it exists, it may require even more correcting terms when trying to telescope and we can get a transcendental problem there.
A: Let me write down here a proof that $\sum f(a,b,c,d)=2$, maybe someone sees how this may be generalized for the second moment. I do not.
We denote the vectors $x=(a,b)$ and $y=(c,d)$ and write $g(x,y)=f(a,b,c,d)$. Denote $S=\sum g(x,y)$, I omit here a routine proof that $S$ is finite. We may fix $x$ and sum up by $y$. For given $x$ the possible $y$ have form $y_0(x)+kx$ for $k=0,1,2,\dots$. The sum by all these $y$  equals
$$
\lim_N \sum_{k=0}^{N-1} \|x\|+\|y_0(x)+kx\|-\|y_0(x)+(k+1)x\|=\\ =\lim_N N\|x\|+\|y_0(x)\|-\|y_0+Nx\|=\|y_0(x)\|-pr_x (y_0(x)),
$$
where $pr_x$ denotes projection onto $x$. So, $S=\sum_x \|y_0(x)\|-pr_x (y_0(x)).$ Analogously we may fix $y$ and sum up by $x$, and get a similar formula $S=\sum_y \|x_0(y)\|-pr_y (x_0(y))$. In the first formula we take a summand corresponding to $x=(1,0)$ and in the second a summand corresponding to $y=(0,1)$. Those summands are equal to 1. Other summands correspond to $x$, resp. $y$, with strictly positive coordinates. Any such vector, which I denote by $z$, is a sum of $x_0(z)$ and $y_0(z)$, thus when we sum up two expressions for $S$ we get 
$$
2S=1+1+\sum_z \|x_0(z)\|+\|y_0(z)\|-\|z\|=2+S,
$$
and $S=2$.
A: I add our explanation and the origin of the problem. 
To obtain the formula we just need to verify the following lemma (by a   straightforward computation). 
Magic Lemma. Let $(a,b),(c,d)$ be as in our summation in the main posting. Draw tangents to the unit circle which are orthogonal to the directions $(a,b),(c,d),(a+c,b+d)$. Then the triangle in the intersection has the area $\frac{f(a,b,c,d)^2}{2}$. 
Now consider the circle given by $x^2+y^2=1$. At every rational point on it we draw the tangent line. Clearly, our disk is the intersection of all the half-planes given by these lines.
So we can approximate the area of the disk by the following sequential cutting. We start with the square $[-1,1]\times [-1,1]$. It has four points of tangency with our circle. We will construct a sequence of polygons converging to the circle. Initially we have sides of directions orthogonal to $(1,0),(0,1),(-1,0),(0,1)$. At every step, we take two vectors $v_1,v_2$(in the same quadrant) which give the basis of $\mathbb Z^2$ and add to our polygon a new side, which is given by the tangent line of the direction orthogonal to $v_1+v_2$.
For example, the first step will be cutting the square by the line $x+y=\sqrt 2$, if we take the vectors $(1,0),(0,1)$. The area of the eliminated triangle is $\frac{(\sqrt2-1)^2}{2}=\frac{f(1,0,0,1)^2}{2}$.
It seems that this procedure can be generalized to higher dimensions — at least a preliminary computation shows that something like the above magic lemma should happen.
