Loading...
Searching...
No Matches
polys.cc
Go to the documentation of this file.
#include "kernel/mod2.h"
#include "misc/options.h"
#include "polys.h"
#include "kernel/ideals.h"
#include "kernel/ideals.h"
#include "polys/clapsing.h"
#include "polys/clapconv.h"
/// Widely used global variable which specifies the current polynomial ring for Singular interpreter and legacy implementations.
/// @Note: one should avoid using it in newer designs, for example due to possible problems in parallelization with threads.
VAR ring currRing = NULL;
VAR coeffs coeffs_BIGINT;
void rChangeCurrRing(ring r)
{
if (currRing!=NULL)
currRing->options = si_opt_1 & TEST_RINGDEP_OPTS;
//------------ set global ring vars --------------------------------
currRing = r;
if( r != NULL )
{
rTest(r);
//------------ global variables related to coefficients ------------
assume( r->cf!= NULL );
nSetChar(r->cf);
//------------ global variables related to polys
p_SetGlobals(r); // also setting TEST_RINGDEP_OPTS
//------------ global variables related to factory -----------------
}
}
poly p_Divide(poly p, poly q, const ring r)
{
assume(q!=NULL);
if (q==NULL)
{
WerrorS("div. by 0");
return NULL;
}
if (p==NULL)
{
p_Delete(&q,r);
return NULL;
}
if ((pNext(q)!=NULL)||rIsPluralRing(r))
{ /* This means that q != 0 consists of at least two terms*/
if(p_GetComp(p,r)==0)
{
if (!rIsNCRing(r))
{
if ((rFieldType(r)==n_Q)
||(rFieldType(r)==n_Zp)
||(rFieldType(r)==n_algExt))
{
poly res=singclap_pdivide(p, q, r);
p_Delete(&p,r);
p_Delete(&q,r);
return res;
}
if((rFieldType(r)==n_transExt)
&&(convSingTrP(p,r))
&&(convSingTrP(q,r)))
{
poly res=singclap_pdivide(p, q, r);
p_Delete(&p,r);
p_Delete(&q,r);
return res;
}
}
// generic division for poly
{
ideal vi=idInit(1,1); vi->m[0]=q;
ideal ui=idInit(1,1); ui->m[0]=p;
ideal R; matrix U;
ring save_ring=currRing;
if (r!=currRing) rChangeCurrRing(r);
BITSET save_opt;
SI_SAVE_OPT1(save_opt);
si_opt_1 &= ~(Sy_bit(OPT_PROT));
ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
SI_RESTORE_OPT1(save_opt);
if (r!=save_ring) rChangeCurrRing(save_ring);
matrix T = id_Module2formatedMatrix(m,1,1,r);
p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
id_Delete((ideal *)&T,r);
id_Delete((ideal *)&U,r);
id_Delete(&R,r);
//vi->m[0]=NULL; ui->m[0]=NULL;
id_Delete(&vi,r);
id_Delete(&ui,r);
return p;
}
}
else
{
int comps=p_MaxComp(p,r);
ideal I=idInit(comps,1);
poly h;
int i;
// conversion to a list of polys:
while (p!=NULL)
{
i=p_GetComp(p,r)-1;
h=pNext(p);
pNext(p)=NULL;
p_SetComp(p,0,r);
I->m[i]=p_Add_q(I->m[i],p,r);
p=h;
}
// division and conversion to vector:
h=NULL;
p=NULL;
for(i=comps-1;i>=0;i--)
{
if (I->m[i]!=NULL)
{
if((rFieldType(r)==n_transExt)
&&(convSingTrP(I->m[i],r))
&&(convSingTrP(q,r))
&&(!rIsNCRing(r)))
{
h=singclap_pdivide(I->m[i],q,r);
}
else if (((rFieldType(r)==n_Q)
||(rFieldType(r)==n_Zp)
||(rFieldType(r)==n_algExt))
&&(!rIsNCRing(r)))
h=singclap_pdivide(I->m[i],q,r);
else
{
ideal vi=idInit(1,1); vi->m[0]=q;
ideal ui=idInit(1,1); ui->m[0]=I->m[i];
ideal R; matrix U;
ring save_ring=currRing;
if (r!=currRing) rChangeCurrRing(r);
BITSET save_opt;
SI_SAVE_OPT1(save_opt);
si_opt_1 &= ~(Sy_bit(OPT_PROT));
ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
SI_RESTORE_OPT1(save_opt);
if (r!=save_ring) rChangeCurrRing(save_ring);
if (idIs0(R))
{
matrix T = id_Module2formatedMatrix(m,1,1,r);
p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
id_Delete((ideal *)&T,r);
}
else p=NULL;
id_Delete((ideal*)&U,r);
id_Delete(&R,r);
vi->m[0]=NULL; ui->m[0]=NULL;
id_Delete(&vi,r);
id_Delete(&ui,r);
}
p_SetCompP(h,i+1,r);
p=p_Add_q(p,h,r);
}
}
id_Delete(&I,r);
p_Delete(&q,r);
return p;
}
}
else
{ /* This means that q != 0 consists of just one term, or LetterPlace */
if (pNext(q)!=NULL)
{
WerrorS("division over a coefficient domain only implemented for terms");
return NULL;
}
return p_DivideM(p,q,r);
}
return NULL;
}
poly pp_Divide(poly p, poly q, const ring r)
{
if (q==NULL)
{
WerrorS("div. by 0");
return NULL;
}
if (p==NULL)
{
return NULL;
}
if ((pNext(q)!=NULL)||rIsPluralRing(r))
{ /* This means that q != 0 consists of at least two terms*/
if(p_GetComp(p,r)==0)
{
if (!rIsNCRing(r))
{
if ((rFieldType(r)==n_Q)
||(rFieldType(r)==n_Zp)
||(rFieldType(r)==n_algExt))
{
poly res=singclap_pdivide(p, q, r);
return res;
}
if((rFieldType(r)==n_transExt)
&&(convSingTrP(p,r))
&&(convSingTrP(q,r)))
{
poly res=singclap_pdivide(p, q, r);
return res;
}
}
// generic division for poly
{
ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
ideal R; matrix U;
ring save_ring=currRing;
if (r!=currRing) rChangeCurrRing(r);
BITSET save_opt;
SI_SAVE_OPT1(save_opt);
si_opt_1 &= ~(Sy_bit(OPT_PROT));
ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
SI_RESTORE_OPT1(save_opt);
if (r!=save_ring) rChangeCurrRing(save_ring);
matrix T = id_Module2formatedMatrix(m,1,1,r);
p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
id_Delete((ideal *)&T,r);
id_Delete((ideal *)&U,r);
id_Delete(&R,r);
//vi->m[0]=NULL; ui->m[0]=NULL;
id_Delete(&vi,r);
id_Delete(&ui,r);
return p;
}
}
else
{
p=p_Copy(p,r);
int comps=p_MaxComp(p,r);
ideal I=idInit(comps,1);
poly h;
int i;
// conversion to a list of polys:
while (p!=NULL)
{
i=p_GetComp(p,r)-1;
h=pNext(p);
pNext(p)=NULL;
p_SetComp(p,0,r);
I->m[i]=p_Add_q(I->m[i],p,r);
p=h;
}
// division and conversion to vector:
h=NULL;
p=NULL;
q=p_Copy(q,r);
for(i=comps-1;i>=0;i--)
{
if (I->m[i]!=NULL)
{
if(!rIsNCRing(r))
{
if((rFieldType(r)==n_transExt)
&&(convSingTrP(I->m[i],r))
&&(convSingTrP(q,r)))
{
h=singclap_pdivide(I->m[i],q,r);
}
else if ((rFieldType(r)==n_Q)
||(rFieldType(r)==n_Zp))
h=singclap_pdivide(I->m[i],q,r);
else
{
ideal vi=idInit(1,1); vi->m[0]=q;
ideal ui=idInit(1,1); ui->m[0]=I->m[i];
ideal R; matrix U;
ring save_ring=currRing;
if (r!=currRing) rChangeCurrRing(r);
BITSET save_opt;
SI_SAVE_OPT1(save_opt);
si_opt_1 &= ~(Sy_bit(OPT_PROT));
ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
SI_RESTORE_OPT1(save_opt);
if (r!=save_ring) rChangeCurrRing(save_ring);
if (idIs0(R))
{
matrix T = id_Module2formatedMatrix(m,1,1,r);
p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
id_Delete((ideal *)&T,r);
}
id_Delete((ideal*)&U,r);
id_Delete(&R,r);
vi->m[0]=NULL; ui->m[0]=NULL;
id_Delete(&vi,r);
id_Delete(&ui,r);
}
}
else
{
ideal vi=idInit(1,1); vi->m[0]=q;
ideal ui=idInit(1,1); ui->m[0]=I->m[i];
ideal R; matrix U;
ring save_ring=currRing;
if (r!=currRing) rChangeCurrRing(r);
BITSET save_opt;
SI_SAVE_OPT1(save_opt);
si_opt_1 &= ~(Sy_bit(OPT_PROT));
ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
SI_RESTORE_OPT1(save_opt);
if (r!=save_ring) rChangeCurrRing(save_ring);
if (idIs0(R))
{
matrix T = id_Module2formatedMatrix(m,1,1,r);
p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
id_Delete((ideal *)&T,r);
}
id_Delete((ideal*)&U,r);
id_Delete(&R,r);
vi->m[0]=NULL; ui->m[0]=NULL;
id_Delete(&vi,r);
id_Delete(&ui,r);
}
p_SetCompP(h,i+1,r);
p=p_Add_q(p,h,r);
}
}
id_Delete(&I,r);
p_Delete(&q,r);
return p;
}
}
else
{ /* This means that q != 0 consists of just one term,
or that r is over a coefficient ring. */
if (pNext(q)!=NULL)
{
WerrorS("division over a coefficient domain only implemented for terms");
return NULL;
}
return pp_DivideM(p,q,r);
}
return NULL;
}
poly p_DivRem(poly p, poly q, poly &rest, const ring r)
{
assume(q!=NULL);
rest=NULL;
if (q==NULL)
{
WerrorS("div. by 0");
return NULL;
}
if (p==NULL)
{
p_Delete(&q,r);
return NULL;
}
if(p_GetComp(p,r)==0)
{
if((rFieldType(r)==n_transExt)
&&(convSingTrP(p,r))
&&(convSingTrP(q,r))
&&(!rIsNCRing(r)))
{
poly res=singclap_pdivide(p, q, r);
rest=singclap_pmod(p,q,r);
p_Delete(&p,r);
p_Delete(&q,r);
return res;
}
else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
&&(!rField_is_Ring(r))
&&(!rIsNCRing(r)))
{
poly res=singclap_pdivide(p, q, r);
rest=singclap_pmod(p,q,r);
p_Delete(&p,r);
p_Delete(&q,r);
return res;
}
else
{
ideal vi=idInit(1,1); vi->m[0]=q;
ideal ui=idInit(1,1); ui->m[0]=p;
ideal R; matrix U;
ring save_ring=currRing;
if (r!=currRing) rChangeCurrRing(r);
BITSET save_opt;
SI_SAVE_OPT1(save_opt);
si_opt_1 &= ~(Sy_bit(OPT_PROT));
ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
SI_RESTORE_OPT1(save_opt);
if (r!=save_ring) rChangeCurrRing(save_ring);
p=m->m[0]; m->m[0]=NULL;
id_Delete(&m,r);
p_SetCompP(p,0,r);
rest=R->m[0]; R->m[0]=NULL;
id_Delete(&R,r);
p_SetCompP(rest,0,r);
id_Delete((ideal *)&U,r);
//vi->m[0]=NULL; ui->m[0]=NULL;
id_Delete(&vi,r);
id_Delete(&ui,r);
return p;
}
}
return NULL;
}
poly singclap_gcd ( poly f, poly g, const ring r )
{
poly res=NULL;
if (f!=NULL)
{
//if (r->cf->has_simple_Inverse) p_Norm(f,r);
if (rField_is_Zp(r)) p_Norm(f,r);
else if (!rField_is_Ring(r)) p_Cleardenom(f, r);
}
if (g!=NULL)
{
//if (r->cf->has_simple_Inverse) p_Norm(g,r);
if (rField_is_Zp(r)) p_Norm(g,r);
else if (!rField_is_Ring(r)) p_Cleardenom(g, r);
}
else return f; // g==0 => gcd=f (but do a p_Cleardenom/pNorm)
if (f==NULL) return g; // f==0 => gcd=g (but do a p_Cleardenom/pNorm)
if(!rField_is_Ring(r)
&& (p_IsConstant(f,r)
||p_IsConstant(g,r)))
{
res=p_One(r);
}
else if (r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
{
res=singclap_gcd_r(f,g,r);
}
else
{
ideal I=idInit(2,1);
I->m[0]=f;
I->m[1]=p_Copy(g,r);
intvec *w=NULL;
ring save_ring=currRing;
if (r!=currRing) rChangeCurrRing(r);
BITSET save_opt;
SI_SAVE_OPT1(save_opt);
si_opt_1 &= ~(Sy_bit(OPT_PROT));
ideal S1=idSyzygies(I,testHomog,&w);
if (w!=NULL) delete w;
// expect S1->m[0]=(-g/gcd,f/gcd)
if (IDELEMS(S1)!=1) WarnS("error in syzygy computation for GCD");
int lp;
p_TakeOutComp(&S1->m[0],1,&res,&lp,r);
p_Delete(&S1->m[0],r);
// GCD is g divided iby (-g/gcd):
res=p_Divide(g,res,r);
// restore, r, opt:
SI_RESTORE_OPT1(save_opt);
if (r!=save_ring) rChangeCurrRing(save_ring);
// clean the result
res=p_Cleardenom(res,r);
if (nCoeff_is_Ring(r->cf)) p_Content(res,r);
return res;
}
p_Delete(&f, r);
p_Delete(&g, r);
return res;
}
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
clapconv.h
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
clapsing.h
intvec
Definition intvec.h:24
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
nSetChar
static FORCE_INLINE void nSetChar(const coeffs r)
initialisations after each ring change
Definition coeffs.h:446
nCoeff_is_Ring
static FORCE_INLINE BOOLEAN nCoeff_is_Ring(const coeffs r)
Definition coeffs.h:732
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:836
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:1111
ideals.h
idIs0
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
Definition simpleideals.cc:960
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
mod2.h
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
options.h
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
p_Add_q
static poly p_Add_q(poly p, poly q, const ring r)
Definition p_polys.h:938
p_SetCompP
static void p_SetCompP(poly p, int i, ring r)
Definition p_polys.h:256
p_SetComp
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition p_polys.h:249
p_IsConstant
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition p_polys.h:1985
p_MaxComp
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
Definition p_polys.h:294
p_Delete
static void p_Delete(poly *p, const ring r)
Definition p_polys.h:903
p_Copy
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition p_polys.h:848
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
rField_is_Zp
static BOOLEAN rField_is_Zp(const ring r)
Definition ring.h:506
rFieldType
static n_coeffType rFieldType(const ring r)
the type of the coefficient filed of r (n_Zp, n_Q, etc)
Definition ring.h:562
rIsPluralRing
static BOOLEAN rIsPluralRing(const ring r)
we must always have this test!
Definition ring.h:406
rIsNCRing
static BOOLEAN rIsNCRing(const ring r)
Definition ring.h:427
rTest
#define rTest(r)
Definition ring.h:794
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:1640
IDELEMS
#define IDELEMS(i)
Definition simpleideals.h:23
R
#define R
Definition sirandom.c:27
testHomog
@ testHomog
Definition structs.h:34