这是另一个问题,最好在深入研究和编写代码之前对其进行检查。
我要说你应该做的第一件事是查看数字网格,不要将它们表示为小数,而是用分数表示。
首先应该清楚的是,您拥有的 的总数只是与原点 的距离的度量。
如果你这样看一个网格,你可以得到所有的分母:
请注意,第一行和第一列并不都是1s - 它们被选择遵循模式,以及适用于所有其他方格的通用公式。
分子有点棘手,但仍然可行。与大多数此类问题一样,答案与组合、阶乘以及一些更复杂的事情有关。这里的典型条目包括Catalan numbers、Stirling's numbers、Pascal's triangle,您几乎总是会看到使用了Hypergeometric functions。
除非您大量数学,否则您不太可能熟悉所有这些,并且有大量的文学作品。所以我有一个更简单的方法来找出你需要的关系,这几乎总是有效的。它是这样的:
- 编写一个简单、低效的算法来获得您想要的序列。
- 将相当多的数字复制到 google 中。
-
希望Online Encyclopedia of Integer Sequences 的结果弹出。
3.b。如果没有,请查看您的序列中的一些差异,或与您的数据相关的其他一些序列。
使用您找到的信息来实现所述序列。
所以,按照这个逻辑,这里是分子:
现在,不幸的是,谷歌搜索没有任何结果。但是,您可以注意到一些关于它们的事情,主要是第一行/列只是 3 的幂,而第二行/列是小于 3 的幂。这种边界与帕斯卡三角形一模一样,还有很多相关的序列。
这是分子和分母之间的差异矩阵:
我们决定 f(0,0) 元素应该遵循相同的模式。这些数字看起来已经简单多了。还要注意——相当有趣的是,这些数字遵循与初始数字相同的规则——除了第一个数字是一个(并且它们被一列和一行偏移)。 T(i,j) = T(i-1,j) + T(i,j-1) + 3*T(i-1,j-1):
1
1 1
1 5 1
1 9 9 1
1 13 33 13 1
1 17 73 73 17 1
1 21 129 245 192 21 1
1 25 201 593 593 201 25 1
这看起来更像你在组合数学中经常看到的序列。
If you google numbers from this matrix, you do get a hit.
然后如果你切断原始数据的链接,你会得到序列A081578,它被描述为“Pascal-(1,3,1) 数组”,这完全有道理 - 如果你旋转矩阵,使0,0元素在顶部,元素组成一个三角形,那么你取左边元素1*,上面元素3*,右边元素1*。
现在的问题是实现用于生成数字的公式。
不幸的是,这往往说起来容易做起来难。比如页面上给出的公式:
T(n,k)=sum{j=0..n, C(k,j-k)*C(n+k-j,k)*3^(j-k)}
是错误的,需要大量阅读 the paper(链接在页面上)才能得出正确的公式。你想要的部分是命题 26,推论 28。在命题 13 之后的表 2 中提到了序列。注意r=4
命题 26 给出了正确的公式,但那里也有一个错字:/。总和中的k=0 应该是j=0:
其中T 是包含系数的三角矩阵。
OEIS 页面确实提供了几个计算数字的实现,但是它们都不是 java 的,而且它们都不能轻易地转录为 java:
有一个数学例子:
Table[ Hypergeometric2F1[-k, k-n, 1, 4], {n, 0, 10}, {k, 0, n}] // Flatten
一如既往地简洁得可笑。还有一个 Haskell 版本,同样简洁:
a081578 n k = a081578_tabl !! n !! k
a081578_row n = a081578_tabl !! n
a081578_tabl = map fst $ iterate
(\(us, vs) -> (vs, zipWith (+) (map (* 3) ([0] ++ us ++ [0])) $
zipWith (+) ([0] ++ vs) (vs ++ [0]))) ([1], [1, 1])
我知道您在 java 中这样做,但我懒得抄写我对 java 的回答(抱歉)。这是一个python实现:
from __future__ import division
import math
#
# Helper functions
#
def cache(function):
cachedResults = {}
def wrapper(*args):
if args in cachedResults:
return cachedResults[args]
else:
result = function(*args)
cachedResults[args] = result
return result
return wrapper
@cache
def fact(n):
return math.factorial(n)
@cache
def binomial(n,k):
if n < k: return 0
return fact(n) / ( fact(k) * fact(n-k) )
def numerator(i,j):
"""
Naive way to calculate numerator
"""
if i == j == 0:
return 0
elif i == 0 or j == 0:
return 3**(max(i,j)-1)
else:
return numerator(i-1,j) + numerator(i,j-1) + 3*numerator(i-1,j-1)
def denominator(i,j):
return 3**(i+j-1)
def A081578(n,k):
"""
http://oeis.org/A081578
"""
total = 0
for j in range(n-k+1):
total += binomial(k, j) * binomial(n-k, j) * 4**(j)
return int(total)
def diff(i,j):
"""
Difference between the numerator, and the denominator.
Answer will then be 1-diff/denom.
"""
if i == j == 0:
return 1/3
elif i==0 or j==0:
return 0
else:
return A081578(j+i-2,i-1)
def answer(i,j):
return 1 - diff(i,j) / denominator(i,j)
# And a little bit at the end to demonstrate it works.
N, M = 10,10
for i in range(N):
row = "%10.5f"*M % tuple([numerator(i,j)/denominator(i,j) for j in range(M)])
print row
print ""
for i in range(N):
row = "%10.5f"*M % tuple([answer(i,j) for j in range(M)])
print row
所以,对于一个封闭的形式:
只是二项式系数。
结果如下:
最后一个加法,如果您希望对大数执行此操作,那么您将需要以不同的方式计算二项式系数,因为您会溢出整数。不过,您的答案是 lal 浮点数,而且由于您显然对大 f(n) = T(n,n) 感兴趣,那么我想您可以使用斯特林近似值之类的。