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$?
 A: 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.
A: Here's how it plays out when $p\equiv1\pmod8$. I think the assumption is needed for then we know that both $-1$ and $2$ are quadratic residues. With $p\equiv5\pmod 8$ there is probably a modification, but I'm down with a flu, so it will have to wait.
Let $Q$ stand for the set of (non-zero) quadratic residues modulo $p$. Let $D\subset Q$ be a set of representative of the quotient group $Q/\langle -1\rangle$. For example, we can choose $D$ to be the set of quadratic residues $\le (p-1)/2$. Let $N$ be the set of quadratic non-residues modulo $p$.
I will need the following fact (proof is not difficult, and very likely done elsewhere on the site):
Fact. If $j$ and $k$ range over $0<j,k\le(p-1)/2$, then the sum $j^2+k^2$ takes the values in $\Bbb{Z}_p^*$ with the following frequencies:
$$
j^2+k^2=\begin{cases}
0,\ \text{$(p-1)/2$ times},\\
q,\ \text{for every element $q\in Q$ exactly $(p-5)/4$ times},\\
n,\ \text{for every element $n\in N$ exactly $(p-1)/4$ times}.
\end{cases}
$$
Let us define
$$
S=\sum_{j=0}^{p-1}e^{2\pi j^2 i/p}=2\sum_{x\in Q}e^{2\pi x i/p}.
$$
The fact implies that in the sum $S^2+S$ every root of unity $e^{2\pi i\ell/p}$
appears exactly $(p-1)/4$ times except the root $1$ coming from $\ell=0$ that occurs $(p-1)/4$ times more often than the others.
By the well known result that the sum of all the roots of unity vanishes, this gives the equation
$$
S^2+S=(p-1)/4.\qquad(*)
$$
Applying the quadratic formula to $(*)$ implies that
$$
S=\frac{-1\pm\sqrt p}2.
$$
The famous Gauss's sign problem tells us that the plus sign always applies.

On with the sums of trig functions. Because $-1\in Q$, the roots of unity gives twice the cosines of $2\pi k/p$, $k\in D$. Next we use the trig identity
$$
\cos^2\alpha=\frac{1+\cos2\alpha}2.
$$
Because also $2\in Q$, we see that $2k$ ranges over $Q$ as $k$ does. Therefore
$$
\begin{aligned}
\sum_{k\in D}\cos^2\frac{2\pi k}p&=\sum_{k\in D}\frac12\left(1+\cos(2\pi 2k/p)\right)\\
&=\frac{p-1}8+\frac S4\\
&=\frac{p-2+\sqrt p}8.
\end{aligned}
$$
The sum of squares of sines can be handled similarly using the formula $\sin^2\alpha=(1-\cos2\alpha)/2$.

A Mathematica-example with $p=17$.
In[1]:=Sort[Table[Mod[k^2,17],{k,1,8}]]
Out[1]:={1, 2, 4, 8, 9, 13, 15, 16}

This implies that we can use $D=\{1,2,4,8\}$.
In[2]:=exp17={1,2,4,8};
In[3]:=N[Sum[Cos[2 Pi exp17[[ k]]/17]^2, {k, 1, 4}]]
Out[3]=2.39039

And the final check
In[4]:=N[(15 + Sqrt[17])/8]
Out[4]=2.39039

More later, if/when needed. I need to rest.
A: The ones that have all even or all odd quadratic residues relative to p, will do the trick.
The proof for this lies is polygon isomorphisms.  The numbers are invariant if taken to an e power (ie an even quadratic), but not an odd power.  So you end up with two numbers a, b whose sum would give p, and the difference is something in sqrt(p),
