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



最近我在几年前写的一些代码中发现了这一点。它被用来通过确定合适的分母来合理化实数(在公差范围内),然后检查原始实数和有理数之间的差是否足够小。

编辑以澄清:我实际上不想转换所有的实数。例如,我可以选择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²阶的公差。以下典型输出中显示的通过计数(其中targetmaxden是通过随机数设置的)表明,在超过一半的时间内立即达到了这样的容差,但在其他情况下,使用了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]

最新更新