一种小数转分数的算法（不限整除）C++

最近需要用到小数转分数算法，便研究了一下。

先看一下最终程序的效果：

说一下数学中有理小数转分数的过程：
有理小数分为有限小数和无限循环小数

1. 有限小数：

有限小数直接去小数点再约分即可。

例：$1.55=\dfrac{155}{100}=\dfrac{31}{20}$

2. 无限循环小数：

先判断循环节长度$n$，原小数$a$乘以$10^n$后再减去$a$可化为有限小数。

例：$1.8123123...$

设 $a=1.8123123...$

循环节为123，一共三位，乘以$10^3$

$10^3a=1812.3123123...$

两式相减可得$999 a=1810.5$

$a=\dfrac{1810.5}{999}=\dfrac{18105}{9990}=\dfrac{1207}{666}$

刚开始想使用上述方法，于是做了一些实验。
计算器算了分数，123/321=0.38317757009345794392523364485981，可见，123、321这两个数并不是很大，但是相除后得到的小数循环节过长，一个double类型尾数长度是52bit，对应十进制只有15~16位的精度，不能保证用算法分析出循环节的长度，上述方法不适用。下面使用连分数做转换。

---

使用连分数做转换

一、连分数的定义：

形如$a_0+\dfrac1{a_1+\dfrac1{a_2+\dfrac1{a_3+\dfrac1{...}}}}$的分式叫做连分式，记为$[a_0;a_1,a_2,a_3,...]$。
有理数（整数、有限小数、无限循环小数）的连分式是有限的。

二、有理数与连分数的转换：

1. 分数或小数转换为连分数

步骤1：令a=这个数，即a←这个数；
步骤2：将a的整数部分⌊a⌋记录下来；
步骤3：令a为a的小数部分，即a←a-⌊a⌋；
步骤4：如果a≠0，跳转到步骤2，否则向下执行步骤5；
步骤5：步骤2中记录的数即为连分数的$[a_0;a_1,a_2,a_3,...]$。
例如将$\dfrac{321}{123}$和$2.6097560975...$转换为连分数，步骤如下表：

$\dfrac{321}{123}$	$2.6097560975...$
$\dfrac{321}{123}=2\dfrac{75}{123}$	$2.6097560975...=2+0.6097560975$ $1/0.6097560975...=1.64$
$\dfrac{123}{75}=1\dfrac{48}{75}$	$1.64=1+0.64$ $1/0.64=1.5625$
$\dfrac{75}{48}=1\dfrac{27}{48}$	$1.5625=1+0.5625$ $1/0.5625=1.7777...$
$\dfrac{48}{27}=1\dfrac{21}{27}$	$1.7777...=1+0.7777...$ $1/0.7777=1.285714285714...$
$\dfrac{27}{21}=1\dfrac6{21}$	$1.285714285714...=1+0.285714285714...$ $1/0.285714285714...=3.5$
$\dfrac{21}6=3\dfrac36=3\dfrac12$	$3.5=3+0.5$ $1/0.5=2$
$\dfrac21=2$	$2=2+0$

可见$\dfrac{321}{123}=[2;1,1,1,1,3,2]=2+\dfrac1{1+\dfrac1{1+\dfrac1{1+\dfrac1{1+{\dfrac1{3+\dfrac1{2}}}}}}}$。

2. 连分数转换为分数或小数

直接通分即可：如$[2;1,1,1,1,3,2]=\frac{107}{41}$。

3. 利用上述方法，将有理数转为连分数再通分，可将$2.6097560975...$转为分数$\frac{107}{41}$。

三、转换算法的C++程序

class ContinuedFraction
{
private:
    //用于保存连分数列表
    std::vector<unsigned> nums;
public:
    //小数转为连分数的构造函数
    ContinuedFraction(double a);
    //连分数转为小数
    double ToDouble();
};

ContinuedFraction::ContinuedFraction(double a)
{
    if (a < 1e-12)
    {
        nums.push_back(0);
        return;
    }
    for (;;)
    {
        unsigned tr = (unsigned)a;
        nums.push_back(tr);
        a -= tr;
        if (a < 1e-12) //小数部分是0, 退出
        {
            return;
        }
        a = 1 / a;
    }
}

double ContinuedFraction::ToDouble()
{
    double a = nums.back();
    for (auto it = nums.rbegin() + 1; it != nums.rend(); it++)
    {
        a = 1 / a;
        a += *it;
    }
    return a;
}

四、计算机转换时的问题

上述转换算法，我们需要做步骤3中的判断（a≠0）才能知道什么时候跳出循环，由于浮点型的误差，程序中a与0作比较应该这样写：abs(a - 0) < 1e-12，其中1e-12为很小的数，也就是说浮点型a与b比较相等，由于误差，需要写成a与b的差的绝对值是否充分小。
程序中的unsigned tr = (unsigned)a; a -= tr;功能是取a中的小数部分，初始浮点型为15位的有效数字，但减去整数部分后会引起浮点型的误差增大。
例如123.00123（8个有效数字）取小数后为0.00123（3个有效数字），有效数字的个数减少了5。
也就是说有效数字减少的个数为 整数的总位数3(123为3位)+小数点后零的个数2(0.00123有2个0)=5。
由于循环的次数，误差将越来越大，所以代码中应该计算实时的误差为多少，才能使a≠0的判断更有效。
程序如下

//初始浮点型的误差
double esp = 1e-15;
for (;;)
{
    unsigned tr = (unsigned)a;
    nums.push_back(tr);
    a -= tr;
    if (a < 1e-12) //小数部分是0, 直接退出
    {
        return;
    }
    //因a减整数部分而引起的浮点型误差增大
    //例如123.123，变为0.123，精度也由6位变为3位
    while (tr > 0)
    {
        tr /= 10;
        esp *= 10;
    }
    //因a取小数引起的浮点型误差增大
    //例如123.00123，变为0.00123，精度也由8位变为3位
    double t = a;
    while (t < 0.1)
    {
        t *= 10;
        esp *= 10;
    }
    a = 1 / a;
}

上面代码已经实时的计算出浮点数的误差esp，理论上在最后一句a = 1 / a;前面加上if (a < esp) return;即可。但是我们计算浮点数的误差过于严谨，可能导致还没有到0就提前退出循环了。所有这里我们把a < esp时候的插入nums的索引indexZheng记录下来，向后再多计算几次一直到2*indexZheng后再退出循环。
转换为分数时，可以分别将nums个数取indexZheng到2*indexZheng对应的double值与给定的double值比较，取误差最小对应的分数。

---

最终的程序如下：

ContinuedFraction.h

#pragma once
#include <vector>

class ContinuedFraction
{
private:
    //用于保存连分数列表
    std::vector<unsigned> nums;
    int indexZheng = -1;
public:
    //小数转为连分数的构造函数
    ContinuedFraction(double a);
    //连分数转为小数
    double ToDouble();
    friend class Fraction;
};

ContinuedFraction.cpp

#include "ContinuedFraction.h"

ContinuedFraction::ContinuedFraction(double a)
{
    if (a < 1e-12)
    {
        nums.push_back(0);
        return;
    }
    //初始浮点型的误差，经测试1e-12效果最好
    double esp = 1e-12;
    for (int i = 0; indexZheng == -1 || i <= 2 * (indexZheng + 5); i++)
    {
        unsigned tr = (unsigned)a;
        nums.push_back(tr);
        a -= tr;
        if (a < 1e-12) //小数部分是0, 直接退出
        {
            if (indexZheng == -1)
            {
                indexZheng = i;
            }
            return;
        }
        //因a减整数部分而引起的浮点型误差增大
        //例如123.123，变为0.123，精度也由6位变为3位
        while (tr > 0)
        {
            tr /= 10;
            esp *= 10;
        }
        因a取小数引起的浮点型误差增大
        例如123.00123，变为0.00123，精度也由8位变为3位
        //double t = a;
        //while (t < 0.1)
        //{
        //  t *= 10;
        //  esp *= 10;
        //}
        if (indexZheng == -1 && a < esp)
        {
            indexZheng = i;
        }
        a = 1 / a;
    }
}

double ContinuedFraction::ToDouble()
{
    double a = nums.back();
    for (auto it = nums.rbegin() + 1; it != nums.rend(); it++)
    {
        a = 1 / a;
        a += *it;
    }
    return a;
}

Fraction.h

#pragma once
#include "ContinuedFraction.h"

class Fraction
{
private:
    //符号位
    bool sign;
    //倒数
    void Invert();
public:
    //分子
    unsigned a;
    //分母
    unsigned b;
    Fraction() {}
    Fraction(double a);
    Fraction(ContinuedFraction confrac);
    void Print();
    double ToDouble();
};

Fraction.cpp

#include "Fraction.h"
#include <cmath>
#include <vector>
#include <iostream>

Fraction::Fraction(ContinuedFraction confrac) : b(1), sign(false)
{
    a = confrac.nums.back();
    for (auto it = confrac.nums.rbegin() + 1; it != confrac.nums.rend(); it++)
    {
        Invert();
        a += b * *it;
    }
}

Fraction::Fraction(double a)
{
    sign = a < 0;
    if (sign) a = -a;

    if (a < 1e-12)
    {
        this->a = 0;
        this->b = 1;
        this->sign = false;
        return;
    }
    ContinuedFraction confrac(a);
    std::vector<unsigned> nums = confrac.nums;

    int firstIndex = confrac.indexZheng - 4;
    firstIndex = firstIndex >= 0 ? firstIndex : 0;
    confrac.nums.resize(firstIndex);

    std::vector<Fraction> vfr;
    std::vector<double> vesp;
    //一次一次增加连分数精度，并计算误差
    for (size_t i = firstIndex; i < nums.size(); i++)
    {
        confrac.nums.push_back(nums[i]);
        Fraction fr(confrac);
        double d = fr.ToDouble();
        vfr.push_back(fr);
        vesp.push_back(abs(d - a));
    }
    //查找误差最小的index
    int mindex = 0;
    double m = 1e100;
    for (size_t i = 0; i < vesp.size(); i++)
    {
        if (vesp[i] < m)
        {
            mindex = i;
            m = vesp[i];
        }
    }
    this->a = vfr[mindex].a;
    this->b = vfr[mindex].b;
    return;
}

void Fraction::Invert()
{
    unsigned t = a;
    a = b;
    b = t;
}

void Fraction::Print()
{
    std::cout << (sign ? "-" : "") << a << "/" << b << std::endl;
}

double Fraction::ToDouble()
{
    double n = (double)a / b;
    return sign ? -n : n;
}

main.cpp

#include "Fraction.h"

int main()
{
    Fraction f1(46531463. / 34554517);
    f1.Print();

    Fraction f2(4.25);
    f2.Print();

    Fraction f3(-1.3333333333333);
    f3.Print();

    return 0;
}