According to Weisstein's World of Mathematics, it was first proved by Gauss in 1823. The reference is volume 4 of Gauss' Werke, which can be read at https://archive.org/details/werkecarlf04gausrich. The relevant pages seem to be 47-49. It seems that Gauss investigated the question and came up with a proof. I don't read Latin, but there is a German summary in the text. Pages 103-104 explain what he did (Edit: I added a rough translation):
Allein da man nicht berechtigt ist, die sichersten Werthe fuer die
wahren Werthe selbst zu halten, so ueberzeugt man sich leicht, dass
man durch dieses Verfahren allemal den wahrscheinlichsten und
mittleren Fehler zu klein finden muss, und daher die gegebenen
Resultaten eine groessere Genauigkeit beilegt, als sie wirklich
besitzen. [But since one is not entitled to treat the most probable values as though they were the actual values, one can easily convince oneself that one must always find that the most probable error and the average error are too small, and that therefore the given results possess a greater accuracy than they really have.]
from which it would seem that it was well-known that the sample variance is a biased estimate of the population variance. The article goes on to say that the difference between the two is usually ignored because it's not important if the sample size is big enough. Then it says:
Der Verfasser hat daher diesen Gegenstand eine besondere Untersuchung
unterworfen, die zu einem sehr Merkwuerdigen hoechst einfachen
Resultate gefuehrt hat. Man braucht nemlich den nach dem angezeigten
fahlerhaften Verfahren gefundenen mittleren Fehler, um ihn in die
richtigen zu verwandeln, nur mit
$$\sqrt{\frac{\pi-\rho}{\pi}}$$
zu multiplicieren, wo $\pi$ die Anzahl der beobachtungen (number of
observations) und $\rho$ die Anzahl der unbekannten groessen (number
of unknowns) bedeutet. [The author has therefore made a special study of this object which has led to a very strange and extremely simple result. Namely, one needs only to multiply the average error found by the above erroneous process by (the given expression) to change it into the right one, where $\pi$ is the number of observations and $\rho$ is the number of unknown quantities.]
So if this is indeed the first time that the correction was found, then it seems that it was found by a clever calculation by Gauss, but people were already aware that some correction was required, so perhaps someone else could have found it empirically before this. Or possibly previous authors didn't care to derive the precise answer because they were working with fairly large data sets anyway.
Summary: manual, but people already knew that $n$ in the denominator wasn't quite right.