1 #include "kernel/mod2.h"
2
3 #include "misc/options.h"
4
5 #include "polys.h"
6 #include "kernel/ideals.h"
7 #include "kernel/ideals.h"
8 #include "polys/clapsing.h"
9 #include "polys/clapconv.h"
10
11 /// Widely used global variable which specifies the current polynomial ring for Singular interpreter and legacy implementatins.
12 /// @Note: one should avoid using it in newer designs, for example due to possible problems in parallelization with threads.
13 VAR ring currRing = NULL;
14
rChangeCurrRing(ring r)15 void rChangeCurrRing(ring r)
16 {
17 //------------ set global ring vars --------------------------------
18 currRing = r;
19 if( r != NULL )
20 {
21 rTest(r);
22 //------------ global variables related to coefficients ------------
23 assume( r->cf!= NULL );
24 nSetChar(r->cf);
25 //------------ global variables related to polys
26 p_SetGlobals(r); // also setting TEST_RINGDEP_OPTS
27 //------------ global variables related to factory -----------------
28 }
29 }
30
p_Divide(poly p,poly q,const ring r)31 poly p_Divide(poly p, poly q, const ring r)
32 {
33 assume(q!=NULL);
34 if (q==NULL)
35 {
36 WerrorS("div. by 0");
37 return NULL;
38 }
39 if (p==NULL)
40 {
41 p_Delete(&q,r);
42 return NULL;
43 }
44 if ((pNext(q)!=NULL)||rIsPluralRing(r))
45 { /* This means that q != 0 consists of at least two terms*/
46 if(p_GetComp(p,r)==0)
47 {
48 if((rFieldType(r)==n_transExt)
49 &&(convSingTrP(p,r))
50 &&(convSingTrP(q,r))
51 &&(!rIsNCRing(r)))
52 {
53 poly res=singclap_pdivide(p, q, r);
54 p_Delete(&p,r);
55 p_Delete(&q,r);
56 return res;
57 }
58 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
59 &&(!rField_is_Ring(r))
60 &&(!rIsNCRing(r)))
61 {
62 poly res=singclap_pdivide(p, q, r);
63 p_Delete(&p,r);
64 p_Delete(&q,r);
65 return res;
66 }
67 else
68 {
69 ideal vi=idInit(1,1); vi->m[0]=q;
70 ideal ui=idInit(1,1); ui->m[0]=p;
71 ideal R; matrix U;
72 ring save_ring=currRing;
73 if (r!=currRing) rChangeCurrRing(r);
74 int save_opt;
75 SI_SAVE_OPT1(save_opt);
76 si_opt_1 &= ~(Sy_bit(OPT_PROT));
77 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
78 SI_RESTORE_OPT1(save_opt);
79 if (r!=save_ring) rChangeCurrRing(save_ring);
80 matrix T = id_Module2formatedMatrix(m,1,1,r);
81 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
82 id_Delete((ideal *)&T,r);
83 id_Delete((ideal *)&U,r);
84 id_Delete(&R,r);
85 //vi->m[0]=NULL; ui->m[0]=NULL;
86 id_Delete(&vi,r);
87 id_Delete(&ui,r);
88 return p;
89 }
90 }
91 else
92 {
93 int comps=p_MaxComp(p,r);
94 ideal I=idInit(comps,1);
95 poly h;
96 int i;
97 // conversion to a list of polys:
98 while (p!=NULL)
99 {
100 i=p_GetComp(p,r)-1;
101 h=pNext(p);
102 pNext(p)=NULL;
103 p_SetComp(p,0,r);
104 I->m[i]=p_Add_q(I->m[i],p,r);
105 p=h;
106 }
107 // division and conversion to vector:
108 h=NULL;
109 p=NULL;
110 for(i=comps-1;i>=0;i--)
111 {
112 if (I->m[i]!=NULL)
113 {
114 if((rFieldType(r)==n_transExt)
115 &&(convSingTrP(I->m[i],r))
116 &&(convSingTrP(q,r))
117 &&(!rIsNCRing(r)))
118 {
119 h=singclap_pdivide(I->m[i],q,r);
120 }
121 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
122 &&(!rField_is_Ring(r))
123 &&(!rIsNCRing(r)))
124 h=singclap_pdivide(I->m[i],q,r);
125 else
126 {
127 ideal vi=idInit(1,1); vi->m[0]=q;
128 ideal ui=idInit(1,1); ui->m[0]=I->m[i];
129 ideal R; matrix U;
130 ring save_ring=currRing;
131 if (r!=currRing) rChangeCurrRing(r);
132 int save_opt;
133 SI_SAVE_OPT1(save_opt);
134 si_opt_1 &= ~(Sy_bit(OPT_PROT));
135 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
136 SI_RESTORE_OPT1(save_opt);
137 if (r!=save_ring) rChangeCurrRing(save_ring);
138 if (idIs0(R))
139 {
140 matrix T = id_Module2formatedMatrix(m,1,1,r);
141 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
142 id_Delete((ideal *)&T,r);
143 }
144 else p=NULL;
145 id_Delete((ideal*)&U,r);
146 id_Delete(&R,r);
147 vi->m[0]=NULL; ui->m[0]=NULL;
148 id_Delete(&vi,r);
149 id_Delete(&ui,r);
150 }
151 p_SetCompP(h,i+1,r);
152 p=p_Add_q(p,h,r);
153 }
154 }
155 id_Delete(&I,r);
156 p_Delete(&q,r);
157 return p;
158 }
159 }
160 else
161 { /* This means that q != 0 consists of just one term, or LetterPlace */
162 #ifdef HAVE_RINGS
163 if (pNext(q)!=NULL)
164 {
165 WerrorS("division over a coefficient domain only implemented for terms");
166 return NULL;
167 }
168 #endif
169 return p_DivideM(p,q,r);
170 }
171 return NULL;
172 }
173
pp_Divide(poly p,poly q,const ring r)174 poly pp_Divide(poly p, poly q, const ring r)
175 {
176 assume(q!=NULL);
177 if (q==NULL)
178 {
179 WerrorS("div. by 0");
180 return NULL;
181 }
182 if (p==NULL)
183 {
184 return NULL;
185 }
186 if ((pNext(q)!=NULL)||rIsPluralRing(r))
187 { /* This means that q != 0 consists of at least two terms*/
188 if(p_GetComp(p,r)==0)
189 {
190 if((rFieldType(r)==n_transExt)
191 &&(convSingTrP(p,r))
192 &&(convSingTrP(q,r))
193 &&(!rIsNCRing(r)))
194 {
195 poly res=singclap_pdivide(p, q, r);
196 return res;
197 }
198 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
199 &&(!rField_is_Ring(r))
200 &&(!rIsNCRing(r)))
201 {
202 poly res=singclap_pdivide(p, q, r);
203 return res;
204 }
205 else
206 {
207 ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
208 ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
209 ideal R; matrix U;
210 ring save_ring=currRing;
211 if (r!=currRing) rChangeCurrRing(r);
212 int save_opt;
213 SI_SAVE_OPT1(save_opt);
214 si_opt_1 &= ~(Sy_bit(OPT_PROT));
215 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
216 SI_RESTORE_OPT1(save_opt);
217 if (r!=save_ring) rChangeCurrRing(save_ring);
218 matrix T = id_Module2formatedMatrix(m,1,1,r);
219 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
220 id_Delete((ideal *)&T,r);
221 id_Delete((ideal *)&U,r);
222 id_Delete(&R,r);
223 //vi->m[0]=NULL; ui->m[0]=NULL;
224 id_Delete(&vi,r);
225 id_Delete(&ui,r);
226 return p;
227 }
228 }
229 else
230 {
231 p=p_Copy(p,r);
232 int comps=p_MaxComp(p,r);
233 ideal I=idInit(comps,1);
234 poly h;
235 int i;
236 // conversion to a list of polys:
237 while (p!=NULL)
238 {
239 i=p_GetComp(p,r)-1;
240 h=pNext(p);
241 pNext(p)=NULL;
242 p_SetComp(p,0,r);
243 I->m[i]=p_Add_q(I->m[i],p,r);
244 p=h;
245 }
246 // division and conversion to vector:
247 h=NULL;
248 p=NULL;
249 q=p_Copy(q,r);
250 for(i=comps-1;i>=0;i--)
251 {
252 if (I->m[i]!=NULL)
253 {
254 if((rFieldType(r)==n_transExt)
255 &&(convSingTrP(I->m[i],r))
256 &&(convSingTrP(q,r))
257 &&(!rIsNCRing(r)))
258 {
259 h=singclap_pdivide(I->m[i],q,r);
260 }
261 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
262 &&(!rField_is_Ring(r))
263 &&(!rIsNCRing(r)))
264 h=singclap_pdivide(I->m[i],q,r);
265 else
266 {
267 ideal vi=idInit(1,1); vi->m[0]=q;
268 ideal ui=idInit(1,1); ui->m[0]=I->m[i];
269 ideal R; matrix U;
270 ring save_ring=currRing;
271 if (r!=currRing) rChangeCurrRing(r);
272 int save_opt;
273 SI_SAVE_OPT1(save_opt);
274 si_opt_1 &= ~(Sy_bit(OPT_PROT));
275 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
276 SI_RESTORE_OPT1(save_opt);
277 if (r!=save_ring) rChangeCurrRing(save_ring);
278 if (idIs0(R))
279 {
280 matrix T = id_Module2formatedMatrix(m,1,1,r);
281 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
282 id_Delete((ideal *)&T,r);
283 }
284 else p=NULL;
285 id_Delete((ideal*)&U,r);
286 id_Delete(&R,r);
287 vi->m[0]=NULL; ui->m[0]=NULL;
288 id_Delete(&vi,r);
289 id_Delete(&ui,r);
290 }
291 p_SetCompP(h,i+1,r);
292 p=p_Add_q(p,h,r);
293 }
294 }
295 id_Delete(&I,r);
296 p_Delete(&q,r);
297 return p;
298 }
299 }
300 else
301 { /* This means that q != 0 consists of just one term,
302 or that r is over a coefficient ring. */
303 #ifdef HAVE_RINGS
304 if (pNext(q)!=NULL)
305 {
306 WerrorS("division over a coefficient domain only implemented for terms");
307 return NULL;
308 }
309 #endif
310 return pp_DivideM(p,q,r);
311 }
312 return NULL;
313 }
314
p_DivRem(poly p,poly q,poly & rest,const ring r)315 poly p_DivRem(poly p, poly q, poly &rest, const ring r)
316 {
317 assume(q!=NULL);
318 rest=NULL;
319 if (q==NULL)
320 {
321 WerrorS("div. by 0");
322 return NULL;
323 }
324 if (p==NULL)
325 {
326 p_Delete(&q,r);
327 return NULL;
328 }
329 if(p_GetComp(p,r)==0)
330 {
331 if((rFieldType(r)==n_transExt)
332 &&(convSingTrP(p,r))
333 &&(convSingTrP(q,r))
334 &&(!rIsNCRing(r)))
335 {
336 poly res=singclap_pdivide(p, q, r);
337 rest=singclap_pmod(p,q,r);
338 p_Delete(&p,r);
339 p_Delete(&q,r);
340 return res;
341 }
342 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
343 &&(!rField_is_Ring(r))
344 &&(!rIsNCRing(r)))
345 {
346 poly res=singclap_pdivide(p, q, r);
347 rest=singclap_pmod(p,q,r);
348 p_Delete(&p,r);
349 p_Delete(&q,r);
350 return res;
351 }
352 else
353 {
354 ideal vi=idInit(1,1); vi->m[0]=q;
355 ideal ui=idInit(1,1); ui->m[0]=p;
356 ideal R; matrix U;
357 ring save_ring=currRing;
358 if (r!=currRing) rChangeCurrRing(r);
359 int save_opt;
360 SI_SAVE_OPT1(save_opt);
361 si_opt_1 &= ~(Sy_bit(OPT_PROT));
362 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
363 SI_RESTORE_OPT1(save_opt);
364 if (r!=save_ring) rChangeCurrRing(save_ring);
365 matrix T = id_Module2formatedMatrix(m,1,1,r);
366 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
367 id_Delete((ideal *)&T,r);
368 T = id_Module2formatedMatrix(R,1,1,r);
369 rest=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
370 id_Delete((ideal *)&T,r);
371 id_Delete((ideal *)&U,r);
372 id_Delete(&R,r);
373 //vi->m[0]=NULL; ui->m[0]=NULL;
374 id_Delete(&vi,r);
375 id_Delete(&ui,r);
376 return p;
377 }
378 }
379 return NULL;
380 }
381
singclap_gcd(poly f,poly g,const ring r)382 poly singclap_gcd ( poly f, poly g, const ring r )
383 {
384 poly res=NULL;
385
386 if (f!=NULL)
387 {
388 //if (r->cf->has_simple_Inverse) p_Norm(f,r);
389 if (rField_is_Zp(r)) p_Norm(f,r);
390 else p_Cleardenom(f, r);
391 }
392 if (g!=NULL)
393 {
394 //if (r->cf->has_simple_Inverse) p_Norm(g,r);
395 if (rField_is_Zp(r)) p_Norm(g,r);
396 else p_Cleardenom(g, r);
397 }
398 else return f; // g==0 => gcd=f (but do a p_Cleardenom/pNorm)
399 if (f==NULL) return g; // f==0 => gcd=g (but do a p_Cleardenom/pNorm)
400 if(!rField_is_Ring(r)
401 && (p_IsConstant(f,r)
402 ||p_IsConstant(g,r)))
403 {
404 res=p_One(r);
405 }
406 else if (r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
407 {
408 res=singclap_gcd_r(f,g,r);
409 }
410 else
411 {
412 ideal I=idInit(2,1);
413 I->m[0]=f;
414 I->m[1]=p_Copy(g,r);
415 intvec *w=NULL;
416 ring save_ring=currRing;
417 if (r!=currRing) rChangeCurrRing(r);
418 int save_opt;
419 SI_SAVE_OPT1(save_opt);
420 si_opt_1 &= ~(Sy_bit(OPT_PROT));
421 ideal S1=idSyzygies(I,testHomog,&w);
422 if (w!=NULL) delete w;
423 // expect S1->m[0]=(-g/gcd,f/gcd)
424 if (IDELEMS(S1)!=1) WarnS("error in syzygy computation for GCD");
425 int lp;
426 p_TakeOutComp(&S1->m[0],1,&res,&lp,r);
427 p_Delete(&S1->m[0],r);
428 // GCD is g divided iby (-g/gcd):
429 res=p_Divide(g,res,r);
430 // restore, r, opt:
431 SI_RESTORE_OPT1(save_opt);
432 if (r!=save_ring) rChangeCurrRing(save_ring);
433 // clean the result
434 res=p_Cleardenom(res,r);
435 p_Content(res,r);
436 return res;
437 }
438 p_Delete(&f, r);
439 p_Delete(&g, r);
440 return res;
441 }
442