取法有几种（整数拆分）

lrh2006

从1--999这999个自然数中，任意取出n个相加（可以重复取），使其和为1000，有几种取法？
这个问题有解吗？或者从编程的角度讲，算法应该怎么进行？
题目是自己瞎想的，请大家指教，谢谢

业余的业余

回复 1# lrh2006

后来查了一下，这个才是所谓的堆垒数论问题，和丢番图方程有某种联系。我只有题和解，没有书，所以也搞不懂堆垒数论怎么玩。

不过这个问题用编程的思路来接，却是必杀，而且其算法可以简单地套用到纸笔求解。

关键词: 动态规划算法。

这个算法的思路我就不介绍了，本质上要递归地表达问题的解。动态规划是用于合并相同的子问题，从而减少计算工作量。不要小看这个减少，很多时候它是把一个指数/阶乘复杂度的问题降低为一个多项式复杂度的问题，从几乎不可解到可轻松求解。

我们把 正整数 $ n $, 分解为 从 $ k $ 开始， 不大于$ 999 $ 的若干个整数项(允许重复项)的和的方法的个数记为 $ P(n,k) $. 那么显然，你的问题就是找 $ P(1000,1) $.

第一步，我们可以先选定不同个数的$1$， 然后要求剩下的数用$2-999$的项来分解，所以
$\begin{align} \begin{aligned}P(1000,1)=&P(999,2)& \text{取了1个1 ，剩下的交给2和以上的去分割}\\
     +&P(998,2)   &  \text{ 取了2个1， 剩下的交给2和以上的去分割}\\
     &\vdots \\
     +&P(2,2)&\text{ 取了998个1， 剩下的交给2和以上的去分割}\\
     +&1&\text{ 全部取1}\end{aligned}\end{align} $

然后 $ P(999,2), P(998,2), \cdots $ 等可以递归的类似求解。用一个小一点的数，比如 $ 5 $, 在纸上比划下就明白了，用大体 $ 5\times 5 $ 的表格，注意求解的顺序。

我不知道用堆垒数论的方法与之相比，更优雅到何种程度，因为不会。这个问题用电脑就是毫秒级问题，时间和空间复杂度都是 $ O(n^2) $, 用纸笔算还费点事，关键很 boring.

---
PS: 这个算法的描述是有问题的。参考后面的程序。具体而言，上面的主要错误有2
1. P(n,k) 定义为 把 n 用 k...n 的任意多个正整数分解的不同式子的个数。(n 包含在 可能的项之中);
2. 要允许选择0个 k 的情形，上面漏掉了。参后面的程序。

realnumber

回复 2# 业余的业余

编一个，学习学习，问题规模小点，修改成1~20内，和是10 好了.

业余的业余

回复 3# realnumber

:)

用动态规划的变体备忘录法，效率稍差，堆栈要求高些，不过对现代的电脑应该不是问题。

1. #include <iostream>
3. // 用了全局变量，满去吧 :)
4. unsigned map[10001][1000];
6. unsigned P(int n, int k)
7. {
8.       unsigned v=0;
10.       if(n<=k)
11.             return 1;
12.       if(map[n][k]!=0)
13.             return map[n][k];
14.       
15.        for(int j=n-1;j>0; j-=k)
16.             v+=P(j,k+1);
17.       
18.         return map[n][k]=v;
19. }
21. int main()
22. {
23.         std::cout<<" 总取法为 "<<P(1000,1)<<std::endl;
24. }

复制代码

没有测试，可能会有小bug.

业余的业余

好多年不玩程序了，现在改玩数学

kuing

理论上用多项式展开也可以，只是我不清楚其时间和空间复杂度有多大，主要我没学过算法……

其原理倒是非常好理解，设
\[F=\left( \sum_{k=0}^{1000}x^k \right)\left( \sum_{k=0}^{\lfloor1000/2\rfloor}x^{2k} \right)\left( \sum_{k=0}^{\lfloor1000/3\rfloor}x^{3k} \right)\cdots\left( \sum_{k=0}^{\lfloor1000/999\rfloor}x^{999k} \right),\]
那么 `F` 展开式的 `x^{1000}` 的系数就是所求。
理由就是：展开时，第 `m` 个括号所取的项就代表拆分中含有 `m` 的个数。
如果你理解了前两天这帖 kuing.cjhb.site/forum.php?mod=redirect&goto=findpost&ptid=5872&pid=29759 （11#）的解析，这个你应该也能明白。

用 Mathematica 来试一下：

2. f[sum_,max_] := Coefficient[Product[Sum[x^(m k), {k, 0, Floor[sum/m]}], {m, 1, max}], x^sum]

复制代码

那么楼主要求的就是 f[1000,999]，这个太大了，来算小一点的：
和为 10，无最大数限制的（即 10=10 也算一种），那就是 f[10, 10]，结果是 42，如果 10=10 不算那就是 41。
我手工拆了一次来验证此结果：

1. 10 = 10
2. 10 = 9+1
3. 10 = 8+2
4. 10 = 8+1+1
5. 10 = 7+3
6. 10 = 7+2+1
7. 10 = 7+1+1+1
8. 10 = 6+4
9. 10 = 6+3+1
10. 10 = 6+2+2
11. 10 = 6+2+1+1
12. 10 = 6+1+1+1+1
13. 10 = 5+5
14. 10 = 5+4+1
15. 10 = 5+3+2
16. 10 = 5+3+1+1
17. 10 = 5+2+2+1
18. 10 = 5+2+1+1+1
19. 10 = 5+1+1+1+1+1
20. 10 = 4+4+2
21. 10 = 4+4+1+1
22. 10 = 4+3+3
23. 10 = 4+3+2+1
24. 10 = 4+3+1+1+1
25. 10 = 4+2+2+2
26. 10 = 4+2+2+1+1
27. 10 = 4+2+1+1+1+1
28. 10 = 4+1+1+1+1+1+1
29. 10 = 3+3+3+1
30. 10 = 3+3+2+2
31. 10 = 3+3+2+1+1
32. 10 = 3+3+1+1+1+1
33. 10 = 3+2+2+2+1
34. 10 = 3+2+2+1+1+1
35. 10 = 3+2+1+1+1+1+1
36. 10 = 3+1+1+1+1+1+1+1
37. 10 = 2+2+2+2+2
38. 10 = 2+2+2+2+1+1
39. 10 = 2+2+2+1+1+1+1
40. 10 = 2+2+1+1+1+1+1+1
41. 10 = 2+1+1+1+1+1+1+1+1
42. 10 = 1+1+1+1+1+1+1+1+1+1

复制代码

和为 50 的，f[50, 50]，运行一分多钟得出 204226。
试了下 100，等好久都不出来，看来复杂度挺大，算了……

业余的业余

回复 6# kuing

业余的业余

回复 7# 业余的业余

增长的速度比我预期的要快得多，可能内建的整数类型不够用，要用自定义的整数类型了。好吧，我测试下。

业余的业余

回复 8# 业余的业余

kuing 的边界不一样，结果不能比较了。我上面的算法有很大的漏洞，测试修改的版本对小的数字测试通过，不过增长速度远低于kuing 的结果，这个似乎不能用边界的细微调整解释。存疑，童鞋们帮查查错吧。

程序(c++)

2. #include <iostream>
5. const int N=5;
6. // 用了全局变量，满去吧 :)
7. unsigned map[N+1][N];
9. unsigned P(int n, int k)
10. {
11.       unsigned v=0;
13.           if(n==0 || n==k)
14.               return 1;
15.       else if(n<k)
16.             return 0;
17.       if(map[n][k]!=0)
18.             return map[n][k];
19.       
20.        for(int j=n;j>=0; j-=k)
21.             v+=P(j,k+1);
22.       
23.        return map[n][k]=v;
24. }
26. int main()
27. {
28.         std::cout<<" 总取法有 "<<P(N,1)<<"种。"<<std::endl;
29.        
30.             // 调试用途
31.                 for(int i=0; i<N+1; ++i)
32.                         for(int j=0; j<N-1; ++j)
33.                         {
34.                                 std::cout<<map[i][j]<<(j==N-2?'\n':' ');
35.                         }
36. }

复制代码

N=1时 得 1
N=2时 得 2，
N=3时 得 3，
N=4时 得 5，
N=5时 得 7，
..
N=10时，得42
N=1000时， 得2853978871

把调试输出关掉，在程序未优化时，耗时不到一秒。

不排除程序有未发现的subtle error 和整数溢出的可能。发现增长的没有kuing 给出的那么快，还是用了 unsigned 类型。

业余的业余

回复 9# 业余的业余

50 的结果也是 204226. 仔细一看，我用的是和kuing 一样的允许 最大项取到1000. 这样和kuing 的两厢验证，程序大体是对的了。而允许1000项时只有唯一的一种，在这个基础上减1，就得到了贴主要求的结果。

kuing

两个帖子都明白了，下次碰到能不能用起来？

@realnumber 为啥不能用？这些又不是新东西……

PS、修改了下标题。

kuing

再丢个链接：oeis.org/A000041

realnumber

回复 11# kuing
我是说刚理解，对自己能熟练运用表示怀疑.

kuing

唉！我早就应该想到，这么经典的问题，Mathematica 应该已经自带相应命令，一查果然如此！！！

它们是 PartitionsP 和 IntegerPartitions，前者就是拆分方法数（同样包含 10=10 这种），后者则是列出具体的拆分。
参考：
reference.wolfram.com/language/ref/PartitionsP.html
reference.wolfram.com/language/ref/IntegerPartitions.html
比如，PartitionsP[10] 得出 42，而 IntegerPartitions[10] 得出：
{{10},{9,1},{8,2},{8,1,1},{7,3},{7,2,1},{7,1,1,1},{6,4},{6,3,1},{6,2,2},{6,2,1,1},{6,1,1,1,1},{5,5},{5,4,1},{5,3,2},{5,3,1,1},{5,2,2,1},{5,2,1,1,1},{5,1,1,1,1,1},{4,4,2},{4,4,1,1},{4,3,3},{4,3,2,1},{4,3,1,1,1},{4,2,2,2},{4,2,2,1,1},{4,2,1,1,1,1},{4,1,1,1,1,1,1},{3,3,3,1},{3,3,2,2},{3,3,2,1,1},{3,3,1,1,1,1},{3,2,2,2,1},{3,2,2,1,1,1},{3,2,1,1,1,1,1},{3,1,1,1,1,1,1,1},{2,2,2,2,2},{2,2,2,2,1,1},{2,2,2,1,1,1,1},{2,2,1,1,1,1,1,1},{2,1,1,1,1,1,1,1,1},{1,1,1,1,1,1,1,1,1,1}}
排列的方式正和我上面手工列的一模一样，都是由大到小。

PartitionsP[50] 的结果是 204226，也和我们得到的一样。

但是 PartitionsP[1000] 得出的是 24061467864032622473692149727991（楼主想求的就是这个数减 1）那就和 9# 的不一样了。

业余的业余

回复 14# kuing

可能是整数溢出，我换成64位整数试试。

64bit 整数也不行。 64bit二进制 大致表达19位十进制数 (64 * log 2), 你的答案有 32 位了。我写个简单的长整数类试试。算法的逻辑本身有问题的可能微乎其微了。

答案一致。

2. #include <iostream>
3. #include <iomanip>
5. // 违章建筑无符号大整数型，仅供解决 N=1000 的情形. 支持36位十进制数
6. class luint
7. {
8. public:
9.        
10.         luint(unsigned long a=0){ data[0]=a; data[1]=data[2]=data[3]=0;}
11.         luint(const luint& rhs){
12.                 data[0]=rhs.data[0];
13.                 data[1]=rhs.data[1];
14.                 data[2]=rhs.data[2];
15.                 data[3]=rhs.data[3];
16.         }
17.        
18.         luint& operator +=(const luint& rhs)
19.         {
20.                 unsigned long carry=0;
21.                 for(int i=0; i<width; ++i)
22.                 {
23.                         carry+=data[i]+rhs.data[i];
24.                         data[i]=carry%1000000000;
25.                         carry/=1000000000;
26.                 }
27.                 if(carry)
28.                         throw "overflow！";
29.                 return *this;
30.         }
31.        
32.         luint& operator =(const luint& rhs)
33.         {
34.                 return *(luint*)new(this)luint(rhs);
35.         }
36.        
37.        
38.         bool operator == (const luint& rhs)const
39.         {
40.                 return data[0]==rhs.data[0] &&
41.                         data[1]==rhs.data[1] &&
42.                         data[2]==rhs.data[2] &&
43.                         data[3]==rhs.data[3];
44.         }
45. friend std::ostream& operator <<(std::ostream& os, const luint& lu)
46. {
47.         int firstNZ; // first non-zero ulong. from 3 to 0
48.         for(firstNZ=3; firstNZ>=0 && lu.data[firstNZ]==0; --firstNZ)
49.                 ;
50.         if(firstNZ<=0)
51.                 return os<<lu.data[0];
52.        
53.         os<<lu.data[firstNZ];
54.         while(--firstNZ>=0)
55.                 os<<std::setw(9)<<std::setfill('0')<<lu.data[firstNZ];
56.         return os;
57. }
58.        
59. private:
60.         static const int width=4;
61.         unsigned long data[width];
62. };
66. const int N=1000;
67. // 用了全局变量，满去吧 :)
68. luint map[N+1][N];
70. luint one(1),zero;
72. luint& P(int n, int k)
73. {
74.       luint v=0;
76.           if(n==0 || n==k)
77.               return one;
78.       else if(n<k)
79.             return zero;
80.       if(!(map[n][k]==zero))
81.             return map[n][k];
82.       
83.        for(int j=n;j>=0; j-=k)
84.             v+=P(j,k+1);
85.       
86.        return map[n][k]=v;
87. }
89. int main()
90. {
93.         std::cout<<" 总取法有 "<<P(N,1)<<"种。"<<std::endl;
94.        
95.             // 调试用途
96. //                for(int i=0; i<N+1; ++i)
97. //                        for(int j=0; j<N-1; ++j)
98. //                        {
99. //                                std::cout<<map[i][j]<<(j==N-2?'\n':' ');
100. //                        }
101. }

复制代码

lrh2006

谢谢各位，谢谢论坛，信息量有点大，容我再琢磨琢磨
什么时候我问个问题，kk回答不了，论坛也回答不了，哈哈哈，让我再想一想

isee

回复 16# lrh2006

这种题太多了，其实

Czhang271828

100 的 $n$-互异拆分, 参考回答 kuing.cjhb.site/forum.php?mod=viewthread&tid=6607&extra=page%3D1