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;
14
15void 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
31poly 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 p=m->m[0]; m->m[0]=NULL;
81 id_Delete(&m,r);
82 p_SetCompP(p,0,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
174poly pp_Divide(poly p, poly q, const ring r)
175{
176 if (q==NULL)
177 {
178 WerrorS("div. by 0");
179 return NULL;
180 }
181 if (p==NULL)
182 {
183 return NULL;
184 }
185 if ((pNext(q)!=NULL)||rIsPluralRing(r))
186 { /* This means that q != 0 consists of at least two terms*/
187 if(p_GetComp(p,r)==0)
188 {
189 if((rFieldType(r)==n_transExt)
190 &&(convSingTrP(p,r))
191 &&(convSingTrP(q,r))
192 &&(!rIsNCRing(r)))
193 {
194 poly res=singclap_pdivide(p, q, r);
195 return res;
196 }
197 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
198 &&(!rField_is_Ring(r))
199 &&(!rIsNCRing(r)))
200 {
201 poly res=singclap_pdivide(p, q, r);
202 return res;
203 }
204 else
205 {
206 ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
207 ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
208 ideal R; matrix U;
209 ring save_ring=currRing;
210 if (r!=currRing) rChangeCurrRing(r);
211 int save_opt;
212 SI_SAVE_OPT1(save_opt);
213 si_opt_1 &= ~(Sy_bit(OPT_PROT));
214 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
215 SI_RESTORE_OPT1(save_opt);
216 if (r!=save_ring) rChangeCurrRing(save_ring);
217 matrix T = id_Module2formatedMatrix(m,1,1,r);
218 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
219 id_Delete((ideal *)&T,r);
220 id_Delete((ideal *)&U,r);
221 id_Delete(&R,r);
222 //vi->m[0]=NULL; ui->m[0]=NULL;
223 id_Delete(&vi,r);
224 id_Delete(&ui,r);
225 return p;
226 }
227 }
228 else
229 {
230 p=p_Copy(p,r);
231 int comps=p_MaxComp(p,r);
232 ideal I=idInit(comps,1);
233 poly h;
234 int i;
235 // conversion to a list of polys:
236 while (p!=NULL)
237 {
238 i=p_GetComp(p,r)-1;
239 h=pNext(p);
240 pNext(p)=NULL;
241 p_SetComp(p,0,r);
242 I->m[i]=p_Add_q(I->m[i],p,r);
243 p=h;
244 }
245 // division and conversion to vector:
246 h=NULL;
247 p=NULL;
248 q=p_Copy(q,r);
249 for(i=comps-1;i>=0;i--)
250 {
251 if (I->m[i]!=NULL)
252 {
253 if((rFieldType(r)==n_transExt)
254 &&(convSingTrP(I->m[i],r))
255 &&(convSingTrP(q,r))
256 &&(!rIsNCRing(r)))
257 {
258 h=singclap_pdivide(I->m[i],q,r);
259 }
260 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
261 &&(!rField_is_Ring(r))
262 &&(!rIsNCRing(r)))
263 h=singclap_pdivide(I->m[i],q,r);
264 else
265 {
266 ideal vi=idInit(1,1); vi->m[0]=q;
267 ideal ui=idInit(1,1); ui->m[0]=I->m[i];
268 ideal R; matrix U;
269 ring save_ring=currRing;
270 if (r!=currRing) rChangeCurrRing(r);
271 int save_opt;
272 SI_SAVE_OPT1(save_opt);
273 si_opt_1 &= ~(Sy_bit(OPT_PROT));
274 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
275 SI_RESTORE_OPT1(save_opt);
276 if (r!=save_ring) rChangeCurrRing(save_ring);
277 if (idIs0(R))
278 {
279 matrix T = id_Module2formatedMatrix(m,1,1,r);
280 p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
281 id_Delete((ideal *)&T,r);
282 }
283 else p=NULL;
284 id_Delete((ideal*)&U,r);
285 id_Delete(&R,r);
286 vi->m[0]=NULL; ui->m[0]=NULL;
287 id_Delete(&vi,r);
288 id_Delete(&ui,r);
289 }
290 p_SetCompP(h,i+1,r);
291 p=p_Add_q(p,h,r);
292 }
293 }
294 id_Delete(&I,r);
295 p_Delete(&q,r);
296 return p;
297 }
298 }
299 else
300 { /* This means that q != 0 consists of just one term,
301 or that r is over a coefficient ring. */
302#ifdef HAVE_RINGS
303 if (pNext(q)!=NULL)
304 {
305 WerrorS("division over a coefficient domain only implemented for terms");
306 return NULL;
307 }
308#endif
309 return pp_DivideM(p,q,r);
310 }
311 return NULL;
312}
313
314poly p_DivRem(poly p, poly q, poly &rest, const ring r)
315{
316 assume(q!=NULL);
317 rest=NULL;
318 if (q==NULL)
319 {
320 WerrorS("div. by 0");
321 return NULL;
322 }
323 if (p==NULL)
324 {
325 p_Delete(&q,r);
326 return NULL;
327 }
328 if(p_GetComp(p,r)==0)
329 {
330 if((rFieldType(r)==n_transExt)
331 &&(convSingTrP(p,r))
332 &&(convSingTrP(q,r))
333 &&(!rIsNCRing(r)))
334 {
335 poly res=singclap_pdivide(p, q, r);
336 rest=singclap_pmod(p,q,r);
337 p_Delete(&p,r);
338 p_Delete(&q,r);
339 return res;
340 }
341 else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
342 &&(!rField_is_Ring(r))
343 &&(!rIsNCRing(r)))
344 {
345 poly res=singclap_pdivide(p, q, r);
346 rest=singclap_pmod(p,q,r);
347 p_Delete(&p,r);
348 p_Delete(&q,r);
349 return res;
350 }
351 else
352 {
353 ideal vi=idInit(1,1); vi->m[0]=q;
354 ideal ui=idInit(1,1); ui->m[0]=p;
355 ideal R; matrix U;
356 ring save_ring=currRing;
357 if (r!=currRing) rChangeCurrRing(r);
358 int save_opt;
359 SI_SAVE_OPT1(save_opt);
360 si_opt_1 &= ~(Sy_bit(OPT_PROT));
361 ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
362 SI_RESTORE_OPT1(save_opt);
363 if (r!=save_ring) rChangeCurrRing(save_ring);
364 p=m->m[0]; m->m[0]=NULL;
365 id_Delete(&m,r);
366 p_SetCompP(p,0,r);
367 rest=R->m[0]; R->m[0]=NULL;
368 id_Delete(&R,r);
369 p_SetCompP(rest,0,r);
370 id_Delete((ideal *)&U,r);
371 //vi->m[0]=NULL; ui->m[0]=NULL;
372 id_Delete(&vi,r);
373 id_Delete(&ui,r);
374 return p;
375 }
376 }
377 return NULL;
378}
379
380poly singclap_gcd ( poly f, poly g, const ring r )
381{
382 poly res=NULL;
383
384 if (f!=NULL)
385 {
386 //if (r->cf->has_simple_Inverse) p_Norm(f,r);
387 if (rField_is_Zp(r)) p_Norm(f,r);
388 else if (!rField_is_Ring(r)) p_Cleardenom(f, r);
389 }
390 if (g!=NULL)
391 {
392 //if (r->cf->has_simple_Inverse) p_Norm(g,r);
393 if (rField_is_Zp(r)) p_Norm(g,r);
394 else if (!rField_is_Ring(r)) p_Cleardenom(g, r);
395 }
396 else return f; // g==0 => gcd=f (but do a p_Cleardenom/pNorm)
397 if (f==NULL) return g; // f==0 => gcd=g (but do a p_Cleardenom/pNorm)
398 if(!rField_is_Ring(r)
399 && (p_IsConstant(f,r)
400 ||p_IsConstant(g,r)))
401 {
402 res=p_One(r);
403 }
404 else if (r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
405 {
406 res=singclap_gcd_r(f,g,r);
407 }
408 else
409 {
410 ideal I=idInit(2,1);
411 I->m[0]=f;
412 I->m[1]=p_Copy(g,r);
413 intvec *w=NULL;
414 ring save_ring=currRing;
415 if (r!=currRing) rChangeCurrRing(r);
416 int save_opt;
417 SI_SAVE_OPT1(save_opt);
418 si_opt_1 &= ~(Sy_bit(OPT_PROT));
419 ideal S1=idSyzygies(I,testHomog,&w);
420 if (w!=NULL) delete w;
421 // expect S1->m[0]=(-g/gcd,f/gcd)
422 if (IDELEMS(S1)!=1) WarnS("error in syzygy computation for GCD");
423 int lp;
424 p_TakeOutComp(&S1->m[0],1,&res,&lp,r);
425 p_Delete(&S1->m[0],r);
426 // GCD is g divided iby (-g/gcd):
427 res=p_Divide(g,res,r);
428 // restore, r, opt:
429 SI_RESTORE_OPT1(save_opt);
430 if (r!=save_ring) rChangeCurrRing(save_ring);
431 // clean the result
432 res=p_Cleardenom(res,r);
433 if (nCoeff_is_Ring(r->cf)) p_Content(res,r);
434 return res;
435 }
436 p_Delete(&f, r);
437 p_Delete(&g, r);
438 return res;
439}
TRUE
#define TRUE
Definition: auxiliary.h:100
FALSE
#define FALSE
Definition: auxiliary.h:96
m
int m
Definition: cfEzgcd.cc:128
i
int i
Definition: cfEzgcd.cc:132
p
int p
Definition: cfModGcd.cc:4078
g
g
Definition: cfModGcd.cc:4090
f
FILE * f
Definition: checklibs.c:9
convSingTrP
BOOLEAN convSingTrP(poly p, const ring r)
Definition: clapconv.cc:352
singclap_pmod
poly singclap_pmod(poly f, poly g, const ring r)
Definition: clapsing.cc:702
singclap_pdivide
poly singclap_pdivide(poly f, poly g, const ring r)
Definition: clapsing.cc:624
singclap_gcd_r
poly singclap_gcd_r(poly f, poly g, const ring r)
Definition: clapsing.cc:68
intvec
Definition: intvec.h:23
ip_smatrix
Definition: matpol.h:15
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
nSetChar
static FORCE_INLINE void nSetChar(const coeffs r)
initialisations after each ring change
Definition: coeffs.h:437
nCoeff_is_Ring
static FORCE_INLINE BOOLEAN nCoeff_is_Ring(const coeffs r)
Definition: coeffs.h:727
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:830
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:1105
idIs0
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
Definition: simpleideals.cc:917
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:3496
MATELEM
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
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
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
p_DivideM
poly p_DivideM(poly a, poly b, const ring r)
Definition: p_polys.cc:1578
p_Content
void p_Content(poly ph, const ring r)
Definition: p_polys.cc:2295
pp_DivideM
poly pp_DivideM(poly a, poly b, const ring r)
Definition: p_polys.cc:1633
p_Norm
void p_Norm(poly p1, const ring r)
Definition: p_polys.cc:3719
p_Cleardenom
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2845
p_One
poly p_One(const ring r)
Definition: p_polys.cc:1313
p_Add_q
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
p_SetCompP
static void p_SetCompP(poly p, int i, ring r)
Definition: p_polys.h:252
p_SetComp
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:245
p_IsConstant
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1962
p_MaxComp
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
Definition: p_polys.h:290
p_Delete
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
p_Copy
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:844
rChangeCurrRing
void rChangeCurrRing(ring r)
Definition: polys.cc:15
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:31
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:174
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:380
p_DivRem
poly p_DivRem(poly p, poly q, poly &rest, const ring r)
Definition: polys.cc:314
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:3415
rField_is_Zp
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:500
rFieldType
static n_coeffType rFieldType(const ring r)
the type of the coefficient filed of r (n_Zp, n_Q, etc)
Definition: ring.h:556
rIsPluralRing
static BOOLEAN rIsPluralRing(const ring r)
we must always have this test!
Definition: ring.h:400
rIsNCRing
static BOOLEAN rIsNCRing(const ring r)
Definition: ring.h:421
rTest
#define rTest(r)
Definition: ring.h:782
rField_is_Ring
#define rField_is_Ring(R)
Definition: ring.h:485
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:1522
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:23
R
#define R
Definition: sirandom.c:27
testHomog
@ testHomog
Definition: structs.h:38
Generated on Mon Feb 27 2023 10:53:49 for My Project by   doxygen 1.9.5 for Singular