免費論壇 繁體 | 簡體
Sclub交友聊天~加入聊天室當版主
分享
返回列表 发帖

[组合] 取法有几种(整数拆分)

从1--999这999个自然数中,任意取出n个相加(可以重复取),使其和为1000,有几种取法?
这个问题有解吗?或者从编程的角度讲,算法应该怎么进行?
题目是自己瞎想的,请大家指教,谢谢
分享到: QQ空间QQ空间 腾讯微博腾讯微博 腾讯朋友腾讯朋友

本帖最后由 业余的业余 于 2019-2-10 01:56 编辑

回复 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 的情形,上面漏掉了。参后面的程序。

TOP

回复 2# 业余的业余

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

TOP

回复 3# realnumber

:)

用动态规划的变体备忘录法,效率稍差,堆栈要求高些,不过对现代的电脑应该不是问题。
  1. #include <iostream>

  2. // 用了全局变量,满去吧 :)
  3. unsigned map[10001][1000];

  4. unsigned P(int n, int k)
  5. {
  6.       unsigned v=0;

  7.       if(n<=k)
  8.             return 1;
  9.       if(map[n][k]!=0)
  10.             return map[n][k];
  11.       
  12.        for(int j=n-1;j>0; j-=k)
  13.             v+=P(j,k+1);
  14.       
  15.         return map[n][k]=v;
  16. }

  17. int main()
  18. {
  19.         std::cout<<" 总取法为 "<<P(1000,1)<<std::endl;
  20. }
复制代码
没有测试,可能会有小bug.
1

评分人数

TOP

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

TOP

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

其原理倒是非常好理解,设
\[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` 的个数。
如果你理解了前两天这帖 http://kuing.orzweb.net/redirect ... =5872&pid=29759 (11#)的解析,这个你应该也能明白。

用 Mathematica 来试一下:
  1. 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,等好久都不出来,看来复杂度挺大,算了……
1

评分人数

    • realnumber: 两个帖子都明白了,下次碰到能不能用起来? ...威望 + 1
$\href{https://kuingggg.github.io/}{\text{About Me}}$

TOP



回复 6# kuing

TOP

回复 7# 业余的业余

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

TOP

回复 8# 业余的业余

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

程序(c++)
  1. #include <iostream>


  2. const int N=5;
  3. // 用了全局变量,满去吧 :)
  4. unsigned map[N+1][N];

  5. unsigned P(int n, int k)
  6. {
  7.       unsigned v=0;

  8.           if(n==0 || n==k)
  9.               return 1;
  10.       else if(n<k)
  11.             return 0;
  12.       if(map[n][k]!=0)
  13.             return map[n][k];
  14.       
  15.        for(int j=n;j>=0; j-=k)
  16.             v+=P(j,k+1);
  17.       
  18.        return map[n][k]=v;
  19. }

  20. int main()
  21. {
  22.         std::cout<<" 总取法有 "<<P(N,1)<<"种。"<<std::endl;
  23.        
  24.             // 调试用途
  25.                 for(int i=0; i<N+1; ++i)
  26.                         for(int j=0; j<N-1; ++j)
  27.                         {
  28.                                 std::cout<<map[i][j]<<(j==N-2?'\n':' ');
  29.                         }
  30. }
复制代码
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 类型。

TOP

回复 9# 业余的业余

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

TOP

两个帖子都明白了,下次碰到能不能用起来?
@realnumber 为啥不能用?这些又不是新东西……

PS、修改了下标题。

TOP

再丢个链接:http://oeis.org/A000041

TOP

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

TOP

唉!我早就应该想到,这么经典的问题,Mathematica 应该已经自带相应命令,一查果然如此!!!
它们是 PartitionsP 和 IntegerPartitions,前者就是拆分方法数(同样包含 10=10 这种),后者则是列出具体的拆分。
参考:
https://reference.wolfram.com/language/ref/PartitionsP.html
https://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# 的不一样了。
$\href{https://kuingggg.github.io/}{\text{About Me}}$

TOP

回复 14# kuing

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

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

答案一致。
  1. #include <iostream>
  2. #include <iomanip>

  3. // 违章建筑无符号大整数型,仅供解决 N=1000 的情形. 支持36位十进制数
  4. class luint
  5. {
  6. public:
  7.        
  8.         luint(unsigned long a=0){ data[0]=a; data[1]=data[2]=data[3]=0;}
  9.         luint(const luint& rhs){
  10.                 data[0]=rhs.data[0];
  11.                 data[1]=rhs.data[1];
  12.                 data[2]=rhs.data[2];
  13.                 data[3]=rhs.data[3];
  14.         }
  15.        
  16.         luint& operator +=(const luint& rhs)
  17.         {
  18.                 unsigned long carry=0;
  19.                 for(int i=0; i<width; ++i)
  20.                 {
  21.                         carry+=data[i]+rhs.data[i];
  22.                         data[i]=carry%1000000000;
  23.                         carry/=1000000000;
  24.                 }
  25.                 if(carry)
  26.                         throw "overflow!";
  27.                 return *this;
  28.         }
  29.        
  30.         luint& operator =(const luint& rhs)
  31.         {
  32.                 return *(luint*)new(this)luint(rhs);
  33.         }
  34.        
  35.        
  36.         bool operator == (const luint& rhs)const
  37.         {
  38.                 return data[0]==rhs.data[0] &&
  39.                         data[1]==rhs.data[1] &&
  40.                         data[2]==rhs.data[2] &&
  41.                         data[3]==rhs.data[3];
  42.         }
  43. friend std::ostream& operator <<(std::ostream& os, const luint& lu)
  44. {
  45.         int firstNZ; // first non-zero ulong. from 3 to 0
  46.         for(firstNZ=3; firstNZ>=0 && lu.data[firstNZ]==0; --firstNZ)
  47.                 ;
  48.         if(firstNZ<=0)
  49.                 return os<<lu.data[0];
  50.        
  51.         os<<lu.data[firstNZ];
  52.         while(--firstNZ>=0)
  53.                 os<<std::setw(9)<<std::setfill('0')<<lu.data[firstNZ];
  54.         return os;
  55. }
  56.        
  57. private:
  58.         static const int width=4;
  59.         unsigned long data[width];
  60. };



  61. const int N=1000;
  62. // 用了全局变量,满去吧 :)
  63. luint map[N+1][N];

  64. luint one(1),zero;

  65. luint& P(int n, int k)
  66. {
  67.       luint v=0;

  68.           if(n==0 || n==k)
  69.               return one;
  70.       else if(n<k)
  71.             return zero;
  72.       if(!(map[n][k]==zero))
  73.             return map[n][k];
  74.       
  75.        for(int j=n;j>=0; j-=k)
  76.             v+=P(j,k+1);
  77.       
  78.        return map[n][k]=v;
  79. }

  80. int main()
  81. {


  82.         std::cout<<" 总取法有 "<<P(N,1)<<"种。"<<std::endl;
  83.        
  84.             // 调试用途
  85. //                for(int i=0; i<N+1; ++i)
  86. //                        for(int j=0; j<N-1; ++j)
  87. //                        {
  88. //                                std::cout<<map[i][j]<<(j==N-2?'\n':' ');
  89. //                        }
  90. }
复制代码

TOP

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

TOP

回复 16# lrh2006

这种题太多了,其实

TOP

100 的 $n$-互异拆分, 参考回答 http://kuing.orzweb.net/viewthre ... &extra=page%3D1
无钱佮歹看、无样佮歹生、无汉草佮无文采、无学历佮无能力、无高度无速度无力度共闲无代记。(闽南话)
口号:珍爱生命,远离内卷。

TOP

返回列表 回复 发帖