On generalizing $\frac{17-\sqrt{17}}8 =\sin^2(t)+\sin^2(2t)+\sin^2(4t)+\sin^2(8t)$? Given $\color{blue}{t = 2\pi/p}$ for the appropriate prime $p=4m+1$.
I. Sine
$$\begin{align}
\frac{5+\sqrt{5}}8 &=\sin^2(t)\\ 
\frac{13+\sqrt{13}}8 &=\sin^2(t)+\sin^2(3t)+\sin^2(4t)\\
\frac{17+\sqrt{17}}8 &=\sin^2(3t)+\sin^2(5t)+\sin^2(6t)+\sin^2(7t)\\
\frac{29+\sqrt{29}}8 &=\sum_{k=1}^7\sin^2(a_k\, t)
\end{align}$$
with the seven $a_k = 1,4,5,6,7,9,13.$ And so on for other prime $p=4m+1.$ For the opposite sign, one uses the remaining integers $b_k \leq \frac{p-1}2.$ For example,
$$\frac{17-\sqrt{17}}8 =\sin^2(t)+\sin^2(2t)+\sin^2(4t)+\sin^2(8t)$$
where the $a_k$ are simply $2^n$. (The next Fermat prime $p=257$ isn't so nice since it has $256/4 = 64$ sine terms.)

II. Cosine
This uses the same set of multipliers $a_k$.
$$\begin{align}
\frac{3-\sqrt{5}}8 &=\cos^2(t)\\ 
\frac{11-\sqrt{13}}8 &=\cos^2(t)+\cos^2(3t)+\cos^2(4t)\\
\frac{15-\sqrt{17}}8 &=\cos^2(3t)+\cos^2(5t)+\cos^2(6t)+\cos^2(7t)\\
\frac{27-\sqrt{29}}8 &=\sum_{k=1}^7\cos^2(a_k\, t)
\end{align}$$
with the same seven $a_k = 1,4,5,6,7,9,13.$ For the opposite sign,
$$\frac{15+\sqrt{17}}8 =\cos^2(t)+\cos^2(2t)+\cos^2(4t)+\cos^2(8t)$$

III. Conclusion
Given prime $p=4m+1$ and $t = 2\pi/p.$ The pattern clearly is,
$$\frac{p\pm\sqrt{p}}8 = \sum_{k=1}^m \sin^2(a_k\, t)$$
$$\frac{(p-2)\mp\sqrt{p}}8 = \sum_{k=1}^m \cos^2(a_k\, t)$$
Question: I used Mathematica's integer relations to find the above examples. But, for any prime $p=4m+1$, what is a clever and faster algorithm to derive the correct set of $a_k$?

I found an algorithm to derive the correct set $H$ of $a_k$ for any prime $\ p=4m+1 \ $ (or its complement).
$H=\{ 1\}$.
For $k$ from $2$ to $2m$:

*

*Compute $\alpha$ such that $\ \ k^2 \equiv \alpha \pmod{4m+1}\ \ $ and $\ \ 1\leqslant \alpha \leqslant 4m \ \ $

*if $\ 1\leqslant \alpha \leqslant 2m \ $ then  $\ \beta =\alpha \ $ else $\ \beta = 4m+1-\alpha $

*$H = H\cup \{\beta \}$
We can use "Quadratic Gauss sums" to show the result.
