ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Tue, 30 Jan 2024 03:44:12 +0100Fastest way to compute the coefficients of this generating function?https://ask.sagemath.org/question/75699/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.)Mon, 29 Jan 2024 21:05:39 +0100https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/Comment by Max Alekseyev for <p>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:</p>
<pre><code>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()
</code></pre>
<p><strong>Question</strong>: 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? </p>
<p>(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.)</p>
https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?comment=75702#post-id-75702It 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.Mon, 29 Jan 2024 22:53:42 +0100https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?comment=75702#post-id-75702Answer by Max Alekseyev for <p>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:</p>
<pre><code>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()
</code></pre>
<p><strong>Question</strong>: 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? </p>
<p>(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.)</p>
https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?answer=75704#post-id-75704Let $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)`.Mon, 29 Jan 2024 23:13:30 +0100https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?answer=75704#post-id-75704Comment by res0001 for <p>Let $c_k$ be the coefficient of $z^k$ in <code>sum( z^n/(1 - 2*z + z^(n+1))^2 for n in range(N))</code>, 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:</p>
<pre><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))
</code></pre>
<p>Then, the coefficient of $z^k$ in <code>(1-z)^2 * sum( z^n/(1 - 2*z + z^(n+1))^2 for n in range(N))</code> can be obtained simply as <code>c(k) - 2*c(k-1) + c(k-2)</code>.</p>
https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?comment=75707#post-id-75707Thanks 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]$.)Tue, 30 Jan 2024 02:15:15 +0100https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?comment=75707#post-id-75707Comment by Max Alekseyev for <p>Let $c_k$ be the coefficient of $z^k$ in <code>sum( z^n/(1 - 2*z + z^(n+1))^2 for n in range(N))</code>, 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:</p>
<pre><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))
</code></pre>
<p>Then, the coefficient of $z^k$ in <code>(1-z)^2 * sum( z^n/(1 - 2*z + z^(n+1))^2 for n in range(N))</code> can be obtained simply as <code>c(k) - 2*c(k-1) + c(k-2)</code>.</p>
https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?comment=75708#post-id-75708Yes, it was a typo. Now fixed.Tue, 30 Jan 2024 03:44:12 +0100https://ask.sagemath.org/question/75699/fastest-way-to-compute-the-coefficients-of-this-generating-function/?comment=75708#post-id-75708