CGAL 5.6.1 - Modular Arithmetic（模板化算法）

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 arithmetic
  CGAL::Protect_FPU_rounding<true> pfr(CGAL_FE_TONEAREST);
  // Use Modular_traits to convert to polynomials with modular coefficients
  typedef 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 degree
  typename 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 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<Polynomial> 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<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 整合的结果。