My Project
Loading...
Searching...
No Matches
polys.cc
Go to the documentation of this file.
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 implementations.
12/// @Note: one should avoid using it in newer designs, for example due to possible problems in parallelization with threads.
13VAR ring currRing = NULL;
14VAR coeffs coeffs_BIGINT;
15
16void rChangeCurrRing(ring r)
17{
18 if (currRing!=NULL)
19 currRing->options = si_opt_1 & TEST_RINGDEP_OPTS;
20 //------------ set global ring vars --------------------------------
21 currRing = r;
22 if( r != NULL )
23 {
24 rTest(r);
25 //------------ global variables related to coefficients ------------
26 assume( r->cf!= NULL );
27 nSetChar(r->cf);
28 //------------ global variables related to polys
29 p_SetGlobals(r); // also setting TEST_RINGDEP_OPTS
30 //------------ global variables related to factory -----------------
31 }
32}
33
34poly p_Divide(poly p, poly q, const ring r)
35{
36 assume(q!=NULL);
37 if (q==NULL)
38 {
39 WerrorS("div. by 0");
40 return NULL;
41 }
42 if (p==NULL)
43 {
44 p_Delete(&q,r);
45 return NULL;
46 }
47 if ((pNext(q)!=NULL)||rIsPluralRing(r))
48 { /* This means that q != 0 consists of at least two terms*/
49 if(p_GetComp(p,r)==0)
50 {
51 if (!rIsNCRing(r))
52 {
53 if ((rFieldType(r)==n_Q)
54 ||(rFieldType(r)==n_Zp)
55 ||(rFieldType(r)==n_algExt))
56 {
57 poly res=singclap_pdivide(p, q, r);
58 p_Delete(&p,r);
59 p_Delete(&q,r);
60 return res;
61 }
62 if((rFieldType(r)==n_transExt)
63 &&(convSingTrP(p,r))
64 &&(convSingTrP(q,r)))
65 {
66 poly res=singclap_pdivide(p, q, r);
67 p_Delete(&p,r);
68 p_Delete(&q,r);
69 return res;
70 }
71 }
72 // generic division for poly
73 {
74 ideal vi=idInit(1,1); vi->m[0]=q;
75 ideal ui=idInit(1,1); ui->m[0]=p;
76 ideal R; matrix U;
77 ring save_ring=currRing;
78 if (r!=currRing) rChangeCurrRing(r);
79 BITSET save_opt;
80 SI_SAVE_OPT1(save_opt);
81 si_opt_1 &= ~(Sy_bit(OPT_PROT));
82 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
83 SI_RESTORE_OPT1(save_opt);
84 if (r!=save_ring) rChangeCurrRing(save_ring);
85 matrix T = id_Module2formatedMatrix(m,1,1,r);
86 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
87 id_Delete((ideal *)&T,r);
88 id_Delete((ideal *)&U,r);
89 id_Delete(&R,r);
90 //vi->m[0]=NULL; ui->m[0]=NULL;
91 id_Delete(&vi,r);
92 id_Delete(&ui,r);
93 return p;
94 }
95 }
96 else
97 {
98 int comps=p_MaxComp(p,r);
99 ideal I=idInit(comps,1);
100 poly h;
101 int i;
102 // conversion to a list of polys:
103 while (p!=NULL)
104 {
105 i=p_GetComp(p,r)-1;
106 h=pNext(p);
107 pNext(p)=NULL;
108 p_SetComp(p,0,r);
109 I->m[i]=p_Add_q(I->m[i],p,r);
110 p=h;
111 }
112 // division and conversion to vector:
113 h=NULL;
114 p=NULL;
115 for(i=comps-1;i>=0;i--)
116 {
117 if (I->m[i]!=NULL)
118 {
119 if((rFieldType(r)==n_transExt)
120 &&(convSingTrP(I->m[i],r))
121 &&(convSingTrP(q,r))
122 &&(!rIsNCRing(r)))
123 {
124 h=singclap_pdivide(I->m[i],q,r);
125 }
126 else if (((rFieldType(r)==n_Q)
127 ||(rFieldType(r)==n_Zp)
128 ||(rFieldType(r)==n_algExt))
129 &&(!rIsNCRing(r)))
130 h=singclap_pdivide(I->m[i],q,r);
131 else
132 {
133 ideal vi=idInit(1,1); vi->m[0]=q;
134 ideal ui=idInit(1,1); ui->m[0]=I->m[i];
135 ideal R; matrix U;
136 ring save_ring=currRing;
137 if (r!=currRing) rChangeCurrRing(r);
138 BITSET save_opt;
139 SI_SAVE_OPT1(save_opt);
140 si_opt_1 &= ~(Sy_bit(OPT_PROT));
141 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
142 SI_RESTORE_OPT1(save_opt);
143 if (r!=save_ring) rChangeCurrRing(save_ring);
144 if (idIs0(R))
145 {
146 matrix T = id_Module2formatedMatrix(m,1,1,r);
147 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
148 id_Delete((ideal *)&T,r);
149 }
150 else p=NULL;
151 id_Delete((ideal*)&U,r);
152 id_Delete(&R,r);
153 vi->m[0]=NULL; ui->m[0]=NULL;
154 id_Delete(&vi,r);
155 id_Delete(&ui,r);
156 }
157 p_SetCompP(h,i+1,r);
158 p=p_Add_q(p,h,r);
159 }
160 }
161 id_Delete(&I,r);
162 p_Delete(&q,r);
163 return p;
164 }
165 }
166 else
167 { /* This means that q != 0 consists of just one term, or LetterPlace */
168 if (pNext(q)!=NULL)
169 {
170 WerrorS("division over a coefficient domain only implemented for terms");
171 return NULL;
172 }
173 return p_DivideM(p,q,r);
174 }
175 return NULL;
176}
177
178poly pp_Divide(poly p, poly q, const ring r)
179{
180 if (q==NULL)
181 {
182 WerrorS("div. by 0");
183 return NULL;
184 }
185 if (p==NULL)
186 {
187 return NULL;
188 }
189 if ((pNext(q)!=NULL)||rIsPluralRing(r))
190 { /* This means that q != 0 consists of at least two terms*/
191 if(p_GetComp(p,r)==0)
192 {
193 if (!rIsNCRing(r))
194 {
195 if ((rFieldType(r)==n_Q)
196 ||(rFieldType(r)==n_Zp)
197 ||(rFieldType(r)==n_algExt))
198 {
199 poly res=singclap_pdivide(p, q, r);
200 return res;
201 }
202 if((rFieldType(r)==n_transExt)
203 &&(convSingTrP(p,r))
204 &&(convSingTrP(q,r)))
205 {
206 poly res=singclap_pdivide(p, q, r);
207 return res;
208 }
209 }
210 // generic division for poly
211 {
212 ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
213 ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
214 ideal R; matrix U;
215 ring save_ring=currRing;
216 if (r!=currRing) rChangeCurrRing(r);
217 BITSET save_opt;
218 SI_SAVE_OPT1(save_opt);
219 si_opt_1 &= ~(Sy_bit(OPT_PROT));
220 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
221 SI_RESTORE_OPT1(save_opt);
222 if (r!=save_ring) rChangeCurrRing(save_ring);
223 matrix T = id_Module2formatedMatrix(m,1,1,r);
224 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
225 id_Delete((ideal *)&T,r);
226 id_Delete((ideal *)&U,r);
227 id_Delete(&R,r);
228 //vi->m[0]=NULL; ui->m[0]=NULL;
229 id_Delete(&vi,r);
230 id_Delete(&ui,r);
231 return p;
232 }
233 }
234 else
235 {
236 p=p_Copy(p,r);
237 int comps=p_MaxComp(p,r);
238 ideal I=idInit(comps,1);
239 poly h;
240 int i;
241 // conversion to a list of polys:
242 while (p!=NULL)
243 {
244 i=p_GetComp(p,r)-1;
245 h=pNext(p);
246 pNext(p)=NULL;
247 p_SetComp(p,0,r);
248 I->m[i]=p_Add_q(I->m[i],p,r);
249 p=h;
250 }
251 // division and conversion to vector:
252 h=NULL;
253 p=NULL;
254 q=p_Copy(q,r);
255 for(i=comps-1;i>=0;i--)
256 {
257 if (I->m[i]!=NULL)
258 {
259 if(!rIsNCRing(r))
260 {
261 if((rFieldType(r)==n_transExt)
262 &&(convSingTrP(I->m[i],r))
263 &&(convSingTrP(q,r)))
264 {
265 h=singclap_pdivide(I->m[i],q,r);
266 }
267 else if ((rFieldType(r)==n_Q)
268 ||(rFieldType(r)==n_Zp))
269 h=singclap_pdivide(I->m[i],q,r);
270 else
271 {
272 ideal vi=idInit(1,1); vi->m[0]=q;
273 ideal ui=idInit(1,1); ui->m[0]=I->m[i];
274 ideal R; matrix U;
275 ring save_ring=currRing;
276 if (r!=currRing) rChangeCurrRing(r);
277 BITSET save_opt;
278 SI_SAVE_OPT1(save_opt);
279 si_opt_1 &= ~(Sy_bit(OPT_PROT));
280 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
281 SI_RESTORE_OPT1(save_opt);
282 if (r!=save_ring) rChangeCurrRing(save_ring);
283 if (idIs0(R))
284 {
285 matrix T = id_Module2formatedMatrix(m,1,1,r);
286 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
287 id_Delete((ideal *)&T,r);
288 }
289 id_Delete((ideal*)&U,r);
290 id_Delete(&R,r);
291 vi->m[0]=NULL; ui->m[0]=NULL;
292 id_Delete(&vi,r);
293 id_Delete(&ui,r);
294 }
295 }
296 else
297 {
298 ideal vi=idInit(1,1); vi->m[0]=q;
299 ideal ui=idInit(1,1); ui->m[0]=I->m[i];
300 ideal R; matrix U;
301 ring save_ring=currRing;
302 if (r!=currRing) rChangeCurrRing(r);
303 BITSET save_opt;
304 SI_SAVE_OPT1(save_opt);
305 si_opt_1 &= ~(Sy_bit(OPT_PROT));
306 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
307 SI_RESTORE_OPT1(save_opt);
308 if (r!=save_ring) rChangeCurrRing(save_ring);
309 if (idIs0(R))
310 {
311 matrix T = id_Module2formatedMatrix(m,1,1,r);
312 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
313 id_Delete((ideal *)&T,r);
314 }
315 id_Delete((ideal*)&U,r);
316 id_Delete(&R,r);
317 vi->m[0]=NULL; ui->m[0]=NULL;
318 id_Delete(&vi,r);
319 id_Delete(&ui,r);
320 }
321 p_SetCompP(h,i+1,r);
322 p=p_Add_q(p,h,r);
323 }
324 }
325 id_Delete(&I,r);
326 p_Delete(&q,r);
327 return p;
328 }
329 }
330 else
331 { /* This means that q != 0 consists of just one term,
332 or that r is over a coefficient ring. */
333 if (pNext(q)!=NULL)
334 {
335 WerrorS("division over a coefficient domain only implemented for terms");
336 return NULL;
337 }
338 return pp_DivideM(p,q,r);
339 }
340 return NULL;
341}
342
343poly p_DivRem(poly p, poly q, poly &rest, const ring r)
344{
345 assume(q!=NULL);
346 rest=NULL;
347 if (q==NULL)
348 {
349 WerrorS("div. by 0");
350 return NULL;
351 }
352 if (p==NULL)
353 {
354 p_Delete(&q,r);
355 return NULL;
356 }
357 if(p_GetComp(p,r)==0)
358 {
359 if((rFieldType(r)==n_transExt)
360 &&(convSingTrP(p,r))
361 &&(convSingTrP(q,r))
362 &&(!rIsNCRing(r)))
363 {
364 poly res=singclap_pdivide(p, q, r);
365 rest=singclap_pmod(p,q,r);
366 p_Delete(&p,r);
367 p_Delete(&q,r);
368 return res;
369 }
370 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
371 &&(!rField_is_Ring(r))
372 &&(!rIsNCRing(r)))
373 {
374 poly res=singclap_pdivide(p, q, r);
375 rest=singclap_pmod(p,q,r);
376 p_Delete(&p,r);
377 p_Delete(&q,r);
378 return res;
379 }
380 else
381 {
382 ideal vi=idInit(1,1); vi->m[0]=q;
383 ideal ui=idInit(1,1); ui->m[0]=p;
384 ideal R; matrix U;
385 ring save_ring=currRing;
386 if (r!=currRing) rChangeCurrRing(r);
387 BITSET save_opt;
388 SI_SAVE_OPT1(save_opt);
389 si_opt_1 &= ~(Sy_bit(OPT_PROT));
390 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
391 SI_RESTORE_OPT1(save_opt);
392 if (r!=save_ring) rChangeCurrRing(save_ring);
393 p=m->m[0]; m->m[0]=NULL;
394 id_Delete(&m,r);
395 p_SetCompP(p,0,r);
396 rest=R->m[0]; R->m[0]=NULL;
397 id_Delete(&R,r);
398 p_SetCompP(rest,0,r);
399 id_Delete((ideal *)&U,r);
400 //vi->m[0]=NULL; ui->m[0]=NULL;
401 id_Delete(&vi,r);
402 id_Delete(&ui,r);
403 return p;
404 }
405 }
406 return NULL;
407}
408
409poly singclap_gcd ( poly f, poly g, const ring r )
410{
411 poly res=NULL;
412
413 if (f!=NULL)
414 {
415 //if (r->cf->has_simple_Inverse) p_Norm(f,r);
416 if (rField_is_Zp(r)) p_Norm(f,r);
417 else if (!rField_is_Ring(r)) p_Cleardenom(f, r);
418 }
419 if (g!=NULL)
420 {
421 //if (r->cf->has_simple_Inverse) p_Norm(g,r);
422 if (rField_is_Zp(r)) p_Norm(g,r);
423 else if (!rField_is_Ring(r)) p_Cleardenom(g, r);
424 }
425 else return f; // g==0 => gcd=f (but do a p_Cleardenom/pNorm)
426 if (f==NULL) return g; // f==0 => gcd=g (but do a p_Cleardenom/pNorm)
427 if(!rField_is_Ring(r)
428 && (p_IsConstant(f,r)
429 ||p_IsConstant(g,r)))
430 {
431 res=p_One(r);
432 }
433 else if (r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
434 {
435 res=singclap_gcd_r(f,g,r);
436 }
437 else
438 {
439 ideal I=idInit(2,1);
440 I->m[0]=f;
441 I->m[1]=p_Copy(g,r);
442 intvec *w=NULL;
443 ring save_ring=currRing;
444 if (r!=currRing) rChangeCurrRing(r);
445 BITSET save_opt;
446 SI_SAVE_OPT1(save_opt);
447 si_opt_1 &= ~(Sy_bit(OPT_PROT));
448 ideal S1=idSyzygies(I,testHomog,&w);
449 if (w!=NULL) delete w;
450 // expect S1->m[0]=(-g/gcd,f/gcd)
451 if (IDELEMS(S1)!=1) WarnS("error in syzygy computation for GCD");
452 int lp;
453 p_TakeOutComp(&S1->m[0],1,&res,&lp,r);
454 p_Delete(&S1->m[0],r);
455 // GCD is g divided iby (-g/gcd):
456 res=p_Divide(g,res,r);
457 // restore, r, opt:
458 SI_RESTORE_OPT1(save_opt);
459 if (r!=save_ring) rChangeCurrRing(save_ring);
460 // clean the result
461 res=p_Cleardenom(res,r);
462 if (nCoeff_is_Ring(r->cf)) p_Content(res,r);
463 return res;
464 }
465 p_Delete(&f, r);
466 p_Delete(&g, r);
467 return res;
468}
BITSET
#define BITSET
Definition auxiliary.h:85
TRUE
#define TRUE
Definition auxiliary.h:101
FALSE
#define FALSE
Definition auxiliary.h:97
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
p
int p
Definition cfModGcd.cc:4086
g
g
Definition cfModGcd.cc:4098
f
FILE * f
Definition checklibs.c:9
convSingTrP
BOOLEAN convSingTrP(poly p, const ring r)
Definition clapconv.cc:375
singclap_pmod
poly singclap_pmod(poly f, poly g, const ring r)
Definition clapsing.cc:732
singclap_pdivide
poly singclap_pdivide(poly f, poly g, const ring r)
Definition clapsing.cc:649
singclap_gcd_r
poly singclap_gcd_r(poly f, poly g, const ring r)
Definition clapsing.cc:70
n_Q
@ n_Q
rational (GMP) numbers
Definition coeffs.h:30
n_algExt
@ n_algExt
used for all algebraic extensions, i.e., the top-most extension in an extension tower is algebraic
Definition coeffs.h:35
n_Zp
@ n_Zp
\F{p < 2^31}
Definition coeffs.h:29
n_transExt
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition coeffs.h:38
WarnS
#define WarnS
Definition emacs.cc:78
res
CanonicalForm res
Definition facAbsFact.cc:60
w
const CanonicalForm & w
Definition facAbsFact.cc:51
WerrorS
void WerrorS(const char *s)
Definition feFopen.cc:24
VAR
#define VAR
Definition globaldefs.h:5
idSyzygies
ideal idSyzygies(ideal h1, tHomog h, intvec **w, BOOLEAN setSyzComp, BOOLEAN setRegularity, int *deg, GbVariant alg)
Definition ideals.cc:834
idLift
ideal idLift(ideal mod, ideal submod, ideal *rest, BOOLEAN goodShape, BOOLEAN isSB, BOOLEAN divide, matrix *unit, GbVariant alg)
represents the generators of submod in terms of the generators of mod (Matrix(SM)*U-Matrix(rest)) = M...
Definition ideals.cc:1109
idIs0
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
Definition simpleideals.cc:1008
T
STATIC_VAR jList * T
Definition janet.cc:30
h
STATIC_VAR Poly * h
Definition janet.cc:971
p_TakeOutComp
void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r)
Definition p_polys.cc:3620
MATELEM
#define MATELEM(mat, i, j)
1-based access to matrix
Definition matpol.h:29
matrix
ip_smatrix * matrix
Definition matpol.h:43
assume
#define assume(x)
Definition mod2.h:389
p_GetComp
#define p_GetComp(p, r)
Definition monomials.h:64
pNext
#define pNext(p)
Definition monomials.h:36
coeffs
The main handler for Singular numbers which are suitable for Singular polynomials.
ndConvSingNFactoryN
CanonicalForm ndConvSingNFactoryN(number, BOOLEAN, const coeffs)
Definition numbers.cc:313
NULL
#define NULL
Definition omList.c:12
si_opt_1
VAR unsigned si_opt_1
Definition options.c:5
OPT_PROT
#define OPT_PROT
Definition options.h:76
SI_SAVE_OPT1
#define SI_SAVE_OPT1(A)
Definition options.h:21
SI_RESTORE_OPT1
#define SI_RESTORE_OPT1(A)
Definition options.h:24
Sy_bit
#define Sy_bit(x)
Definition options.h:31
TEST_RINGDEP_OPTS
#define TEST_RINGDEP_OPTS
Definition options.h:101
p_DivideM
poly p_DivideM(poly a, poly b, const ring r)
Definition p_polys.cc:1582
p_Content
void p_Content(poly ph, const ring r)
Definition p_polys.cc:2343
pp_DivideM
poly pp_DivideM(poly a, poly b, const ring r)
Definition p_polys.cc:1637
p_Norm
void p_Norm(poly p1, const ring r)
Definition p_polys.cc:3844
p_Cleardenom
poly p_Cleardenom(poly p, const ring r)
Definition p_polys.cc:2893
p_One
poly p_One(const ring r)
Definition p_polys.cc:1314
rChangeCurrRing
void rChangeCurrRing(ring r)
Definition polys.cc:16
coeffs_BIGINT
VAR coeffs coeffs_BIGINT
Definition polys.cc:14
p_Divide
poly p_Divide(poly p, poly q, const ring r)
polynomial division a/b, ignoring the rest via singclap_pdivide resp. idLift destroys a,...
Definition polys.cc:34
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
pp_Divide
poly pp_Divide(poly p, poly q, const ring r)
polynomial division a/b, ignoring the rest via singclap_pdivide resp. idLift does not destroy a,...
Definition polys.cc:178
singclap_gcd
poly singclap_gcd(poly f, poly g, const ring r)
polynomial gcd via singclap_gcd_r resp. idSyzygies destroys f and g
Definition polys.cc:409
p_DivRem
poly p_DivRem(poly p, poly q, poly &rest, const ring r)
Definition polys.cc:343
polys.h
Compatibility layer for legacy polynomial operations (over currRing).
p_SetGlobals
void p_SetGlobals(const ring r, BOOLEAN complete)
set all properties of a new ring - also called by rComplete
Definition ring.cc:3494
rTest
#define rTest(r)
Definition ring.h:799
rField_is_Ring
#define rField_is_Ring(R)
Definition ring.h:491
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition simpleideals.cc:35
id_Delete
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
Definition simpleideals.cc:123
id_Module2formatedMatrix
matrix id_Module2formatedMatrix(ideal mod, int rows, int cols, const ring R)
Definition simpleideals.cc:1721
IDELEMS
#define IDELEMS(i)
Definition simpleideals.h:23
R
#define R
Definition sirandom.c:27
testHomog
@ testHomog
Definition structs.h:34
Generated on for My Project by   doxygen 1.18.0 for Singular