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

edit retag close merge delete

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.

( 2024-01-29 22:53:42 +0200 )edit

Sort by ยป oldest newest most voted

Let $c_k$ be the coefficient of $z^k$ in sum( z^n/(1 - 2*z + z^(n+1))^2 for n in range(N)), then for $0 < k < N$ we have $$c_k = [z^k] \sum_{n\geq 1} z^{n} (1-2z+z^{n+1})^{-2} = [z^k] \sum_{n\geq1} \sum_{m\geq0} z^{n+m} \binom{-2}{m} (-2+z^n)^m$$ $$=\sum_{n=1}^k \sum_{m=0\atop m\equiv k\pmod{n}}^{k-n} \binom{-2}{m} \binom{m}{\frac{k-n-m}n}(-2)^{m-\frac{k-n-m}n}.$$ Here is a corresponding code:

def c(k):
return 1 if k==0 else sum(binomial(-2,m) * binomial(m,(t:=(k-m-n)//n)) * (-2)^(m-t) for n in range(1,k+1) for m in range(k%n,k-n+1,n))


Then, the coefficient of $z^k$ in (1-z)^2 * sum( z^n/(1 - 2*z + z^(n+1))^2 for n in range(N)) can be obtained simply as c(k) - 2*c(k-1) + c(k-2).

more

Thanks for this answer! I hadn't realized one might find such a formula. Using your c(k) in a program to compute the final coefficients for a selected range of k-values, I find that it takes about 7 minutes to compute the ten of them for k=30000..30009. (I'll have to study the method you used to derive that formula. BTW, I'm sure you intended to write $[z^k]$.)

( 2024-01-30 02:15:15 +0200 )edit

Yes, it was a typo. Now fixed.

( 2024-01-30 03:44:12 +0200 )edit