File: [local] / OpenXM / src / kan96xx / Kan / red.c (download)
Revision 1.5, Wed Jul 30 09:00:52 2003 UTC (21 years, 2 months ago) by takayama
Branch: MAIN
Changes since 1.4: +2 -2
lines
Fixed a bug of groebner_gen().
ecartd.gb : groebner basis computation with ecart division and dehomogenization.
|
/* $OpenXM: OpenXM/src/kan96xx/Kan/red.c,v 1.5 2003/07/30 09:00:52 takayama Exp $ */
#include <stdio.h>
#include "datatype.h"
#include "extern2.h"
#include "gradedset.h"
#define mymax(p,q) (p>q?p:q)
int DebugReductionRed = 0;
extern int Sugar;
struct spValue sp_gen(f,g)
POLY f;
POLY g;
/* the results may be rewritten. */
{
struct spValue r;
POLY a;
POLY b;
MONOMIAL tf,tg;
MONOMIAL ta,tb;
int n,i;
struct coeff *ac;
struct coeff *bc;
int amodify,bmodify,c,ell;
if ((f ISZERO) || (g ISZERO)) {
warningGradedSet("sp_gen(): f and g must not be zero. Returns zero.");
r.a = ZERO;
r.b = ZERO;
return(r);
}
checkRingSp(f,g,r);
tf = f->m; tg = g->m;
n = tf->ringp->n;
ta = newMonomial(tf->ringp);
tb = newMonomial(tf->ringp);
for (i=0; i<n; i++) {
ta->e[i].x = mymax(tf->e[i].x,tg->e[i].x) - tf->e[i].x;
ta->e[i].D = mymax(tf->e[i].D,tg->e[i].D) - tf->e[i].D;
tb->e[i].x = mymax(tf->e[i].x,tg->e[i].x) - tg->e[i].x;
tb->e[i].D = mymax(tf->e[i].D,tg->e[i].D) - tg->e[i].D;
}
switch (f->coeffp->tag) {
case INTEGER:
a = newCell(coeffCopy(g->coeffp),ta);
b = newCell(coeffCopy(f->coeffp),tb);
b->coeffp->val.i = -(b->coeffp->val.i);
break;
case MP_INTEGER:
{
MP_INT *gcd,*ac;
gcd = newMP_INT();
mpz_gcd(gcd,f->coeffp->val.bigp,g->coeffp->val.bigp);
ac = newMP_INT();
mpz_mdiv(ac,g->coeffp->val.bigp,gcd);
mpz_mdiv(gcd,f->coeffp->val.bigp,gcd);
mpz_neg(gcd,gcd);
a = newCell(mpintToCoeff(ac,tf->ringp),ta);
b = newCell(mpintToCoeff(gcd,tf->ringp),tb);
}
break;
case POLY_COEFF:
c = f->m->ringp->c; ell = f->m->ringp->l;
if (ell-c > 0) { /* q-case */
amodify = 0;
for (i=c; i<ell; i++) {
amodify += (tb->e[i].D)*(tg->e[i].x);
}
bmodify = 0;
for (i=c; i<ell; i++) {
bmodify += (ta->e[i].D)*(tf->e[i].x);
}
if (amodify > bmodify) {
amodify = amodify-bmodify;
bmodify = 0;
}else{
bmodify = bmodify - amodify;
amodify = 0;
}
}
if (amodify) {
/* a ---> q^amodify a, b ---> - b */
ac = polyToCoeff(cxx(1,0,amodify,f->m->ringp->next),f->m->ringp);
Cmult(ac,ac,g->coeffp);
a = newCell(ac,ta);
bc = coeffNeg(f->coeffp,f->m->ringp);
b = newCell(bc,tb);
}else{
/* a ---> a, b ---> -q^bmodify b */
a = newCell(coeffCopy(g->coeffp),ta);
bc = coeffNeg(f->coeffp,f->m->ringp);
Cmult(bc,polyToCoeff(cxx(1,0,bmodify,f->m->ringp->next),f->m->ringp),bc);
b = newCell(bc,tb);
}
break;
default:
errorGradedSet("sp_gen(): Unsupported tag.");
break;
}
/* r.sp = ppAddv(ppMult(a,f),ppMult(b,g)); */
r.a = a;
r.b = b;
return(r);
}
POLY reduction1_gen_debug(f,g,needSyz,c,h)
POLY f;
POLY g;
int needSyz;
POLY *c; /* set */
POLY *h; /* set */
/* f must be reducible by g. r = c*f + h*g */
{
extern struct ring *CurrentRingp;
struct ring *rp;
struct spValue sv;
POLY f2;
int showLength = 0;
/*DebugReductionRed = 1;*/
if (needSyz) {
if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
*c = cxx(1,0,0,rp);
*h = ZERO;
}
sv = (*sp)(f,g);
f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
if (showLength) {printf(" [%d] ",pLength(f)); fflush(stdout);}
if (!isOrdered(f2) || !isOrdered(f)) {
if (!isOrdered(f)) {
printf("\nf is not ordered polynomial.");
}else if (!isOrdered(f)) {
printf("\nf2 is not ordered polynomial.");
}
printf("f=%s,\nf2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
getchar();
getchar();
}
if ((*mmLarger)(f2,f) >= 1) {
printf("error in reduction1.");
printf("f=%s --> f2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
printf("f2 = (%s)*f+(%s)*(%s)\n",POLYToString(sv.a,'*',1),
POLYToString(sv.b,'*',1),POLYToString(g,'*',1));
getchar();
getchar();
}
if (DebugReductionRed) {
printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
}
f = f2;
if (needSyz) {
*c = ppMult(sv.a,*c);
*h = ppAdd(ppMult(sv.a,*h),sv.b);
}
while ((*isReducible)(f,g)) {
sv = (*sp)(f,g);
f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
if (DebugReductionRed) {
printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
}
if (showLength) {printf(" [%d] ",pLength(f)); fflush(stdout);}
if (!isOrdered(f2) || !isOrdered(f)) {
if (!isOrdered(f)) {
printf("\nf is not ordered polynomial.");
}else if (!isOrdered(f)) {
printf("\nf2 is not ordered polynomial.");
}
printf("f=%s,\nf2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
getchar();
getchar();
}
if ((*mmLarger)(f2,f) >= 1) {
printf("error in reduction1.");
printf("f=%s --> f2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
printf("f2 = (%s)*f+(%s)*(%s)\n",POLYToString(sv.a,'*',1),
POLYToString(sv.b,'*',1),
POLYToString(g,'*',1));
getchar();
getchar();
}
f = f2;
if (needSyz) {
*c = ppMult(sv.a,*c);
*h = ppAdd(ppMult(sv.a,*h),sv.b);
}
}
return(f);
}
POLY reduction1_gen(f,g,needSyz,c,h)
POLY f;
POLY g;
int needSyz;
POLY *c; /* set */
POLY *h; /* set */
/* f must be reducible by g. r = c*f + h*g */
{
extern struct ring *CurrentRingp;
struct ring *rp;
struct spValue sv;
POLY f2;
if (needSyz) {
if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
*c = cxx(1,0,0,rp);
*h = ZERO;
}
sv = (*sp)(f,g);
f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
if (DebugReductionRed) {
printf("c=%s, d=%s, g=%s: f --> c*f + d*g.\n",
POLYToString(sv.a,'*',1),POLYToString(sv.b,'*',1),POLYToString(g,'*',1));
printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
}
f = f2;
if (needSyz) {
*c = ppMult(sv.a,*c);
*h = ppAdd(ppMult(sv.a,*h),sv.b);
}
while ((*isReducible)(f,g)) {
sv = (*sp)(f,g);
f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
if (DebugReductionRed) {
printf("! c=%s, d=%s, g=%s: f --> c*f + d*g.\n",
POLYToString(sv.a,'*',1),POLYToString(sv.b,'*',1),POLYToString(g,'*',1));
printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
}
f = f2;
if (needSyz) {
*c = ppMult(sv.a,*c);
*h = ppAdd(ppMult(sv.a,*h),sv.b);
}
}
return(f);
}
POLY reduction1Cdr_gen(f,fs,g,needSyz,c,h)
POLY f;
POLY fs;
POLY g;
int needSyz;
POLY *c; /* set */
POLY *h; /* set */
/* f must be reducible by g. r = c*f + h*g */
{
extern struct ring *CurrentRingp;
struct ring *rp;
struct spValue sv;
POLY f2;
if (needSyz) {
if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
*c = cxx(1,0,0,rp);
*h = ZERO;
}
sv = (*sp)(fs,g);
f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
if (DebugReductionRed) {
printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
}
f = f2;
if (needSyz) {
*c = ppMult(sv.a,*c);
*h = ppAdd(ppMult(sv.a,*h),sv.b);
}
while ((fs = (*isCdrReducible)(f,g)) != ZERO) {
sv = (*sp)(fs,g);
f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
if (DebugReductionRed) {
printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
}
f = f2;
if (needSyz) {
*c = ppMult(sv.a,*c);
*h = ppAdd(ppMult(sv.a,*h),sv.b);
}
}
return(f);
}
/* for debug */
int isOrdered(f)
POLY f;
{ POLY g;
if (f ISZERO) return(1);
g = f->next;
while (g != POLYNULL) {
if ((*mmLarger)(g,f)>=1) return(0);
f = g;
g = f->next;
}
return(1);
}
POLY reduction_gen(f,gset,needSyz,syzp)
POLY f;
struct gradedPolySet *gset;
int needSyz;
struct syz0 *syzp; /* set */
{
int reduced,reduced1,reduced2;
int grd;
struct polySet *set;
POLY cf,syz;
int i;
POLY cc,cg;
extern struct ring *CurrentRingp;
struct ring *rp;
if (needSyz) {
if (f ISZERO) { rp = CurrentRingp; } else { rp = f->m->ringp; }
cf = cxx(1,0,0,rp);
syz = ZERO;
}
reduced = 0; /* no */
do {
reduced1 = 0; /* no */
grd = 0;
while (grd < gset->maxGrade) {
if (!Sugar) {
if (grd > (*grade)(f)) break;
}
set = gset->polys[grd];
do {
reduced2 = 0; /* no */
for (i=0; i<set->size; i++) {
if (f ISZERO) goto ss;
if ((*isReducible)(f,set->g[i])) {
f = (*reduction1)(f,set->g[i],needSyz,&cc,&cg);
if (needSyz) {
cf = ppMult(cc,cf);
syz = cpMult(toSyzCoeff(cc),syz);
syz = ppAddv(syz,toSyzPoly(cg,grd,i));
}
reduced = reduced1 = reduced2 = 1; /* yes */
}
}
} while (reduced2 != 0);
grd++;
}
}while (reduced1 != 0);
ss: ;
if (needSyz) {
syzp->cf = cf; /* cf is in the CurrentRingp */
syzp->syz = syz; /* syz is in the SyzRingp */
}
return(f);
}
POLY reduction_gen_rev(f,gset,needSyz,syzp)
POLY f;
struct gradedPolySet *gset;
int needSyz;
struct syz0 *syzp; /* set */
{
int reduced,reduced1,reduced2;
int grd;
struct polySet *set;
POLY cf,syz;
int i;
POLY cc,cg;
extern struct ring *CurrentRingp;
struct ring *rp;
if (needSyz) {
if (f ISZERO) { rp = CurrentRingp; } else { rp = f->m->ringp; }
cf = cxx(1,0,0,rp);
syz = ZERO;
}
reduced = 0; /* no */
do {
reduced1 = 0; /* no */
grd = ((*grade)(f) < gset->maxGrade-1?(*grade)(f): gset->maxGrade-1);
while (grd >= 0) { /* reverse order for reduction */
set = gset->polys[grd];
do {
reduced2 = 0; /* no */
for (i=0; i<set->size; i++) {
if (f ISZERO) goto ss;
if ((*isReducible)(f,set->g[i])) {
f = (*reduction1)(f,set->g[i],needSyz,&cc,&cg);
if (needSyz) {
cf = ppMult(cc,cf);
syz = cpMult(toSyzCoeff(cc),syz);
syz = ppAddv(syz,toSyzPoly(cg,grd,i));
}
reduced = reduced1 = reduced2 = 1; /* yes */
}
}
} while (reduced2 != 0);
grd--;
}
}while (reduced1 != 0);
ss: ;
if (needSyz) {
syzp->cf = cf; /* cf is in the CurrentRingp */
syzp->syz = syz; /* syz is in the SyzRingp */
}
return(f);
}
POLY reductionCdr_gen(f,gset,needSyz,syzp)
POLY f;
struct gradedPolySet *gset;
int needSyz;
struct syz0 *syzp; /* set */
{
int reduced,reduced1,reduced2;
int grd;
struct polySet *set;
POLY cf,syz;
int i;
POLY cc,cg;
POLY fs;
extern struct ring *CurrentRingp;
struct ring *rp;
if (needSyz) {
if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
cf = cxx(1,0,0,rp);
syz = ZERO;
}
reduced = 0; /* no */
do {
reduced1 = 0; /* no */
grd = 0;
while (grd < gset->maxGrade) {
if (!Sugar) {
if (grd > (*grade)(f)) break;
}
set = gset->polys[grd];
do {
reduced2 = 0; /* no */
for (i=0; i<set->size; i++) {
if (f ISZERO) goto ss;
if ((fs =(*isCdrReducible)(f,set->g[i])) != ZERO) {
f = (*reduction1Cdr)(f,fs,set->g[i],needSyz,&cc,&cg);
if (needSyz) {
cf = ppMult(cc,cf);
syz = cpMult(toSyzCoeff(cc),syz);
syz = ppAddv(syz,toSyzPoly(cg,grd,i));
}
reduced = reduced1 = reduced2 = 1; /* yes */
}
}
} while (reduced2 != 0);
grd++;
}
}while (reduced1 != 0);
ss: ;
if (needSyz) {
syzp->cf = cf;
syzp->syz = syz;
}
return(f);
}
int isReducible_gen(f,g)
POLY f;
POLY g;
{
int n,i;
MONOMIAL tf;
MONOMIAL tg;
if (f ISZERO) return(0);
if (g ISZERO) return(0);
checkRingIsR(f,g);
tf = f->m; tg = g->m;
n = tf->ringp->n;
for (i=0; i<n; i++) {
if (tf->e[i].x < tg->e[i].x) return(0);
if (tf->e[i].D < tg->e[i].D) return(0);
}
return(1);
}
POLY isCdrReducible_gen(f,g)
POLY f;
POLY g;
{
while (f != POLYNULL) {
if ((*isReducible)(f,g)) {
return(f);
}
f = f->next;
}
return(ZERO);
}
POLY lcm_gen(f,g)
POLY f;
POLY g;
{
MONOMIAL tf,tg;
MONOMIAL lcm;
int n;
int i;
tf = f->m; tg = g->m;
n = tf->ringp->n;
lcm = newMonomial(tf->ringp);
for (i=0; i<n; i++) {
lcm->e[i].x = mymax(tf->e[i].x,tg->e[i].x);
lcm->e[i].D = mymax(tf->e[i].D,tg->e[i].D);
}
return(newCell(intToCoeff(1,tf->ringp),lcm));
}
int grade_gen(f)
POLY f;
{
int r;
int i,n;
MONOMIAL tf;
if (f ISZERO) return(-1);
tf = f->m;
n = tf->ringp->n;
r = 0;
for (i=0; i<n; i++) {
r += tf->e[i].x;
r += tf->e[i].D;
}
return(r);
}
/* constructors */
POLY toSyzPoly(cg,grd,index)
POLY cg;
int grd;
int index;
/* the result is read only. */
{
extern struct ring *SyzRingp;
POLY r;
r = newCell(toSyzCoeff(cg),newMonomial(SyzRingp));
r->m->e[0].x = grd;
r->m->e[0].D = index;
return(r);
}
struct coeff *toSyzCoeff(f)
POLY f;
{
extern struct ring *SyzRingp;
struct coeff *c;
c = newCoeff();
c->tag = POLY_COEFF;
c->val.f = f;
c->p = SyzRingp->p;
return(c);
}
void initSyzRingp() {
extern struct ring *SyzRingp;
extern struct ring *CurrentRingp;
static char *x[]={"grade"};
static char *d[]={"index"};
static int order[]={1,0,
0,1};
static int outputOrderForSyzRing[] = {0,1};
static int ringSerial = 0;
char *ringName = NULL;
SyzRingp = (struct ring *)sGC_malloc(sizeof(struct ring));
if (SyzRingp == (struct ring *)NULL)
errorGradedSet("initSyzRingp(); No more memory");
SyzRingp->p = CurrentRingp->p;
SyzRingp->n = 1; SyzRingp->m = 1; SyzRingp->l = 1; SyzRingp->c = 1;
SyzRingp->nn = 1; SyzRingp->mm = 1; SyzRingp->ll = 1;
SyzRingp->cc = 1;
SyzRingp->x = x;
SyzRingp->D = d;
SyzRingp->order = order;
SyzRingp->orderMatrixSize = 2;
setFromTo(SyzRingp);
SyzRingp->next = CurrentRingp;
SyzRingp->multiplication = mpMult_poly; /* Multiplication is not used.*/
SyzRingp->schreyer = 0;
SyzRingp->gbListTower = NULL;
SyzRingp->outputOrder = outputOrderForSyzRing;
ringName = (char *)sGC_malloc(20);
if (ringName == NULL) errorGradedSet("No more memory.");
sprintf(ringName,"syzring%05d",ringSerial);
SyzRingp->name = ringName;
}
POLY reductionCdr_except_grd_i(POLY f,struct gradedPolySet *gset,
int needSyz,struct syz0 *syzp,
int skipGrd,int skipi, int *reducedp)
{
int reduced,reduced1,reduced2;
int grd;
struct polySet *set;
POLY cf,syz;
int i;
POLY cc,cg;
POLY fs;
extern struct ring *CurrentRingp;
struct ring *rp;
*reducedp = 0;
if (needSyz) {
if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
cf = cxx(1,0,0,rp);
syz = ZERO;
}
reduced = 0; /* no */
do {
reduced1 = 0; /* no */
grd = 0;
while (grd < gset->maxGrade) {
/*
if (!Sugar) {
if (grd > (*grade)(f)) break;
}
*/
set = gset->polys[grd];
do {
reduced2 = 0; /* no */
for (i=0; i<set->size; i++) {
if (f ISZERO) goto ss;
if ((!((grd == skipGrd) && (i == skipi))) && (set->del[i]==0)) {
/* Do not use deleted element.*/
if ((fs =(*isCdrReducible)(f,set->g[i])) != ZERO) {
f = (*reduction1Cdr)(f,fs,set->g[i],needSyz,&cc,&cg);
/* What is cg? */
if (needSyz) {
cf = ppMult(cc,cf);
syz = cpMult(toSyzCoeff(cc),syz);
syz = ppAddv(syz,toSyzPoly(cg,grd,i));
}
*reducedp = reduced = reduced1 = reduced2 = 1; /* yes */
}
}
}
} while (reduced2 != 0);
grd++;
}
}while (reduced1 != 0);
ss: ;
if (needSyz) {
syzp->cf = cf;
syzp->syz = syz;
}
return(f);
}