A fancy approach to Lagrange's identity: differentiation.
$$\frac{\partial}{\partial a}(a^2+b^2)(c^2+d^2) = 2a(c^2+d^2),$$
$$\frac{\partial}{\partial a}\left((ac-bd)^2+(ad+bc)^2\right)=2c(ac-bd)+2d(ad+bc) $$
they match, and the same happens for $\frac{\partial}{\partial b},\frac{\partial}{\partial c},\frac{\partial}{\partial d}$, so $(a^2+b^2)(c^2+d^2)$ and $(ac-bd)^2+(ad+bc)^2$ differ by a constant, that is clearly zero by evaluating both functions at $(1,1,1,1)$. The usual approach: the norm on $\mathbb{Z}[i]$ (the Gaussian integers) defined through
$$ N(z) = N(a+ib) = z\cdot\overline{z} = (a+ib)(a-ib) = a^2+b^2$$
is multiplicative:
$$ N(zw) = zw\cdot\overline{zw} = z\overline{z}\cdot w\overline{w} = N(z)N(w)$$
hence by taking $z=a+ib,w=b+id$ Lagrange's identity easily follows.
Lagrange's identity gives that the set $S=\{n:n=a^2+b^2,a\in\mathbb{Z},b\in\mathbb{Z}\}$ is a semigroup, i.e. if $s_1,s_2\in S$, then $s_1\cdot s_2\in S$. Using Fermat's descent, it is not difficult to show that if $m\in\mathbb{N}$ belongs to $S$ and $m=uv$ with $\gcd(u,v)=1$, then both $u$ and $v$ belong to $S$. Moreover, the only prime numbers belonging to $S$ are $2$ and the primes of the form $4k+1$. Since
$$25988=2^2\cdot 73\cdot 89 $$
and both $73$ and $89$ belong to $S$, $25988\in S$, too. To be clear, from
$$ 73=8^2+3^2,\qquad 89=8^2+5^2 $$
it follows that:
$$ 73\cdot 89 = (64-15)^2+(24+40)^2 = 49^2+64^2, $$
hence $25988 = 98^2+128^2$.