=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/dalg.c,v retrieving revision 1.17 retrieving revision 1.18 diff -u -p -r1.17 -r1.18 --- OpenXM_contrib2/asir2000/engine/dalg.c 2017/08/31 02:36:21 1.17 +++ OpenXM_contrib2/asir2000/engine/dalg.c 2018/03/29 01:32:51 1.18 @@ -1,5 +1,5 @@ /* - * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.16 2013/11/05 02:55:03 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.17 2017/08/31 02:36:21 noro Exp $ */ #include "ca.h" @@ -15,979 +15,979 @@ void dalgtoalg(DAlg da,Alg *r); NumberField get_numberfield() { - return current_numberfield; + return current_numberfield; } void setfield_dalg(NODE alist) { - NumberField nf; - VL vl,vl1,vl2; - int n,i,dim; - Alg *gen; - P *defpoly; - P p; - Q c,iq,two; - DP *ps,*mb; - DP one; - NODE t,b,b1,b2,hlist,mblist; - struct order_spec *current_spec; + NumberField nf; + VL vl,vl1,vl2; + int n,i,dim; + Alg *gen; + P *defpoly; + P p; + Q c,iq,two; + DP *ps,*mb; + DP one; + NODE t,b,b1,b2,hlist,mblist; + struct order_spec *current_spec; - nf = (NumberField)MALLOC(sizeof(struct oNumberField)); - current_numberfield = nf; - vl = 0; - for ( t = alist; t; t = NEXT(t) ) { - clctalg((P)BDY((Alg)BDY(t)),&vl1); - mergev(ALG,vl,vl1,&vl2); vl = vl2; - } - for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ ); - nf->n = n; - nf->vl = vl; - nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P)); - nf->ps = ps = (DP *)MALLOC(n*sizeof(DP)); - current_spec = dp_current_spec; - STOQ(2,two); - create_order_spec(0,(Obj)two,&nf->spec); - initd(nf->spec); - for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) { - ptozp(vl1->v->attr,1,&c,&defpoly[i]); - ptod(ALG,vl,defpoly[i],&ps[i]); - STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1; - MKNODE(b2,(pointer)ps[i],hlist); hlist = b2; - } - ptod(ALG,vl,(P)ONE,&one); - MKDAlg(one,ONE,nf->one); - nf->ind = b; - dp_mbase(hlist,&mblist); - initd(current_spec); - nf->dim = dim = length(mblist); - nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP)); - for ( i = 0, t = mblist; t; t = NEXT(t), i++ ) - mb[dim-i-1] = (DP)BDY(t); + nf = (NumberField)MALLOC(sizeof(struct oNumberField)); + current_numberfield = nf; + vl = 0; + for ( t = alist; t; t = NEXT(t) ) { + clctalg((P)BDY((Alg)BDY(t)),&vl1); + mergev(ALG,vl,vl1,&vl2); vl = vl2; + } + for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ ); + nf->n = n; + nf->vl = vl; + nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P)); + nf->ps = ps = (DP *)MALLOC(n*sizeof(DP)); + current_spec = dp_current_spec; + STOQ(2,two); + create_order_spec(0,(Obj)two,&nf->spec); + initd(nf->spec); + for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) { + ptozp(vl1->v->attr,1,&c,&defpoly[i]); + ptod(ALG,vl,defpoly[i],&ps[i]); + STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1; + MKNODE(b2,(pointer)ps[i],hlist); hlist = b2; + } + ptod(ALG,vl,(P)ONE,&one); + MKDAlg(one,ONE,nf->one); + nf->ind = b; + dp_mbase(hlist,&mblist); + initd(current_spec); + nf->dim = dim = length(mblist); + nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP)); + for ( i = 0, t = mblist; t; t = NEXT(t), i++ ) + mb[dim-i-1] = (DP)BDY(t); } void setfield_gb(NODE gb,VL vl,struct order_spec *spec) { - NumberField nf; - VL vl1,vl2; - int n,i,dim; - Alg *gen; - P *defpoly; - P p; - Q c,iq,two; - DP *ps,*mb; - DP one; - NODE t,b,b1,b2,hlist,mblist; - struct order_spec *current_spec; + NumberField nf; + VL vl1,vl2; + int n,i,dim; + Alg *gen; + P *defpoly; + P p; + Q c,iq,two; + DP *ps,*mb; + DP one; + NODE t,b,b1,b2,hlist,mblist; + struct order_spec *current_spec; - nf = (NumberField)MALLOC(sizeof(struct oNumberField)); - current_numberfield = nf; - for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ ); - nf->n = n; - nf->psn = length(gb); - nf->vl = vl; - nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P)); - nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP)); - current_spec = dp_current_spec; - nf->spec = spec; - initd(nf->spec); - for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) { - ptozp((P)BDY(t),1,&c,&defpoly[i]); - ptod(CO,vl,defpoly[i],&ps[i]); - STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1; - MKNODE(b2,(pointer)ps[i],hlist); hlist = b2; - } - ptod(ALG,vl,(P)ONE,&one); - MKDAlg(one,ONE,nf->one); - nf->ind = b; - dp_mbase(hlist,&mblist); - initd(current_spec); - nf->dim = dim = length(mblist); - nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP)); - for ( i = 0, t = mblist; t; t = NEXT(t), i++ ) - mb[dim-i-1] = (DP)BDY(t); + nf = (NumberField)MALLOC(sizeof(struct oNumberField)); + current_numberfield = nf; + for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ ); + nf->n = n; + nf->psn = length(gb); + nf->vl = vl; + nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P)); + nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP)); + current_spec = dp_current_spec; + nf->spec = spec; + initd(nf->spec); + for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) { + ptozp((P)BDY(t),1,&c,&defpoly[i]); + ptod(CO,vl,defpoly[i],&ps[i]); + STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1; + MKNODE(b2,(pointer)ps[i],hlist); hlist = b2; + } + ptod(ALG,vl,(P)ONE,&one); + MKDAlg(one,ONE,nf->one); + nf->ind = b; + dp_mbase(hlist,&mblist); + initd(current_spec); + nf->dim = dim = length(mblist); + nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP)); + for ( i = 0, t = mblist; t; t = NEXT(t), i++ ) + mb[dim-i-1] = (DP)BDY(t); } void qtodalg(Q q,DAlg *r) { - NumberField nf; - Q t; - DP nm; + NumberField nf; + Q t; + DP nm; - if ( !(nf=current_numberfield) ) - error("qtodalg : current_numberfield is not set"); - if ( !q ) - *r = 0; - else if ( NID(q) == N_DA ) - *r = (DAlg)q; - else if ( NID(q) == N_Q ) { - if ( INT(q) ) { - muldc(CO,nf->one->nm,(Obj)q,&nm); - MKDAlg(nm,ONE,*r); - } else { - NTOQ(NM(q),SGN(q),t); - muldc(CO,nf->one->nm,(Obj)t,&nm); - NTOQ(DN(q),1,t); - MKDAlg(nm,t,*r); - } - } else - error("qtodalg : invalid argument"); + if ( !(nf=current_numberfield) ) + error("qtodalg : current_numberfield is not set"); + if ( !q ) + *r = 0; + else if ( NID(q) == N_DA ) + *r = (DAlg)q; + else if ( NID(q) == N_Q ) { + if ( INT(q) ) { + muldc(CO,nf->one->nm,(Obj)q,&nm); + MKDAlg(nm,ONE,*r); + } else { + NTOQ(NM(q),SGN(q),t); + muldc(CO,nf->one->nm,(Obj)t,&nm); + NTOQ(DN(q),1,t); + MKDAlg(nm,t,*r); + } + } else + error("qtodalg : invalid argument"); } void obj_algtodalg(Obj obj,Obj *r) { - DAlg d; - DCP dc,dcr0,dcr; - P c,p; - Obj t; - Obj nm,dn; - NODE b,s,s0; - R rat; - VECT v; - MAT mat; - LIST list; - pointer *a; - pointer **m; - int len,row,col,i,j,l; + DAlg d; + DCP dc,dcr0,dcr; + P c,p; + Obj t; + Obj nm,dn; + NODE b,s,s0; + R rat; + VECT v; + MAT mat; + LIST list; + pointer *a; + pointer **m; + int len,row,col,i,j,l; - if ( !obj ) { - *r = 0; - return; - } - switch ( OID(obj) ) { - case O_N: - algtodalg((Alg)obj,&d); *r = (Obj)d; - break; - case O_P: - for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) { - obj_algtodalg((Obj)COEF(dc),&t); - if ( t ) { - NEXTDC(dcr0,dcr); - COEF(dcr) = (P)t; - DEG(dcr) = DEG(dc); - } - } - if ( dcr0 ) { - MKP(VR((P)obj),dcr0,p); - *r = (Obj)p; - } else - *r = 0; - break; - case O_R: - obj_algtodalg((Obj)NM((R)obj),&nm); - obj_algtodalg((Obj)DN((R)obj),&dn); - if ( !dn ) - error("obj_algtodalg : division by 0"); - if ( !nm ) - *r = 0; - else { - MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat; - } - break; - case O_LIST: - s0 = 0; - for ( b = BDY((LIST)obj); b; b = NEXT(b) ) { - NEXTNODE(s0,s); - obj_algtodalg((Obj)BDY(b),&t); - BDY(s) = (pointer)t; - } - NEXT(s) = 0; - MKLIST(list,s0); - *r = (Obj)list; - break; - case O_VECT: - l = ((VECT)obj)->len; - a = BDY((VECT)obj); - MKVECT(v,l); - for ( i = 0; i < l; i++ ) { - obj_algtodalg((Obj)a[i],&t); - BDY(v)[i] = (pointer)t; - } - *r = (Obj)v; - break; - case O_MAT: - row = ((MAT)obj)->row; col = ((MAT)obj)->col; - m = BDY((MAT)obj); - MKMAT(mat,row,col); - for ( i = 0; i < row; i++ ) - for ( j = 0; j < col; j++ ) { - obj_algtodalg((Obj)m[i][j],&t); - BDY(mat)[i][j] = (pointer)t; - } - *r = (Obj)mat; - break; - default: - *r = obj; - break; - } + if ( !obj ) { + *r = 0; + return; + } + switch ( OID(obj) ) { + case O_N: + algtodalg((Alg)obj,&d); *r = (Obj)d; + break; + case O_P: + for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) { + obj_algtodalg((Obj)COEF(dc),&t); + if ( t ) { + NEXTDC(dcr0,dcr); + COEF(dcr) = (P)t; + DEG(dcr) = DEG(dc); + } + } + if ( dcr0 ) { + MKP(VR((P)obj),dcr0,p); + *r = (Obj)p; + } else + *r = 0; + break; + case O_R: + obj_algtodalg((Obj)NM((R)obj),&nm); + obj_algtodalg((Obj)DN((R)obj),&dn); + if ( !dn ) + error("obj_algtodalg : division by 0"); + if ( !nm ) + *r = 0; + else { + MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat; + } + break; + case O_LIST: + s0 = 0; + for ( b = BDY((LIST)obj); b; b = NEXT(b) ) { + NEXTNODE(s0,s); + obj_algtodalg((Obj)BDY(b),&t); + BDY(s) = (pointer)t; + } + NEXT(s) = 0; + MKLIST(list,s0); + *r = (Obj)list; + break; + case O_VECT: + l = ((VECT)obj)->len; + a = BDY((VECT)obj); + MKVECT(v,l); + for ( i = 0; i < l; i++ ) { + obj_algtodalg((Obj)a[i],&t); + BDY(v)[i] = (pointer)t; + } + *r = (Obj)v; + break; + case O_MAT: + row = ((MAT)obj)->row; col = ((MAT)obj)->col; + m = BDY((MAT)obj); + MKMAT(mat,row,col); + for ( i = 0; i < row; i++ ) + for ( j = 0; j < col; j++ ) { + obj_algtodalg((Obj)m[i][j],&t); + BDY(mat)[i][j] = (pointer)t; + } + *r = (Obj)mat; + break; + default: + *r = obj; + break; + } } void obj_dalgtoalg(Obj obj,Obj *r) { - Alg d; - DCP dc,dcr0,dcr; - P c,p; - Obj t; - Obj nm,dn; - NODE b,s,s0; - R rat; - VECT v; - MAT mat; - LIST list; - pointer *a; - pointer **m; - int len,row,col,i,j,l; + Alg d; + DCP dc,dcr0,dcr; + P c,p; + Obj t; + Obj nm,dn; + NODE b,s,s0; + R rat; + VECT v; + MAT mat; + LIST list; + pointer *a; + pointer **m; + int len,row,col,i,j,l; - if ( !obj ) { - *r = 0; - return; - } - switch ( OID(obj) ) { - case O_N: - dalgtoalg((DAlg)obj,&d); *r = (Obj)d; - break; - case O_P: - for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) { - obj_dalgtoalg((Obj)COEF(dc),&t); - if ( t ) { - NEXTDC(dcr0,dcr); - COEF(dcr) = (P)t; - DEG(dcr) = DEG(dc); - } - } - if ( dcr0 ) { - MKP(VR((P)obj),dcr0,p); - *r = (Obj)p; - } else - *r = 0; - break; - case O_R: - obj_dalgtoalg((Obj)NM((R)obj),&nm); - obj_dalgtoalg((Obj)DN((R)obj),&dn); - if ( !dn ) - error("obj_dalgtoalg : division by 0"); - if ( !nm ) - *r = 0; - else { - MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat; - } - break; - case O_LIST: - s0 = 0; - for ( b = BDY((LIST)obj); b; b = NEXT(b) ) { - NEXTNODE(s0,s); - obj_dalgtoalg((Obj)BDY(b),&t); - BDY(s) = (pointer)t; - } - NEXT(s) = 0; - MKLIST(list,s0); - *r = (Obj)list; - break; - case O_VECT: - l = ((VECT)obj)->len; - a = BDY((VECT)obj); - MKVECT(v,l); - for ( i = 0; i < l; i++ ) { - obj_dalgtoalg((Obj)a[i],&t); - BDY(v)[i] = (pointer)t; - } - *r = (Obj)v; - break; - case O_MAT: - row = ((MAT)obj)->row; col = ((MAT)obj)->col; - m = BDY((MAT)obj); - MKMAT(mat,row,col); - for ( i = 0; i < row; i++ ) - for ( j = 0; j < col; j++ ) { - obj_dalgtoalg((Obj)m[i][j],&t); - BDY(mat)[i][j] = (pointer)t; - } - *r = (Obj)mat; - break; - default: - *r = obj; - break; - } + if ( !obj ) { + *r = 0; + return; + } + switch ( OID(obj) ) { + case O_N: + dalgtoalg((DAlg)obj,&d); *r = (Obj)d; + break; + case O_P: + for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) { + obj_dalgtoalg((Obj)COEF(dc),&t); + if ( t ) { + NEXTDC(dcr0,dcr); + COEF(dcr) = (P)t; + DEG(dcr) = DEG(dc); + } + } + if ( dcr0 ) { + MKP(VR((P)obj),dcr0,p); + *r = (Obj)p; + } else + *r = 0; + break; + case O_R: + obj_dalgtoalg((Obj)NM((R)obj),&nm); + obj_dalgtoalg((Obj)DN((R)obj),&dn); + if ( !dn ) + error("obj_dalgtoalg : division by 0"); + if ( !nm ) + *r = 0; + else { + MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat; + } + break; + case O_LIST: + s0 = 0; + for ( b = BDY((LIST)obj); b; b = NEXT(b) ) { + NEXTNODE(s0,s); + obj_dalgtoalg((Obj)BDY(b),&t); + BDY(s) = (pointer)t; + } + NEXT(s) = 0; + MKLIST(list,s0); + *r = (Obj)list; + break; + case O_VECT: + l = ((VECT)obj)->len; + a = BDY((VECT)obj); + MKVECT(v,l); + for ( i = 0; i < l; i++ ) { + obj_dalgtoalg((Obj)a[i],&t); + BDY(v)[i] = (pointer)t; + } + *r = (Obj)v; + break; + case O_MAT: + row = ((MAT)obj)->row; col = ((MAT)obj)->col; + m = BDY((MAT)obj); + MKMAT(mat,row,col); + for ( i = 0; i < row; i++ ) + for ( j = 0; j < col; j++ ) { + obj_dalgtoalg((Obj)m[i][j],&t); + BDY(mat)[i][j] = (pointer)t; + } + *r = (Obj)mat; + break; + default: + *r = obj; + break; + } } void algtodalg(Alg a,DAlg *r) { - P ap,p,p1; - Q c,c1,d1,dn,nm; - DP dp; - DAlg da; - NumberField nf; - struct order_spec *current_spec; - VL vl,tvl,svl; - V v; + P ap,p,p1; + Q c,c1,d1,dn,nm; + DP dp; + DAlg da; + NumberField nf; + struct order_spec *current_spec; + VL vl,tvl,svl; + V v; - if ( !(nf=current_numberfield) ) - error("algtodalg : current_numberfield is not set"); - if ( !a ) { - *r = 0; - return; - } - switch (NID((Num)a) ) { - case N_Q: - c = (Q)a; - if ( INT(c) ) { - muldc(CO,nf->one->nm,(Obj)c,&dp); - MKDAlg(dp,ONE,*r); - } else { - NTOQ(NM(c),SGN(c),c1); - NTOQ(DN(c),1,d1); - muldc(CO,nf->one->nm,(Obj)c1,&dp); - MKDAlg(dp,d1,*r); - } - break; - case N_A: - ap = (P)BDY(a); - ptozp(ap,1,&c,&p); - if ( INT(c) ) { - p = ap; - dn = ONE; - } else { - NTOQ(NM(c),SGN(c),nm); - NTOQ(DN(c),1,dn); - mulpq(p,(P)nm,&p1); p = p1; - } - current_spec = dp_current_spec; initd(nf->spec); - get_vars((Obj)p,&vl); - for ( tvl = vl; tvl; tvl = NEXT(tvl) ) { - v = tvl->v; - for ( svl = nf->vl; svl; svl = NEXT(svl) ) - if ( v == svl->v ) - break; - if ( !svl ) - error("algtodalg : incompatible numberfield"); - } - ptod(ALG,nf->vl,p,&dp); - MKDAlg(dp,dn,da); - simpdalg(da,r); - break; - default: - error("algtodalg : invalid argument"); - break; - } + if ( !(nf=current_numberfield) ) + error("algtodalg : current_numberfield is not set"); + if ( !a ) { + *r = 0; + return; + } + switch (NID((Num)a) ) { + case N_Q: + c = (Q)a; + if ( INT(c) ) { + muldc(CO,nf->one->nm,(Obj)c,&dp); + MKDAlg(dp,ONE,*r); + } else { + NTOQ(NM(c),SGN(c),c1); + NTOQ(DN(c),1,d1); + muldc(CO,nf->one->nm,(Obj)c1,&dp); + MKDAlg(dp,d1,*r); + } + break; + case N_A: + ap = (P)BDY(a); + ptozp(ap,1,&c,&p); + if ( INT(c) ) { + p = ap; + dn = ONE; + } else { + NTOQ(NM(c),SGN(c),nm); + NTOQ(DN(c),1,dn); + mulpq(p,(P)nm,&p1); p = p1; + } + current_spec = dp_current_spec; initd(nf->spec); + get_vars((Obj)p,&vl); + for ( tvl = vl; tvl; tvl = NEXT(tvl) ) { + v = tvl->v; + for ( svl = nf->vl; svl; svl = NEXT(svl) ) + if ( v == svl->v ) + break; + if ( !svl ) + error("algtodalg : incompatible numberfield"); + } + ptod(ALG,nf->vl,p,&dp); + MKDAlg(dp,dn,da); + simpdalg(da,r); + break; + default: + error("algtodalg : invalid argument"); + break; + } } void dalgtoalg(DAlg da,Alg *r) { - NumberField nf; - P p,p1; - Q inv; + NumberField nf; + P p,p1; + Q inv; - if ( !(nf=current_numberfield) ) - error("dalgtoalg : current_numberfield is not set"); - if ( !da ) *r = 0; - else { - dtop(ALG,nf->vl,da->nm,(Obj *)&p); - invq(da->dn,&inv); - mulpq(p,(P)inv,&p1); - MKAlg(p1,*r); - } + if ( !(nf=current_numberfield) ) + error("dalgtoalg : current_numberfield is not set"); + if ( !da ) *r = 0; + else { + dtop(ALG,nf->vl,da->nm,(Obj *)&p); + invq(da->dn,&inv); + mulpq(p,(P)inv,&p1); + MKAlg(p1,*r); + } } void simpdalg(DAlg da,DAlg *r) { - NumberField nf; - DP nm; - DAlg d; - Q dn,dn1; - struct order_spec *current_spec; + NumberField nf; + DP nm; + DAlg d; + Q dn,dn1; + struct order_spec *current_spec; - if ( !(nf=current_numberfield) ) - error("simpdalg : current_numberfield is not set"); - if ( !da ) { - *r = 0; - return; - } - current_spec = dp_current_spec; initd(nf->spec); - dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,(P *)&dn); - if ( !nm ) *r = 0; - else { - initd(current_spec); - mulq(da->dn,dn,&dn1); - MKDAlg(nm,dn1,d); - rmcontdalg(d,r); - } + if ( !(nf=current_numberfield) ) + error("simpdalg : current_numberfield is not set"); + if ( !da ) { + *r = 0; + return; + } + current_spec = dp_current_spec; initd(nf->spec); + dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,(P *)&dn); + if ( !nm ) *r = 0; + else { + initd(current_spec); + mulq(da->dn,dn,&dn1); + MKDAlg(nm,dn1,d); + rmcontdalg(d,r); + } } void adddalg(DAlg a,DAlg b,DAlg *c) { - NumberField nf; - Q dna,dnb,a1,b1,dn,g; - N an,bn,gn; - DAlg t; - DP ta,tb,nm; - struct order_spec *current_spec; + NumberField nf; + Q dna,dnb,a1,b1,dn,g; + N an,bn,gn; + DAlg t; + DP ta,tb,nm; + struct order_spec *current_spec; - if ( !(nf=current_numberfield) ) - error("adddalg : current_numberfield is not set"); - if ( !a ) - *c = b; - else if ( !b ) - *c = a; - else { - qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; - dna = a->dn; - dnb = b->dn; - gcdn(NM(dna),NM(dnb),&gn); - divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn); - NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1); - /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */ - muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb); - current_spec = dp_current_spec; initd(nf->spec); - addd(CO,ta,tb,&nm); - initd(current_spec); - if ( !nm ) - *c = 0; - else { - mulq(dna,b1,&dn); - MKDAlg(nm,dn,*c); - } - } + if ( !(nf=current_numberfield) ) + error("adddalg : current_numberfield is not set"); + if ( !a ) + *c = b; + else if ( !b ) + *c = a; + else { + qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; + dna = a->dn; + dnb = b->dn; + gcdn(NM(dna),NM(dnb),&gn); + divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn); + NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1); + /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */ + muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb); + current_spec = dp_current_spec; initd(nf->spec); + addd(CO,ta,tb,&nm); + initd(current_spec); + if ( !nm ) + *c = 0; + else { + mulq(dna,b1,&dn); + MKDAlg(nm,dn,*c); + } + } } void subdalg(DAlg a,DAlg b,DAlg *c) { - NumberField nf; - Q dna,dnb,a1,b1,dn,g; - N an,bn,gn; - DP ta,tb,nm; - DAlg t; - struct order_spec *current_spec; + NumberField nf; + Q dna,dnb,a1,b1,dn,g; + N an,bn,gn; + DP ta,tb,nm; + DAlg t; + struct order_spec *current_spec; - if ( !(nf=current_numberfield) ) - error("subdalg : current_numberfield is not set"); - if ( !a ) - *c = b; - else if ( !b ) - *c = a; - else { - qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; - dna = a->dn; - dnb = b->dn; - gcdn(NM(dna),NM(dnb),&gn); - divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn); - NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1); - /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */ - muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb); - current_spec = dp_current_spec; initd(nf->spec); - subd(CO,ta,tb,&nm); - initd(current_spec); - if ( !nm ) - *c = 0; - else { - mulq(dna,b1,&dn); - MKDAlg(nm,dn,*c); - } - } + if ( !(nf=current_numberfield) ) + error("subdalg : current_numberfield is not set"); + if ( !a ) + *c = b; + else if ( !b ) + *c = a; + else { + qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; + dna = a->dn; + dnb = b->dn; + gcdn(NM(dna),NM(dnb),&gn); + divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn); + NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1); + /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */ + muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb); + current_spec = dp_current_spec; initd(nf->spec); + subd(CO,ta,tb,&nm); + initd(current_spec); + if ( !nm ) + *c = 0; + else { + mulq(dna,b1,&dn); + MKDAlg(nm,dn,*c); + } + } } void muldalg(DAlg a,DAlg b,DAlg *c) { - NumberField nf; - DP nm; - Q dn; - DAlg t; - struct order_spec *current_spec; + NumberField nf; + DP nm; + Q dn; + DAlg t; + struct order_spec *current_spec; - if ( !(nf=current_numberfield) ) - error("muldalg : current_numberfield is not set"); - if ( !a || !b ) - *c = 0; - else { - qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; - current_spec = dp_current_spec; initd(nf->spec); - muld(CO,a->nm,b->nm,&nm); - initd(current_spec); - mulq(a->dn,b->dn,&dn); - MKDAlg(nm,dn,t); - simpdalg(t,c); - } + if ( !(nf=current_numberfield) ) + error("muldalg : current_numberfield is not set"); + if ( !a || !b ) + *c = 0; + else { + qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; + current_spec = dp_current_spec; initd(nf->spec); + muld(CO,a->nm,b->nm,&nm); + initd(current_spec); + mulq(a->dn,b->dn,&dn); + MKDAlg(nm,dn,t); + simpdalg(t,c); + } } void divdalg(DAlg a,DAlg b,DAlg *c) { - DAlg inv,t; - int ret; + DAlg inv,t; + int ret; - if ( !current_numberfield ) - error("divdalg : current_numberfield is not set"); - if ( !b ) - error("divdalg : division by 0"); - if ( !a ) - c = 0; - else { - qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; - ret = invdalg(b,&inv); - if ( !ret ) { - error("divdalg : the denominator is not invertible"); - } - muldalg(a,inv,c); - } + if ( !current_numberfield ) + error("divdalg : current_numberfield is not set"); + if ( !b ) + error("divdalg : division by 0"); + if ( !a ) + c = 0; + else { + qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t; + ret = invdalg(b,&inv); + if ( !ret ) { + error("divdalg : the denominator is not invertible"); + } + muldalg(a,inv,c); + } } void rmcontdalg(DAlg a, DAlg *r) { - DP u,u1; - Q cont,c,d; - N gn,cn,dn; + DP u,u1; + Q cont,c,d; + N gn,cn,dn; - if ( !a ) - *r = a; - else { - dp_ptozp(a->nm,&u); - divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont); - gcdn(NM(cont),NM(a->dn),&gn); - divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c); - divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d); - muldc(CO,u,(Obj)c,&u1); - MKDAlg(u1,d,*r); - } + if ( !a ) + *r = a; + else { + dp_ptozp(a->nm,&u); + divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont); + gcdn(NM(cont),NM(a->dn),&gn); + divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c); + divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d); + muldc(CO,u,(Obj)c,&u1); + MKDAlg(u1,d,*r); + } } int invdalg(DAlg a,DAlg *c) { - NumberField nf; - int dim,n,i,j,k,l; - DP *mb; - DP m,d,u; - N ln,gn,qn; - DAlg *simp; - DAlg t,a0,r; - Q dn,dnsol,mul,nmc,dn1; - MAT mobj,sol; - Q **mat,**solmat; - MP mp0,mp; - int *rinfo,*cinfo; - int rank,nparam; - NODE nd0,nd,ndt; - struct order_spec *current_spec; - struct oEGT eg0,eg1; - extern struct oEGT eg_le; + NumberField nf; + int dim,n,i,j,k,l; + DP *mb; + DP m,d,u; + N ln,gn,qn; + DAlg *simp; + DAlg t,a0,r; + Q dn,dnsol,mul,nmc,dn1; + MAT mobj,sol; + Q **mat,**solmat; + MP mp0,mp; + int *rinfo,*cinfo; + int rank,nparam; + NODE nd0,nd,ndt; + struct order_spec *current_spec; + struct oEGT eg0,eg1; + extern struct oEGT eg_le; - if ( !(nf=current_numberfield) ) - error("invdalg : current_numberfield is not set"); - if ( !a ) - error("invdalg : division by 0"); - else if ( NID(a) == N_Q ) { - invq((Q)a,&dn); *c = (DAlg)dn; - return 1; - } - dim = nf->dim; - mb = nf->mb; - n = nf->n; - ln = ONEN; - dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc); - MKDAlg(u,ONE,a0); - simp = (DAlg *)ALLOCA(dim*sizeof(DAlg)); - current_spec = dp_current_spec; initd(nf->spec); - for ( i = 0; i < dim; i++ ) { - m = mb[i]; - for ( j = i-1; j >= 0; j-- ) - if ( dp_redble(m,mb[j]) ) - break; - if ( j >= 0 ) { - dp_subd(m,mb[j],&d); - muld(CO,d,simp[j]->nm,&u); - MKDAlg(u,simp[j]->dn,t); - simpdalg(t,&simp[i]); - } else { - MKDAlg(m,ONE,t); - muldalg(t,a0,&simp[i]); - } - gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn); - muln(NM(simp[i]->dn),qn,&ln); - } - initd(current_spec); - NTOQ(ln,1,dn); - MKMAT(mobj,dim,dim+1); - mat = (Q **)BDY(mobj); - mulq(dn,a->dn,&mat[0][dim]); - for ( j = 0; j < dim; j++ ) { - divq(dn,simp[j]->dn,&mul); - for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- ) - if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) { - mulq(mul,(Q)mp->c,&mat[i][j]); - mp = NEXT(mp); - } - } - get_eg(&eg0); - rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo); - get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1); - if ( cinfo[0] == dim ) { - /* the input is invertible */ - solmat = (Q **)BDY(sol); - for ( i = dim-1, mp0 = 0; i >= 0; i-- ) - if ( solmat[i][0] ) { - NEXTMP(mp0,mp); - mp->c = (Obj)solmat[i][0]; - mp->dl = BDY(mb[i])->dl; - } - NEXT(mp) = 0; MKDP(n,mp0,u); - mulq(dnsol,nmc,&dn1); - MKDAlg(u,dn1,r); - rmcontdalg(r,c); - return 1; - } else - return 0; + if ( !(nf=current_numberfield) ) + error("invdalg : current_numberfield is not set"); + if ( !a ) + error("invdalg : division by 0"); + else if ( NID(a) == N_Q ) { + invq((Q)a,&dn); *c = (DAlg)dn; + return 1; + } + dim = nf->dim; + mb = nf->mb; + n = nf->n; + ln = ONEN; + dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc); + MKDAlg(u,ONE,a0); + simp = (DAlg *)ALLOCA(dim*sizeof(DAlg)); + current_spec = dp_current_spec; initd(nf->spec); + for ( i = 0; i < dim; i++ ) { + m = mb[i]; + for ( j = i-1; j >= 0; j-- ) + if ( dp_redble(m,mb[j]) ) + break; + if ( j >= 0 ) { + dp_subd(m,mb[j],&d); + muld(CO,d,simp[j]->nm,&u); + MKDAlg(u,simp[j]->dn,t); + simpdalg(t,&simp[i]); + } else { + MKDAlg(m,ONE,t); + muldalg(t,a0,&simp[i]); + } + gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn); + muln(NM(simp[i]->dn),qn,&ln); + } + initd(current_spec); + NTOQ(ln,1,dn); + MKMAT(mobj,dim,dim+1); + mat = (Q **)BDY(mobj); + mulq(dn,a->dn,&mat[0][dim]); + for ( j = 0; j < dim; j++ ) { + divq(dn,simp[j]->dn,&mul); + for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- ) + if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) { + mulq(mul,(Q)mp->c,&mat[i][j]); + mp = NEXT(mp); + } + } + get_eg(&eg0); + rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo); + get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1); + if ( cinfo[0] == dim ) { + /* the input is invertible */ + solmat = (Q **)BDY(sol); + for ( i = dim-1, mp0 = 0; i >= 0; i-- ) + if ( solmat[i][0] ) { + NEXTMP(mp0,mp); + mp->c = (Obj)solmat[i][0]; + mp->dl = BDY(mb[i])->dl; + } + NEXT(mp) = 0; MKDP(n,mp0,u); + mulq(dnsol,nmc,&dn1); + MKDAlg(u,dn1,r); + rmcontdalg(r,c); + return 1; + } else + return 0; } NODE inv_or_split_dalg(DAlg a,DAlg *c) { - NumberField nf; - int dim,n,i,j,k,l; - DP *mb; - DP m,d,u; - N ln,gn,qn; - DAlg *simp; - DAlg t,a0,r; - Q dn,dnsol,mul,nmc,dn1; - MAT mobj,sol; - Q **mat,**solmat; - MP mp0,mp; - int *rinfo,*cinfo; - int rank,nparam; - NODE nd0,nd,ndt; - struct order_spec *current_spec; - struct oEGT eg0,eg1; - extern struct oEGT eg_le; - extern int DP_Print; + NumberField nf; + int dim,n,i,j,k,l; + DP *mb; + DP m,d,u; + N ln,gn,qn; + DAlg *simp; + DAlg t,a0,r; + Q dn,dnsol,mul,nmc,dn1; + MAT mobj,sol; + Q **mat,**solmat; + MP mp0,mp; + int *rinfo,*cinfo; + int rank,nparam; + NODE nd0,nd,ndt; + struct order_spec *current_spec; + struct oEGT eg0,eg1; + extern struct oEGT eg_le; + extern int DP_Print; - if ( !(nf=current_numberfield) ) - error("invdalg : current_numberfield is not set"); - if ( !a ) - error("invdalg : division by 0"); - else if ( NID(a) == N_Q ) { - invq((Q)a,&dn); *c = (DAlg)dn; - return 0; - } - dim = nf->dim; - mb = nf->mb; - n = nf->n; - ln = ONEN; - dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc); - MKDAlg(u,ONE,a0); - simp = (DAlg *)MALLOC(dim*sizeof(DAlg)); - current_spec = dp_current_spec; initd(nf->spec); - for ( i = 0; i < dim; i++ ) { - if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); } - m = mb[i]; - for ( j = i-1; j >= 0; j-- ) - if ( dp_redble(m,mb[j]) ) - break; - if ( j >= 0 ) { - dp_subd(m,mb[j],&d); - if ( simp[j] ) { - muld(CO,d,simp[j]->nm,&u); - MKDAlg(u,simp[j]->dn,t); - simpdalg(t,&simp[i]); - } else - simp[i] = 0; - } else { - MKDAlg(m,ONE,t); - muldalg(t,a0,&simp[i]); - } - if ( simp[i] ) { - gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn); - muln(NM(simp[i]->dn),qn,&ln); - } - } - initd(current_spec); - NTOQ(ln,1,dn); - MKMAT(mobj,dim,dim+1); - mat = (Q **)BDY(mobj); - mulq(dn,a->dn,&mat[0][dim]); - for ( j = 0; j < dim; j++ ) { - if ( simp[j] ) { - divq(dn,simp[j]->dn,&mul); - for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- ) - if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) { - mulq(mul,(Q)mp->c,&mat[i][j]); - mp = NEXT(mp); - } - } - } - get_eg(&eg0); - rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo); - get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1); - if ( cinfo[0] == dim ) { - /* the input is invertible */ - solmat = (Q **)BDY(sol); - for ( i = dim-1, mp0 = 0; i >= 0; i-- ) - if ( solmat[i][0] ) { - NEXTMP(mp0,mp); - mp->c = (Obj)solmat[i][0]; - mp->dl = BDY(mb[i])->dl; - } - NEXT(mp) = 0; MKDP(n,mp0,u); - mulq(dnsol,nmc,&dn1); - MKDAlg(u,dn1,r); - rmcontdalg(r,c); - return 0; - } else { - /* the input is not invertible */ - nparam = sol->col; - solmat = (Q **)BDY(sol); - nd0 = 0; - for ( k = 0; k < nparam; k++ ) { - /* construct a new basis element */ - m = mb[cinfo[k]]; - mp0 = 0; - NEXTMP(mp0,mp); - chsgnq(dnsol,&dn1); mp->c = (Obj)dn1; - mp->dl = BDY(m)->dl; - /* skip the last parameter */ - for ( l = rank-2; l >= 0; l-- ) { - if ( solmat[l][k] ) { - NEXTMP(mp0,mp); - mp->c = (Obj)solmat[l][k]; - mp->dl = BDY(mb[rinfo[l]])->dl; - } - } - NEXT(mp) = 0; MKDP(n,mp0,u); - NEXTNODE(nd0,nd); - BDY(nd) = (pointer)u; - NEXT(nd) = 0; - } - NEXT(nd) = 0; - return nd0; - } + if ( !(nf=current_numberfield) ) + error("invdalg : current_numberfield is not set"); + if ( !a ) + error("invdalg : division by 0"); + else if ( NID(a) == N_Q ) { + invq((Q)a,&dn); *c = (DAlg)dn; + return 0; + } + dim = nf->dim; + mb = nf->mb; + n = nf->n; + ln = ONEN; + dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc); + MKDAlg(u,ONE,a0); + simp = (DAlg *)MALLOC(dim*sizeof(DAlg)); + current_spec = dp_current_spec; initd(nf->spec); + for ( i = 0; i < dim; i++ ) { + if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); } + m = mb[i]; + for ( j = i-1; j >= 0; j-- ) + if ( dp_redble(m,mb[j]) ) + break; + if ( j >= 0 ) { + dp_subd(m,mb[j],&d); + if ( simp[j] ) { + muld(CO,d,simp[j]->nm,&u); + MKDAlg(u,simp[j]->dn,t); + simpdalg(t,&simp[i]); + } else + simp[i] = 0; + } else { + MKDAlg(m,ONE,t); + muldalg(t,a0,&simp[i]); + } + if ( simp[i] ) { + gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn); + muln(NM(simp[i]->dn),qn,&ln); + } + } + initd(current_spec); + NTOQ(ln,1,dn); + MKMAT(mobj,dim,dim+1); + mat = (Q **)BDY(mobj); + mulq(dn,a->dn,&mat[0][dim]); + for ( j = 0; j < dim; j++ ) { + if ( simp[j] ) { + divq(dn,simp[j]->dn,&mul); + for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- ) + if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) { + mulq(mul,(Q)mp->c,&mat[i][j]); + mp = NEXT(mp); + } + } + } + get_eg(&eg0); + rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo); + get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1); + if ( cinfo[0] == dim ) { + /* the input is invertible */ + solmat = (Q **)BDY(sol); + for ( i = dim-1, mp0 = 0; i >= 0; i-- ) + if ( solmat[i][0] ) { + NEXTMP(mp0,mp); + mp->c = (Obj)solmat[i][0]; + mp->dl = BDY(mb[i])->dl; + } + NEXT(mp) = 0; MKDP(n,mp0,u); + mulq(dnsol,nmc,&dn1); + MKDAlg(u,dn1,r); + rmcontdalg(r,c); + return 0; + } else { + /* the input is not invertible */ + nparam = sol->col; + solmat = (Q **)BDY(sol); + nd0 = 0; + for ( k = 0; k < nparam; k++ ) { + /* construct a new basis element */ + m = mb[cinfo[k]]; + mp0 = 0; + NEXTMP(mp0,mp); + chsgnq(dnsol,&dn1); mp->c = (Obj)dn1; + mp->dl = BDY(m)->dl; + /* skip the last parameter */ + for ( l = rank-2; l >= 0; l-- ) { + if ( solmat[l][k] ) { + NEXTMP(mp0,mp); + mp->c = (Obj)solmat[l][k]; + mp->dl = BDY(mb[rinfo[l]])->dl; + } + } + NEXT(mp) = 0; MKDP(n,mp0,u); + NEXTNODE(nd0,nd); + BDY(nd) = (pointer)u; + NEXT(nd) = 0; + } + NEXT(nd) = 0; + return nd0; + } } NODE dp_inv_or_split(NODE gb,DP f,struct order_spec *spec, DP *inv) { - int dim,n,i,j,k,l,nv; - DP *mb,*ps; - DP m,d,u,nm; - N ln,gn,qn; - DAlg *simp; - DAlg a0,r; - Q dn,dnsol,mul,nmc,dn1,iq; - MAT mobj,sol; - Q **mat,**solmat; - MP mp0,mp; - int *rinfo,*cinfo; - int rank,nparam; - NODE nd0,nd,ndt,ind,indt,t,mblist; - struct oEGT eg0,eg1; - extern struct oEGT eg_le; - extern int DP_Print; - initd(spec); - dp_ptozp(f,&u); f = u; + int dim,n,i,j,k,l,nv; + DP *mb,*ps; + DP m,d,u,nm; + N ln,gn,qn; + DAlg *simp; + DAlg a0,r; + Q dn,dnsol,mul,nmc,dn1,iq; + MAT mobj,sol; + Q **mat,**solmat; + MP mp0,mp; + int *rinfo,*cinfo; + int rank,nparam; + NODE nd0,nd,ndt,ind,indt,t,mblist; + struct oEGT eg0,eg1; + extern struct oEGT eg_le; + extern int DP_Print; + initd(spec); + dp_ptozp(f,&u); f = u; - n = length(gb); - ps = (DP *)MALLOC(n*sizeof(DP)); - for ( ind = 0, i = 0, t = gb; i < n; i++, t = NEXT(t) ) { - ps[i] = (DP)BDY(t); - NEXTNODE(ind,indt); - STOQ(i,iq); BDY(indt) = iq; - } - if ( ind ) NEXT(indt) = 0; - dp_true_nf(ind,f,ps,1,&nm,(P *)&dn); - if ( !nm ) error("dp_inv_or_split : input is 0"); - f = nm; + n = length(gb); + ps = (DP *)MALLOC(n*sizeof(DP)); + for ( ind = 0, i = 0, t = gb; i < n; i++, t = NEXT(t) ) { + ps[i] = (DP)BDY(t); + NEXTNODE(ind,indt); + STOQ(i,iq); BDY(indt) = iq; + } + if ( ind ) NEXT(indt) = 0; + dp_true_nf(ind,f,ps,1,&nm,(P *)&dn); + if ( !nm ) error("dp_inv_or_split : input is 0"); + f = nm; - dp_mbase(gb,&mblist); - dim = length(mblist); - mb = (DP *)MALLOC(dim*sizeof(DP)); - for ( i = 0, t = mblist; i < dim; i++, t = NEXT(t) ) - mb[dim-i-1] = (DP)BDY(t); - nv = mb[0]->nv; - ln = ONEN; - simp = (DAlg *)MALLOC(dim*sizeof(DAlg)); - for ( i = 0; i < dim; i++ ) { - if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); } - m = mb[i]; - for ( j = i-1; j >= 0; j-- ) - if ( dp_redble(m,mb[j]) ) - break; - if ( j >= 0 ) { - dp_subd(m,mb[j],&d); - if ( simp[j] ) { - muld(CO,d,simp[j]->nm,&u); - dp_true_nf(ind,u,ps,1,&nm,(P *)&dn); - mulq(simp[j]->dn,dn,&dn1); - MKDAlg(nm,dn1,simp[i]); - } else - simp[i] = 0; - } else { - dp_true_nf(ind,f,ps,1,&nm,(P *)&dn); - MKDAlg(nm,dn,simp[i]); - } - if ( simp[i] ) { - gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn); - muln(NM(simp[i]->dn),qn,&ln); - } - } - NTOQ(ln,1,dn); - MKMAT(mobj,dim,dim+1); - mat = (Q **)BDY(mobj); - mat[0][dim] = dn; - for ( j = 0; j < dim; j++ ) { - if ( simp[j] ) { - divq(dn,simp[j]->dn,&mul); - for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- ) - if ( dl_equal(nv,BDY(mb[i])->dl,mp->dl) ) { - mulq(mul,(Q)mp->c,&mat[i][j]); - mp = NEXT(mp); - } - } - } - get_eg(&eg0); - rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo); - get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1); - if ( cinfo[0] == dim ) { - /* the input is invertible */ - solmat = (Q **)BDY(sol); - for ( i = dim-1, mp0 = 0; i >= 0; i-- ) - if ( solmat[i][0] ) { - NEXTMP(mp0,mp); - mp->c = (Obj)solmat[i][0]; - mp->dl = BDY(mb[i])->dl; - } - NEXT(mp) = 0; MKDP(nv,mp0,*inv); - return 0; - } else { - /* the input is not invertible */ - nparam = sol->col; - solmat = (Q **)BDY(sol); - nd0 = 0; - for ( k = 0; k < nparam; k++ ) { - /* construct a new basis element */ - m = mb[cinfo[k]]; - mp0 = 0; - NEXTMP(mp0,mp); - chsgnq(dnsol,&dn1); mp->c = (Obj)dn1; - mp->dl = BDY(m)->dl; - /* skip the last parameter */ - for ( l = rank-2; l >= 0; l-- ) { - if ( solmat[l][k] ) { - NEXTMP(mp0,mp); - mp->c = (Obj)solmat[l][k]; - mp->dl = BDY(mb[rinfo[l]])->dl; - } - } - NEXT(mp) = 0; MKDP(nv,mp0,u); - NEXTNODE(nd0,nd); - BDY(nd) = (pointer)u; - NEXT(nd) = 0; - } - NEXT(nd) = 0; - return nd0; - } + dp_mbase(gb,&mblist); + dim = length(mblist); + mb = (DP *)MALLOC(dim*sizeof(DP)); + for ( i = 0, t = mblist; i < dim; i++, t = NEXT(t) ) + mb[dim-i-1] = (DP)BDY(t); + nv = mb[0]->nv; + ln = ONEN; + simp = (DAlg *)MALLOC(dim*sizeof(DAlg)); + for ( i = 0; i < dim; i++ ) { + if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); } + m = mb[i]; + for ( j = i-1; j >= 0; j-- ) + if ( dp_redble(m,mb[j]) ) + break; + if ( j >= 0 ) { + dp_subd(m,mb[j],&d); + if ( simp[j] ) { + muld(CO,d,simp[j]->nm,&u); + dp_true_nf(ind,u,ps,1,&nm,(P *)&dn); + mulq(simp[j]->dn,dn,&dn1); + MKDAlg(nm,dn1,simp[i]); + } else + simp[i] = 0; + } else { + dp_true_nf(ind,f,ps,1,&nm,(P *)&dn); + MKDAlg(nm,dn,simp[i]); + } + if ( simp[i] ) { + gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn); + muln(NM(simp[i]->dn),qn,&ln); + } + } + NTOQ(ln,1,dn); + MKMAT(mobj,dim,dim+1); + mat = (Q **)BDY(mobj); + mat[0][dim] = dn; + for ( j = 0; j < dim; j++ ) { + if ( simp[j] ) { + divq(dn,simp[j]->dn,&mul); + for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- ) + if ( dl_equal(nv,BDY(mb[i])->dl,mp->dl) ) { + mulq(mul,(Q)mp->c,&mat[i][j]); + mp = NEXT(mp); + } + } + } + get_eg(&eg0); + rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo); + get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1); + if ( cinfo[0] == dim ) { + /* the input is invertible */ + solmat = (Q **)BDY(sol); + for ( i = dim-1, mp0 = 0; i >= 0; i-- ) + if ( solmat[i][0] ) { + NEXTMP(mp0,mp); + mp->c = (Obj)solmat[i][0]; + mp->dl = BDY(mb[i])->dl; + } + NEXT(mp) = 0; MKDP(nv,mp0,*inv); + return 0; + } else { + /* the input is not invertible */ + nparam = sol->col; + solmat = (Q **)BDY(sol); + nd0 = 0; + for ( k = 0; k < nparam; k++ ) { + /* construct a new basis element */ + m = mb[cinfo[k]]; + mp0 = 0; + NEXTMP(mp0,mp); + chsgnq(dnsol,&dn1); mp->c = (Obj)dn1; + mp->dl = BDY(m)->dl; + /* skip the last parameter */ + for ( l = rank-2; l >= 0; l-- ) { + if ( solmat[l][k] ) { + NEXTMP(mp0,mp); + mp->c = (Obj)solmat[l][k]; + mp->dl = BDY(mb[rinfo[l]])->dl; + } + } + NEXT(mp) = 0; MKDP(nv,mp0,u); + NEXTNODE(nd0,nd); + BDY(nd) = (pointer)u; + NEXT(nd) = 0; + } + NEXT(nd) = 0; + return nd0; + } } void chsgndalg(DAlg a,DAlg *c) { - DP nm; - Q t; + DP nm; + Q t; - if ( !a ) *c = 0; - else if ( NID(a) == N_Q ) { - chsgnq((Q)a,&t); *c = (DAlg)t; - } else { - chsgnd(a->nm,&nm); - MKDAlg(nm,a->dn,*c); - } + if ( !a ) *c = 0; + else if ( NID(a) == N_Q ) { + chsgnq((Q)a,&t); *c = (DAlg)t; + } else { + chsgnd(a->nm,&nm); + MKDAlg(nm,a->dn,*c); + } } void pwrdalg(DAlg a,Q e,DAlg *c) { - NumberField nf; - DAlg t,z,y; - Q q; - N en,qn; - int r; - int ret; + NumberField nf; + DAlg t,z,y; + Q q; + N en,qn; + int r; + int ret; - if ( !(nf=current_numberfield) ) - error("pwrdalg : current_numberfield is not set"); - if ( !a ) - *c = !e ? (DAlg)ONE : 0; - else if ( NID(a) == N_Q ) { - pwrq((Q)a,e,&q); *c = (DAlg)q; - } else if ( !e ) - *c = nf->one; - else if ( UNIQ(e) ) - *c = a; - else { - if ( SGN(e) < 0 ) { - ret = invdalg(a,&t); - if ( !ret ) - error("pwrdalg : the denominator is not invertible"); - a = t; - } - en = NM(e); - y = nf->one; - z = a; - while ( 1 ) { - r = divin(en,2,&qn); en = qn; - if ( r ) { - muldalg(z,y,&t); y = t; - if ( !en ) { - *c = y; - return; - } - } - muldalg(z,z,&t); z = t; - } - } + if ( !(nf=current_numberfield) ) + error("pwrdalg : current_numberfield is not set"); + if ( !a ) + *c = !e ? (DAlg)ONE : 0; + else if ( NID(a) == N_Q ) { + pwrq((Q)a,e,&q); *c = (DAlg)q; + } else if ( !e ) + *c = nf->one; + else if ( UNIQ(e) ) + *c = a; + else { + if ( SGN(e) < 0 ) { + ret = invdalg(a,&t); + if ( !ret ) + error("pwrdalg : the denominator is not invertible"); + a = t; + } + en = NM(e); + y = nf->one; + z = a; + while ( 1 ) { + r = divin(en,2,&qn); en = qn; + if ( r ) { + muldalg(z,y,&t); y = t; + if ( !en ) { + *c = y; + return; + } + } + muldalg(z,z,&t); z = t; + } + } } int cmpdalg(DAlg a,DAlg b) { - DAlg c; + DAlg c; - subdalg(a,b,&c); - if ( !c ) return 0; - else - return SGN((Q)BDY(c->nm)->c); + subdalg(a,b,&c); + if ( !c ) return 0; + else + return SGN((Q)BDY(c->nm)->c); } /* convert da to a univariate poly; return the position of variable */ int dalgtoup(DAlg da,P *up,Q *dn) { - int nv,i,hi,current_d; - DCP dc0,dc; - MP h,mp0,mp,t; - DL hd,d; - DP c; - DAlg cc; - P v; + int nv,i,hi,current_d; + DCP dc0,dc; + MP h,mp0,mp,t; + DL hd,d; + DP c; + DAlg cc; + P v; - nv = da->nm->nv; - h = BDY(da->nm); - *dn = da->dn; - hd = h->dl; - for ( i = 0; i < nv; i++ ) - if ( hd->d[i] ) break; - hi = i; - current_d = hd->d[i]; - dc0 = 0; - mp0 = 0; - for ( t = h; t; t = NEXT(t) ) { - NEWDL(d,nv); - for ( i = 0; i <= hi; i++ ) d->d[i] = 0; - for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i]; - d->td = t->dl->td - t->dl->d[hi]; - if ( t->dl->d[hi] != current_d ) { - NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc); - NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc; - current_d = t->dl->d[hi]; - mp0 = 0; - } - NEXTMP(mp0,mp); - mp->c = t->c; mp->dl = d; - } - NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc); - NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc; - NEXT(dc) = 0; - makevar("x",&v); - MKP(VR(v),dc0,*up); - return hi; + nv = da->nm->nv; + h = BDY(da->nm); + *dn = da->dn; + hd = h->dl; + for ( i = 0; i < nv; i++ ) + if ( hd->d[i] ) break; + hi = i; + current_d = hd->d[i]; + dc0 = 0; + mp0 = 0; + for ( t = h; t; t = NEXT(t) ) { + NEWDL(d,nv); + for ( i = 0; i <= hi; i++ ) d->d[i] = 0; + for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i]; + d->td = t->dl->td - t->dl->d[hi]; + if ( t->dl->d[hi] != current_d ) { + NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc); + NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc; + current_d = t->dl->d[hi]; + mp0 = 0; + } + NEXTMP(mp0,mp); + mp->c = t->c; mp->dl = d; + } + NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc); + NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc; + NEXT(dc) = 0; + makevar("x",&v); + MKP(VR(v),dc0,*up); + return hi; }