Fastest way to compute the coefficients of this generating function?
I need the first N coefficients of a certain generating function, with N around 30000. Using the following code, about 2.8 hrs were needed for N=10000, and 20 hrs for N=20000; unfortunately, it looks like N=30000 would require days:
N = 30000
R.<z> = PowerSeriesRing(ZZ, default_prec=N+1)
gf = (1-z)^2 * sum( z^n/(1 - 2*z + z^(n+1))^2 for n in range(N))
coeffs = gf.coefficients()
Question: Is there a better (faster!) way to do this in SageMath? E.g., some alternative to PowerSeriesRing, or some improvement to how I'm using it?
(Aside: The coefficients of this gf are theoretically expected to exhibit a certain phenomenon that doesn't show up until N is around 30000, and I'm just interested in confirming it.)
It is possible to derive a formula for a particular coefficient of
gf
. However, getting them all unavoidably requires ~ $N^2$ operations, which would be time-consuming.