# Efficient computation of number of partitions into powers of 2

Consider the function $$a(n)$$ defined as the number of partitions of $$n$$ into powers of $$2$$. The sequence is given at OEIS.

I am trying to calculate $$a(n)$$ modulo some fixed prime, for large $$n$$ and formulate algorithm that uses $$O(\log n)$$ basic arithmetic operations, or equivalently only $$O(\log n)$$ evaluations of the function $$a(\cdot)$$.

$$a(n)$$ has following recursive relations:

Recursive relation 1: \begin{align} a(2k+1)&=a(2k)\\ a(2k)&=a(2k-1)+a(k) \end{align}

Recursive relation 2: \begin{align} a(n)&= \begin{cases} 1, & n=0 \\ \sum_{i=0}^{\lfloor n/2\rfloor} a(i), & n>0 \end{cases} \end{align}

Both of these relations need $$O(n)$$ computations. We can modify the second recurrence and can obtain another recurrence relation for $$a(16k+r)~\forall ~ r\in[0,15]$$. For examples: \begin{align} a(16k) &= a(8k) + a(8k+2) + (4k-2) a(4k+2) + \sum_{i=1}^{k-1} (8k-2-8i) a(2k+2i)\\ a(16k+2) &= a(8k+2) + 4k a(4k+2) + \sum_{i=1}^{k} (8k+2-8i) a(2k+2i) \end{align}

This relation also requires $$O(n)$$ computations. We can find similar recurrence relations of the form $$a(2^mk+r)$$, but still that will be $$O(n)$$ operations. Is it really possible to do better computationally?

• Have you checked the links & references & code in the OEIS entry (and for oeis.org/A000123)? If there is no fast algorithm listed there, then odds are that no fast algorithm is known.
– D.W.
Commented Jul 9 at 6:30
• I think what @Somos has described in the PROG section of the OEIS sequence, is similar to what I have done below, but I don't know PARI so can't be sure.
– EnEm
Commented Jul 10 at 10:37
• I'm guessing this is intended to solve Project Euler Problem 890. It's against the guidelines of Project Euler for solutions to the problems to be openly published, so I think it would be better if this question could somehow be removed. It might be too late now, though. Commented Jul 11 at 0:11

Below is an $$\mathcal{O}(\log^4 N)$$ time algorithm to compute $$a(N)$$ modulo a fixed prime $$p$$.

You can also improve the algorithm to be in $$\mathcal{O} \left( \left( \log N \cdot \log \log N \right)^2 \right)$$ time using NTT to do the evaluations and interpolation (algorithms to do that described here). An extra condition required here is that $$p$$ is NTT-able for polynomials of degree $$\approx 2\lfloor \log_2 N \rfloor + 2$$, i.e., if $$p-1$$ is factorised into $$2^rs$$ where $$s$$ is odd, then $$2^r \ge 4\lfloor \log_2 N \rfloor + 6$$. An example of a prime $$\approx10^9$$ which is NTT-able for very large polynimials is $$p=998244353$$, where $$p-1=2^{23}\times119$$.

#### Algorithm description

First I define a transformation operation on a given finite sequence $$b = (b[0], b[1], \dots b[n])$$ as $$f(b) = \left(\sum_{i=0}^{n} b[i], \sum_{i=2}^{n} b[i], \sum_{i=4}^{n} b[i], \dots \sum_{i = 2 \left\lfloor \frac{n}{2} \right\rfloor}^{n} b[i] \right)$$

Using the recurrence relation $$2$$ in the original post, its easy to see that $$\begin{split} \sum_{i=0}^{n} b[i] \times a(i) &= \sum_{i=0}^n \left(b[i] \times \sum_{j=0}^{\left\lfloor \frac{i}{2} \right\rfloor} a(j) \right) \\ &= \sum_{i=0}^n \sum_{j=0}^{\left\lfloor \frac{i}{2} \right\rfloor} \left(b[i] \times a(j)\right) \\ &= \sum_{j=0}^{\left\lfloor \frac{n}{2} \right\rfloor} \sum_{i=2j}^n \left(b[i] \times a(j) \right) \\ &= \sum_{j=0}^{\left\lfloor \frac{n}{2} \right\rfloor} \left(\left(\sum_{i=2j}^n b[i]\right) \times a(j) \right) \\ \end{split}$$ or $$\sum_{i=0}^{n} b[i] \times a(i) = \sum_{i=0}^{\left\lfloor \frac{n}{2} \right\rfloor} f(b)[i] \times a(i)$$

This means to evaluate $$a(N)$$, we can evaluate the series of sequences $$b_0, b_1, \dots b_{\lfloor \log_2 N \rfloor + 1}$$, where $$b_0 = (0, 0, \dots 0, 1)$$ of size $$N+1$$, and $$b_{i+1} = f(b_i)$$ for each $$i \in \{0, 1, \dots \lfloor \log_2 N \rfloor\}$$. Then because $$a(N) = \sum b_0[i]\cdot a(i) = \sum b_1[i]\cdot a(i) \dots = \sum b_{\lfloor \log_2 N \rfloor + 1}[i]\cdot a(i)$$, and $$b_{\lfloor \log_2 N \rfloor + 1}$$ is just of size $$1$$, we get $$A(N) = b_{\lfloor \log_2 N \rfloor + 1}[0]$$.

Showing an example here, for $$N = 19$$, we get the following series of sequences:

• $$b_0 = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1)$$
• $$b_1 = (1, 1, 1, 1, 1, 1, 1, 1, 1, 1)$$
• $$b_2 = (10, 8, 6, 4, 2)$$
• $$b_3 = (30, 12, 2)$$
• $$b_4 = (44, 2)$$
• $$b_5 = (46)$$

So we get $$a(19) = 46$$.

Now the key to optimising in time, is to notice that for each $$b_i$$ for $$i>0$$, we can describe it as the evaluations of a $$i-1$$ degree polynomial, i.e., there exists a polynomial $$p_i(x)$$ of degree $$i-1$$, s.t., $$b_i[j] = p_i(j)$$ for all valid $$j$$. For the above example $$b_2$$ has the corresponding polynomial $$p_2(x) = 10-2x$$.

We can prove this via induction for all $$N$$: base step $$b_1$$ has polynomial $$p_1(x) = 1$$, and if $$p_i(x)$$ is expressed as $$\sum c_j x^j$$, then for any transformation $$\begin{split} p_{i+1}(x) &= \sum b_i[j] - \sum_{y=0}^{2x-1} p_i(y) \\ &= C - \sum \left( c_j \sum_{y=0}^{2x-1} y^j \right) \\ \end{split}$$ And we know $$\sum_{y=0}^{x} y^k$$ has a degree of $$k+1$$.

So we don't actually have to evaluate the whole sequences $$b_i$$, but can work by only storing the $$i$$ coefficients of $$p_i$$. Then the transform will work by evaluating the last $$2(i+1)$$ values of $$b_i$$ (using $$p_i$$), using them to generate the last $$i+1$$ values of $$b_{i+1}$$, and then interpolating these values to form $$p_{i+1}$$.

Evaluation takes $$\mathcal{O}(n^2)$$ time, and interpolation $$\mathcal{O}(n^3)$$ time for a $$n$$-degree polynomial, which can be both improved to $$\mathcal{O}(n \log^2 n)$$ times using NTT. As our polynomials grow to a max size of $$\log_2 N$$ degree, and we do the transformations a maximum of $$\log_2 N$$ times, we get the time complexities described above.

#### Implementation

I implemented the above idea in C++, to verify my thinking. It seems to work for all the listed values in OEIS, but do let me know if it starts to crack above that.

#include <iostream>
#include <vector>
using namespace std;

class NumberOf2PowerPartitionsModuloP {
private:
int p;

int pow(int x, int y) {
int t = 1;
while (y) {
if (y & 1) t = 1LL * x * t % p;
x = 1LL * x * x % p;
y >>= 1;
}
return t;
}

int inverse(int x) { return pow(x, p - 2); }

void checkPrime() {
for (int x = 2; x * x <= p; x++) {
if (p % x == 0) {
cout << p << " is not prime!!";
exit(1);
}
}
}

public:
NumberOf2PowerPartitionsModuloP(int p) {
this->p = p;
checkPrime();
}

int compute(int n) {
if (n <= 1) return 1;

int currentSize = (n >> 1), currentSelected = 0;
vector<int> polynomial[2] = {vector<int>({1}), vector<int>()};
vector<int> evaluations;

while (currentSize > 0) {
int nv = polynomial[currentSelected].size();

for (int x = evaluations.size(); x <= 2 * nv + (currentSize & 1); x++) {
int t = 1, sm = 0;
for (int i = 0; i < nv; i++) {
sm = (sm + 1LL * polynomial[currentSelected][i] * t) % p;
t = 1LL * x * t % p;
}
evaluations.push_back(sm);
}

int y = 1;
if (currentSize & 1) {
evaluations[0] = (evaluations[0] + evaluations[1]) % p;
y = 2;
}
for (int x = 1; x <= nv; x++) {
evaluations[x] =
(0LL + evaluations[x - 1] + evaluations[y] + evaluations[y + 1]) %
p;
y += 2;
}
evaluations.resize(nv + 1);

currentSelected ^= 1;
currentSize >>= 1;

polynomial[currentSelected].clear();
polynomial[currentSelected].resize(nv + 1, 0);
for (int i = 0; i <= nv; i++) {
vector<int> tmpPoly({1});
for (int j = 0; j <= nv; j++) {
if (i == j) continue;
tmpPoly.push_back(0);
for (int k = tmpPoly.size() - 1; k >= 1; k--) {
tmpPoly[k] = (tmpPoly[k - 1] + p - (1LL * tmpPoly[k] * j % p)) % p;
}
tmpPoly[0] = (p - (1LL * tmpPoly[0] * j % p)) % p;
}

int coefficient = evaluations[i];
for (int j = 0; j <= nv; j++) {
if (i == j) continue;
coefficient = 1LL * coefficient * inverse((i + p - j) % p) % p;
}

for (int j = 0; j <= nv; j++) {
polynomial[currentSelected][j] =
(polynomial[currentSelected][j] +
1LL * tmpPoly[j] * coefficient % p) %
p;
}
}
}

return evaluations[0];
}
};

int main() {
int n, p;
cin >> n >> p;
cout << NumberOf2PowerPartitionsModuloP(p).compute(n);
}