so you want
##\det\big(\mathbf H\big) =
\begin{vmatrix}
1 & 2 & 3 & ... & n \\
2 & 3 & 4 & ... & 1 \\
3 & 4 & 5 & ... & 2 \\
{\vdots}& {\vdots}& {\vdots} & {\vdots} & {\vdots} \\
n & 1 & 2 & ... & n-1 \notag
\end{vmatrix} ##
step 0.) Do the bonus question first and get comfortable with the companion matrix. In particular get comfortable with using it for the polynomial ##x^n -1##
step 1.) multiply your Hankel matrix ##\mathbf H## by the reflection matrix ##\mathbf J## to 'convert' your Hankel matrix into a Toeplitz that is in fact a circulant matrix. i.e. ##\mathbf H \mathbf J = \mathbf S## or
##\mathbf H = \mathbf S\mathbf J^T = \mathbf S\mathbf J## since ##\mathbf J## is a permutation matrix and symmetric (or involutive)
See "Relation between Hankel and Toeplitz matrices" on wikipedia page
https://en.wikipedia.org/wiki/Hankel_matrix
Note determinants multiply and ##\det\big(\mathbf J\big)## is easy to calculate
(a) as a hint it's involutive so you could just count the number of -1 eigenvalues which is implied by taking the trace... and you know that ##\text{trace}\big(\mathbf J\big) = 0## when n is even and ##\text{trace}\big(\mathbf J\big) = 1## when n is odd...
(b) just view ##\mathbf J## as a composition of disjoint pairwise swaps (transpositions), each of which has a permutation matrix with determinant -1 associated with it, so just count the number of these swaps...
step 2.) read this post from a couple years ago
https://www.physicsforums.com/threads/eigenvalues-of-circulant-matrices.927173/
step 3.) consider ##\sum_{k=1}^n k x^{k-1}##
for ##x = 1## this is ##\binom{n+1}{2}## as stated in my earlier post. For other ##x##, and in particular nth roots of unity other than 1, you can get this by differentiating the finite geometric series, or asking wolfram
https://www.wolframalpha.com/input/?i=sum+from+k+=+1+to+n+k+*+x^(k-1)
when you specialize to the remaining (n-1) distinct nth roots of unity (other than 1) this simplifies to
##\sum_{k=1}^n k x^{k-1} = \frac{-n}{1-x}= \frac{-n}{1-\omega^r}## for ##r \in\{1,2,...,n-1\}##
step 4.) The determinant of your circulant matrix is the product of eigenvalues, so
##\det\big(\mathbf S\big) = \binom{n+1}{2}\prod_{r=1}^{n-1}\frac{-n}{1-\omega^r} = (-1)^{n-1} n^{n-1}\binom{n+1}{2}\prod_{r=1}^{n-1}\frac{1}{1-\omega^r}##
step 5.)
*edited* to get the result via easier means of factoring instead of calculus
our polynomial
##p(x)= x^n - 1= q(x) \cdot (x-1) = (x-\omega^1)(x-\omega^2)...(x-\omega^{n-1}) \cdot (x-1)##
but synthetic division, or direct calculation tells us
##q(x)\cdot (x-1) = \big(x^{n-1} + x^{n-2} +... + x + 1\big)\cdot \big(x-1\big) = x^n - 1##
so
##q(x) =\big(x^{n-1} + x^{n-2} +... + x + 1\big) = \big(x-\omega^1\big)\big(x-\omega^2\big)...\big(x-\omega^{n-1}\big)##
and in particular
##q(1) = n = \big(1-\omega^1\big)\big(1-\omega^2\big)...\big(1-\omega^{n-1}\big) = \prod_{r=1}^{n-1} \big(1-\omega^r\big) ##
or
##\prod_{r=1}^{n-1}\frac{1}{1-\omega^r} = \frac{1}{n}##
so plugging this back in
##\det\big(\mathbf S\big) = (-1)^{n-1} n^{n-1}\binom{n+1}{2}\prod_{r=1}^{n-1}\frac{1}{1-\omega^r}= (-1)^{n-1} n^{n-2}\binom{n+1}{2}##
step 6.) fine tune the role of the sign function and the determinant of ##\mathbf J## (which is going to be +/- 1)
edit:
there's a very nasty but very small bug in the above related to zero vs 1 indexing. If you look at the implied circulant matrix ##\mathbf S## it has ##n## as the component on its diagonal, not ##1##.
i.e. with the list
##[1,2,3,...,n]##
the writeup above contemplates (see link from 2 years ago)
##\mathbf S = \sum_{k=0}^{n-1} s_k \mathbf P^k = 1 \mathbf I + 2 \mathbf P^1 + ... + n \mathbf P^{n-1}##
where ##s_0## would be the zeroth component in that list (i.e. value of 1), then ##s_1## the next component (value of 2) and so on,
but what actually occurs in the problem is
##\mathbf S = \sum_{k=1}^{n} s_k \mathbf P^k = 1 \mathbf P + 2 \mathbf P^2 + ... + n \mathbf P^{n}= \mathbf P\big(1 \mathbf I + 2 \mathbf P^1 + ... + n \mathbf P^{n-1}\big)##
i.e. grabbing components, where we start counting at one, so ##s_1 = 1##, ##s_2 = 2## and so on, observing that ##\mathbf P^n = \mathbf I##
This kind of thing is to be expected when I toggle 0 vs 1 indexing I suppose.
The good news is we are still dealing with nth roots of unity, ##\omega##, with just one extra multiplication by them for each respective eigenvalue. This obviously has no impact to the eigenvalue of 1. For the other nth roots, the problem becomes estimating
##\prod_{r=1}^{n-1}\frac{\omega^r}{1-\omega^r}= \big(\prod_{r=1}^{n-1}\frac{1}{1-\omega^r}\big)\big(\prod_{r=1}^{n-1}\omega^r\big) ##
so we can use our before calculation of ##\big(\prod_{r=1}^{n-1}\frac{1}{1-\omega^r}\big)##
and now we need an estimate of ##\big(\prod_{r=1}^{n-1}\omega^r\big) ##
for odd n, ##\big(\prod_{r=1}^{n-1}\omega^r\big) =1 ## since -1 is not an odd nth root of unity i.e. ##(-1)^n = -1 \neq 1## and every term in here is accordingly non-real, coming in conjugate pairs, each product of which is 1, thus the total product is 1.
for even n ##\big(\prod_{r=1}^{n-1}\omega^r\big) = -1##, by the same reasoning above, except the purely real -1 term shows up in there once and ##-1 \cdot 1 = - 1##.
It's a small sign error that got buried in some gorey details. Apologies.
This means for odd n
##\det\big(\mathbf S\big) = n^{n-2}\binom{n+1}{2}##
as before, and for even n it now should be
##\det\big(\mathbf S\big) = n^{n-2}\binom{n+1}{2}##
circling back, and using earlier hints
##\det\big(\mathbf J\big) = (-1)^{\frac{n}{2}}## for even n
##\det\big(\mathbf J\big) = (-1)^{\frac{n-1}{2}}## for odd n
but this is the same as the official answer... i.e. ##(-1)^\frac{n(n-1)}{2}## is just a symmetrized version of this.
One way to check is, for even n:
##\frac{n}{2}\%2 = \frac{n}{2}(n-1)\%2##
because n-1 is odd and equal to 1, mod 2
and for odd n
##\frac{n-1}{2}\%2= \frac{(n-1)}{2}n\%2##
because ##n## is odd and equal to 1, mod 2
so with a rather ugly bandaid at the end, this completes the derivation
