数字重复不超过k次的n位正整数有多少个？

原创已于 2024-11-13 17:34:21 修改 · 8.6k 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #NP问题 #Fortran #组合数学 #优化算法 #ChatGPT

于 2023-02-04 18:23:36 首次发布

文章讲述了通过与ChatGPT交流解决一个关于30位正整数排列组合的问题，指出ChatGPT的初始错误答案以及NP完全问题的复杂性。作者使用Fortran编写了200行代码，采用枚举和排列组合算法，成功计算出正确结果，强调在某些复杂问题上人类思维仍具优势。

昨晚请美国同学帮忙，与ChatGPT聊了一个数学问题：
30位正整数，其中0到9每一种数字最多可以重复17次，一共有多少个这样的正整数？

第一次回答： 有18³⁰个。
这个结果超过了30位整数的总量，显然是错误的。

第二次回答： 没有固定的办法来求解这个问题，因为它是一个NP完全问题，你可以用枚举排列组合的方法来获得近似值。
问题性质判断是对的，没有直接的数学计算方法，只能枚举。问题是，枚举对于30位整数的规模是不可能完成的，因为每秒枚举1亿亿个方案，也需要300万年。

第三次回答： 可以用母函数法；同学提示多位整数首位不能为0，它贴心反馈，母函数运用要有所变化，应扣除0在首位的情形。
方法已经很对路子了，但是母函数法仍然有超大计算量，计算机也无法快速实现。

用古老的fortran写了200行代码。
算法要点：
先枚举数字的组合形式，再计算每一组合对应的排列数，并结合首位预选非0数字，将枚举量从10³⁰量级降到10⁴量级。

测试了一下，枚举组合2823个，耗时一秒钟。这样的30位正整数一共有：899999999763722375523242495085个。

示例是一个规模更小的同类问题，可以直观描述出算法和过程：
12位正整数，每种数字最多可重复7次。
运用枚举组合+计算排列的算法，枚举组合48个，耗时不到十分之一秒。
计算结果，有899969278230个这样的12位正整数。
你甚至可以手工完成计算。

在NP问题上，人脑可能还是略胜一筹

附：n = 12，k = 7 的运行示例

12位正整数，数字最多重复7次，计算有多少个。

在首位不选0的前提下，按照组合分类计算后续位数的排列方案数：
00000001,xabcdefghijj：C(10,1)*C(9,9)*A(11,11)/A(2,2) = 199584000
00000002,xabcdefghiii：C(10,1)*C(9,8)*A(11,11)/A(3,3) = 598752000
00000003,xabcdefghhii：C(10,2)*C(8,7)*A(11,11)/A(2,2)^2 = 3592512000
00000004,xabcdefghhhh：C(10,1)*C(9,7)*A(11,11)/A(4,4) = 598752000
00000005,xabcdefgghhh：C(10,1)*C(9,1)*C(8,6)*A(11,11)/A(3,3)/A(2,2) = 8382528000
00000006,xabcdefggggg：C(10,1)*C(9,6)*A(11,11)/A(5,5) = 279417600
00000007,xabcdeffgghh：C(10,3)*C(7,5)*A(11,11)/A(2,2)^3 = 12573792000
00000008,xabcdeffgggg：C(10,1)*C(9,1)*C(8,5)*A(11,11)/A(4,4)/A(2,2) = 4191264000
00000009,xabcdefffggg：C(10,2)*C(8,5)*A(11,11)/A(3,3)^2 = 2794176000
00000010,xabcdeffffff：C(10,1)*C(9,5)*A(11,11)/A(6,6) = 69854400
00000011,xabcdeeffggg：C(10,1)*C(9,2)*C(7,4)*A(11,11)/A(3,3)/A(2,2)^2 = 20956320000
00000012,xabcdeefffff：C(10,1)*C(9,1)*C(8,4)*A(11,11)/A(5,5)/A(2,2) = 1047816000
00000013,xabcdeeeffff：C(10,1)*C(9,1)*C(8,4)*A(11,11)/A(4,4)/A(3,3) = 1746360000
00000014,xabcdeeeeeee：C(9,1)*C(9,4)*A(11,11)/A(7,7) = 8981280
00000015,xabcddeeffgg：C(10,4)*C(6,3)*A(11,11)/A(2,2)^4 = 10478160000
00000016,xabcddeeffff：C(10,1)*C(9,2)*C(7,3)*A(11,11)/A(4,4)/A(2,2)^2 = 5239080000
00000017,xabcddeeefff：C(10,2)*C(8,1)*C(7,3)*A(11,11)/A(3,3)^2/A(2,2) = 6985440000
00000018,xabcddeeeeee：C(10,1)*C(9,1)*C(8,3)*A(11,11)/A(6,6)/A(2,2) = 139708800
00000019,xabcdddeeeee：C(10,1)*C(9,1)*C(8,3)*A(11,11)/A(5,5)/A(3,3) = 279417600
00000020,xabcddddeeee：C(10,2)*C(8,3)*A(11,11)/A(4,4)^2 = 174636000
00000021,xabccddeefff：C(10,1)*C(9,3)*C(6,2)*A(11,11)/A(3,3)/A(2,2)^3 = 10478160000
00000022,xabccddeeeee：C(10,1)*C(9,2)*C(7,2)*A(11,11)/A(5,5)/A(2,2)^2 = 628689600
00000023,xabccdddeeee：C(10,1)*C(9,1)*C(8,1)*C(7,2)*A(11,11)/A(4,4)/A(3,3)/A(2,2) = 2095632000
00000024,xabccddddddd：C(9,1)*C(9,1)*C(8,2)*A(11,11)/A(7,7)/A(2,2) = 8981280
00000025,xabcccdddeee：C(10,3)*C(7,2)*A(11,11)/A(3,3)^3 = 465696000
00000026,xabcccdddddd：C(10,1)*C(9,1)*C(8,2)*A(11,11)/A(6,6)/A(3,3) = 23284800
00000027,xabccccddddd：C(10,1)*C(9,1)*C(8,2)*A(11,11)/A(5,5)/A(4,4) = 34927200
00000028,xabbccddeeff：C(10,5)*C(5,1)*A(11,11)/A(2,2)^5 = 1571724000
00000029,xabbccddeeee：C(10,1)*C(9,3)*C(6,1)*A(11,11)/A(4,4)/A(2,2)^3 = 1047816000
00000030,xabbccdddeee：C(10,2)*C(8,2)*C(6,1)*A(11,11)/A(3,3)^2/A(2,2)^2 = 2095632000
00000031,xabbccdddddd：C(10,1)*C(9,2)*C(7,1)*A(11,11)/A(6,6)/A(2,2)^2 = 34927200
00000032,xabbcccddddd：C(10,1)*C(9,1)*C(8,1)*C(7,1)*A(11,11)/A(5,5)/A(3,3)/A(2,2) = 139708800
00000033,xabbccccdddd：C(10,2)*C(8,1)*C(7,1)*A(11,11)/A(4,4)^2/A(2,2) = 87318000
00000034,xabbbcccdddd：C(10,1)*C(9,2)*C(7,1)*A(11,11)/A(4,4)/A(3,3)^2 = 116424000
00000035,xabbbccccccc：C(9,1)*C(9,1)*C(8,1)*A(11,11)/A(7,7)/A(3,3) = 855360
00000036,xabbbbcccccc：C(10,1)*C(9,1)*C(8,1)*A(11,11)/A(6,6)/A(4,4) = 1663200
00000037,xabbbbbccccc：C(10,2)*C(8,1)*A(11,11)/A(5,5)^2 = 997920
00000038,xaabbccddeee：C(10,1)*C(9,4)*A(11,11)/A(3,3)/A(2,2)^4 = 523908000
00000039,xaabbccddddd：C(10,1)*C(9,3)*A(11,11)/A(5,5)/A(2,2)^3 = 34927200
00000040,xaabbcccdddd：C(10,1)*C(9,1)*C(8,2)*A(11,11)/A(4,4)/A(3,3)/A(2,2)^2 = 174636000
00000041,xaabbccccccc：C(9,1)*C(9,2)*A(11,11)/A(7,7)/A(2,2)^2 = 641520
00000042,xaabbbcccddd：C(10,3)*C(7,1)*A(11,11)/A(3,3)^3/A(2,2) = 77616000
00000043,xaabbbcccccc：C(10,1)*C(9,1)*C(8,1)*A(11,11)/A(6,6)/A(3,3)/A(2,2) = 3326400
00000044,xaabbbbccccc：C(10,1)*C(9,1)*C(8,1)*A(11,11)/A(5,5)/A(4,4)/A(2,2) = 4989600
00000045,xaaabbbccccc：C(10,1)*C(9,2)*A(11,11)/A(5,5)/A(3,3)^2 = 3326400
00000046,xaaabbbbcccc：C(10,2)*C(8,1)*A(11,11)/A(4,4)^2/A(3,3) = 4158000
00000047,xaaaabbbbbbb：C(9,1)*C(9,1)*A(11,11)/A(7,7)/A(4,4) = 26730
00000048,xaaaaabbbbbb：C(10,1)*C(9,1)*A(11,11)/A(6,6)/A(5,5) = 41580

上述合计 = 99996586470

枚举组合数 = 48

符合条件的方案数 = 9×99996586470 = 899969278230

附：fortran代码

$freeform

! nn位正整数，各数字最多重复ni次，一共有多少种？
! 2023-02-04
! 作者：szw_sh@163.com

! 算法要点：

! x为首位数字，不能为0，9选1；

! 计算其后(nn-1)位所有组合及其排列方案总数；
! 对于重复次数为ni的数字，只能从(10-1=9)个数字中选取n个；
! 重复数字少于ni的其它数字，按照(9-n)数字中选取m，依次进行；
! 叠加(nn-1)的全排列，再除以重复数字排列方案数；

! (nn-1)位的方案数，乘以首位的9种选择，得到总方案数；

! 为使得输出计算式逻辑更加清晰，公式顺序与程序计算顺序不同；
! 排列和除法项记录于ha，组合记录于hc，最后格式化输出；

! 删去计算排列组合子程序，用数组aa和cc加速计算；
! m大于100时，选择输出，减少时耗。


module abc                                                 ! 公共变量
   use big_integers_module
   parameter (nn=30,ni=17,nf=max(nn,10),na=min(10,nn))     ! nn位数，ni数字可重复次数，nf阶乘，na递归枚举数组
   character h(10)                                         ! 代指10个数字的符号，用于表达组合方案
   data h/'a','b','c','d','e','f','g','h','i','j'/
   character*150 ha,hc                                     ! 保存计算式的变量
   type(big_integer) f(0:nf)                               ! f用于加快排列组合计算的阶乘数组
   type(big_integer) mm                                    ! mm累加方案数(总方案数)
   type(big_integer) aa(ni,10)                             ! aa保存A(i,i)**j，减少重复计算
   integer m                                               ! m组合计数
   integer cc(10,10)                                       ! cc保存C(i,j)，减少重复计算
end module



use abc                                                    ! 主程序
integer a(na)

call init                                                  ! 初始化f，aa，cc

write(*,'(/i<(nn+40)/50+1>,a,i<(ni+40)/50+1>,a)') nn,'位正整数，数字最多重复',ni,'次，计算有多少个。'

write(*,'(/a)') '在首位不选0的前提下，按照组合分类计算后续位数的排列方案数：'

a=0
call pp(a,1,0)                                             ! 以首位a的不同个数调用pp，枚举计算(nn-1)位组合的排列方案数

write(*,'(/2a)')  '上述合计 = ',trim(char(mm))
write(*,'(/2a)')  '枚举组合数 = ',trim(char(m))
write(*,'(/2a\)') '符合条件的方案数 = 9×',trim(char(mm))
mm=mm*9
write(*,'(2a)')   ' = ',trim(char(mm))

end



recursive subroutine pp(a,n,b)                             ! 递归子程序，枚举组合，计算排列

use abc                                                    ! b为sum(a(1:n-1))
type(big_integer) t
integer a(na),b,s(ni)
logical key

np=1
if(n.gt.1) np=a(n-1)

do i=np,ni

   a(n)=i
   k=b+i

   if(k.gt.nn-1) then                                      ! 组合位数超出nn，退出递归
     exit

   else if(k+(na-n)*ni.lt.nn-1) then                       ! 后续组合位数达不到nn，避免f下标负数，跳过
     cycle

   else if(k.eq.nn-1) then                                 ! 符合nn位的组合，计数，并计算排列数

      m=m+1
      key=mod(m,321).eq.1
      if(m.le.100) key=.true.

      if(key) write(*,'(i8.8,a,\)') m,','

      s=0                                                  ! 统计各重复次数对应的数字有几个
      if(key) write(*,'(a\)') 'x'

      do j=1,n                                             ! 输出枚举到的组合样式
         s(a(j))=s(a(j))+1
         do j1=1,a(j)
            if(key) write(*,'(a\)') h(j)
         end do
      end do

      ha=' '                                               ! ha,hc 置初值，写入计算式
      ka=0
      hc=' '
      kc=0

      ii=10                                                ! 计算排列数量，并记字符串
      t=f(nn-1)                                            ! 此处不能使用aa，aa最大只能到ni的阶乘
      if(key) write(ha(ka+1:ka+15),'(a,i2,a,i2,a)') '*A(',nn-1,',',nn-1,')'
      ka=ka+15

      do j=ni,1,-1                                         ! 由重复数最多的ni开始，逐个选取并计算C和A的值

         if(s(j).eq.0) cycle

         if(j.eq.ni) then                                  ! 重复数为ni的数字的选择和计算

            t=t*cc(9,s(j))/aa(j,s(j))
            if(key) write(hc(kc+1:kc+10),'(a,i2,a)') '*C(9,',s(j),')'
            kc=kc+15

            if(j.gt.1.and.key) then
               write(ha(ka+1:ka+15),'(a,i2,a,i2,a)') '/A(',j,',',j,')'
               if(s(j).gt.1) write(ha(ka+11:ka+15),'(a,i2)') '^',s(j)
               ka=ka+15
            end if

         else                                              ! 其它情形的计算

            t=t*cc(ii,s(j))/aa(j,s(j))
            if(key) write(hc(kc+1:kc+10),'(a,i2,a,i2,a)') '*C(',ii,',',s(j),')'
            kc=kc+15

            if(j.gt.1.and.key) then
               write(ha(ka+1:ka+15),'(a,i2,a,i2,a)') '/A(',j,',',j,')'
               if(s(j).gt.1) write(ha(ka+11:ka+15),'(a,i2)') '^',s(j)
               ka=ka+15
            end if

         end if

         ii=ii-s(j)                                        ! 调整可选数量ii

      end do

      hc(1:1)=' '                                          ! hc头部的*是为统一写入格式造成，需要清除
      if(key) call fmt_ha_hc                               ! 格式化 ha hc，输出计算式和结果
      if(key) write(*,'(5a)') '：',trim(hc),trim(ha),' = ',trim(char(t))
      mm=mm+t                                              ! 累加方案数量

   else
      call pp(a,n+1,k)                                     ! 未达成组合的，继续递归枚举

   end if

end do

end



subroutine fmt_ha_hc                                       ! 子程序，格式化输出公式字符串

use abc

ii=0
do j=1,len_trim(hc)
   if(hc(j:j).eq.' ') cycle
   ii=ii+1
   hc(ii:ii)=hc(j:j)
end do
hc(ii+1:)=' '

ii=0
do j=1,len_trim(ha)
   if(ha(j:j).eq.' ') cycle
   ii=ii+1
   ha(ii:ii)=ha(j:j)
end do
ha(ii+1:)=' '

end



subroutine init                                            ! 子程序，初始化f aa cc

use abc

f(0)=1
do i=1,nf
   f(i)=f(i-1)*i
end do

do i=1,10
   do j=1,i
      cc(i,j)=f(i)/f(j)/f(i-j)
   end do
end do

do i=1,ni
   do j=1,10
      aa(i,j)=f(i)**j
   end do
end do

end