从抽象谈面向对象与泛型编程 part.2 – 泛型编程的抽象

1. 前言

上一篇浅讨了OOP中抽象常用的方法，而这一篇的内容，主要是对《数学与泛型编程》，后文简称MGP，中部分内容的小结。这本著作有相当多的篇幅，是对数论基础的讨论。其中的一条主线，就是最大公约数(GCD)算法的探讨。GCD是数论的基础，可以说绝大部分数论的定理都是基于GCD推导出来的，而整个数论的发展史就是一部推广GCD算法的历史。GCD算法推广的过程，也演绎出了很多数学中常用的抽象方法。而GP的抽象，从一个具体的算法或者业务出发，抽象出一个适用于一个类型集合的外观，这个过程与数学定理的推广过程很类似。所以这一篇，我们将以最大公约数算法的推广为主线，来探讨泛型编程中的抽象方法。

2. 最大公约数算法

GCD早在古希腊，毕达哥斯拉学派中，就有了很深的研究。古希腊时的人们还没有抽象出整数这个数学工具，甚至连零的概念都没有。那个时候人类研究的都是具体的问题，例如使用单位长度的绳子进行测量。GCD就是其中一个很经典的问题，求出两条线段的公尺量。下面的代码，实现的是最经典的GCD算法（ 这里假设a的长度大于b，即使是b大于a，我们也可以通过交换位置得到相同的条件，因此后面的GCD版本都假设a大于b ）

line_segment gcd(line_segment a, line_segment b) {
    while(a != b) {
        a = remainder(a, b);
        std::swap(a, b);
    }
    return a;
}

GCD的原始版本可以求出两条线段的公尺量，就是求出一条能够同时测量两条线段的单位线段。公尺量作为单位1，建立起了一个类似正整数的模型，因此零没有引入到系统中来。对于非负整数集合而言，零是必须要考虑的情况，下面的代码就是计算两个32位无符号整形的GCD代码：

uint32_t gcd(uint32_t a, uint32_t b) {
    while(b != 0) {
        a = a % b;
        std::swap(a, b);
    }
    return a;
}

这里有一个不严谨的地方，就是无符号整形无法映射到整个非负整数的集合，因为前者是一个有限集合。但是对于gcd算法而言，是不会出现算法溢出的，原因是说给定GCD的两个入参a和b，是不会计算出比a和b更大值的情况。因此我们可以在这里忽略计算机无符号整形不能与非负整数同构。

回到正题，算法很简单，这是一个辗转求余的过程。数论的研究是要把它推广到适应更广泛的对象，而对于程序的实现而言，就是要推广它能够适应更多的类型。添加一个支持64位整数的版本，利用函数重载就很容易实现：

uint64_t gcd(uint64_t a, uint64_t b) {
    while(b != 0) {
        a = a % b;
        std::swap(a, b);
    }
    return a;
}

3. 支持有符号的整数

当我们把GCD从非负整数推广到整个整数集合的时候，其实就是确定如何对负整数求余的操作的。如果把全体整数列在一条直线上，把每个整数i看作是一个从0出发，以1位单位，长度为i的有向线段。那么负整数的绝对值，就是它的长度了。我们把GCD推广到整数的方法，就是把对数求余变成对数的度量进行求余。

int gcd(int a, int b) {
    a = std::abs(a);
    b = std::abs(b);
    while(b != 0) {
        a = a % b;
        std::swap(a, b);
    }
    return a;
}

当然，我们目前所使用的计算机体系，是可以使用对负数取模操作来计算余数的，实际上我们可以一行代码都不需要修改。

int gcd(int a, int b) {
    while(b != 0) {
        a = a % b;
        std::swap(a, b);
    }
    return a;
}

4. 支持更多的整数类型

之前版本的GCD就可以处理所有全部的整形了，但它也只能处理整形。原因就是我们使用了取模操作来计算余数，这个只能对整形数据求余，但对其他整数类型的数据无效，例如高斯整数。

高斯整数是形如x = m + in的整数，它有实部和虚部两个部分。我们可以使用标准库已经实现的complex来表示高斯整数：

using gint32_t = std::complex<int32_t>;

推广GCD到高斯整数，我们就要实现高斯整数的带余除法的求余部分：

gint32_t remainder(gint32_t a, gint32_t b) {
    guint32_t quotient = (a * std::conj(b)) / std::norm(b);
    guint32_t remainder = a - b * quotient;
    return remainder;
}

有了求余的函数，就可以实现高斯整数的GCD算法了：

gint32_t gcd(gint32_t a, gint32_t b) {
    while(b != gint32_t{ 0 }) {
        a = remainder(a, b);
        std::swap(a, b);
    }
    return a;
}

6. 支持多项式

现在我们的GCD算法可以处理整形和高斯整数了，但是GCD还可以处理更多的问题。例如实数域上的一元多项式。先给出多项式的一个简单外观：

template <typename T>
class polynomial {
public:
    using value_type = std::vector<T>;
    using size_type = std::vector<T>::size_type;
    using element_type = T;
    // f(0) = coef
    polynomial(element_type coef) : coefs_({coef}) {} 
    // construct x^order
    polynomial(element_type coef, size_type order); 
    // construct from { c1, c2, ..., cn }
    polynomial(std::initializer_list<element_type> coefs); 
    size_type degree() const;   // the degree of x^3 + x is 3
    // other interfaces ...
private:
    value_type coefs_;
};
 
// operator + - *
template <typename T>
inline static polynomial<T> operator+ (
    polynomial<T> const& lhs, polynomial<T> const& rhs);
template <typename T>
inline static polynomial<T> operator- (
    polynomial<T> const& lhs, polynomial<T> const& rhs);
template <typename T>
inline static polynomial<T> operator* (
    polynomial<T> const& lhs, polynomial<T> const& rhs);

这是一个简单polynomial的实现，还有一些copy assign的接口和实现等，受限于篇幅，这里就不给出详细列出了。为了让GCD算法在这里可以正常工作，我们需要实现一个多项式除法，来求解余式。这里选择长除法作为一个简单的示例。

inline static polynomial<double> remainder(
    polynomial<double> const& a, polynomial<double> const& b) {
    auto quotient = polynomial<double>{ 0 };
    auto remainder = a;
    while(remainder != 0 && remainder.degree() > b.degree()) {
        auto coef = a.front() / b.front();
        auto order = a.degree() - b.degree();
        auto t = T{ coef, order };
        quotient += t;
        remainder -= t * b;
    }
    return remainder;
}

有了求余式的算法，针对多项式的GCD算法，就可以实现了：

polynomial<double> gcd(
    polynomial<double> a, polynomial<double> b) {
    while(b != polynomial<double>{ 0 }) {
        a = remainder(a, b);
        std::swap(a, b);
    }
    return a;
}

6. 抽象统一的外观

代码写到这里，可以发现我们已经实现的几个GCD的版本，算法的步骤几乎都是相同的，都是辗转余，并以b为零作为终止条件。不同的只是类型。我们可以从这里开始使用泛型编程的技术进行抽象了。抽象的统一外观，可以如同下列代码：

template <typename T>
T gcd(T a, T b) {
    while( b != T{ 0 }) {
        a = remainder(a, b);
        std::swap(a, b);
    }
    return a;
}

算法的核心在于remainder能够正确操作，所以还需要对remainder进一步的完善，使得它可以正确处理全部的整形，高斯整数和实数域的一元多项式：

template <typename T>
auto remainder(T a, T b) ->
    std::enable_if_t<std::is_integral_v<T>, T> {
    return a % b;
}
 
template <typename T>
auto remainder(T a, T b) ->
    std::enable_if_t<is_gauss_integral_v<T>, T> {
    T quotient = (a * std::conj(b)) / std::norm(b);
    T remainder = a - b * quotient;
    return remainder;
}
 
template <typename T>
auto remainder(T const& a, T const& b) -> 
    std::enable_if_t<is_real_polynomial_v<T>, T> {
    // ... implementation omitted here
    return remainder;
}

这里使用了C++模板的SFINAE特性，使用std::enable_if来控制不同静态分支的生成。也就是说对于普通整形，上面的版本才会被实例化；而对于高斯整数，中间的版本才会被实例化；对于实数域的多项式，下面的版本才会生成实际的代码。剩下一个判断T类型是否是高斯整数或者多项式的元函数了。

// meta function to check if T is a type of gauss integer
template <typename T>
inline constexpr bool is_gauss_integral_v = 
    std::conjunction<
        is_instantiate_of<T, std::complex>,
        std::is_signed<safe_value_type_t<T>>::value;
// meta function to check if T is a type of real polynomial
template <typename T>
inline constexpr bool is_real_polynomial_v = 
    std::conjunction_v<
        is_instantiate_of<T, polynomial>,
        std::is_floating_point<safe_element_type_t<T>>
    >;

元函数是编译期对类型和编译期常量计算的函数。例如std::is_integral可以返回一个类型是否是整形。而我们这里判断一个类型是否是高斯整数，可以先判断这个类型是否是std::complex的一个实例化类型，并且std::complex的模板参数T必须是一个有符号整形。std::conjunction是编译期的逻辑与元函数。关于模板元编程的话题，我们会在其他系列的文章详细讨论。

7. 抽象能推广到多远？

到目前为止，我们抽象出了一个GCD算法的外观，那么这个抽象可以推广到多远？MGP上告诉我们，GCD可以推广到整个欧几里得整环。那么，我们可以得到一个目前为止，最通用的GCD算法的外观。

template <EuclideanDomain T>
constexpr T gcd(T a, T b) {
    while(b != T{ 0 }) {
        a = remainder(a, b);
        std::swap(a, b);
    }
    return a;
}

constexpr是C++11引入的编译期求值的关键字，remainder和quotient都应该加上constexpr修饰符，并且GCD需要处理的所有类型，都应该支持constexpr构造与operator. 除了我们目前的多项式的实现需要做比较大的更改来支持constexpr，数学中的大部分数值类型添加编译期的支持，对目前C++17标准不是一件难事。

template< EuclideanDomain T >是C++将要到来的新标准Concepts. EuclideanDomain对类型T有一系列的约束要求，有了这个强大的工具，我们可以很容易地对一个类型集合做抽象了。这里不详细展开Concepts的细节，不过可以给出将来可能的EuclideanDomain Concepts的实现。先看看欧几里得整环的定义：

1. 是一个整环；

2. 定义了带余除法：a == quotient(a, b) * b + remainder(a, b)；

3. 存在范数norm(a), 并且满足以下条件:

* if norm(a) == 0, thus a == T{0}，

* if b != 0, thus norm(a * b) >= norm(b)，

* 余式的范数小于除式的范数: norm(remainder(a, b)) < norm(b).

由公理可知，remainder的范数一定小于除式的范数，所以GCD算法的迭代肯定会终止。受篇幅所限，这里假设已经有了整环IntegralDomain的Concepts实现，那欧几里得整环可能的实现如下：

template <typename T>
concept bool EuclideanDomain = IntegralDomain<T> &&
    requires(T a, T b) {
        { quotient(a, b) } -> T;
        { remainder(a, b) } -> T;
        requries (norm(T{0}) == 0);
        requires (a != 0 && b != 0);
        requries norm(a * b) > norm(b);
        requries norm(remainder(a, b)) < norm(b);
    };

目前，concepts标准还处于TS状态，所以还只能用C++的SFINAE进行简单的模拟。concepts相对于SFINAE，除了可以对类型所具有的语法层面上的约束之外，还有语义的约束，语义的约束又可以称为契约编程。这种机制可以极大的增强GP的抽象能力。

8. 小结

鄙文通过简单实现了整数，高斯整数和多项式的GCD算法，从中抽象出了通用GCD算法的外观，并将其适用到了欧几里得整环。可以看出GP的抽象过程，与数学理论与模型的推广过程非常相似。篇幅所限，对GP和OOP抽象的分析与比较，放在下一篇中论述。

Post Views: 1