4
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – D.W.
    Commented Jul 9 at 6:30
  • $\begingroup$ 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. $\endgroup$
    – EnEm
    Commented Jul 10 at 10:37
  • 2
    $\begingroup$ 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. $\endgroup$
    – Ege Erdil
    Commented Jul 11 at 0:11

1 Answer 1

1
$\begingroup$

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);
}
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .