### 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.)