浮点数到分数的转换（一）

---

     暂停一下subspace.m的连载。歇一歇，换个话题。

     在Matlab的Command Window中输入命令
format rat
则所有的结果都会以分数的形式显示出来，如：
>> sqrt(8)

ans =

    3363/1189 
在Python中，我们也可以用module fractions来实现浮点数到分数的转换，如：
In [47]: Fraction.from_float(1.1)
Out[47]: Fraction(2476979795053773, 2251799813685248)

     这种转换是如何实现的呢？在Matlab中大概是用C实现的，在Python中是用float类型的as_integer_ratio实现的，而这个as_integer_ratio估计也是用C实现的，目前手头没有这些C代码，所以只好自己想办法了。还是先看看浮点数的存储机制。下面都是以双精度浮点数为例。双精度浮点数长度为64 bits，最高位为符号位S，后面接着11b的指数E，最后是52b的尾数M，即
SEEE EEEE EEEE MMMM ... MMMM
它对应的数值为
(-1)^S + (1.MMMM ... MMMM)₂ * 2^{(EEEEEEEEEEE - 1023)}
当EEEEEEEEEEE - 1023≥52时，此浮点数没有小数部分，因此直接以整数的形式显示即可，如
>> 2^53

ans =

9007199254740992 
如何将指数部分大于等于52的双精度浮点数转换为对应的整数？依然要利用编写eps函数时定义的DoubleBits（参见《关于Matlab的eps函数（七）——“又续”的续》，稍微修改了一些）：

import ctypes as ct

# --xialulee--
# --2010.09.19--

class DoubleBits(ct.Union):
    class Bits(ct.Structure):
        _fields_ = [('ML',  ct.c_uint, 32),
                    ('MH',  ct.c_uint, 20),
                    ('E',   ct.c_uint, 11),
                    ('S',   ct.c_uint, 1)
                    ]
    _fields_ = [('value',   ct.c_double),
                ('bits',    Bits)
                ]
    @property
    def E(self):
        return self.bits.E - 1023

    @property
    def M(self):
        return self.bits.ML + (self.bits.MH << 32) + (1 << 52)

    @property
    def S(self):
        return (-1) ** self.bits.S



     之所以把52位的尾数部分拆分成ML和MH，是因为位段结构体的一个域不能超过CPU字边界。想必在64位机器上就不用这么麻烦了。有了DoubleBits，我们可以轻松操作双精度浮点数的各个域，因此将这样的数字转化为整数非常简单，在Python中还有任意长整数类型的支持，操作起来非常轻松，如
In [123]: r = 2.0 ** 53

In [124]: d = DoubleBits(r)

In [125]: d.M * 2**(d.E - 52)
Out[125]: 9007199254740992L

和Matlab的结果一致。这样做的原理是：
1.MMMM…MMMM × 2^E == 1MMMM…MMMM × 2^-52× 2^E == 1MMMM…MMMM × 2^(E-52)

     对于指数部分小于52的双精度浮点数，它会含有小数部分。回想一下十进制小数一般是怎样转换为分数的，比如3.1416，我们将其小数点右移4位（即相当于乘以10⁴）作为分子，然后以10⁴作为分母，约一下分，就好了。浮点数转分数也是一样，只是IEEE 754标准的浮点数都是以2为基的，而不是10。

     首先我们忽略浮点数的指数部分，仅仅看尾数。由于尾数有52b，因此将尾数转换为分数就是：
1MMMM...MMMM / 2⁵²

现在在把指数部分考虑进去，整个浮点数就应该是：
     1MMMM...MMMM / 2⁵²× 2^{EE...E - 1023}
== 1MMMM...MMMM / 2^{52 - (EE...E - 1023)}
试一试：
In [126]: r = 1.1

In [127]: d = DoubleBits(r)

In [128]: numerator = d.M

In [129]: denom = 2 ** (52 - d.E)

In [130]: numerator
Out[130]: 4953959590107546L

In [131]: denom
Out[131]: 4503599627370496L

     有了分子和分母，还不算完，约分是必须的，在module fractions中，有名为gcd（囧）的函数，用于寻找两个整数的最大公因子（greatest common divisor），其原理很并不复杂，千年之前就已经出现（辗转相除法），代码也很简单，大家看5秒钟就能明白了，这个函数在fractions.py的开头，很容易找到：
In [132]: from fractions import gcd

In [133]: gcd(numerator, denom)
Out[133]: 2L

将分子和分母同时除以公因子，就得到了约分过后的分数了。

In [134]: numerator / 2
Out[134]: 2476979795053773L

In [135]: denom / 2
Out[135]: 2251799813685248L

     这个结果与Fraction.from_float是一致的。但是大家可能心存疑虑，1.1转换成分数应该是11/10，结果怎么转成了：
Fraction(2476979795053773, 2251799813685248)
这是因为IEEE 754是以2为基表示小数的。也就是说IEEE754的小数转换成分数的话，其分母必然是2的整数次幂。但是十进制小数1.1转换成分数的话其分母为10，10含有因子5，因此0.1这样的数不可能用有限长度的二进制小数精确表示。

     Matlab的默认类型也是基于IEEE 754的双精度浮点类型的，但是为什么Matlab可以呢：
>> format rat
>> 1.1

ans =

      11/10

这是因为Matlab在转换的过程中对精度进行了调整。如果我们在Python中也这样做的话，会得到同样的效果：
In [65]: Fraction.from_float(1.1).limit_denominator(1000)
Out[65]: Fraction(11, 10)
至于limit_denominator的算法，下次吧，这次就不戏说了。