00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX__POLYNOMIAL_RING_DICHO__HPP
00014 #define __MMX__POLYNOMIAL_RING_DICHO__HPP
00015 #include <algebramix/polynomial_ring_naive.hpp>
00016 #include <algebramix/polynomial_dicho.hpp>
00017 #include <algebramix/matrix.hpp>
00018
00019 namespace mmx {
00020 #define TMPL template<typename C>
00021 #define TMPLP template <typename Polynomial>
00022
00023
00024
00025
00026
00027 template<typename V>
00028 struct polynomial_ring_dicho_inc: public V {
00029 typedef typename V::Vec Vec;
00030 typedef typename V::Naive Naive;
00031 typedef typename V::Positive Positive;
00032 typedef polynomial_ring_dicho_inc<typename V::No_simd> No_simd;
00033 typedef polynomial_ring_dicho_inc<typename V::No_thread> No_thread;
00034 typedef polynomial_ring_dicho_inc<typename V::No_scaled> No_scaled;
00035 };
00036
00037 template<typename F, typename V, typename W>
00038 struct implementation<F,V,polynomial_ring_dicho_inc<W> >:
00039 public implementation<F,V,W> {};
00040
00041 DEFINE_VARIANT_1 (typename V, V,
00042 polynomial_ring_dicho,
00043 polynomial_ring_dicho_inc<
00044 polynomial_ring_naive<
00045 polynomial_dicho<V> > >)
00046
00047
00048
00049
00050
00051 template<typename V>
00052 struct polynomial_gcd_ring_dicho_inc:
00053 public V {};
00054
00055 template<typename F, typename V, typename W>
00056 struct implementation<F,V,polynomial_gcd_ring_dicho_inc<W> >:
00057 public implementation<F,V,W> {};
00058
00059 DEFINE_VARIANT_1 (typename V, V,
00060 polynomial_gcd_ring_dicho,
00061 polynomial_gcd_ring_dicho_inc<
00062 polynomial_gcd_ring_naive<
00063 polynomial_ring_dicho<V> > >)
00064
00065
00066
00067
00068
00069 template<typename V, typename W>
00070 struct implementation<polynomial_divide,V,polynomial_ring_dicho_inc<W> >:
00071 public implementation<polynomial_multiply,V>
00072 {
00073 typedef polynomial_divide_threshold<polynomial_ring_dicho<W> > Th;
00074 typedef implementation<polynomial_multiply,V> Pol;
00075 typedef implementation<polynomial_divide,polynomial_naive> Fallback;
00076
00077 TMPL static void
00078 quo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
00079
00080 if (n1 < n2);
00081 else if (n1 < Threshold(C,Th))
00082 Fallback::quo_rem (dest, s1, s2, n1, n2);
00083 else {
00084 nat l= n1-n2+1;
00085 if (n2 > l) {
00086 nat d= n2-l;
00087 quo_rem (dest, s1+d, s2+d, n1-d, l);
00088 nat buf_size= aligned_size<C,V> (l+d-1);
00089 C* buf= mmx_new<C> (buf_size);
00090 Pol::mul (buf, dest, s2, l, d);
00091 Pol::sub (s1, buf, l+d-1);
00092 mmx_delete<C> (buf, buf_size);
00093 }
00094 else {
00095 nat h= (n2+1) >> 1;
00096 nat d= l-h;
00097 quo_rem (dest+d, s1+d, s2, n1-d, n2);
00098 quo_rem (dest, s1, s2, n1-h, n2);
00099 }
00100 }
00101 }
00102
00103 TMPL static void
00104 tquo_rem (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
00105
00106
00107
00108
00109
00110 if (n1 < n2)
00111 Pol::copy (dest, s1, n1);
00112 else if (n1 < Threshold(C,Th))
00113 Fallback::tquo_rem (dest, s1, s2, n1, n2);
00114 else {
00115 nat l= n1-n2+1;
00116 if (n2 > l) {
00117 nat d= n2-l;
00118 nat buf_size= aligned_size<C,V> (l+d-1);
00119 C* buf= mmx_new<C> (buf_size);
00120 nat tmp_size= aligned_size<C,V> (n1);
00121 C* tmp= mmx_new<C> (tmp_size);
00122 Pol::neg (buf , s1 , l+d-1);
00123 Pol::tmul (tmp+n2-1, s2 , buf, d, l);
00124 Pol::add (tmp+n2-1, s1+n2-1, l);
00125 Pol::copy (tmp , s1 , n2-1);
00126 Pol::tquo_rem (dest+d , tmp+d , s2+d, n1-d, l);
00127 Pol::copy (dest , s1 , d);
00128 mmx_delete<C> (buf, buf_size);
00129 mmx_delete<C> (tmp, tmp_size);
00130 }
00131 else {
00132 nat h= (n2+1) >> 1;
00133 nat d= l-h;
00134 nat tmp_size= aligned_size<C,V> (n1-d);
00135 C* tmp= mmx_new<C> (tmp_size);
00136 tquo_rem (dest , s1 , s2, n1-h, n2);
00137 Pol::copy (tmp , dest + d, n2-1);
00138 Pol::copy (tmp+n2-1, s1+n1-h , h);
00139 tquo_rem (dest+d , tmp , s2, n1-d, n2);
00140 mmx_delete<C> (tmp, tmp_size);
00141 }
00142 }
00143 }
00144
00145 TMPL static void
00146 pquo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
00147 if (n1 < n2);
00148 else if (n1 < Threshold(C,Th))
00149 Fallback::pquo_rem (dest, s1, s2, n1, n2);
00150 else {
00151 nat l= n1-n2+1;
00152 if (n2 > l) {
00153 nat d= n2-l;
00154 pquo_rem (dest, s1+d, s2+d, n1-d, l);
00155 nat buf_size= aligned_size<C,V> (l+d-1);
00156 C* buf= mmx_new<C> (buf_size);
00157 Pol::mul (buf, dest, s2, l, d);
00158 Pol::mul_sc (s1, binpow (s2[n2-1], l), d);
00159 Pol::sub (s1, buf, l+d-1);
00160 mmx_delete<C> (buf, buf_size);
00161 }
00162 else {
00163 nat h= (n2+1) >> 1;
00164 nat d= l-h;
00165 pquo_rem (dest+d, s1+d, s2, n1-d, n2);
00166 Pol::mul_sc (s1, binpow (s2[n2-1], h), d);
00167 pquo_rem (dest, s1, s2, n1-h, n2);
00168 Pol::mul_sc (dest+d, binpow (s2[n2-1], d), h);
00169 }
00170 }
00171 }
00172
00173 };
00174
00175
00176
00177
00178
00179 template<typename V, typename BV>
00180 struct implementation<polynomial_subresultant_base,V,
00181 polynomial_gcd_ring_dicho_inc<BV> > {
00182
00183 template<typename T> static inline T
00184 defected_pow (const T& x, const T& y, int n) {
00185
00186 VERIFY (n >= 0, "bug");
00187 if (n == 0) return y;
00188 if (n == 1) return x;
00189 T z= square_op::op (defected_pow (x, y, n >> 1)) / y;
00190 return ((n & 1) == 0) ? z : ((x * z) / y);
00191 }
00192
00193 #define C typename scalar_type_helper<Polynomial>::val
00194
00195 TMPLP static void
00196 subresultant_compose (matrix<Polynomial>& num_R,
00197 const matrix<Polynomial>& num_S,
00198 const C& den_S, nat k= 0) {
00199 if (k == 0)
00200 num_R= num_S * num_R;
00201 else {
00202 num_R(0,0)= range (range (num_R(0,0), 0, k) * range (num_S(0,0), 0, k) +
00203 range (num_R(0,1), 0, k) * range (num_S(1,0), 0, k),
00204 0, k);
00205 num_R(0,1)= range (range (num_R(0,0), 0, k) * range (num_S(0,1), 0, k) +
00206 range (num_R(0,1), 0, k) * range (num_S(1,1), 0, k),
00207 0, k);
00208 num_R(1,0)= range (range (num_R(1,0), 0, k) * range (num_S(0,0), 0, k) +
00209 range (num_R(1,1), 0, k) * range (num_S(1,0), 0, k),
00210 0, k);
00211 num_R(1,1)= range (range (num_R(1,0), 0, k) * range (num_S(0,1), 0, k) +
00212 range (num_R(1,1), 0, k) * range (num_S(1,1), 0, k),
00213 0, k);
00214 }
00215 num_R(0,0) /= den_S; num_R(0,1) /= den_S;
00216 num_R(1,0) /= den_S; num_R(1,1) /= den_S; }
00217
00218 TMPLP static void
00219 subresultant_compose (matrix<Polynomial>& num_R,
00220 const vector<matrix<Polynomial> >& num_S,
00221 const vector<C>& den_S, nat k=0) {
00222
00223
00224 nat n= N(num_S);
00225 if (n == 0) return;
00226 if (n == 1) {
00227 subresultant_compose (num_R, num_S[0], den_S[0], k);
00228 return;
00229 }
00230 nat m= n >> 1;
00231 matrix<Polynomial> num_T (Polynomial(C(0)), 2, 2);
00232 C det_num_R= num_R(0,0)[0] * num_R(1,1)[0] - num_R(1,0)[0] * num_R(0,1)[0];
00233 num_T(0,0)= num_T(1,1)= det_num_R;
00234 subresultant_compose (num_T, range (num_S, 0, m), range (den_S, 0, m));
00235 subresultant_compose (num_R, num_T, det_num_R);
00236 det_num_R= num_R(0,0)[0] * num_R(1,1)[0] - num_R(1,0)[0] * num_R(0,1)[0];
00237 num_T(0,0)= num_T(1,1)= det_num_R;
00238 num_T(0,1)= num_T(1,0)= 0;
00239 subresultant_compose (num_T, range (num_S, m, n), range (den_S, m, n));
00240 subresultant_compose (num_R, num_T, det_num_R); }
00241
00242 TMPLP static void
00243 half_subresultant_rec (const Polynomial& r0, const Polynomial& r1,
00244 bool defected, const C& num,
00245 matrix<Polynomial>& num_R, C& den_R, nat& k_out,
00246 vector<matrix<Polynomial> >& Q, vector<C>& rho,
00247 nat k) {
00248
00249
00250
00251 nat n0= N(r0), n1= N(r1), l, l0, l1, h1, k2= (k << 1) - 1;
00252 VERIFY (n0 > n1, "bad input sizes");
00253 VERIFY (k <= n0, "index k out of range");
00254 if (k < n0 - n1 + 1) {
00255 num_R= matrix<Polynomial> (Polynomial (C(0)), 2, 2);
00256 num_R(0,0)= num; num_R(1,1)= num; den_R= num;
00257 k_out= 1;
00258 return;
00259 }
00260 half_subresultant_rec (r0, r1, defected, num,
00261 num_R, den_R, k_out, Q, rho, k >> 1);
00262 VERIFY (k_out <= (k >> 1), "bug");
00263
00264 Polynomial s0, s1, t0, t1, q;
00265 if (k2 >= n0) {
00266 t0= r0; t1= r1;
00267 }
00268 else {
00269 t0= range (r0, n0 - k2, n0); t1= range (r1, n0 - k2, n1);
00270 defected= true;
00271 }
00272 s1= num_R(1,0) * t0 + num_R(1,1) * t1;
00273 if (s1 == 0) { return; }
00274 s0= num_R(0,0) * t0 + num_R(0,1) * t1;
00275 nat m0= N(s0), m1= N(s1);
00276 if (defected) {
00277 s1= range (s1, k_out - 1, m1);
00278 if (s1 == 0) { return; }
00279 s0= range (s0, k_out - 1, m0);
00280 m0 -= k_out - 1; m1 -= k_out - 1;
00281 }
00282 s0 /= den_R; s1 /= den_R;
00283 if (k < k_out + m0 - m1) return;
00284 C c0= s0[m0-1], c1= s1[m1-1], den_S= square_op::op (c0);
00285 C a= defected_pow (c1, c0, m0 - m1 - 1);
00286 t0= (m0 == m1 + 1) ? s1 : ((a * s1) / c0);
00287 C b0= t0[m1-1];
00288 t1= rem (b0 * c1 * s0, s1, q);
00289 l0= N(t0), l1= N(t1), l= N(q);
00290 matrix<Polynomial> num_S (Polynomial (C(0)), 2, 2);
00291 num_S(0,0)= 0 ; num_S(0,1)= a * c0;
00292 if ((m0-m1-1) & 1) { num_S(1,0)= (-b0) * c1; num_S(1,1)= q; }
00293 else { num_S(1,0)= b0 * c1; num_S(1,1)= -q; }
00294 Q << num_S; rho << den_S;
00295 subresultant_compose (num_R, num_S, den_S); k_out += l - 1;
00296 if (defected && l0 >= l-1 && l1 >= l-1) {
00297 t0= range (t0, l - 1, l0); t1= range (t1, l - 1, l1);
00298 l0 -= l - 1; l1 -= l - 1;
00299 }
00300 t1 /= ((m0-m1-1) & 1) ? (-den_S) : den_S;
00301 C det_num_R= num_R(0,0)[0] * num_R(1,1)[0] - num_R(1,0)[0] * num_R(0,1)[0];
00302 if (t1 == 0) { return; }
00303 half_subresultant_rec (t0, t1, defected, det_num_R,
00304 num_S, den_S, h1, Q, rho, k - k_out + 1);
00305 VERIFY (h1 <= k - k_out + 1, "bug");
00306 subresultant_compose (num_R, num_S, den_S); k_out += h1 - 1; }
00307
00308 TMPLP static void
00309 half_subresultant (const Polynomial& f, const Polynomial& g,
00310 matrix<Polynomial>& R,
00311 vector<matrix<Polynomial> >& Q, vector<C>& rho, nat k) {
00312
00313 nat n= N(f), m= N(g), k_out, k_rem;
00314 VERIFY (m > 1 && n >= m, "bad input sizes");
00315 VERIFY (k <= n, "index k out of range");
00316 C d= binpow (g[m-1], n-m);
00317 Polynomial q;
00318 Polynomial s0= g, s1= prem (f, g, q); k_rem = k - N(q) + 1;
00319 nat m0= N(s0), m1= N(s1);
00320 R= matrix<Polynomial> (Polynomial (C(0)), 2, 2);
00321 R(0,0)= 0 ; R(0,1)= 1;
00322 R(1,0)= g[m-1] * d; R(1,1)= - q;
00323 Q= vec (R); rho= vec (C(1));
00324 if (m1 == 0 || k < n - m1 + 1) return;
00325 C c0= s0[m0-1], c1= s1[m1-1];
00326 C a= defected_pow (c1, d, m0-m1-1);
00327 Polynomial t0= (m0 == m1+1) ? s1 : ((a * s1) / d);
00328 C b0= t0[m1-1];
00329 Polynomial t1= rem (b0 * c1 * s0, s1, q); k_rem -= N(q) - 1;
00330 matrix<Polynomial> num_S (Polynomial (C(0)), 2, 2); C den_S (d * c0);
00331 t1 /= ((m0-m1-1) & 1) ? -den_S : den_S;
00332 num_S(0,0)= 0 ; num_S(0,1)= a * c0;
00333 if ((m0-m1-1) & 1) { num_S(1,0)= -(b0 * c1); num_S(1,1)= q; }
00334 else { num_S(1,0)= b0 * c1 ; num_S(1,1)= -q; }
00335 Q << num_S; rho << den_S;
00336 subresultant_compose (R, num_S, den_S);
00337 if (N(t1) == 0 || k < n - N(t1) + 1) return;
00338 C det_R= R(0,0)[0] * R(1,1)[0] - R(1,0)[0] * R(0,1)[0];
00339 half_subresultant_rec (t0, t1, false, det_R,
00340 num_S, den_S, k_out, Q, rho, k_rem);
00341 subresultant_compose (R, num_S, den_S); }
00342
00343 TMPLP static Polynomial
00344 _half_gcd (const Polynomial& f, const Polynomial& g,
00345 const Polynomial& first_prem) {
00346
00347
00348
00349 nat n= N(f), m= N(g);
00350 VERIFY (m > 1 && n >= m, "bad input sizes");
00351
00352 C d= binpow (g[m-1], n-m);
00353 Polynomial q, s0= g, s1= first_prem;
00354 nat m0= N(s0), m1= N(s1);
00355 if (m1 == 0) return g;
00356 if (m1 == 1) return s1;
00357
00358 C c0= s0[m0-1], c1= s1[m1-1];
00359 C a= defected_pow (c1, d, m0-m1-1);
00360 Polynomial t0= (m0 == m1+1) ? s1 : ((a * s1) / d);
00361 C b0= t0[m1-1];
00362 Polynomial t1= rem (b0 * c1 * s0, s1, q);
00363 if (N(t1) == 0) return s1;
00364 if (N(t1) == 1) return t1;
00365 matrix<Polynomial> num_S (Polynomial (C(0)), 2, 2); C den_S (d * c0);
00366 t1 /= ((m0-m1-1) & 1) ? -den_S : den_S;
00367 num_S(0,0)= 0 ; num_S(0,1)= a * c0;
00368 if ((m0-m1-1) & 1) { num_S(1,0)= -(b0 * c1); num_S(1,1)= q; }
00369 else { num_S(1,0)= b0 * c1 ; num_S(1,1)= -q; }
00370
00371 C den= (- g[m-1] * d * (num_S(0,0)[0] * num_S(1,1)[0]
00372 - num_S(1,0)[0] * num_S(0,1)[0])) / den_S;
00373 vector<matrix<Polynomial> > Q; vector<C> rho;
00374 matrix<Polynomial> num_T (Polynomial (C(0)), 2, 2); C den_T; nat k_out;
00375 half_subresultant_rec (t0, t1, false, den,
00376 num_T, den_T, k_out, Q, rho, N(t1));
00377 Polynomial h;
00378 nat k= N(t1) - degree (num_T(1,0));
00379 h= range (range (num_T(1,0), 0, k) * range (t0, 0, k) +
00380 range (num_T(1,1), 0, k) * range (t1, 0, k), 0, k);
00381 if (h != 0) return (den * h) / den_T;
00382 k= N(t1) - degree (num_T(0,0));
00383 h= range (range (num_T(0,0), 0, k) * range (t0, 0, k) +
00384 range (num_T(0,1), 0, k) * range (t1, 0, k), 0, k);
00385 VERIFY (h != 0, "bug");
00386 return (den * h) / den_T; }
00387
00388 public:
00389
00390 TMPLP static void
00391 subresultant_sequence (const Polynomial& f, const Polynomial& g,
00392 vector<Polynomial>& res,
00393 vector<Polynomial>& cof, vector<Polynomial>& cog,
00394 Polynomial& d1, Polynomial& u1, Polynomial& v1,
00395 Polynomial& d2, Polynomial& u2, Polynomial& v2,
00396 nat dest_deg) {
00397
00398 VERIFY (degree (f) >= degree (g) && degree (g) > 0, "incorrect degrees");
00399 nat n= degree (f), m= degree (g), k;
00400 nat min_deg_res= (d1 != 0 || d2 != 0) ? dest_deg : m;
00401 for (k= min_deg_res; k+1 > 0; k--)
00402 if (k < N(res) && res[k] != 0) min_deg_res= k;
00403 nat min_deg_cof= (u1 != 0 || u2 != 0) ? dest_deg : m;
00404 for (k= min_deg_cof; k+1 > 0; k--)
00405 if (k < N(cof) && cof[k] != 0) min_deg_cof= k;
00406 nat min_deg_cog= (v1 != 0 || v2 != 0) ? dest_deg : m;
00407 for (k= min_deg_cog; k+1 > 0; k--)
00408 if (k < N(cog) && cog[k] != 0) min_deg_cog= k;
00409 nat min_deg_co= min (min_deg_cof, min_deg_cog);
00410 bool opt_res= res != 0 || d1 != 0 || d2 != 0;
00411 bool opt_cof= cof != 0 || u1 != 0 || u2 != 0;
00412 bool opt_cog= cog != 0 || v1 != 0 || v2 != 0;
00413 if (!(opt_res || opt_cof || opt_cog)) return;
00414 vector<matrix<Polynomial> > Q; vector<C> rho; matrix<Polynomial> R;
00415 half_subresultant (f, g, R, Q, rho, n + 1 - dest_deg);
00416 if (u1 != 0) u1= R(0,0); if (v1 != 0) v1= R(0,1);
00417 if (u2 != 0) u2= R(1,0); if (v2 != 0) v2= R(1,1);
00418 if (d1 != 0) {
00419 k= N(g) - degree (R(0,0));
00420 d1= range (range (R(0,0), 0, k) * range (f, 0, k) +
00421 range (R(0,1), 0, k) * range (g, 0, k), 0, k);
00422 }
00423 if (d2 != 0) {
00424 k= N(g) - degree (R(1,0));
00425 d2= range (range (R(1,0), 0, k) * range (f, 0, k) +
00426 range (R(1,1), 0, k) * range (g, 0, k), 0, k);
00427 }
00428 R= identity_matrix<Polynomial> (2);
00429 vector<Polynomial> r= vec<Polynomial> (f, g);
00430 vector<nat> nn= vec<nat> (N(g));
00431
00432 for (nat i= 1; i < N(Q); i++) nn << (nn[i-1] - degree (Q[i](1,1)));
00433 nn << 0; nat l0= 0, l1; k= m;
00434 while (k > 0) {
00435 while ((k-1 >= N(res) || res[k-1] == 0) &&
00436 (k-1 >= N(cof) || cof[k-1] == 0) &&
00437 (k-1 >= N(cog) || cog[k-1] == 0) && k > 0) k--;
00438 if (k == 0) break;
00439 l1= find (nn, k);
00440
00441 if (l1 >= N(nn)) {
00442 l1= find (nn, k+1);
00443
00444
00445 if (l1 >= N(nn)) {
00446 if (k-1 < N(res)) res[k-1] = 0;
00447 if (k-1 < N(cof)) cof[k-1] = 0;
00448 if (k-1 < N(cog)) cog[k-1] = 0;
00449 k--; continue;
00450 }
00451 k++;
00452 }
00453 subresultant_compose (R, range (Q, l0, l1+1), range (rho, l0, l1+1),
00454 k-1 >= min_deg_co ? 0 : k);
00455 if (k-1 < N(res) && res[k-1] != 0)
00456 res[k-1]= range (range (R(0,0), 0, k) * range (f, 0, k) +
00457 range (R(0,1), 0, k) * range (g, 0, k), 0, k);
00458 if (k-1 < N(cof) && cof[k-1] != 0) cof[k-1]= R(0,0);
00459 if (k-1 < N(cog) && cog[k-1] != 0) cog[k-1]= R(0,1);
00460 k--; if (k == 0) break;
00461 if (k-1 < N(res) && res[k-1] != 0)
00462 res[k-1]= range (range (R(1,0), 0, k) * range (f, 0, k) +
00463 range (R(1,1), 0, k) * range (g, 0, k), 0, k);
00464 if (k-1 < N(cof) && cof[k-1] != 0) cof[k-1]= R(1,0);
00465 if (k-1 < N(cog) && cog[k-1] != 0) cog[k-1]= R(1,1);
00466 k--; l0= l1+1;
00467 }
00468 }
00469
00470 #undef C
00471 };
00472
00473
00474
00475
00476
00477 template<typename V, typename W>
00478 struct implementation<polynomial_evaluate,V,polynomial_gcd_ring_dicho_inc<W> >:
00479 public implementation<polynomial_evaluate,V,W>
00480 {
00481 typedef implementation<polynomial_subresultant_base,V,
00482 polynomial_gcd_ring_dicho_inc<W> > Pol;
00483
00484 #define C typename scalar_type_helper<Polynomial>::val
00485 private:
00486 TMPLP static Polynomial
00487 gcd_with_first_prem (const Polynomial& P1, const Polynomial& P2,
00488 const Polynomial& R) {
00489
00490
00491
00492 nat n1= N(P1), n2= N(P2);
00493 if (n1 < n2) return gcd_op::op (P2, P1);
00494 if (n2 == 0) { return P1; }
00495 if (n1 == 0) { return P2; }
00496 if (R == 0) {
00497 C c2; Polynomial T= primitive_part (P2, c2);
00498 return T * gcd_op::op (contents (P1), c2);
00499 }
00500 Polynomial G= Pol::_half_gcd (P1, P2, R);
00501 return primitive_part (G) * gcd_op::op (contents (P1), contents (P2));
00502 }
00503 #undef C
00504
00505 public:
00506 TMPLP static vector<Polynomial>
00507 multi_gcd (const Polynomial& P, const vector<Polynomial>& Q) {
00508 vector<Polynomial> R= prem (P, Q);
00509 for (nat i = 0; i < N(Q); i++)
00510 R[i]= gcd_with_first_prem (P, Q[i], R[i]);
00511 return R;
00512 }
00513
00514 };
00515
00516 #undef TMPL
00517 #undef TMPLP
00518 }
00519 #endif //__MMX__POLYNOMIAL_RING_DICHO__HPP
