algorithm - 将实数转换为有理数时检查终止

Question

最近我在几年前写的一些代码中发现了这一点。它用于通过确定合适的分母，然后检查原始实数和有理数之间的差异是否足够小，来合理化一个实数（在容差内）。

编辑澄清：我实际上不想转换所有实际值。例如，我可以选择 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。

score 1 · Accepted Answer

这可以被认为是错误的 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

如果初始猜测足够接近，则目标函数应该是凸的。

正如 Jacob 所建议的，通过最小化以下误差函数可以更好地解决这个问题：

ArgMin ( p * r - q ), where r is real, q and p are integers

这是线性规划，可以通过任何 ILP（整数线性规划）求解器有效求解。 GD 适用于非线性情况，但在线性问题中缺乏效率。

初始猜测和停止条件可以与上述类似。对于求解器的个人选择可以获得更好的选择。

我建议您仍然应该在局部最小值附近假设凸性，这可以大大降低成本。您也可以尝试 Simplex 方法，该方法非常适合线性规划问题。

在这点上，我把功劳归功于 Jacob。

score 0 · Accepted Answer

与此类似的问题在大约开始的近似部分中得到解决。Bill Gosper 的连 分数算术文档的第 28 页。（参考：后记文件；另见文本版本，从 1984 行开始。）一般的想法是计算低端和高端范围限制数的连分数近似值，直到两个分数不同，然后选择一个值在这两个近似值的范围内。这保证给出一个最简单的分数，使用 Gosper 的术语。

下面的 python 代码（程序“simpleden”）实现了类似的过程。（它可能不如 Gosper 建议的实现好，但足以让您看到该方法产生什么样的结果。）完成的工作量类似于欧几里德算法的工作量，即O(n) 与n 位，所以程序相当快。一些示例测试用例（即程序的输出）显示在代码本身之后。请注意，simpleratio(vlo, vhi)如果 vhi 小于 vlo，则此处显示的函数返回 -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（等）的引用，即删除通过计数代码。

该例程smallden被赋予目标值和允许分母的最大值。给定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]

algorithm - 将实数转换为有理数时检查终止

2 回答 2

Related

Reference