免費論壇 繁體 | 簡體
Sclub交友聊天~加入聊天室當版主
分享
返回列表 发帖
本帖最后由 业余的业余 于 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

本帖最后由 业余的业余 于 2019-2-9 23:39 编辑

回复 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



回复 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

本帖最后由 业余的业余 于 2019-2-11 00:47 编辑

回复 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

返回列表 回复 发帖