1 导言
模块化运算是现代代数系统的基本工具。结合中国余数定理,它是计算 gcd、结果等算法的主力。此外,它还可以作为一种非常有效的过滤器,因为通常只需计算一个素数的模数对应值,就可以排除某个值为零的可能性。
2 留数和模板化
先,本软件包引入了一个 Residue 类型。质数 p 保存在一个静态成员变量中。该类提供静态成员函数来更改该值。
改变质数会使已存在的该类型对象失效。但是,已经存在的对象不会失去其相对于旧质数的价值,并且在恢复旧质数后可以重新使用。由于该类型是基于双运算的,因此质数的值只能小于 226。p 的初始值为 67108859。
此外,软件包还引入了可模块化概念。如果存在从 T 到基于残差类型的代数结构的映射,那么代数结构 T 就被认为是可模块化的。对于标量类型(如整数),这种映射只是进入由 Residue 表示的 Z/pZ 的典型同态映射。对于复合类型(如多项式),该映射适用于复合类型的系数。该映射由类 Modular_traits
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#include #ifdef CGAL_USE_GMP #include #include // 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 arithmetic CGAL::Protect_FPU_rounding pfr(CGAL_FE_TONEAREST); // Use Modular_traits to convert to polynomials with modular coefficients typedef CGAL::Modular_traits 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 degree typename CGAL::Polynomial_traits_d ::Degree degree; typename CGAL::Polynomial_traits_d ::Degree mdegree; if ( degree(p1) != mdegree(mp1)) return true; if ( degree(p2) != mdegree(mp2)) return true; // compute gcd for modular images MPolynomial mg = CGAL::gcd(mp1,mp2); // if the modular gcd is not trivial: return true if ( 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 MT; typedef typename MT::Is_modularizable Is_modularizable; // Try to avoid actual gcd computation if (may_have_common_factor(p1,p2, Is_modularizable())){ // Compute gcd, since the filter indicates a common factor return CGAL::gcd(p1,p2); }else{ typename CGAL::Polynomial_traits_d ::Univariate_content content; typename CGAL::Polynomial_traits_d ::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 Poly; CGAL::Polynomial_traits_d ::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 整合的结果。