1 导言
模块化运算是现代代数系统的基本工具。结合中国余数定理,它是计算 gcd、结果等算法的主力。此外,它还可以作为一种非常有效的过滤器,因为通常只需计算一个素数的模数对应值,就可以排除某个值为零的可能性。
2 留数和模板化
先,本软件包引入了一个 Residue 类型。质数 p 保存在一个静态成员变量中。该类提供静态成员函数来更改该值。
改变质数会使已存在的该类型对象失效。但是,已经存在的对象不会失去其相对于旧质数的价值,并且在恢复旧质数后可以重新使用。由于该类型是基于双运算的,因此质数的值只能小于 226。p 的初始值为 67108859。
此外,软件包还引入了可模块化概念。如果存在从 T 到基于残差类型的代数结构的映射,那么代数结构 T 就被认为是可模块化的。对于标量类型(如整数),这种映射只是进入由 Residue 表示的 Z/pZ 的典型同态映射。对于复合类型(如多项式),该映射适用于复合类型的系数。该映射由类 Modular_traits<T> 提供。Modular_traits<T> 类的设计使 Modularizable 概念可被视为可选概念,也就是说,Modular_traits<T> 提供了一个可用于调度的标记。
2.1 示例
在下面的例子中,为了避免对多项式进行不必要的 gcd 计算,使用了模块化算术作为过滤器。在欧几里得算法中,由于系数的增长,gcd 计算的成本会非常高。
一般的思路是,首先只计算一个素数的 gcd。如果这个模块 gcd 是常数,我们(在大多数情况下)就可以得出结论,实际的 gcd 也是常数。
为此,示例引入了函数 may_have_common_factor()。请注意,这个函数有两个版本,分别适用于系数类型可模块化和不可模块化的情况。如果类型不可模块化,则不应用过滤器,函数返回 true。
需要进一步注意的是,类 Residue 的实现要求尾数精度符合 IEEE 浮点运算标准(IEEE 754)。然而,在某些处理器上,传统 FPU 使用的是扩展精度。因此,在执行任何算术运算之前,必须执行适当的尾数长度。此外,数字必须四舍五入到最接近的数值。可以使用带有 CGAL_FE_TONEAREST 的 Protect_FPU_rounding 来确保这一点,它的副作用还包括强制执行所需的精度。
File Modular_arithmetic/modular_filter.cpp
#include <CGAL/config.h>
#include <iostream>
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpz.h>
#include <CGAL/Polynomial.h>
// Function in case Polynomial is Modularizable
template< typename Polynomial >
bool may_have_common_factor(const Polynomial& p1, const Polynomial& p2, CGAL::Tag_true){std::cout<< "The type is modularizable" << std::endl;// Enforce IEEE double precision and rounding mode to nearest// before using modular arithmeticCGAL::Protect_FPU_rounding<true> pfr(CGAL_FE_TONEAREST);// Use Modular_traits to convert to polynomials with modular coefficientstypedef CGAL::Modular_traits<Polynomial> MT;typedef typename MT::Residue_type MPolynomial;typedef typename MT::Modular_image Modular_image;MPolynomial mp1 = Modular_image()(p1);MPolynomial mp2 = Modular_image()(p2);// check for unlucky primes, the polynomials should not lose a degreetypename CGAL::Polynomial_traits_d<Polynomial>::Degree degree;typename CGAL::Polynomial_traits_d<MPolynomial>::Degree mdegree;if ( degree(p1) != mdegree(mp1)) return true;if ( degree(p2) != mdegree(mp2)) return true;// compute gcd for modular imagesMPolynomial mg = CGAL::gcd(mp1,mp2);// if the modular gcd is not trivial: return trueif ( mdegree(mg) > 0 ){std::cout << "The gcd may be non trivial" << std::endl;return true;}else{std::cout << "The gcd is trivial" << std::endl;return false;}
}
// This function returns true, since the filter is not applicable
template< typename Polynomial >
bool may_have_common_factor(const Polynomial&, const Polynomial&, CGAL::Tag_false){std::cout<< "The type is not modularizable" << std::endl;return true;
}
template< typename Polynomial >
Polynomial modular_filtered_gcd(const Polynomial& p1, const Polynomial& p2){typedef CGAL::Modular_traits<Polynomial> MT;typedef typename MT::Is_modularizable Is_modularizable;// Try to avoid actual gcd computationif (may_have_common_factor(p1,p2, Is_modularizable())){// Compute gcd, since the filter indicates a common factorreturn CGAL::gcd(p1,p2);}else{typename CGAL::Polynomial_traits_d<Polynomial>::Univariate_content content;typename CGAL::Polynomial_traits_d<Polynomial>::Construct_polynomial construct;return construct(CGAL::gcd(content(p1),content(p2))); // return trivial gcd}
}
int main(){CGAL::IO::set_pretty_mode(std::cout);typedef CGAL::Gmpz NT;typedef CGAL::Polynomial<NT> Poly;CGAL::Polynomial_traits_d<Poly>::Construct_polynomial construct;Poly f1=construct(NT(2), NT(6), NT(4));Poly f2=construct(NT(12), NT(4), NT(8));Poly f3=construct(NT(3), NT(4));std::cout << "f1 : " << f1 << std::endl;std::cout << "f2 : " << f2 << std::endl;std::cout << "compute modular filtered gcd(f1,f2): " << std::endl;Poly g1 = modular_filtered_gcd(f1,f2);std::cout << "gcd(f1,f2): " << g1 << std::endl;std::cout << std::endl;Poly p1 = f1*f3;Poly p2 = f2*f3;std::cout << "f3 : " << f3 << std::endl;std::cout << "p1=f1*f3 : " << p1 << std::endl;std::cout << "p2=f2*f3 : " << p2 << std::endl;std::cout << "compute modular filtered gcd(p1,p2): " << std::endl;Poly g2 = modular_filtered_gcd(p1,p2);std::cout << "gcd(p1,p2): " << g2 << std::endl;
}
#else
int main (){std::cout << " This example needs GMP! " << std::endl;
}
#endif
3 历史
Residue 类基于 Sylvain Pion 等人在 [2] 中提出的 C 代码。
软件包的其余部分是将 Exacus [1] 的 NumeriX 库与 CGAL 整合的结果。