numpy 多项式系数的极大递推

mpbci0fu 于 2023-03-30 发布在其他

关注(0)|答案(2)|浏览(104)

最近，我在YouTube上遇到了一个有趣的计数问题，由3Blue1Brown-OlympiadLevelCounting提出。这个问题是找到集合{1，2，...，2000}中元素之和能被5整除的子集的数量。Grant Sanderson提出了一个漂亮的解决方案，但他也讨论了如果在计算机上实现以有效地解决这个问题，算法会变得多么繁琐。
然而，由于欧拉的工作，使用整数分区来解决这个问题是很好的。解决方案是所有小于或等于2000的5的倍数的distinct partitions的总和，这与我多年前在MathstackExchange上发布的内容类似，其中使用“dxiv”的注解是此代码的关键。
毫无疑问，Grant Sanderson的解决方案是优雅的，因为它不需要计算机，并且可以在没有冗长计算的情况下找到。
我的代码是：

import numpy as np

def p(k):
    if k == 1:
        return np.array([1, 1])
    else:
        q = p(k - 1)
    return add(q, np.append(np.zeros(k, int), q))

def add(u, v):
    u = np.append(u, np.zeros(len(v) - len(u), int))
    return np.add(u, v)

def test(n):
    a = 1 << n
    b = 1 << (n//5)
    b *= 4
    return (a+b)//5

n = 1000      
generating_function = sum(p(n)[0:-1:5])+1
grant_sanderson = test(n)
print(generating_function)
print(grant_sanderson)
print(generating_function == grant_sanderson)

我的困难是当n是一个很大的数字时。而且，这些数字太大了，我必须使用gmpy2这是一个图表，显示数组关于n//2对称：

Traceback (most recent call last):
  File "c:\Users\patil\OneDrive\Documents\GitHub\Theorie-der-Zahlen\Distinct partitions\distinctParitions-gmp2\main.py", line 15, in <module>
    generating_function = sum(p(n)[0:-1:5])+1
                              ^^^^
  File "c:\Users\patil\OneDrive\Documents\GitHub\Theorie-der-Zahlen\Distinct partitions\distinctParitions-gmp2\generatingfunction.py", line 8, in p
    q = p(k - 1)
        ^^^^^^^^
  File "c:\Users\patil\OneDrive\Documents\GitHub\Theorie-der-Zahlen\Distinct partitions\distinctParitions-gmp2\generatingfunction.py", line 8, in p
    q = p(k - 1)
        ^^^^^^^^
  File "c:\Users\patil\OneDrive\Documents\GitHub\Theorie-der-Zahlen\Distinct partitions\distinctParitions-gmp2\generatingfunction.py", line 8, in p
    q = p(k - 1)
        ^^^^^^^^
  [Previous line repeated 996 more times]
RecursionError: maximum recursion depth exceeded

这段代码可以说明事情是如何迅速失控的：

import gmpy2

def q(n):
    Q = [gmpy2.mpz(0)] * (n+1)
    Q[0] = gmpy2.mpz(1)
    for k in range(1, n):
        for i in range(n, k-1, -1):
            Q[i] += Q[i-k]
    return Q[n] + 1

numpy

来源：https://stackoverflow.com/questions/75835358/the-max-recursion-for-the-coefficient-of-a-polynomial

2条答案

按热度按时间

bxgwgixi1#

另一种方法是计算达到5个模值（0，1，2，3，4）的总和的组合。将这些计数保存在列表中并逐渐增加它们。最后，将零取模的计数计算为总和为5的倍数的集合的数量：

def sumDiv5Count(numbers):
    modCounts = [0]*5
    for n in numbers:
        modCounts = [c+modCounts[(i-n)%5] for i,c in enumerate(modCounts)]
        modCounts[n%5] += 1        
    return modCounts[0]

每个数的模5与5个可能的模组合以添加更多的方式来产生5个不同的结果模值。
例如13%5 -〉3，

生成0的方法被添加到生成0+3的方法--〉插槽3
将产生1的路添加到产生1+3 --〉时隙4的路
生产2的方法被添加到生产2+3的方法--〉插槽0（2+3）%5 = 0
生产3的方法添加到生产3+3的方法--〉插槽1（3+3）%5 = 1
生产4的方法被添加到生产4+3的方法--〉插槽2（4+3）%5 = 2
还有一种方法可以自己生产3

输出：

print(sumDiv5Count(range(1,10)))
# 103

print(sumDiv5Count(range(1,50)))
# 112589990684671

print(sumDiv5Count(range(1,100)))
# 126765060022822940149670739967

print(sumDiv5Count(range(1,500)))
# 327339060789614187001318969682759915221664204604306478948329136809613379640467455488327009232590415715088668412756007101428785894679831065931534041087

print(sumDiv5Count(range(1,1000)))
# 1071508607186267320948425049060001810561404811705533607443750388370351051124936122493198378815695858127594672917553146825187145285692314043598457757469857480393456777482423098542107460506237114187795418215304647498358194126739876755916554395250481509160715757985439053702508024174143636196837700927487

print(sumDiv5Count(range(1,2001)))
# 22962613905485090484656664023553639680446354041773904009552854736515325227847406277133189726330125398368919292779749255468942379217261106628518627123333063707825997829062456000137755829648008974285785398012697248956323092729277672789463405208093270794180999311632479761788925921124662329907232844394066536268833781796891701120475896961582811780186955300085800543341325166104401626447256258352253576663441319799079283625404355971680808431970636650308177886780418384110991556717934409897816293912852988275811422719154702569434391547265221166310540389294622648560061463880851178273858239474974548427800575

这种方法的一个优点是，它可以很容易地变成一个生成器，当它通过数字时，它会逐渐产生计数。它也不依赖于数字是否在序列中，甚至是不同的。考虑到它是一个O（n）过程，性能很好（几毫秒）。

赞(0）回复(0）举报 2023-03-30

u1ehiz5o2#

这不会抛出递归错误！对原始代码所做的更改：

sys.setrecursionlimit(2500)
不知道为什么，但是gmpy2抛出了错误的数字。所以使用了decimal.Decimal

from decimal import Decimal, getcontext
import sys
sys.setrecursionlimit(2500)

def p(k):
    if k == 1:
        return [Decimal(1), Decimal(1)]
    else:
        q = p(k - 1)
        return add(q, [Decimal(0)] * k + q)

def add(u, v):
    u = u + [Decimal(0)] * (len(v) - len(u))
    return [x + y for x, y in zip(u, v)]



def test(n):
    getcontext().prec = n+10  
    a = Decimal(1 << n)
    b = Decimal(1 << (n//5))
    b *= Decimal(4)
    result = (a+b)/Decimal(5)
    return Decimal(result)  

n = 2000
generating_function = sum(p(n)[0:-1:5])+1
print(generating_function)

赞(0）回复(0）举报 2023-03-30

我来回答

numpy 多项式系数的极大递推

2条答案

相关问题

热门标签

最新问答