从抽象谈面向对象与泛型编程 part.2 - 泛型编程的抽象
22 Jul 20181. 前言
上一篇浅讨了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的实现。先看看欧几里得整环的定义:
- 是一个整环;
- 定义了带余除法:a == quotient(a, b) * b + remainder(a, b);
- 存在范数norm(a), 并且满足以下条件:
- let norm(a) == 0, thus a == T{0},
- let 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抽象的分析与比较,放在下一篇中论述。