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$?
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$.
This implies that we can use $D=\{1,2,4,8\}$.
And the final check
More later, if/when needed. I need to rest.