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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
|
// Copyright (c) 2002-2008 Max-Planck-Institute Saarbruecken (Germany)
//
// This file is part of CGAL (www.cgal.org)
//
// $URL: https://github.com/CGAL/cgal/blob/v6.1.1/Polynomial/include/CGAL/Polynomial/modular_gcd_utcf_algorithm_M.h $
// $Id: include/CGAL/Polynomial/modular_gcd_utcf_algorithm_M.h 08b27d3db14 $
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Dominik Huelse <dominik.huelse@gmx.de>
// Michael Hemmer <mhemmer@uni-mainz.de>
//
// ============================================================================
/*! \file CGAL/Polynomial/modular_gcd_utcf_algorithm_M.h
provides gcd for Polynomials, based on Modular arithmetic.
*/
#ifndef CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_ALGORITHM_M_H
#define CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_ALGORITHM_M_H 1
#include <CGAL/basic.h>
#include <CGAL/Residue.h>
#include <CGAL/Polynomial/modular_gcd.h>
#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial/Cached_extended_euclidean_algorithm.h>
#include <CGAL/Scalar_factor_traits.h>
#include <CGAL/Chinese_remainder_traits.h>
#include <CGAL/Cache.h>
#include <CGAL/Real_timer.h>
// algorithm M for integer polynomials, without denominator bound
namespace CGAL {
namespace internal{
template <class NT> Polynomial<NT> gcd_utcf_UFD(Polynomial<NT>,Polynomial<NT>);
template <class NT>
Polynomial< Polynomial<NT> > modular_gcd_utcf_algorithm_M(
const Polynomial< Polynomial<NT> >& FF1 ,
const Polynomial< Polynomial<NT> >& FF2 ){
return gcd_utcf_UFD(FF1, FF2);
}
template <class NT>
Polynomial<NT> modular_gcd_utcf_algorithm_M(
const Polynomial<NT>& FF1 ,
const Polynomial<NT>& FF2 ){
// Enforce IEEE double precision and to nearest before using modular arithmetic
CGAL::Protect_FPU_rounding<true> pfr(CGAL_FE_TONEAREST);
// std::cout << "start modular_gcd_utcf_algorithm_M " << std::endl;
#ifdef CGAL_MODULAR_GCD_TIMER
timer_init.start();
#endif
typedef Polynomial<NT> Poly;
// will play the role of content
typedef typename CGAL::Scalar_factor_traits<Poly>::Scalar Scalar;
typedef typename CGAL::Modular_traits<Poly>::Residue_type MPoly;
typedef typename CGAL::Modular_traits<Scalar>::Residue_type MScalar;
typedef Chinese_remainder_traits<Poly> CRT;
typename CRT::Chinese_remainder chinese_remainder;
CGAL::Real_timer timer;
if(FF1.is_zero()){
if(FF2.is_zero()){
return Poly(1);// TODO: return 0 for CGAL
}
else{
// std::cout<<"\nFF1 is zero"<<std::endl;
return CGAL::canonicalize(FF2);
}
}
if(FF2.is_zero()){
return CGAL::canonicalize(FF1);
}
if(FF1.degree() == 0 || FF2.degree() == 0){
Poly result;
result = Poly(CGAL::gcd(FF1.content(),FF2.content()));
return CGAL::canonicalize(result);
}
Poly F1 = CGAL::canonicalize(FF1);
Poly F2 = CGAL::canonicalize(FF2);
Scalar f1 = scalar_factor(F1.lcoeff()); // ilcoeff(F1)
Scalar f2 = scalar_factor(F2.lcoeff()); // ilcoeff(F2)
Scalar g_ = scalar_factor(f1,f2);
//std::cout <<" g_ : "<< g_ << std::endl;
bool solved = false;
int prime_index = -1;
int n = 0; // number of lucky primes
int degree_F1 = F1.degree();
int degree_F2 = F2.degree();
int degree_e = (std::min)(degree_F1,degree_F2);
MScalar mg_;
MPoly mF1,mF2,mG_;
typename CRT::Scalar_type p(0),q(0),pq,s,t;
Poly Gs,H1s,H2s, Gs_old; // s =^ star
#ifdef CGAL_MODULAR_GCD_TIMER
timer_init.stop();
#endif
while(!solved){
do{
//---------------------------------------
//choose prime not dividing f1 or f2
MScalar tmp1, tmp2;
do{
int current_prime = -1;
prime_index++;
if(prime_index >= 2000){
std::cerr<<"primes in the array exhausted"<<std::endl;
current_prime = internal::get_next_lower_prime(current_prime);
}
else{
current_prime = internal::primes[prime_index];
}
CGAL::Residue::set_current_prime(current_prime);
#ifdef CGAL_MODULAR_GCD_TIMER
timer_image.start();
#endif
tmp1 = CGAL::modular_image(f1);
tmp2 = CGAL::modular_image(f2);
#ifdef CGAL_MODULAR_GCD_TIMER
timer_image.stop();
#endif
}
while(!(( tmp1 != 0 ) && ( tmp2 != 0 )));
// --------------------------------------
// invoke gcd for current prime
#ifdef CGAL_MODULAR_GCD_TIMER
timer_image.start();
#endif
mg_ = CGAL::modular_image(g_);
mF1 = CGAL::modular_image(F1);
mF2 = CGAL::modular_image(F2);
#ifdef CGAL_MODULAR_GCD_TIMER
timer_image.stop();
// replace mG_ = CGAL::gcd (mF1,mF2)*MPoly(mg_); for multivariat
timer_gcd.start();
#endif
mG_ = CGAL::gcd(mF1,mF2)*MPoly(mg_);
#ifdef CGAL_MODULAR_GCD_TIMER
timer_gcd.stop();
#endif
//mH1 = CGAL::integral_div(mF1,mG_);
//mH2 = CGAL::integral_div(mF2,mG_);
//---------------------------------------
// return if G is constant
if (mG_ == MPoly(1)) return Poly(1);
// --------------------------------------
}// repeat until mG_ degree is less equal the known bound
// check prime
while( mG_.degree() > degree_e);
if( mG_.degree() < degree_e ){
if( n != 0 ) std::cout << "UNLUCKY PRIME !!"<< std::endl;
// restart chinese remainder
// ignore previous unlucky primes
n=1;
degree_e= mG_.degree();
}else{
CGAL_postcondition( mG_.degree() == degree_e);
n++; // increase number of lucky primes
}
// --------------------------------------
// try chinese remainder
// std::cout <<" chinese remainder round :" << n << std::endl;
typename CGAL::Modular_traits<Poly>::Modular_image_representative inv_map;
if(n == 1){
// init chinese remainder
q = CGAL::Residue::get_current_prime(); // implicit !
Gs_old = Gs = inv_map(mG_);
//H1s_old = H1s = inv_map(mH1);
//H2s_old = H2s = inv_map(mH2);
}else{
// continue chinese remainder
p = CGAL::Residue::get_current_prime(); // implicit!
Gs_old = Gs ;
//H1s_old = H1s ;
//H2s_old = H2s ;
#ifdef CGAL_MODULAR_GCD_TIMER
timer_CR.start();
#endif
// chinese_remainder(q,Gs ,p,inv_map(mG_),pq,Gs);
// cached_extended_euclidean_algorithm(q,p,s,t);
internal::Cached_extended_euclidean_algorithm
<typename CRT::Scalar_type, 1> ceea;
ceea(q,p,s,t);
pq =p*q;
chinese_remainder(q,p,pq,s,t,Gs,inv_map(mG_),Gs);
#ifdef CGAL_MODULAR_GCD_TIMER
timer_CR.stop();
#endif
q=pq;
}
try{
if( n != 1 && Gs == Gs_old ){
Poly r1,r2;
#ifdef CGAL_MODULAR_GCD_TIMER
timer_division.start();
#endif
typedef CGAL::Algebraic_structure_traits< Poly > ASTE_Poly;
typename ASTE_Poly::Divides divides;
bool div1=divides(Gs,g_*F1,H1s);
bool div2=divides(Gs,g_*F2,H2s);
if (div1 && div2){
solved = true;
}
// this is the old code
// NT dummy;
// Poly::euclidean_division(g_*F1,Gs,H1s,r1);
// Poly::euclidean_division(g_*F2,Gs,H2s,r2);
// if (r1.is_zero() && r2.is_zero())
// solved = true;
#ifdef CGAL_MODULAR_GCD_TIMER
timer_division.stop();
#endif
// std::cout << "number of primes used : "<< n << std::endl;
} // end while
}catch(...){}
}
//TODO CGAL: change this to multivariat content
// Scalar scalar_content_f1 = scalar_factor(FF1);
// Scalar scalar_content_f2 = scalar_factor(FF2);
// Scalar scalar_content_gcd = CGAL::gcd(scalar_content_f1,scalar_content_f2);
// Poly result = CGAL::canonicalize(Gs)*Poly(scalar_content_gcd);
// return result;
return CGAL::canonicalize(Gs);
}
} // namespace internal
} // namespace CGAL
#endif // CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_ALGORITHM_M_H
|
