Determining coefficient of x^m term in (x^2 + x + 1)^n is even or odd

Asked 29/4, 2017 at 9:46 Answered 29/4, 2017 at 16:57

Solved algorithm math algebra number-theory modulo

For a given integers n and m, determine that coefficient of x^m term in (x^2+x+1)^n is even or odd?

For example, if n=3 and m=4, (x^2+x+1)^3 = x^6 + 3x^5 + [[6x^4]] + 7x^3 + 6x^2 + 3x + 1, so coefficient of x^4 term is 6 (=even).

n and m is as large as 10^12 and I want to calculate in a few seconds, so you can't calculate in linear time.

Do you have any efficient algorithm?

Dismember answered 29/4, 2017 at 9:46 Comment(4)

What do you mean with x^m coefficient of? – Periodical 29/4, 2017 at 9:48

@Periodical Added example. – Dismember 29/4, 2017 at 9:53

@VincentvanderWeele Really? It is not everything is mathematics. It's a programming question because I want to calculate it FASTER, not 10 or 20 minutes in one pair of value. It is an algorithm question. – Dismember 29/4, 2017 at 9:56

If the question is asked like that, you can be pretty sure that you don't have to calculate the actual coefficient. I'd probably start with converting everything to binary and look for which terms I can eliminate straight away. – Medlock 29/4, 2017 at 10:46

Yes, linear time in the number of bits in the input.

The coefficients in question are trinomial coefficients T(n, m). For binomial coefficients, we would use Lucas's theorem; let's work out the trinomial analog for p = 2.

Working mod 2 and following the proof of Nathan Fine,

(1 + x + x^2)^{2^i} = 1 + x^{2^i} + x^{2^{2 i}}

(1 + x + x^2)^n
    = prod_i ((1 + x + x^2)^{2^i n_i})
        where n = sum_i n_i 2^i and n_i in {0, 1} for all i
        (i.e., n_i is the binary representation of n
    = prod_i (1 + x^{2^i n_i} + x^{2^i 2 n_i})
    = prod_i sum_{m_i = 0}^{2 n_i} x^{2^i}
    = sum_{(m_i)} prod_i x^{2^i m_i}
        taken over sequences (m_i) where 0 ≤ m_i ≤ 2 n_i.

In the binomial case, the next step is to observe that, for the coefficient of x^m, there's at most one choice of (m_i) whose x^{2^i m_i} factors have the right product, i.e., the binary representation of m.

In the trinomial case, we have to consider binary pseudo-representations (m_i) of m where pseudo-bits can be zero, one, or two. There is a contribution to the sum if and only if for all i such that n_i = 0, we have m_i = 0.

We can write an automaton that scans n and m bit by bit. State a is initial and accepting.

a (0:0:nm') -> a nm'    [emit 0]
a (1:0:nm') -> a nm'    [emit 0]
            -> b nm'    [emit 2]
a (1:1:nm') -> a nm'    [emit 1]

b (0:1:nm') -> a nm'    [emit 0]
b (1:0:nm') -> b nm'    [emit 1]
b (1:1:nm') -> a nm'    [emit 0]
            -> b nm'    [emit 2]

We can use dynamic programming to count the paths. In code form:

def trinomial_mod_two(n, m):
    a, b = 1, 0
    while m:
        n1, n = n & 1, n >> 1
        m1, m = m & 1, m >> 1
        if n1:
            if m1:
                a, b = a ^ b, b
            else:
                a, b = a, a ^ b
        elif m1:
            a, b = b, 0
        else:
            a, b = a, 0
    return a

Branchless version for giggles:

def trinomial_mod_two_branchless(n, m):
    a, b = 1, 0
    while m:
        n1, n = n & 1, n >> 1
        m1, m = m & 1, m >> 1
        a, b = ((n1 | ~m1) & a) ^ (m1 & b), ((n1 & ~m1) & a) ^ (n1 & b)
    return a

Accomplished answered 29/4, 2017 at 16:57 Comment(2)

Very clever, and I wish I'd thought of this approach. I think there must be a clearer way to explain it though. Something about counting ways to construct m from adding some combination of the bits of n (or those bits shifted left once). – Acuna 29/4, 2017 at 19:48

Wow! Very nice solution! – Dismember 30/4, 2017 at 7:10

First note, that if one is only interested in whether the coefficient of x^m is odd or even, one can consider the coefficients of the polynomial to be elements of the finite field F2.

Note (1+x+x^2)^2 = (1+x^2+x^4) mod 2 because cross terms are all even. In fact, if n is a power of 2, then (1+x+x^2)^n = (1 + x^n + x^2n) mod 2.

For general n, write it as a sum of powers of 2 (that is, in binary).

n = (2^a1 + 2^a2 + 2^a3 + ... + 2^ak).

Then multiply together the powers corresponding to each power of 2:

(1+x+x^2)^n = (1+x^(2^a1)+x^(2^(a1+1))) * ((1+x^(2^a2)+x^(2^(a2+1))) * ...

Each of the terms in this product now has only 3 factors, and there's at most 35 or 36 of them if n is bounded by 10^12. So it's easy to multiply them together.

Here's some Python code that does this:

# poly_times computes the product of polynomials
# p and q over the field F2. They are each
# represented by a set of non-zero coefficients.
# For example set([0, 2, 5]) corresponds to x^0 + x^2 + x^5.
def poly_times(p, q):
    result = set()
    for i in p:
        for j in q:
            if i+j in result:
                result.remove(i+j)
            else:
                result.add(i+j)
    return result

# Return the coefficient of x^m in (1+x+x^2)^n.
def coeff(n, m):
    prod = set([0])
    i = 0
    while n:
        if n % 2:
            prod = poly_times(prod, [0, 2**i, 2**(i+1)])
        i += 1
        n >>= 1
    return 1 if m in prod else 0

for m in xrange(10):
    print m, coeff(3, m)

print coeff(1947248745, 1947248745034)

Optimization

For n with a large number of bits set, this gets too slow as n approached 10^12.

But, one can speed things up greatly by splitting the polynomial power into two parts, and then finding the coefficient of m in the last step not by doing a full polynomial multiplication but instead by counting pairs of coefficient in each part which sum to m. Here's the optimized coeff:

# poly_times computes the product of polynomials
# p and q over the field F2. They are each
# represented by a set of non-zero coefficients.
# For example set([0, 2, 5]) corresponds to x^0 + x^2 + x^5.
# Coefficients larger than m are discarded.
def poly_times(p, q, m):
    result = set()
    for i in p:
        for j in q:
            if i + j > m:
                continue
            if i+j in result:
                result.remove(i+j)
            else:
                result.add(i+j)
    return result

# Return the coefficient of x^m in (1+x+x^2)^n.
def coeff(n, m):
    if m > 2*n: return 0
    prod = [set([0]), set([0])]
    i = 0
    while n:
        if n % 2:
            prod[i//20] = poly_times(prod[i//20], [0, 2**i, 2**(i+1)], m)
        i += 1
        n >>= 1
    s = 0
    for x in prod[0]:
        s += m-x in prod[1]
    return s % 2

for m in xrange(10):
    print m, coeff(3, m)

print coeff(0xffffffffff, 0xffffffffff)

Note that this can compute coeff(0xffffffffff, 0xffffffffff) in a few seconds, and 0xffffffffff is larger than 10**12.

Acuna answered 29/4, 2017 at 11:37 Comment(5)

Hack case: n=2^35-1, m=2^35-1 (I think the number of odd coefficient will be much) – Dismember 29/4, 2017 at 11:48

It's definitely linear in the worst case, but it's pretty difficult to calculate since if n has lots of bits set then there's lots of common terms. But I thought the question was about finding a solution that ran in a few seconds rather than about computational complexity? I guessed from the question that it was an online judge type question. If it's not fast enough, there's plenty of scope for optimizing this solution further even beyond the couple of suggestions I made in the question. – Acuna 29/4, 2017 at 11:58

I think you're right that it's too slow when n has many bits set. – Acuna 29/4, 2017 at 12:10

I have an optimization which makes it plenty fast. – Acuna 29/4, 2017 at 12:12

OK, now it takes 4 seconds for n=2^40-1, m=2^40-1. – Acuna 29/4, 2017 at 12:23

Yes, linear time in the number of bits in the input.

The coefficients in question are trinomial coefficients T(n, m). For binomial coefficients, we would use Lucas's theorem; let's work out the trinomial analog for p = 2.

Working mod 2 and following the proof of Nathan Fine,

(1 + x + x^2)^{2^i} = 1 + x^{2^i} + x^{2^{2 i}}

(1 + x + x^2)^n
    = prod_i ((1 + x + x^2)^{2^i n_i})
        where n = sum_i n_i 2^i and n_i in {0, 1} for all i
        (i.e., n_i is the binary representation of n
    = prod_i (1 + x^{2^i n_i} + x^{2^i 2 n_i})
    = prod_i sum_{m_i = 0}^{2 n_i} x^{2^i}
    = sum_{(m_i)} prod_i x^{2^i m_i}
        taken over sequences (m_i) where 0 ≤ m_i ≤ 2 n_i.

We can write an automaton that scans n and m bit by bit. State a is initial and accepting.

a (0:0:nm') -> a nm'    [emit 0]
a (1:0:nm') -> a nm'    [emit 0]
            -> b nm'    [emit 2]
a (1:1:nm') -> a nm'    [emit 1]

b (0:1:nm') -> a nm'    [emit 0]
b (1:0:nm') -> b nm'    [emit 1]
b (1:1:nm') -> a nm'    [emit 0]
            -> b nm'    [emit 2]

We can use dynamic programming to count the paths. In code form:

def trinomial_mod_two(n, m):
    a, b = 1, 0
    while m:
        n1, n = n & 1, n >> 1
        m1, m = m & 1, m >> 1
        if n1:
            if m1:
                a, b = a ^ b, b
            else:
                a, b = a, a ^ b
        elif m1:
            a, b = b, 0
        else:
            a, b = a, 0
    return a

Branchless version for giggles:

def trinomial_mod_two_branchless(n, m):
    a, b = 1, 0
    while m:
        n1, n = n & 1, n >> 1
        m1, m = m & 1, m >> 1
        a, b = ((n1 | ~m1) & a) ^ (m1 & b), ((n1 & ~m1) & a) ^ (n1 & b)
    return a

Accomplished answered 29/4, 2017 at 16:57 Comment(2)

Wow! Very nice solution! – Dismember 30/4, 2017 at 7:10

The coefficient of interest depends on the number of ways n terms can be selected from x² + x + 1 so that the sum of the powers of the selected terms is m. These ways can be grouped in groups that have the same number of selected x² terms and x terms (the number of times 1 is selected follows from that).

Let a be the number of x² terms, and b the number of x terms, and c the number of 1 terms in a particular group.

Then the following equalities hold:

2a + b = m
a + b + c = n

Obviously there are in general several groups with different values for a, b, c. Once a is determined, the values for b and c are also determined. One needs thus to only iterate over the possible values for a to get all groups.

If you were to write a brute force algorithm to get the coefficient itself, it would look like this in Python:

def combi(n, k): # number of ways to take k elements from n elements
    import math
    f = math.factorial
    return f(n) // f(k) // f(n-k)    

def get_coeff(n, m):
    if m > n * 2 or n < 0 or m < 0: # basic argument check
        return None
    if m > n: # mirrored situation is the same
        m = 2*n - m            
    coeff = 0
    for a in range(0, m//2+1):
        b = m - 2*a
        coeff += combi(n, a) * combi(n-a, b)
    return coeff

The function combi(n, k) would return the number of ways to take k elements from n elements, i.e. the binomial coefficient.

The product of the two calls to combi answers the following question:

In how many ways can I take a times the x² term and b times the x term? Note that the number of ways that the constant term can be taken is 1 once the other 2 choices have been made.

Now, since we only need to know whether the final coefficient is odd or even, we also only need to know whether the binomial coefficient is odd or even. As explained on math.stackexchange.com this can be determined with a simple bit operation:

def combi_parity(n, k):
    return (k & (n-k)) == 0

def get_coeff_parity(n, m):
    if m > n * 2 or n < 0 or m < 0: # basic argument check
        return None
    if m > n:
        m = 2*n - m # mirrored situation is the same
    coeff_odd = 0
    for a in range(0, m//2+1):
        b = m - 2*a
        # A product is odd only when both factors are odd (so we perform an "and")
        # A sum changes parity whenever the added term is odd (so we perform a "xor")
        coeff_odd ^= combi_parity(n, a) and combi_parity(n-a, b) 
    return coeff_odd

Periodical answered 29/4, 2017 at 14:0 Comment(6)

Upvoted. I tried this approach too, but couldn't get it fast enough for large m,n. I even optimized a little so that a ranges over only numbers with combi_parity(n, a)==1. Here's my C code, which is still running on my laptop: gist.github.com/paulhankin/5f5aa05231ce8659b6b040d849fd800e . – Acuna 29/4, 2017 at 14:8

Your formulation of c(n, k) is equivalent but better than the one I used: n&k==k. I'm going to steal it for my program! – Acuna 29/4, 2017 at 14:13

You're not using large numbers in your test case: 2^40-1 should be 2**40-1. Sorry, I mixed up python and math notation. 2^40-1 is 37 ;( – Acuna 29/4, 2017 at 14:17

Ahhhhhh yes, thanks! I updated the sample call in the referenced link. – Periodical 29/4, 2017 at 14:19

My nested loop is actually an (attempted) optimization that steps a (which is called k in my code) over only numbers for which combi_parity(n, a) is True. – Acuna 29/4, 2017 at 14:19

OK, I missed that. – Periodical 29/4, 2017 at 14:19

Okay, I just thought of a solution. Here goes:

Think of the equation as written n times, (a.x^2 + b.x^1 + c).(a.x^2 + b.x^1 + c)...n times. a, b and c are constants I assumed in general.
Now, we have to pick a term out of each so that the result of multiplication of all such terms result in x^m
I can say now that, we have to find solutions of the equation, t1.2 + t2 = m where t1 is no of occurrence of x^2 and t2 of x. This is because t1 and t2 then would make the term to be of the form k.x^m(k is constant). This is finding integral Diophantine solutions of this equation, that is finding all satisfying pairs of {t1, t2}
Now, we have to apply a bit of permutation for finding coefficient here. Assuming you have one of the solution {1, 2} for previous step then for this diad, coefficient would be (1^1.nC1).(2^2.(n-1)C2) which would be one of the coefficient constituent. If you sum all such coefficient terms corresponding to all the Diophantine solutions you would get the coefficient.

Implementing the above algorithmically would take some time for me so I have posted the steps.

Note: I searched a bit and there are various algorithms for Diophantine solutions. Here is one post related to that: How to solve Linear Diophantine equations in programming?

EDIT: As asked for an example,

Lets say, we have the equation, (x^2 + x^1 + x^1)^3 and we have to find coefficient of x^3. Therefore, we have m = 3.
Separately writing the equation is for visually seeing the step. It is,

(x^2 + x^1 + x^1).(x^2 + x^1 + x^1).(x^2 + x^1 + x^1)
Now, we want to pick either of {x^2, x^1, 1} from each of these. There will be several ways to pick it out to make the multiplication of the form, x^3
To solve this, we can write the equation, 2.a + b = 3 where a is no of times x^2 is picked and b is no of times x is picked. The solutions of this equation are, {0, 3} and {1, 1}. Now because, we have to also account for the order in which we pick them, we shall apply combinatorics in the next step.
The coefficient would be, 2^0.3C0.3^3.3C3 + 2^1.3C1.3^1.2C1. Now here, in the first term, 2^0.3C0.3^3.3C3, 3C0 means picking 0 occurences of x^2 with 3C3 means 3 occurrences of x which would give x^3 but we also multiply with 2^0 because 2 is the coefficient of x^2 in the equation and similarly, 3^3 because 3 is the coefficient of x. Similarly, you go about the second term corresponding to {1, 1}
This adds up to 63 which you can verify by manually multiplying and you will get 63.

Hope I am clear.

Neom answered 29/4, 2017 at 10:41 Comment(5)

I don't understand. Can you give me an example? – Dismember 29/4, 2017 at 11:16

@Dismember : I am preparing the example. Give me a few minutes. – Neom 29/4, 2017 at 11:47

@Dismember : See the example. – Neom 29/4, 2017 at 12:26

But if n and m is large you can't brute force {a, b}. – Dismember 29/4, 2017 at 13:39

Summarizing this method: compute sum(c(n, k)c(n-k, m-2k) for k=0..m/2). You can make it pretty fast by computing it mod 2, and noting that c(a, b)=1 mod 2 only if a&b==b (so you only have to iterate over k with bits a subset of the bits of n). But I couldn't get a C version to run fast enough when n and m are large. If you want to play with my code: gist.github.com/paulhankin/5f5aa05231ce8659b6b040d849fd800e – Acuna 29/4, 2017 at 13:54

Hot tags

Godot Unity Godot Help Programming Godot 4.X GUI GDScript 3D 2D Physics CSharp Godot 3.X VR XR Projects C++

Optimization

Recommended topics

Hot tags