昨晚请美国同学帮忙，与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
2562
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
