确定 (x^2 + x + 1)^n 中 x^m 项的系数是偶数还是奇数

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

对于给定的整数nm,确定(x^2+x+1)^nx^m项的系数是偶数还是奇数?

例如,如果 n=3 且 m=4,(x^2+x+1)^3 = x^6 + 3x^5 + [[6x^4]] + 7x^3 + 6x^2 + 3x + 1,则 x^4 项的系数为 6(=偶数)。

nm 有10^12那么大,我想在几秒内计算,所以你不能用线性时间计算。

你有什么高效的算法吗?

好吧,我刚刚想到了解决办法。这里是:

  1. 将等式想象成 n 次,(a.x^2 + b.x^1 + c).(a.x^2 + b.x^1 + c)...n 次。 a、b 和 c 是我一般假设的常数。
  2. 现在,我们必须从每一项中选出一项,以便所有这些项相乘的结果为 x^m
  3. 我现在可以说,我们必须找到方程的解,t1.2 + t2 = m 其中 t1 是没有出现 x^2t2x。这是因为 t1t2 会使项的形式为 k.x^m(k 是常数)。这是找到这个方程的积分丢番图解,也就是找到所有令人满意的 {t1, t2}
  4. 现在,我们必须在这里应用一些排列来找到系数。假设您有上一步的解决方案 {1, 2} 之一,那么对于这个二元组,系数将是 (1^1.nC1).(2^2.(n-1)C2),这将是系数成分之一。如果对所有丢番图解对应的所有此类系数项求和,您将得到系数。

通过算法实现上述算法对我来说需要一些时间,所以我已经 post 编辑了这些步骤。

注意:我搜索了一下,丢番图解有多种算法。这是一个与此相关的post:How to solve Linear Diophantine equations in programming?

编辑:举个例子,

  1. 比方说,我们有方程 (x^2 + x^1 + x^1)^3,我们必须找到 x^3 的系数。因此,我们有 m = 3.
  2. 单独写方程是为了直观看步。它是,

    (x^2 + x^1 + x^1).(x^2 + x^1 + x^1).(x^2 + x^1 + x^1)

  3. 现在,我们要从其中每一个中选择一个 {x^2, x^1, 1}。有几种方法可以挑出来做乘法的形式,x^3

  4. 为了解决这个问题,我们可以写出等式,2.a + b = 3 其中 a 是 x^2 被选取的次数,b 是 x 被选取的次数。这个方程的解是 {0, 3}{1, 1}。现在因为,我们还必须考虑我们选择它们的顺序,我们将在下一步中应用组合学。

  5. 系数为 2^0.3C0.3^3.3C3 + 2^1.3C1.3^1.2C1。现在在这里, 在第一项中,2^0.3C0.3^3.3C33C0 表示选择 0 次出现的 x^2,而 3C3 表示出现 3 次 x,这将给出 x^3 但我们还乘以2^0因为2是方程中x^2的系数,同理,3^3因为3是x的系数。类似地,你去处理对应于 {1, 1}

  6. 的第二项
  7. 这加起来是63,你可以通过手动相乘来验证,你会得到63。

希望我清楚。

首先注意,如果只对x^m的系数是奇数还是偶数感兴趣,可以将多项式的系数视为有限域F2的元素。

注意 (1+x+x^2)^2 = (1+x^2+x^4) mod 2 因为交叉项都是偶数。事实上,如果n是2的幂,那么(1+x+x^2)^n = (1 + x^n + x^2n) mod 2.

对于一般的n,写成2的幂和(即二进制)。

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

然后将2的各个次方对应的次方相乘:

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

此乘积中的每个项现在只有 3 个因子,如果 n 以 10^12 为界,则最多有 35 或 36 个因子。所以很容易将它们相乘。

下面是一些 Python 执行此操作的代码:

# 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)    

优化

对于设置了大量位的 n,当 n 接近 10^12 时,这会变得太慢。

但是,可以通过将多项式幂分成两部分来大大加快速度,然后在最后一步中找到 m 的系数,而不是通过进行完整的多项式乘法,而是通过计算对每个部分的系数总和为 m。这是优化后的 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)

请注意,这可以在几秒钟内计算出 coeff(0xffffffffff, 0xffffffffff),并且 0xffffffffff 大于 10**12。

感兴趣的系数取决于 n 项可以从 x² + x + 1 中选择的方式的数量,以便总和所选项的幂是 m。这些方式可以分组为具有相同数量的选定 项和 x 项(次数 1 被选中)。

a项的数量,bx 项,以及 c 特定组中 1 项的数量。

则下列等式成立:

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

显然,通常有几个组 a、b、c 具有不同的值。一旦 a 被确定,bc 的值也被确定。因此,只需迭代 a 的可能值即可获得所有组。

如果你要编写一个强力算法来获取系数本身,它在 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

函数combi(n, k)将return从n个元素中取出k个元素的方法数,即binomial coefficient.

两次调用 combi 的结果回答了以下问题:

我可以用多少种方法来计算 a 项和 bx 项?请注意,一旦做出其他 2 种选择,常数项的取法数量为 1。

现在,既然我们只需要知道最后的系数是奇数还是偶数,那么我们也只需要知道二项式系数是奇数还是偶数。 正如 math.stackexchange.com 中所解释的,这可以通过简单的位操作来确定:

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

repl.it 上查看 运行。

是的,输入位数的线性时间。

所讨论的系数是三项式系数T(n, m)。对于二项式系数,我们将使用 Lucas's theorem;让我们计算出 p = 2.

的三项式模拟

正在工作 mod 2 并遵循 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.

在二项式情况下,下一步是观察,对于x^m的系数,最多只有一个(m_i)的选择,其x^{2^i m_i}个因子有正确的乘积,即 m.

的二进制表示

在三项式情况下,我们必须考虑 m 的二进制伪表示 (m_i),其中伪位可以是零、一或二。当且仅当对于所有 i 使得 n_i = 0,我们有 m_i = 0.

时,对总和有贡献

我们可以编写一个自动机来逐位扫描nm。状态 a 初始并接受。

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]

我们可以使用动态规划来计算路径。代码形式:

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

咯咯笑的无分支版本:

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