Esta no es una novela, sino más bien una alternativa de derivación de @Start wearing purple's resultado. Traté de mantener todo simple y se explica, lo que resultó en un poco detallado de la solución. Así que si usted no está interesado en los detalles, usted puede seguir sólo el etiquetado ecuaciones hasta el Paso 3.
Paso 1 (Reducción de la Integral). Comenzamos con la siguiente fórmula:
$$ \mathcal{S}_m
= \int_{0}^{1} \frac{\log(1-x^m)}{1+x} \, dx
= \sum_{\omega \ : \ \omega^m = 1} \int_{0}^{1} \frac{\log(1 - \omega x)}{1+x} \, dx. \etiqueta{1} $$
Con el fin de trabajar con los RHS, la consideramos como una función de $\omega \in \Bbb{C} \setminus (1, \infty)$. El uso de la diferenciación bajo el signo integral técnica, nos encontramos con que
\begin{align*}
\int_{0}^{1} \frac{\log(1 - \omega x)}{1+x} \, dx
&= \int_{0}^{\omega} \left( \frac{d}{dz} \int_{0}^{1} \frac{\log(1 - zx)}{1+x} \, dx \right) \, dz\\
&= \int_{0}^{\omega} \left( \int_{0}^{1} \left( \frac{1}{1+x} - \frac{1}{1-zx}\right) dx \right) \frac{dz}{1+z} \\
&= \int_{0}^{\omega} \left( \frac{\log(1-z)}{z} - \frac{\log\left(\frac{1-z}{2}\right)}{1+z} \right) \, dz \\
&= -\operatorname{Li}_2(\omega) + I(\omega),
\end{align*}
donde la integral w.r.t. $z$ es llevada al interior de la región de $\Bbb{C}\setminus(1, \infty)$ y la función de $I(\omega)$ se define en $\Bbb{C}\setminus(1,\infty)$
$$ I(\omega) := -\int_{0}^{\omega} \frac{\log\big(\frac{1-z}{2}\big)}{1+z} \, dz. \tag{*} $$
Paso 2 (Propiedades de $I(\omega)$). ¿Por qué consideramos que esta integral es el que satisface las siguientes dos propiedades:
El uso de la identidad de $\text{(1)}$ y la relación que hemos desarrollado, $\mathcal{S}_m$ es escrito como
$$
\mathcal{S}_m
=\sum_{\omega \ : \ \omega^m = 1} (I(\omega) - \operatorname{Li}_2(\omega))
= -\frac{\zeta(2)}{m} + \sum_{\omega \ : \ \omega^m = 1} I(\omega). \etiqueta{2}
$$
Esto se deduce del teorema de la multiplicación para polylogarithm.
Lo que es más importante, para $\omega \in \Bbb{C}$ fuera de $(-\infty, 1) \cup (1, \infty)$, tenemos
$$
I(-\omega) + I(\omega) = \log^2 2 - \log\big(\tfrac{1+\omega}{2}\big)\log\big(\tfrac{1-\omega}{2}\big). \etiqueta{3}
$$
Aquí, tomamos un convenio que el registro de la parte es $0$ al $\omega = \pm 1$. Esto es consistente tanto con el límite de $\omega \to 1$ o $\omega \to -1$ y con el real conocidos los valores de $I(1)+I(-1)$. Observe también que $\text{(3)}$ fácilmente de la siguiente manera a partir de la integración por partes:
$$ I(-\omega) = \left[ -\log\big(\tfrac{1+z}{2}\big)\log\big(\tfrac{1-z}{2}\big) \right]_{0}^{-\omega} - \int_{0}^{-\omega} \frac{\log\big(\frac{1+z}{2}\big)}{1-z} \, dz. $$
Paso 3 (Fórmula para $\mathcal{S}_{2p}$). Al $m = 2p$ @Start wearing purple observado, $\omega^m = 1$ implica $(-\omega)^m = 1$. Por lo tanto la combinación de $\text{(2)}$$\text{(3)}$, tenemos
\begin{align*}
\mathcal{S}_{2p}
&= -\frac{\pi^2}{12p} + \frac{1}{2} \sum_{\omega \ : \ \omega^{2p} = 1} ( I(\omega) + I(-\omega)) \\
&= -\frac{\pi^2}{12p} + \frac{1}{2} \sum_{\omega \ : \ \omega^{2p} = 1} \left( \log^2 2 - \log\left(\frac{1+\omega}{2}\right)\log\left(\frac{1-\omega}{2}\right) \right).
\end{align*}
(Aquí, como se señaló en el Paso 2, se considera $\log(\frac{1+\omega}{2})\log(\frac{1-\omega}{2}) = 0$ al $\omega = \pm 1$.) Finalmente, podemos simplificar la última suma teniendo en cuenta solamente a $\omega$ $\Im(\omega) > 0$ como sigue:
\begin{align*}
\mathcal{S}_{2p}
&= p\log^2 2 - \frac{\pi^2}{12p} - \Re \sum_{\substack{\omega \ : \ \omega^{2p} = 1 \\ \Im(\omega) > 0}} \log\left(\frac{1+\omega}{2}\right)\log\left(\frac{1-\omega}{2}\right) \\
&= p\log^2 2 - \frac{\pi^2}{12p} - \Re \sum_{k=1}^{p-1} \log\left(\frac{1+e^{ik\pi/p}}{2}\right)\log\left(\frac{1+e^{i(k-p)\pi/p}}{2}\right) \\
&= p\log^2 2 - \frac{\pi^2}{12p} - \sum_{k=1}^{p-1} \left( \log\left(\cos\frac{\pi k}{2p}\right)\log\left(\sin\frac{\pi k}{2p}\right) + \frac{\pi^2}{4p^2}k(p-k) \right) \\
&= p\log^2 2 - \frac{(p^2+1)\pi^2}{24p} - \sum_{k=1}^{p-1} \log\left(\cos\frac{\pi k}{2p}\right)\log\left(\sin\frac{\pi k}{2p}\right).
\end{align*}
Esta solución se basa en la simetría entre las $I(\omega)$$I(-\omega)$, por lo que dudo de que esto funcionará para los impares $m$.