最近我在几年前写的一些代码中发现了这一点。它被用来通过确定合适的分母来合理化实数(在公差范围内),然后检查原始实数和有理数之间的差是否足够小。
编辑以澄清:我实际上不想转换所有的实数。例如,我可以选择14的最大分母,而等于7/15的实数将保持原样。这并不清楚,因为它是我在这里写的算法中的一个外部变量。
得到分母的算法是这样的(伪代码):
denominator(x)
frac = fractional part of x
recip = 1/frac
if (frac < tol)
return 1
else
return recip * denominator(recip)
end
end
似乎是基于连续分数的,尽管再次看它就清楚了,它是错误的。(它对我来说很有效,因为它最终会吐出无穷大,这是我在外面处理的,但它通常会非常慢。)tol的值实际上没有任何作用,除非是在终止或数字接近的情况下。我不认为这与真正理性转换的容忍度有关。
我用迭代版本取代了它,它不仅更快,而且我很确定它理论上不会失败(一开始d=1,分数部分返回正,所以recip总是>=1):
denom_iter(x d)
return d if d > maxd
frac = fractional part of x
recip = 1/frac
if (frac = 0)
return d
else
return denom_iter(recip d*recip)
end
end
我想知道的是,是否有一种方法可以选择maxd,以确保它转换给定公差下可能的所有值。我假设是1/tol,但不想错过什么。我还想知道在这种方法中是否有一种方法可以真正限制分母大小——这允许一些分母大于maxd。
这可以被认为是错误的2D最小化问题:
ArgMin ( r - q / p ), where r is real, q and p are integers
我建议使用 Gradient Descent algorithm
。该目标函数中的梯度为:
f'(q, p) = (-1/p, q/p^2)
初始猜测 r_o
可以是q是最接近r的整数,p是1。
停止条件可以是误差的阈值。
GD的伪代码可以在wiki中找到:http://en.wikipedia.org/wiki/Gradient_descent
如果最初的猜测足够接近,那么目标函数应该是凸的。
正如雅各布所建议的,通过最小化以下错误函数可以更好地解决这个问题:
ArgMin ( p * r - q ), where r is real, q and p are integers
这就是线性规划,它可以通过任何ILP(整数线性规划)求解器有效求解。 GD
处理非线性问题,但在线性问题上缺乏效率。
初始猜测和停止条件可以类似于上面所述。单独选择求解器可以获得更好的选择。
我建议你仍然应该假设凸性接近局部最小值,这可以大大降低成本。您也可以尝试Simplex方法,它在线性规划问题上很好。
我在这方面赞扬雅各布。
类似的问题在Bill Gosper的连分式算术文档的第28页开始的Approximations部分中得到了解决。(参考:postscript文件;另请参阅1984行的文本版本。)总体思路是计算低端和高端范围限制数的连续分数近似值,直到两个分数不同,然后在这两个近似值的范围内选择一个值。使用Gosper的术语,这保证给出一个最简单的分数。
下面的python代码(程序"simpleden")实现了类似的过程。(它可能没有Gosper建议的实现好,但足够好,你可以看到该方法产生了什么样的结果。)所做的工作量与欧几里得算法的工作量相似,ieO(n)用于n位的数字,因此该程序相当快。一些示例测试用例(即程序的输出)显示在代码本身之后。注意,如果vhi小于vlo,这里显示的函数simpleratio(vlo, vhi)
返回-1。
#!/usr/bin/env python
def simpleratio(vlo, vhi):
rlo, rhi, eps = vlo, vhi, 0.0000001
if vhi < vlo: return -1
num = denp = 1
nump = den = 0
while 1:
klo, khi = int(rlo), int(rhi)
if klo != khi or rlo-klo < eps or rhi-khi < eps:
tlo = denp + klo * den
thi = denp + khi * den
if tlo < thi:
return tlo + (rlo-klo > eps)*den
elif thi < tlo:
return thi + (rhi-khi > eps)*den
else:
return tlo
nump, num = num, nump + klo * num
denp, den = den, denp + klo * den
rlo, rhi = 1/(rlo-klo), 1/(rhi-khi)
def test(vlo, vhi):
den = simpleratio(vlo, vhi);
fden = float(den)
ilo, ihi = int(vlo*den), int(vhi*den)
rlo, rhi = ilo/fden, ihi/fden;
izok = 'ok' if rlo <= vlo <= rhi <= vhi else 'wrong'
print '{:4d}/{:4d} = {:0.8f} vlo:{:0.8f} {:4d}/{:4d} = {:0.8f} vhi:{:0.8f} {}'.format(ilo,den,rlo,vlo, ihi,den,rhi,vhi, izok)
test (0.685, 0.695)
test (0.685, 0.7)
test (0.685, 0.71)
test (0.685, 0.75)
test (0.685, 0.76)
test (0.75, 0.76)
test (2.173, 2.177)
test (2.373, 2.377)
test (3.484, 3.487)
test (4.0, 4.87)
test (4.0, 8.0)
test (5.5, 5.6)
test (5.5, 6.5)
test (7.5, 7.3)
test (7.5, 7.5)
test (8.534537, 8.534538)
test (9.343221, 9.343222)
程序输出:
> ./simpleden
8/ 13 = 0.61538462 vlo:0.68500000 9/ 13 = 0.69230769 vhi:0.69500000 ok
6/ 10 = 0.60000000 vlo:0.68500000 7/ 10 = 0.70000000 vhi:0.70000000 ok
6/ 10 = 0.60000000 vlo:0.68500000 7/ 10 = 0.70000000 vhi:0.71000000 ok
2/ 4 = 0.50000000 vlo:0.68500000 3/ 4 = 0.75000000 vhi:0.75000000 ok
2/ 4 = 0.50000000 vlo:0.68500000 3/ 4 = 0.75000000 vhi:0.76000000 ok
3/ 4 = 0.75000000 vlo:0.75000000 3/ 4 = 0.75000000 vhi:0.76000000 ok
36/ 17 = 2.11764706 vlo:2.17300000 37/ 17 = 2.17647059 vhi:2.17700000 ok
18/ 8 = 2.25000000 vlo:2.37300000 19/ 8 = 2.37500000 vhi:2.37700000 ok
114/ 33 = 3.45454545 vlo:3.48400000 115/ 33 = 3.48484848 vhi:3.48700000 ok
4/ 1 = 4.00000000 vlo:4.00000000 4/ 1 = 4.00000000 vhi:4.87000000 ok
4/ 1 = 4.00000000 vlo:4.00000000 8/ 1 = 8.00000000 vhi:8.00000000 ok
11/ 2 = 5.50000000 vlo:5.50000000 11/ 2 = 5.50000000 vhi:5.60000000 ok
5/ 1 = 5.00000000 vlo:5.50000000 6/ 1 = 6.00000000 vhi:6.50000000 ok
-7/ -1 = 7.00000000 vlo:7.50000000 -7/ -1 = 7.00000000 vhi:7.30000000 wrong
15/ 2 = 7.50000000 vlo:7.50000000 15/ 2 = 7.50000000 vhi:7.50000000 ok
8030/ 941 = 8.53347503 vlo:8.53453700 8031/ 941 = 8.53453773 vhi:8.53453800 ok
24880/2663 = 9.34284641 vlo:9.34322100 24881/2663 = 9.34322193 vhi:9.34322200 ok
如果不是一个范围内最简单的分数,而是在给定分母大小上限的情况下寻求最佳近似,请考虑以下代码,它将替换def test(vlo, vhi)
向前的所有代码。
def smallden(target, maxden):
global pas
pas = 0
tol = 1/float(maxden)**2
while 1:
den = simpleratio(target-tol, target+tol);
if den <= maxden: return den
tol *= 2
pas += 1
# Test driver for smallden(target, maxden) routine
import random
totalpass, trials, passes = 0, 20, [0 for i in range(20)]
print 'Maxden Num Den Num/Den Target Error Passes'
for i in range(trials):
target = random.random()
maxden = 10 + round(10000*random.random())
den = smallden(target, maxden)
num = int(round(target*den))
got = float(num)/den
print '{:4d} {:4d}/{:4d} = {:10.8f} = {:10.8f} + {:12.9f} {:2}'.format(
int(maxden), num, den, got, target, got - target, pas)
totalpass += pas
passes[pas-1] += 1
print 'Average pass count: {:0.3}nPass histo: {}'.format(
float(totalpass)/trials, passes)
在生产代码中,去掉对pas
(等)的所有引用,即去掉通过计数代码。
例程CCD_ 7被赋予允许分母的目标值和最大值。给定maxden
分母的可能选择,可以合理地假设可以实现1/maxden²
阶的公差。以下典型输出中显示的通过计数(其中target
和maxden
是通过随机数设置的)表明,在超过一半的时间内立即达到了这样的容差,但在其他情况下,使用了2倍、4倍或8倍大的容差,需要对simpleratio
进行额外调用。请注意,10000个数字的测试运行的最后两行输出显示在20个数字测试运行的完整输出之后。
Maxden Num Den Num/Den Target Error Passes
1198 32/ 509 = 0.06286837 = 0.06286798 + 0.000000392 1
2136 115/ 427 = 0.26932084 = 0.26932103 + -0.000000185 1
4257 839/2670 = 0.31423221 = 0.31423223 + -0.000000025 1
2680 449/ 509 = 0.88212181 = 0.88212132 + 0.000000486 3
2935 440/1853 = 0.23745278 = 0.23745287 + -0.000000095 1
6128 347/1285 = 0.27003891 = 0.27003899 + -0.000000077 3
8041 1780/4243 = 0.41951449 = 0.41951447 + 0.000000020 2
7637 3926/7127 = 0.55086292 = 0.55086293 + -0.000000010 1
3422 27/ 469 = 0.05756930 = 0.05756918 + 0.000000113 2
1616 168/1507 = 0.11147976 = 0.11147982 + -0.000000061 1
260 62/ 123 = 0.50406504 = 0.50406378 + 0.000001264 1
3775 52/3327 = 0.01562970 = 0.01562750 + 0.000002195 6
233 6/ 13 = 0.46153846 = 0.46172772 + -0.000189254 5
3650 3151/3514 = 0.89669892 = 0.89669890 + 0.000000020 1
9307 2943/7528 = 0.39094049 = 0.39094048 + 0.000000013 2
962 206/ 225 = 0.91555556 = 0.91555496 + 0.000000594 1
2080 564/1975 = 0.28556962 = 0.28556943 + 0.000000190 1
6505 1971/2347 = 0.83979548 = 0.83979551 + -0.000000022 1
1944 472/ 833 = 0.56662665 = 0.56662696 + -0.000000305 2
3244 291/1447 = 0.20110574 = 0.20110579 + -0.000000051 1
Average pass count: 1.85
Pass histo: [12, 4, 2, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
10000个测试运行的最后两行输出:
Average pass count: 1.77
Pass histo: [56659, 25227, 10020, 4146, 2072, 931, 497, 233, 125, 39, 33, 17, 1, 0, 0, 0, 0, 0, 0, 0]