パイこね変換を分数演算で"正しく"変換する

programming python

『渋滞学』を最後まで読みきらずにいたのを思い出して読んでみた。
『渋滞学』によるとパイこね変換というアルゴリズムがあり、”コンピュータは原理的にこの計算はできない”らしい。

0 <= x <= 1
while True:
  if x < 0.5: x = 2x
  else: x = 2 - 2x
  print x

この計算は浮動小数点の誤差をどんどん大きくさせるのでだんだん誤差がひどくなりまともな結果にならなくなるらしい。
Pythonのインタプリタで確認した。

x = 0.1
>>> for _ in xrange(100):
...  if x < 0.5: x = 2*x
...  else: x = 2 - 2*x
...  print x
...
0.2
0.4
0.8
0.4
0.8
0.4
0.8
0.4
0.8
0.4
0.8
0.4
0.8
0.4
0.8
0.4
0.799999999999
0.400000000001
0.800000000003
0.399999999994
0.799999999988
0.400000000023
0.800000000047
0.399999999907
0.799999999814
0.400000000373
0.800000000745
0.39999999851
0.79999999702
0.40000000596
0.800000011921
0.399999976158
0.799999952316
0.400000095367
0.800000190735
0.39999961853
0.799999237061
0.400001525879
0.800003051758
0.399993896484
0.799987792969
0.400024414062
0.800048828125
0.39990234375
0.7998046875
0.400390625
0.80078125
0.3984375
0.796875
0.40625
0.8125
0.375
0.75
0.5
1.0
0.0
0.0
0.0
...

たしかにー。本ではこの問題をあげて

コンピュータが必ず計算ミスをするような例題を簡単に作れるのでコンピュータを過信してはいけない。

とある。そーかも。でもこれって浮動小数点の誤差の問題だから整数演算すれば計算できるだろ、とか思って粗末な分数クラスで計算してみた。

# fraction.py
def gcd(x, y):
    while y!=0:
        x,y = y, x%y
    return x
class Fraction:
    def __init__(self, n, d):
        gcdnd = gcd(n, d)
        self.num = n / gcdnd
        self.den = d / gcdnd
    def __repr__(self):return self.__str__()
    def __str__(self):
        if self.den == 1:
            return "%d" % (self.num)
        else:
            return "%d/%d (%f)" % (self.num, self.den, float(self.num)/self.den)
    def __mul__(self, n):
        return Fraction(self.num*n, self.den)
    def __sub__(self, f):
        return Fraction(self.num * f.den - self.den * f.num, self.den * f.den)
    def __lt__(self, x):
        a = float(self.num) / self.den
        b = float(x.num)/ x.den
        return a < b
    def __eq__(self, x):
        return self.num == x.num and self.den == x.den
    def __gt__(self, x):
        return not self.__eq__(x) and not self.__lt__(x)

def main():
    import sys
    if sys.argv.__len__() != 4:
        print "usage: %s numerator denominator iterate" % __file__
        sys.exit()
    x = Fraction(*[int(x) for x in sys.argv[1:3]])
    iterate = int(sys.argv[3])

    onesecond = Fraction(1, 2)
    two = Fraction(2, 1)
    print x
    for _ in xrange(iterate):
        if x < onesecond:
            x = x * 2
        else:
            x = two - x * 2
        print x

if __name__ == '__main__':
    main()

初めて分数演算のコードを見たのはSICPだったのだけれど、それと比べて上のコードはすごく汚く見える。コーディングが悪いのだろうか。それともSchemeがそういう問題に向いているのだろうか。

$ python fraction.py
usage: fraction.py numerator denominator iterate
$ python fraction.py 1 10 100
1/10 (0.100000)
1/5 (0.200000)
2/5 (0.400000)
4/5 (0.800000)
...(略)
4/5 (0.800000)
2/5 (0.400000)
$ python fraction.py 1 10 1000000
1/10 (0.100000)
1/5 (0.200000)
2/5 (0.400000)
4/5 (0.800000)
...(略)
4/5 (0.800000)
2/5 (0.400000)

できる。これなら1万回だろうが100万回だろうがちゃんと変換できる。

面白そうなページを見つけたのであとで読む。
http://gascon.cocolog-nifty.com/blog/2008/09/post_e202.html