找到所有可能的加权和等于一定整数的正半整数
我遇到了fortran 90代码的问题。所以,我有n+1
正半定整数,k2
,k3
,...,kn
和n
。现在对于给定的n
,我需要找到,k2
,k3
,kn
的所有可能的组合,使得k1+2*k2+3*k3+...+n*kn = n
。有人可能会想到使用一个n级嵌套循环,每个ki
它从0
到n
运行,但实际上我会把这个代码放在一个子程序和n
(也就是,k
的数量)将是一个输入这个子程序。因此,如果我要使用嵌套循环,那么嵌套层次的数量应该是自动的,我认为使用Fortran很困难(如果不是不可能的话)。有没有更好的方法来解决这个问题?找到所有可能的加权和等于一定整数的正半整数
如果n
比较小(比如5),我认为编写多维环路以获得满足k1 + k2 * 2的所需k向量(k1,k2,...,kn)会更简单+ ... + kn * n = n。否则,它可能是一个使用递归的选项。示例代码可能是这样的(注意,我们需要recursive
关键字):
module recur_mod
implicit none
integer, parameter :: nvecmax = 1000
contains
subroutine recur_main(n, kveclist, nvec)
integer, intent(in) :: n
integer, allocatable, intent(out) :: kveclist(:,:) !! list of k-vectors
integer, intent(out) :: nvec !! number of k-vectors
integer kvec(n), ind
allocate(kveclist(n, nvecmax))
kveclist(:,:) = 0
nvec = 0
ind = 1
kvec(:) = 0
call recur_sub(n, ind, kvec, kveclist, nvec) !! now entering recursion...
endsubroutine
recursive subroutine recur_sub(n, ind, kvec, kveclist, nvec)
integer, intent(in) :: n, ind
integer, intent(inout) :: kvec(:), kveclist(:,:), nvec
integer k, ksum, t, ind_next
do k = 0, n
kvec(ind) = k
ksum = sum([(kvec(t) * t, t = 1, ind)]) !! k1 + k2*2 + ... + ki*i
if (ksum > n) cycle !! early rejection
if (ind < n) then
ind_next = ind + 1
call recur_sub(n, ind_next, kvec, kveclist, nvec) !! go to the next index
endif
if (ind == n .and. ksum == n) then
nvec = nvec + 1
if (nvec > nvecmax) stop "nvecmax too small"
kveclist(:, nvec) = kvec(:) !! save k-vectors
endif
enddo
endsubroutine
end
program main
use recur_mod
implicit none
integer, allocatable :: kveclist(:,:)
integer :: n, nvec, ivec
do n = 1, 4
call recur_main(n, kveclist, nvec)
print *
print *, "For n = ", n
do ivec = 1, nvec
print *, "kvec = ", kveclist(:, ivec)
enddo
enddo
end
这给(与gfortran 7.1)
For n = 1
kvec = 1
For n = 2
kvec = 0 1
kvec = 2 0
For n = 3
kvec = 0 0 1
kvec = 1 1 0
kvec = 3 0 0
For n = 4
kvec = 0 0 0 1
kvec = 0 2 0 0
kvec = 1 0 1 0
kvec = 2 1 0 0
kvec = 4 0 0 0
在这里,我们可以看到,例如,kvec = [k1,k2,k3,k4] = [2,1,0,0]
对于n = 4满足k1 + k2*2 + k3*3 + k4*4 = 2 + 1*2 + 0 + 0 = 4
。使用这些k向量,我们可以评估f(g(x))的n阶导数(如OP所述,在page之后)。
要看到递归是如何工作的,它往往是有用的插入很多print
的 语句和监视变量如何变化。例如,代码的“冗长”的版本可能看起来像这样(在这里,我已经删除了早期排斥事情简单):
recursive subroutine recur_sub(n, ind, kvec, kveclist, nvec)
integer, intent(in) :: n, ind
integer, intent(inout) :: kvec(:), kveclist(:,:), nvec
integer k, ksum, t, ind_next
print *, "Top of recur_sub: ind = ", ind
do k = 0, n
kvec(ind) = k
print *, "ind = ", ind, " k = ", k, "kvec = ", kvec
if (ind < n) then
ind_next = ind + 1
print *, " > now going to the next level"
call recur_sub(n, ind_next, kvec, kveclist, nvec)
print *, " > returned to the current level"
endif
if (ind == n) then
ksum = sum([(kvec(t) * t, t = 1, n)])
if (ksum == n) then
nvec = nvec + 1
if (nvec > nvecmax) stop "nvecmax too small"
kveclist(:, nvec) = kvec(:)
endif
endif
enddo
print *, "Exiting recur_sub"
endsubroutine
这给了n = 2
:
Top of recur_sub: ind = 1
ind = 1 k = 0 kvec = 0 0
> now going to the next level
Top of recur_sub: ind = 2
ind = 2 k = 0 kvec = 0 0
ind = 2 k = 1 kvec = 0 1
ind = 2 k = 2 kvec = 0 2
Exiting recur_sub
> returned to the current level
ind = 1 k = 1 kvec = 1 2
> now going to the next level
Top of recur_sub: ind = 2
ind = 2 k = 0 kvec = 1 0
ind = 2 k = 1 kvec = 1 1
ind = 2 k = 2 kvec = 1 2
Exiting recur_sub
> returned to the current level
ind = 1 k = 2 kvec = 2 2
> now going to the next level
Top of recur_sub: ind = 2
ind = 2 k = 0 kvec = 2 0
ind = 2 k = 1 kvec = 2 1
ind = 2 k = 2 kvec = 2 2
Exiting recur_sub
> returned to the current level
Exiting recur_sub
For n = 2
kvec = 0 1
kvec = 2 0
请注意,在递归调用的进入和退出时,数组kvec
和kveclist
的值总是保留。特别是,我们覆盖kvec
的各个维度的元素,以获得可能模式的完整列表(基本上等同于多维循环)。我们还注意到,kvec
的值仅在最终递归级别(即,ind = n
)有意义。
另一种方法来获取工作是改写为等价的,非递归的那些特定n
递归调用。例如,n = 2
的情况可以被重写如下。这不依赖于递归,而是执行与上述代码相同的操作(对于n = 2
)。如果我们想象将recur_sub2
内嵌到recur_sub
中,则很明显这些例程相当于和k2
上的二维循环。
!! routine specific for n = 2
subroutine recur_sub(n, ind, kvec, kveclist, nvec)
integer, intent(in) :: n, ind
integer, intent(inout) :: kvec(:), kveclist(:,:), nvec
integer k1
!! now ind == 1
do k1 = 0, n
kvec(ind) = k1
call recur_sub2(n, ind + 1, kvec, kveclist, nvec) !! go to next index
enddo
endsubroutine
!! routine specific for n = 2
subroutine recur_sub2(n, ind, kvec, kveclist, nvec)
integer, intent(in) :: n, ind
integer, intent(inout) :: kvec(:), kveclist(:,:), nvec
integer k2, ksum, t
!! now ind == n == 2
do k2 = 0, n
kvec(ind) = k2
ksum = sum([(kvec(t) * t, t = 1, n)])
if (ksum == 2) then
nvec = nvec + 1
if (nvec > nvecmax) stop "nvecmax too small"
kveclist(:, nvec) = kvec(:) !! save k-vectors
endif
enddo
endsubroutine
附录
以下是对于f的(G(X))(只是为了好玩)第n衍生物的 “相当打印” 例程(在朱)。如果必要的话,请用Fortran相应的程序(但请小心大N!)
function recur_main(n)
kvec = zeros(Int, n) # [k1,k2,...,kn] (work vector)
kveclist = Vector{Int}[] # list of k-vectors (output)
recur_sub(n, 1, kvec, kveclist) # now entering recursion over {ki}...
return kveclist
end
function recur_sub(n, i, kvec, kveclist)
for ki = 0 : n
kvec[ i ] = ki
ksum = sum(kvec[ t ] * t for t = 1:i) # k1 + k2*2 + ... + ki*i
ksum > n && continue # early rejection
if i < n
recur_sub(n, i + 1, kvec, kveclist) # go to the next index
end
if i == n && ksum == n
push!(kveclist, copy(kvec)) # save k-vectors
end
end
end
function showderiv(n)
kveclist = recur_main(n)
println()
println("(f(g))_$(n) = ")
fact(k) = factorial(big(k))
for (term, kvec) in enumerate(kveclist)
fac1 = div(fact(n), prod(fact(kvec[ i ]) for i = 1:n))
fac2 = prod(fact(i)^kvec[ i ] for i = 1:n)
coeff = div(fac1, fac2)
term == 1 ? print(" ") : print(" + ")
coeff == 1 ? print(" "^15) : @printf("%15i", coeff)
print(" (f$(sum(kvec)))")
for i = 1 : length(kvec)
ki = kvec[ i ]
if ki > 0
print("(g$(i))")
ki > 1 && print("^$(ki)")
end
end
println()
end
end
#--- test ---
if false
for n = 1 : 4
kveclist = recur_main(n)
println("\nFor n = ", n)
for kvec in kveclist
println("kvec = ", kvec)
end
end
end
showderiv(1)
showderiv(2)
showderiv(3)
showderiv(4)
showderiv(5)
showderiv(8)
showderiv(10)
结果(其中fm
和gm
分别指的f
和g
,第m个衍生物):
(f(g))_1 =
(f1)(g1)
(f(g))_2 =
(f1)(g2)
+ (f2)(g1)^2
(f(g))_3 =
(f1)(g3)
+ 3 (f2)(g1)(g2)
+ (f3)(g1)^3
(f(g))_4 =
(f1)(g4)
+ 3 (f2)(g2)^2
+ 4 (f2)(g1)(g3)
+ 6 (f3)(g1)^2(g2)
+ (f4)(g1)^4
(f(g))_5 =
(f1)(g5)
+ 10 (f2)(g2)(g3)
+ 5 (f2)(g1)(g4)
+ 15 (f3)(g1)(g2)^2
+ 10 (f3)(g1)^2(g3)
+ 10 (f4)(g1)^3(g2)
+ (f5)(g1)^5
(f(g))_8 =
(f1)(g8)
+ 35 (f2)(g4)^2
+ 56 (f2)(g3)(g5)
+ 28 (f2)(g2)(g6)
+ 280 (f3)(g2)(g3)^2
+ 210 (f3)(g2)^2(g4)
+ 105 (f4)(g2)^4
+ 8 (f2)(g1)(g7)
+ 280 (f3)(g1)(g3)(g4)
+ 168 (f3)(g1)(g2)(g5)
+ 840 (f4)(g1)(g2)^2(g3)
+ 28 (f3)(g1)^2(g6)
+ 280 (f4)(g1)^2(g3)^2
+ 420 (f4)(g1)^2(g2)(g4)
+ 420 (f5)(g1)^2(g2)^3
+ 56 (f4)(g1)^3(g5)
+ 560 (f5)(g1)^3(g2)(g3)
+ 70 (f5)(g1)^4(g4)
+ 210 (f6)(g1)^4(g2)^2
+ 56 (f6)(g1)^5(g3)
+ 28 (f7)(g1)^6(g2)
+ (f8)(g1)^8
(f(g))_10 =
(f1)(g10)
+ 126 (f2)(g5)^2
+ 210 (f2)(g4)(g6)
+ 120 (f2)(g3)(g7)
+ 2100 (f3)(g3)^2(g4)
+ 45 (f2)(g2)(g8)
+ 1575 (f3)(g2)(g4)^2
+ 2520 (f3)(g2)(g3)(g5)
+ 630 (f3)(g2)^2(g6)
+ 6300 (f4)(g2)^2(g3)^2
+ 3150 (f4)(g2)^3(g4)
+ 945 (f5)(g2)^5
+ 10 (f2)(g1)(g9)
+ 1260 (f3)(g1)(g4)(g5)
+ 840 (f3)(g1)(g3)(g6)
+ 2800 (f4)(g1)(g3)^3
+ 360 (f3)(g1)(g2)(g7)
+ 12600 (f4)(g1)(g2)(g3)(g4)
+ 3780 (f4)(g1)(g2)^2(g5)
+ 12600 (f5)(g1)(g2)^3(g3)
+ 45 (f3)(g1)^2(g8)
+ 1575 (f4)(g1)^2(g4)^2
+ 2520 (f4)(g1)^2(g3)(g5)
+ 1260 (f4)(g1)^2(g2)(g6)
+ 12600 (f5)(g1)^2(g2)(g3)^2
+ 9450 (f5)(g1)^2(g2)^2(g4)
+ 4725 (f6)(g1)^2(g2)^4
+ 120 (f4)(g1)^3(g7)
+ 4200 (f5)(g1)^3(g3)(g4)
+ 2520 (f5)(g1)^3(g2)(g5)
+ 12600 (f6)(g1)^3(g2)^2(g3)
+ 210 (f5)(g1)^4(g6)
+ 2100 (f6)(g1)^4(g3)^2
+ 3150 (f6)(g1)^4(g2)(g4)
+ 3150 (f7)(g1)^4(g2)^3
+ 252 (f6)(g1)^5(g5)
+ 2520 (f7)(g1)^5(g2)(g3)
+ 210 (f7)(g1)^6(g4)
+ 630 (f8)(g1)^6(g2)^2
+ 120 (f8)(g1)^7(g3)
+ 45 (f9)(g1)^8(g2)
+ (f10)(g1)^10
它也可以用这种方式表示,但主要问题是让Fortran找到满足上述要求的'ki'的所有可能组合。 – nougako
这是拖钓或家庭作业。如果不是,你能解释你正在做什么,为什么? – starmole
这既不是驯服也不是功课。我实际上试图用所谓的Faa di Bruno公式来计算函数f(g(x))相对于x的任意阶的导数。如果你看[链接](http://mathworld.wolfram.com/FaadiBrunosFormula.html)中的等式(1),你会明白为什么我需要做我在我的文章中解释的内容。 – nougako