Does anyone know how to prove (or have a reference for) the fact that the inversion $$I : \Bbb R^n \setminus \{0\} \rightarrow \Bbb R^n \setminus \{0\}\\ x \mapsto \frac x {\|x\|^2}$$ sends circles to circles or straight lines ?
Since $I = L^{-1} \circ I \circ L$ on $\Bbb R^n \setminus \{0\}$ for any linear invertible isometry $L : \Bbb R^n \rightarrow \Bbb R^n$, we can assume that the circle $C$ we consider is contained in an affine $2$-plane, say $\Bbb R^2 \times \{h\}$ for some $h \in \Bbb R^{n-2}$.
If $h = 0$, the proof is straightforward. However, I don't know how to prove the case $h \neq 0$.
But maybe this is not the way to go to prove this statement.
I saw here https://mathoverflow.net/questions/25251/properties-of-the-n-dimensional-stereographic-projection someone saying something about circles being the intersection of several $(n-1)$-spheres (because it's quite easy to prove that $(n-1)$-spheres are sent to $(n-1)$-spheres or hyperplanes), but I don't really want to use this claim.
The inverse of a vector $x$ is defined by
$$\frac1x = \frac{x}{x^2} = \Big(\frac1{x\cdot x}\Big)x$$
If the circle's centre is $c$, its radius is $\lVert a\rVert=\lVert b\rVert$, and its plane is spanned by vectors $a$ and $b$ (with $a\cdot b=0$), then it can be parametrized by
$$x = c + a\cos\theta + b\sin\theta$$
Without loss of generality, $a$ and $b$ can be rotated until $b$ is orthogonal to $c$. So now we have
$$a^2 = b^2,\qquad a\cdot b = b\cdot c = 0$$
and
$$\frac1x = \frac{c+a\cos\theta+b\sin\theta}{c^2+2c\cdot a\cos\theta+a^2}$$
This equation doesn't look like a circle, but of course it must be. The new circle's centre $c'$ should be the midpoint of the two extremal points at $\theta=0,\;\theta=\pi$ :
$$c' = \frac12\Big(\frac1{c+a} + \frac1{c-a}\Big) = \overset{\text{algebra}}{\cdots} = \frac{(c^2+a^2)c-2(c\cdot a)a}{(c+a)^2(c-a)^2}$$
And the new (directed) radius $a'$ should be half the displacement between these points:
$$a' = \frac12\Big(\frac1{c+a} - \frac1{c-a}\Big) = \cdots = \frac{(c^2+a^2)a-2(a\cdot c)c}{(c+a)^2(c-a)^2}$$
and $b' = \lVert a'\rVert\frac{b}{\lVert b\rVert}$ .
(These are undefined if $c+a=0$ or $c-a=0$, but then the result is a straight line, and the problem is 2-dimensional, which is easily solved.)
If these guesses are correct, then the new point must be in this plane, and have constant radius:
(The wedge product $\wedge$ is a measure of linear independence; any vector $a$ has $a\wedge a=0$.)
To verify these equations, first note that the denominator is
$$(c+a)^2(c-a)^2 = (c^2+a^2+2c\cdot a)(c^2+a^2-2c\cdot a) = (c^2+a^2)^2-4(c\cdot a)^2$$
thus
$$c'^2 = \bigg(\frac{(c^2+a^2)c-2(c\cdot a)a}{(c+a)^2(c-a)^2}\bigg)^2 = \frac{(c^2+a^2)^2c^2-4(c^2+a^2)(c\cdot a)^2+4(c\cdot a)^2a^2}{(c+a)^4(c-a)^4}$$
$$= \frac{\big((c^2+a^2)^2-4(c\cdot a)^2\big)c^2}{(c+a)^4(c-a)^4} = \frac{c^2}{(c+a)^2(c-a)^2}$$
$$a'^2 = \cdots = \frac{a^2}{(c+a)^2(c-a)^2}$$
and
$$c'\wedge a' = \bigg(\frac{(c^2+a^2)c-2(c\cdot a)a}{(c+a)^2(c-a)^2}\bigg)\wedge\bigg(\frac{(c^2+a^2)a-2(a\cdot c)c}{(c+a)^2(c-a)^2}\bigg)$$
$$= \frac{(c^2+a^2)^2(c\wedge a)-2(c^2+a^2)(c\cdot a)(c\wedge c)-2(c^2+a^2)(c\cdot a)(a\wedge a)+4(c\cdot a)^2(a\wedge c)}{(c+a)^4(c-a)^4}$$
$$= \frac{\big((c^2+a^2)^2-4(c\cdot a)^2\big)(c\wedge a)}{(c+a)^4(c-a)^4} = \frac{c\wedge a}{(c+a)^2(c-a)^2}$$
Now we can calculate the two expressions:
$$\Big(\frac1x - c'\Big)\wedge a' = \frac{c+a\cos\theta+b\sin\theta}{c^2+2c\cdot a\cos\theta+a^2}\wedge\frac{(c^2+a^2)a-2(a\cdot c)c}{(c+a)^2(c-a)^2}-c'\wedge a'$$
$$= \frac{(c^2+a^2)(c\wedge a)-0+0-2(c\cdot a\cos\theta)(a\wedge c)}{(c^2+2c\cdot a\cos\theta+a^2)(c+a)^2(c-a)^2}+\frac{b\sin\theta}{x^2}\wedge a'-c'\wedge a'$$
$$= \frac{(c^2+2c\cdot a\cos\theta+a^2)(c\wedge a)}{(c^2+2c\cdot a\cos\theta+a^2)(c+a)^2(c-a)^2}+\frac{b\sin\theta}{x^2}\wedge a'-c'\wedge a'$$
$$= \frac{c\wedge a}{(c+a)^2(c-a)^2}-c'\wedge a'+\frac{b\sin\theta}{x^2}\wedge a'$$
$$= 0+\frac{b\sin\theta}{x^2}\wedge a'$$
$$\Big(\frac1x - c'\Big)\wedge a'\wedge b' = \frac{b\sin\theta}{x^2}\wedge a'\wedge\frac{b\lVert a'\rVert}{\lVert b\rVert} = 0$$
and
$$\Big(\frac1x - c'\Big)^2 = \frac1{x^2} - 2\frac{x\cdot c'}{x^2} + c'^2$$
$$= \frac1{x^2}-\frac2{x^2}\frac{(c+a\cos\theta+b\sin\theta)\cdot\big((c^2+a^2)c-2(c\cdot a)a\big)}{(c+a)^2(c-a)^2}+c'^2$$
$$= \frac{(c+a)^2(c-a)^2}{x^2(c+a)^2(c-a)^2}-2\frac{(c^2+a^2)c^2-2(c\cdot a)^2+(c^2+a^2)(c\cdot a\cos\theta)-2a^2(c\cdot a\cos\theta)+0-0}{x^2(c+a)^2(c-a)^2}+c'^2$$
$$= \frac{(c^2+a^2)^2-4(c\cdot a)^2}{x^2(c+a)^2(c-a)^2}-2\frac{(c^2+a^2)c^2-2(c\cdot a)^2+(c^2-a^2)(c\cdot a\cos\theta)}{x^2(c+a)^2(c-a)^2}+c'^2$$
$$= \frac{(c^2+a^2)(c^2+a^2-2c^2)-2(c^2-a^2)(c\cdot a\cos\theta)}{x^2(c+a)^2(c-a)^2}+c'^2$$
$$= \frac{(c^2+a^2+2c\cdot a\cos\theta)(a^2-c^2)}{(c^2+2c\cdot a\cos\theta+a^2)(c+a)^2(c-a)^2}+c'^2$$
$$= \frac{a^2-c^2}{(c+a)^2(c-a)^2}+\frac{c^2}{(c+a)^2(c-a)^2}$$
$$= \frac{a^2}{(c+a)^2(c-a)^2}$$
$$=a'^2$$