I am looking for a formula for the maximal order of an element in the group $\operatorname{GL}\left(n,p\right)$, where $ p$ is prime.
I recall seeing such a formula in a paper from the mid- or early 20th century, but could not find again this reference. I will be grateful for any hint.
Well, by Cayley--Hamilton, each matrix $A\in {\rm GL}(n,p)$ generates an at most $n$-dimensional subalgebra ${\mathbb F}_p[A]\subseteq M(n,p)$ thus containing at most $p^n-1$ nonzero elements. Hence the order of $A$ cannot exceed $p^n-1$.
On the other hand, consider a degree $n$ monic polynomial $P_n$ whose root is a generator $\xi$ of ${\mathbb F}_{p^n}^*$. Then a matrix with $P_n$ as its characteristic polynomial has order at least $p^n-1$ since $\xi$ is its eigenvalue.
ADDENDUM. if you wish the order to be the power of $p$, then the answer is $d=p^{\lceil \log_p n\rceil}$. Since the order of $A$ is divisible by the multiplicative orders of its eigenvalues, all the eigenvalues should be $1$. Hence the characteristic polynomial is $(x-1)^n$, so $A^d-I=(A-I)^d=0$.
On the other hand, if $A=I+J$ is the Jordan cell of size $n$ (with eigenvalue 1), then $A^{d/p}=I^{d/p}+J^{d/p}\neq I$, but $A^d=I+J^d=I$.
NB. The subgroup of all (upper-)unitriangular matrices is a Sylow $p$-subgroup in ${\rm GL}(n,p)$. So you may concentrate on it when looking at the elements of this kind.
ADDENDUM-2 (much later). This is to answer the question in the comments about the maximal order of an element $f\in AGL(n,q)$, where $q$ is a power of $p$. Write $f(x)=Ax+b$.
If $1$ is not an eigenvalue of $A$, then $f$ has a fixed point (the equation $f(x)=x$ has a solution), so we may regard it as an element of $GL(n,q)$, and the maximal order of $f$ is again $q^n-1$.
So we are concerned with the case when the minimal polynomial $\mu(x)$ of $A$ vanishes at $1$, say $\mu(x)=(x-1)^k\nu(x)$, where $\nu(1)\neq 0$. Then $A$ is similar to a block-diagonal matrix with blocks having minimal polynomials $(x-1)^k$ and $\nu(x)$ (in view of $\mathbb F_q^n=\mathop{\mathrm {Ker}}(A-I)^k\oplus\mathop{\rm Ker}\nu(A)$). So the order $d$ of $A$ does not exceed $p^{\lceil\log_p k\rceil}(q^{n-k}-1)$ if $n<k$, and $p^{\lceil\log_p n\rceil}$ otherwise. If $n>3$ (or $n=3$ and $q>2$), one may easily see that this bound does not exceed $q^{n-1}-1$. So $f^d$ is a translation, which yields that $f^{pd}=\mathord{\rm id}$, and $pd\leq p(q^{n-1}-1)<q^n-1$. Thus the maximal order of $f$ in these cases is still $q^n-1$.
We are left with the cases $n=1$, $n=2$, or $n=3$, $p=2$. When $n=1$, the answer is obviously $\max(p,q-1)$. When $n=2$, the only case left is $d=p=q$ achieved when $A$ is similar to the Jordan cell $\begin{pmatrix}1&1\\0&1\end{pmatrix}$. In this case, $f^p(x)=x+(A^{p-1}+\dots+I)b=x$ unless $p=2$, when $f^2(x)=x+\begin{pmatrix}0&1\\0&0\end{pmatrix}b$. So the answer is still $q^2-1$ when $q>2$, and $4$ otherwise.
Finally, if $n=3$ and $q=2$, then the order of $A$ having eigenvalue $1$ exceeds 3 only if $A$ is similar to $3\times 3$ Jordan cell (then $d=4$); but in this case $f^4=\mathord{\rm id}$. So this is not an exception.
Summing up, the only cases when the order may be greater than $q^n-1$ are: (1) $n=1$, $q=p$ (the maximal order is $p$), and (2) $n=2$, $q=2$ (the maximal order is $4$).
