悠闲数学娱乐论坛(第2版)'s Archiver

lrh2006 发表于 2019-2-9 21:07

取法有几种(整数拆分)

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

业余的业余 发表于 2019-2-9 21:41

[i=s] 本帖最后由 业余的业余 于 2019-2-10 01:56 编辑 [/i]

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29779&ptid=5883]1#[/url] [i]lrh2006[/i] [/b]

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

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

关键词: 动态规划算法。

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

我们把 正整数 $ 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 发表于 2019-2-9 22:31

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29780&ptid=5883]2#[/url] [i]业余的业余[/i] [/b]

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

业余的业余 发表于 2019-2-9 23:31

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29783&ptid=5883]3#[/url] [i]realnumber[/i] [/b]

:)

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

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

unsigned P(int n, int k)
{
      unsigned v=0;

      if(n<=k)
            return 1;
      if(map[n][k]!=0)
            return map[n][k];
      
       for(int j=n-1;j>0; j-=k)
            v+=P(j,k+1);
      
        return map[n][k]=v;
}

int main()
{
        std::cout<<" 总取法为 "<<P(1000,1)<<std::endl;
}
[/code]没有测试,可能会有小bug.

业余的业余 发表于 2019-2-9 23:37

好多年不玩程序了,现在改玩数学{:biggrin:}

kuing 发表于 2019-2-9 23:38

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

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

用 Mathematica 来试一下:[code]
f[sum_,max_] := Coefficient[Product[Sum[x^(m k), {k, 0, Floor[sum/m]}], {m, 1, max}], x^sum]
[/code]那么楼主要求的就是 f[1000,999],这个太大了,来算小一点的:
和为 10,无最大数限制的(即 10=10 也算一种),那就是 f[10, 10],结果是 42,如果 10=10 不算那就是 41。
我手工拆了一次来验证此结果:[code]10 = 10
10 = 9+1
10 = 8+2
10 = 8+1+1
10 = 7+3
10 = 7+2+1
10 = 7+1+1+1
10 = 6+4
10 = 6+3+1
10 = 6+2+2
10 = 6+2+1+1
10 = 6+1+1+1+1
10 = 5+5
10 = 5+4+1
10 = 5+3+2
10 = 5+3+1+1
10 = 5+2+2+1
10 = 5+2+1+1+1
10 = 5+1+1+1+1+1
10 = 4+4+2
10 = 4+4+1+1
10 = 4+3+3
10 = 4+3+2+1
10 = 4+3+1+1+1
10 = 4+2+2+2
10 = 4+2+2+1+1
10 = 4+2+1+1+1+1
10 = 4+1+1+1+1+1+1
10 = 3+3+3+1
10 = 3+3+2+2
10 = 3+3+2+1+1
10 = 3+3+1+1+1+1
10 = 3+2+2+2+1
10 = 3+2+2+1+1+1
10 = 3+2+1+1+1+1+1
10 = 3+1+1+1+1+1+1+1
10 = 2+2+2+2+2
10 = 2+2+2+2+1+1
10 = 2+2+2+1+1+1+1
10 = 2+2+1+1+1+1+1+1
10 = 2+1+1+1+1+1+1+1+1
10 = 1+1+1+1+1+1+1+1+1+1 [/code]和为 50 的,f[50, 50],运行一分多钟得出 204226。
试了下 100,等好久都不出来,看来复杂度挺大,算了……

业余的业余 发表于 2019-2-9 23:41

{:victory:}

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29787&ptid=5883]6#[/url] [i]kuing[/i] [/b]

业余的业余 发表于 2019-2-9 23:42

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29788&ptid=5883]7#[/url] [i]业余的业余[/i] [/b]

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

业余的业余 发表于 2019-2-10 00:22

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29789&ptid=5883]8#[/url] [i]业余的业余[/i] [/b]

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

程序(c++)[code]
#include <iostream>


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

unsigned P(int n, int k)
{
      unsigned v=0;

          if(n==0 || n==k)
              return 1;
      else if(n<k)
            return 0;
      if(map[n][k]!=0)
            return map[n][k];
      
       for(int j=n;j>=0; j-=k)
            v+=P(j,k+1);
      
       return map[n][k]=v;
}

int main()
{
        std::cout<<" 总取法有 "<<P(N,1)<<"种。"<<std::endl;
       
            // 调试用途
                for(int i=0; i<N+1; ++i)
                        for(int j=0; j<N-1; ++j)
                        {
                                std::cout<<map[i][j]<<(j==N-2?'\n':' ');
                        }
}[/code]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 类型。

业余的业余 发表于 2019-2-10 00:30

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29790&ptid=5883]9#[/url] [i]业余的业余[/i] [/b]

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

kuing 发表于 2019-2-10 01:05

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

PS、修改了下标题。

kuing 发表于 2019-2-10 01:11

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

realnumber 发表于 2019-2-10 09:23

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29792&ptid=5883]11#[/url] [i]kuing[/i] [/b]
我是说刚理解,对自己能熟练运用表示怀疑.

kuing 发表于 2019-2-10 14:36

唉!我早就应该想到,这么经典的问题,Mathematica 应该已经自带相应命令,一查果然如此!!![img]https://qzonestyle.gtimg.cn/qzone/em/e400824.gif[/img][img]https://qzonestyle.gtimg.cn/qzone/em/e400824.gif[/img][img]https://qzonestyle.gtimg.cn/qzone/em/e400824.gif[/img]
它们是 PartitionsP 和 IntegerPartitions,前者就是拆分方法数(同样包含 10=10 这种),后者则是列出具体的拆分。
参考:
[url]https://reference.wolfram.com/language/ref/PartitionsP.html[/url]
[url]https://reference.wolfram.com/language/ref/IntegerPartitions.html[/url]
比如,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# 的不一样了。

业余的业余 发表于 2019-2-10 22:54

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=29795&ptid=5883]14#[/url] [i]kuing[/i] [/b]

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

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

答案一致。[code]
#include <iostream>
#include <iomanip>

// 违章建筑无符号大整数型,仅供解决 N=1000 的情形. 支持36位十进制数
class luint
{
public:
       
        luint(unsigned long a=0){ data[0]=a; data[1]=data[2]=data[3]=0;}
        luint(const luint& rhs){
                data[0]=rhs.data[0];
                data[1]=rhs.data[1];
                data[2]=rhs.data[2];
                data[3]=rhs.data[3];
        }
       
        luint& operator +=(const luint& rhs)
        {
                unsigned long carry=0;
                for(int i=0; i<width; ++i)
                {
                        carry+=data[i]+rhs.data[i];
                        data[i]=carry%1000000000;
                        carry/=1000000000;
                }
                if(carry)
                        throw "overflow!";
                return *this;
        }
       
        luint& operator =(const luint& rhs)
        {
                return *(luint*)new(this)luint(rhs);
        }
       
       
        bool operator == (const luint& rhs)const
        {
                return data[0]==rhs.data[0] &&
                        data[1]==rhs.data[1] &&
                        data[2]==rhs.data[2] &&
                        data[3]==rhs.data[3];
        }
friend std::ostream& operator <<(std::ostream& os, const luint& lu)
{
        int firstNZ; // first non-zero ulong. from 3 to 0
        for(firstNZ=3; firstNZ>=0 && lu.data[firstNZ]==0; --firstNZ)
                ;
        if(firstNZ<=0)
                return os<<lu.data[0];
       
        os<<lu.data[firstNZ];
        while(--firstNZ>=0)
                os<<std::setw(9)<<std::setfill('0')<<lu.data[firstNZ];
        return os;
}
       
private:
        static const int width=4;
        unsigned long data[width];
};



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

luint one(1),zero;

luint& P(int n, int k)
{
      luint v=0;

          if(n==0 || n==k)
              return one;
      else if(n<k)
            return zero;
      if(!(map[n][k]==zero))
            return map[n][k];
      
       for(int j=n;j>=0; j-=k)
            v+=P(j,k+1);
      
       return map[n][k]=v;
}

int main()
{


        std::cout<<" 总取法有 "<<P(N,1)<<"种。"<<std::endl;
       
            // 调试用途
//                for(int i=0; i<N+1; ++i)
//                        for(int j=0; j<N-1; ++j)
//                        {
//                                std::cout<<map[i][j]<<(j==N-2?'\n':' ');
//                        }
}
[/code]

lrh2006 发表于 2019-3-9 21:23

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

isee 发表于 2019-3-9 23:19

[b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=30226&ptid=5883]16#[/url] [i]lrh2006[/i] [/b]

这种题太多了,其实

Czhang271828 发表于 2022-2-12 14:49

100 的 $n$-互异拆分, 参考回答 [url]http://kuing.orzweb.net/viewthread.php?tid=6607&extra=page%3D1[/url]

页: [1]

Powered by Discuz! Archiver 7.2  © 2001-2009 Comsenz Inc.