El método estándar para resolverlos consiste en filtro de raíces de la unidad (o análisis de Fourier finito para los pretenciosos).
La idea es mirar el polinomio f(x)=x^{d-a_0}(1+x)^n . Esto tiene la propiedad de que el coeficiente de x^{dk} es \binom{n}{a_0+dk} . Así que si pudiéramos extraer de alguna manera sólo los coeficientes de x^{dk} y sumarlos, habríamos terminado.
El truco está en tener en cuenta \frac{f(1)+f(\omega)+\dots+f\left(\omega^{d-1}\right)}{d}, donde \omega=\exp(2i\pi/d) es el d -enésima raíz de la unidad. ¿Puedes ver por qué esto extrae sólo los coeficientes de x^{dk} en la suma?
Para concretar, veamos el siguiente ejemplo: \sum_{k\geq0}\binom{2021}{3k+1}. Tomemos el polinomio f(x)=x^2(1+x)^{2021} por lo que escribir \omega=\exp(2i\pi/3) la suma que queremos es \frac{(1+1)^{2021}+\omega^2(1+\omega)^{2021}+\omega(1+\omega^2)^{2021}}{3}. Ahora, en general, en este punto se acaba de convertir a la forma polar 1+\omega^m=re^{i\theta} y golpearlo, pero para n=3 puede aprovechar el ingenioso truco que 1+\omega+\omega^2=0 para facilitar el cálculo: \begin{align*} \frac{(1+1)^{2021}+\omega^2(1+\omega)^{2021}+\omega(1+\omega^2)^{2021}}{3} &= \frac{1}{3}\left[2^{2021}+\omega^2(-\omega^2)^{2021}+\omega(-\omega)^{2021}\right] \\ &= \frac{1}{3}\left[2^{2021}-\omega^{4044}-\omega^{2022}\right] \\ &=\frac{1}{3}\left[2^{2021}-2\right]. \end{align*}