# Calculating a really large binomial sum

Suppose I want to calculate the coefficient of a certain term (modulo a prime $$p$$) in the series expansion of $$f(x) = \frac{1}{1 - x^u - x^v}$$ where $$u,v$$ are positive integers. We know that $$\frac{1}{1 - x} = \sum_{k = 0}^\infty x^k$$ Therefore $$f(x) = \sum_{k = 0}^\infty \left(x^u + x^v\right)^k$$ The expansion of $$(x^u + x^v)^k$$ contains a term of degree $$m$$, if and only if there exist integers $$r,s$$ satisfying \left\{\begin{aligned}r + s &= k\\ru + sv &= m\end{aligned}\right. In this case the coefficient of that term is $$\binom{k}{r}$$.

In general all integral solutions to the equation $$ru + sv = m$$ can be represented as \left\{\begin{aligned}r &= a + ck\\s &= b + dk\end{aligned}\right. where $$a,b,c,d$$ are integer constants and $$k$$ is an integer variable.

Therefore what I want to calculate is $$\sum_{k = k_1}^{k_2} \binom{r + s}{r} \quad (\mathrm{mod}\; p)$$ where $$k_1,k_2$$ are necessary bounds on $$k$$ to ensure both $$r,s$$ are positive. If $$k_2 - k_1$$ is small I could just use Lucas's theorem to calculate the residue for each term. However in my problem $$k_2 - k_1$$ is really large (around $$10^{11}$$). Is there a better algorithm to calculate this sum?

There is: note that your coefficients $$c_n$$ satisfy a fibonnaci-style recurrence relation $$c_n=c_{n-u}+c_{n-v}$$ (this can be proven by classic generating function methods - note that $$(1-x^u-x^v)f(x)=1$$ and expand the series term-by-term on both sides; initial conditions should be easy to figure out). From there, there are plenty of methods that can be used, but the best for your purposes (large-$$n$$ coefficients modulo a prime) are probably the matrix multiplication methods - you can create a vector of coefficients $$\mathfrak{C_n}=\langle c_n-1, c_{n-1}, \ldots, c_{n-v}\rangle$$ and then the coefficients of $$\mathfrak{C_{n+1}}$$ are derivable by just taking $$\mathfrak{C_{n+1}}=M\mathfrak{C_n}$$ for some matrix $$M$$. This means that $$\mathfrak{C}_{n+k}=M^k\mathfrak{C}_n$$, and so by computing large powers of $$M$$ (for which the usual binary method can be used), large values of $$c_n$$ can be computed. If your $$v$$ is large then the matrices may be large, but since all your work can be done mod $$p$$, then none of the coefficients should get out of hand.
• @user2249675 What are the specific parameters on your problem? You mentioned values of $k_2-k_1$ but it'd be good to have that directly translated back into (rough) $u$ / $v$ and $n$ ranges. – Steven Stadnicki Mar 13 at 20:50
• Thank you. I have solved the original problem (ProjectEuler 258). In that problem $u = 1999,v = 2000,n = 10^{18}$. The trick is to reduce matrix multiplication via Cayley-Hamilton theorem. – user2249675 Mar 13 at 20:56
• It seems at first blush that using Cayley-Hamilton is basically just equivalent to solving the original problem, though it's very possible I'm missing something there. In any case, though, an $n$ of $10^{18}$ means about 60 squaring operations of the matrix, and kerrywong.com/2009/03/07/matrix-multiplication-performance-in-c suggests that a $2000\times2000$ multiply in naive code takes a bit less than a minute, so even the relatively naive approach here is pretty feasible in the grand scheme of things. – Steven Stadnicki Mar 13 at 21:47