1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#include <CGAL/basic.h>
#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 useing 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::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 examples needs GMP! " << std::endl;
}
#endif
|
