=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/algnum.c,v retrieving revision 1.15 retrieving revision 1.16 diff -u -p -r1.15 -r1.16 --- OpenXM_contrib2/asir2000/builtin/algnum.c 2017/08/31 02:36:20 1.15 +++ OpenXM_contrib2/asir2000/builtin/algnum.c 2018/03/29 01:32:50 1.16 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.14 2013/11/17 17:34:59 ohara Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.15 2017/08/31 02:36:20 noro Exp $ */ #include "ca.h" #include "parse.h" @@ -72,841 +72,841 @@ void Pdalgtodp(); void Pdptodalg(); struct ftab alg_tab[] = { - {"set_field",Pset_field,-3}, - {"get_field_defpoly",Pget_field_defpoly,1}, - {"get_field_generator",Pget_field_generator,1}, - {"algtodalg",Palgtodalg,1}, - {"dalgtoalg",Pdalgtoalg,1}, - {"dalgtodp",Pdalgtodp,1}, - {"dalgtoup",Pdalgtoup,1}, - {"dptodalg",Pdptodalg,1}, - {"inv_or_split_dalg",Pinv_or_split_dalg,1}, - {"invalg_chrem",Pinvalg_chrem,2}, - {"invalg_le",Pinvalg_le,1}, - {"defpoly",Pdefpoly,1}, - {"newalg",Pnewalg,1}, - {"mainalg",Pmainalg,1}, - {"algtorat",Palgtorat,1}, - {"rattoalg",Prattoalg,1}, - {"getalg",Pgetalg,1}, - {"getalgtree",Pgetalgtree,1}, - {"alg",Palg,1}, - {"algv",Palgv,1}, - {0,0,0}, + {"set_field",Pset_field,-3}, + {"get_field_defpoly",Pget_field_defpoly,1}, + {"get_field_generator",Pget_field_generator,1}, + {"algtodalg",Palgtodalg,1}, + {"dalgtoalg",Pdalgtoalg,1}, + {"dalgtodp",Pdalgtodp,1}, + {"dalgtoup",Pdalgtoup,1}, + {"dptodalg",Pdptodalg,1}, + {"inv_or_split_dalg",Pinv_or_split_dalg,1}, + {"invalg_chrem",Pinvalg_chrem,2}, + {"invalg_le",Pinvalg_le,1}, + {"defpoly",Pdefpoly,1}, + {"newalg",Pnewalg,1}, + {"mainalg",Pmainalg,1}, + {"algtorat",Palgtorat,1}, + {"rattoalg",Prattoalg,1}, + {"getalg",Pgetalg,1}, + {"getalgtree",Pgetalgtree,1}, + {"alg",Palg,1}, + {"algv",Palgv,1}, + {0,0,0}, }; static int UCN,ACNT; void Pset_field(NODE arg,Q *rp) { - int ac; - NODE a0,a1; - VL vl0,vl; - struct order_spec *spec; + int ac; + NODE a0,a1; + VL vl0,vl; + struct order_spec *spec; - if ( (ac = argc(arg)) == 1 ) - setfield_dalg(BDY((LIST)ARG0(arg))); - else if ( ac == 3 ) { - a0 = BDY((LIST)ARG0(arg)); - a1 = BDY((LIST)ARG1(arg)); - for ( vl0 = 0; a1; a1 = NEXT(a1) ) { - NEXTVL(vl0,vl); - vl->v = VR((P)BDY(a1)); - } - if ( vl0 ) NEXT(vl) = 0; - create_order_spec(0,ARG2(arg),&spec); - setfield_gb(a0,vl0,spec); - } - *rp = 0; + if ( (ac = argc(arg)) == 1 ) + setfield_dalg(BDY((LIST)ARG0(arg))); + else if ( ac == 3 ) { + a0 = BDY((LIST)ARG0(arg)); + a1 = BDY((LIST)ARG1(arg)); + for ( vl0 = 0; a1; a1 = NEXT(a1) ) { + NEXTVL(vl0,vl); + vl->v = VR((P)BDY(a1)); + } + if ( vl0 ) NEXT(vl) = 0; + create_order_spec(0,ARG2(arg),&spec); + setfield_gb(a0,vl0,spec); + } + *rp = 0; } void Palgtodalg(NODE arg,DAlg *rp) { - algtodalg((Alg)ARG0(arg),rp); + algtodalg((Alg)ARG0(arg),rp); } void Pdalgtoalg(NODE arg,Alg *rp) { - dalgtoalg((DAlg)ARG0(arg),rp); + dalgtoalg((DAlg)ARG0(arg),rp); } void Pdalgtodp(NODE arg,LIST *r) { - NODE b; - DP nm; - Q dn; - DAlg da; + NODE b; + DP nm; + Q dn; + DAlg da; - da = (DAlg)ARG0(arg); - nm = da->nm; - dn = da->dn; - b = mknode(2,nm,dn); - MKLIST(*r,b); + da = (DAlg)ARG0(arg); + nm = da->nm; + dn = da->dn; + b = mknode(2,nm,dn); + MKLIST(*r,b); } void Pdptodalg(NODE arg,DAlg *r) { - DP d,nm,nm1; - MP m; - Q c,a; - DAlg t; - - d = (DP)ARG0(arg); - if ( !d ) *r = 0; - else { - for ( m = BDY(d); m; m = NEXT(m) ) - if ( !INT((Q)m->c) ) break; - if ( !m ) { - MKDAlg(d,(Q)ONE,t); - } else { - dp_ptozp(d,&nm); - divq((Q)BDY(d)->c,(Q)BDY(nm)->c,&c); - NTOQ(NM(c),SGN(c),a); - muldc(CO,nm,(Obj)a,&nm1); - NTOQ(DN(c),1,a); - MKDAlg(nm1,a,t); - } - simpdalg(t,r); - } + DP d,nm,nm1; + MP m; + Q c,a; + DAlg t; + + d = (DP)ARG0(arg); + if ( !d ) *r = 0; + else { + for ( m = BDY(d); m; m = NEXT(m) ) + if ( !INT((Q)m->c) ) break; + if ( !m ) { + MKDAlg(d,(Q)ONE,t); + } else { + dp_ptozp(d,&nm); + divq((Q)BDY(d)->c,(Q)BDY(nm)->c,&c); + NTOQ(NM(c),SGN(c),a); + muldc(CO,nm,(Obj)a,&nm1); + NTOQ(DN(c),1,a); + MKDAlg(nm1,a,t); + } + simpdalg(t,r); + } } void Pdalgtoup(NODE arg,LIST *r) { - NODE b; - int pos; - P up; - DP nm; - Q dn,q; + NODE b; + int pos; + P up; + DP nm; + Q dn,q; - pos = dalgtoup((DAlg)ARG0(arg),&up,&dn); - STOQ(pos,q); - b = mknode(3,up,dn,q); - MKLIST(*r,b); + pos = dalgtoup((DAlg)ARG0(arg),&up,&dn); + STOQ(pos,q); + b = mknode(3,up,dn,q); + MKLIST(*r,b); } NODE inv_or_split_dalg(DAlg,DAlg *); -NumberField get_numberfield(); +NumberField get_numberfield(); void Pget_field_defpoly(NODE arg,DAlg *r) { - NumberField nf; - DP d; + NumberField nf; + DP d; - nf = get_numberfield(); - d = nf->ps[QTOS((Q)ARG0(arg))]; - MKDAlg(d,ONE,*r); + nf = get_numberfield(); + d = nf->ps[QTOS((Q)ARG0(arg))]; + MKDAlg(d,ONE,*r); } void Pget_field_generator(NODE arg,DAlg *r) { - int index,n,i; - DL dl; - MP m; - DP d; + int index,n,i; + DL dl; + MP m; + DP d; - index = QTOS((Q)ARG0(arg)); - n = get_numberfield()->n; - NEWDL(dl,n); - for ( i = 0; i < n; i++ ) dl->d[i] = 0; - dl->d[index] = 1; dl->td = 1; - NEWMP(m); m->dl = dl; m->c = (Obj)ONE; NEXT(m) = 0; - MKDP(n,m,d); - MKDAlg(d,ONE,*r); + index = QTOS((Q)ARG0(arg)); + n = get_numberfield()->n; + NEWDL(dl,n); + for ( i = 0; i < n; i++ ) dl->d[i] = 0; + dl->d[index] = 1; dl->td = 1; + NEWMP(m); m->dl = dl; m->c = (Obj)ONE; NEXT(m) = 0; + MKDP(n,m,d); + MKDAlg(d,ONE,*r); } void Pinv_or_split_dalg(NODE arg,Obj *rp) { - NODE gen,t,nd0,nd; - LIST list; - int l,i,j,n; - DP *ps,*ps1,*psw; - NumberField nf; - DAlg inv; - extern struct order_spec *dp_current_spec; - struct order_spec *current_spec; + NODE gen,t,nd0,nd; + LIST list; + int l,i,j,n; + DP *ps,*ps1,*psw; + NumberField nf; + DAlg inv; + extern struct order_spec *dp_current_spec; + struct order_spec *current_spec; - gen = inv_or_split_dalg((DAlg)ARG0(arg),&inv); - if ( !gen ) - *rp = (Obj)inv; - else { - nf = get_numberfield(); - current_spec = dp_current_spec; initd(nf->spec); - l = length(gen); - n = nf->n; - ps = nf->ps; - psw = (DP *)ALLOCA((n+l)*sizeof(DP)); - for ( i = j = 0; i < n; i++ ) { - for ( t = gen; t; t = NEXT(t) ) - if ( dp_redble(ps[i],(DP)BDY(t)) ) break; - if ( !t ) - psw[j++] = ps[i]; - } - nd0 = 0; - /* gen[0] < gen[1] < ... */ - /* psw[0] > psw[1] > ... */ - for ( i = j-1, t = gen; i >= 0 && t; ) { - NEXTNODE(nd0,nd); - if ( compd(CO,psw[i],(DP)BDY(t)) > 0 ) { - BDY(nd) = BDY(t); t = NEXT(t); - } else - BDY(nd) = (pointer)psw[i--]; - } - for ( ; i >= 0; i-- ) { - NEXTNODE(nd0,nd); BDY(nd) = (pointer)psw[i]; - } - for ( ; t; t = NEXT(t) ) { - NEXTNODE(nd0,nd); BDY(nd) = BDY(t); - } - NEXT(nd) = 0; - MKLIST(list,nd0); - initd(current_spec); - *rp = (Obj)list; - } + gen = inv_or_split_dalg((DAlg)ARG0(arg),&inv); + if ( !gen ) + *rp = (Obj)inv; + else { + nf = get_numberfield(); + current_spec = dp_current_spec; initd(nf->spec); + l = length(gen); + n = nf->n; + ps = nf->ps; + psw = (DP *)ALLOCA((n+l)*sizeof(DP)); + for ( i = j = 0; i < n; i++ ) { + for ( t = gen; t; t = NEXT(t) ) + if ( dp_redble(ps[i],(DP)BDY(t)) ) break; + if ( !t ) + psw[j++] = ps[i]; + } + nd0 = 0; + /* gen[0] < gen[1] < ... */ + /* psw[0] > psw[1] > ... */ + for ( i = j-1, t = gen; i >= 0 && t; ) { + NEXTNODE(nd0,nd); + if ( compd(CO,psw[i],(DP)BDY(t)) > 0 ) { + BDY(nd) = BDY(t); t = NEXT(t); + } else + BDY(nd) = (pointer)psw[i--]; + } + for ( ; i >= 0; i-- ) { + NEXTNODE(nd0,nd); BDY(nd) = (pointer)psw[i]; + } + for ( ; t; t = NEXT(t) ) { + NEXTNODE(nd0,nd); BDY(nd) = BDY(t); + } + NEXT(nd) = 0; + MKLIST(list,nd0); + initd(current_spec); + *rp = (Obj)list; + } } void Pnewalg(arg,rp) NODE arg; Alg *rp; { - P p; - VL vl; - P c; + P p; + VL vl; + P c; - p = (P)ARG0(arg); - if ( !p || OID(p) != O_P ) - error("newalg : invalid argument"); - clctv(CO,p,&vl); - if ( NEXT(vl) ) - error("newalg : invalid argument"); - c = COEF(DC(p)); - if ( !NUM(c) || !RATN(c) ) - error("newalg : invalid argument"); - mkalg(p,rp); + p = (P)ARG0(arg); + if ( !p || OID(p) != O_P ) + error("newalg : invalid argument"); + clctv(CO,p,&vl); + if ( NEXT(vl) ) + error("newalg : invalid argument"); + c = COEF(DC(p)); + if ( !NUM(c) || !RATN(c) ) + error("newalg : invalid argument"); + mkalg(p,rp); } void mkalg(p,r) P p; Alg *r; { - VL vl,mvl,nvl; - V a,tv; - char buf[BUFSIZ]; - char *name; - P x,t,s; - Num c; - DCP dc,dcr,dcr0; + VL vl,mvl,nvl; + V a,tv; + char buf[BUFSIZ]; + char *name; + P x,t,s; + Num c; + DCP dc,dcr,dcr0; - for ( vl = ALG; vl; vl = NEXT(vl) ) - if ( !cmpalgp(p,(P)vl->v->attr) ) { - a = vl->v; break; - } - if ( !vl ) { - NEWVL(vl); NEXT(vl) = ALG; ALG = vl; - NEWV(a); vl->v = a; - sprintf(buf,"#%d",ACNT++); - name = (char *)MALLOC(strlen(buf)+1); - strcpy(name,buf); NAME(a) = name; + for ( vl = ALG; vl; vl = NEXT(vl) ) + if ( !cmpalgp(p,(P)vl->v->attr) ) { + a = vl->v; break; + } + if ( !vl ) { + NEWVL(vl); NEXT(vl) = ALG; ALG = vl; + NEWV(a); vl->v = a; + sprintf(buf,"#%d",ACNT++); + name = (char *)MALLOC(strlen(buf)+1); + strcpy(name,buf); NAME(a) = name; - for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { - NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); c = (Num)COEF(dc); - if ( NID(c) != N_A ) - COEF(dcr) = (P)c; - else - COEF(dcr) = (P)BDY(((Alg)c)); - } - NEXT(dcr) = 0; MKP(a,dcr0,t); a->attr = (pointer)t; - - sprintf(buf,"t%s",name); makevar(buf,&s); + for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { + NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); c = (Num)COEF(dc); + if ( NID(c) != N_A ) + COEF(dcr) = (P)c; + else + COEF(dcr) = (P)BDY(((Alg)c)); + } + NEXT(dcr) = 0; MKP(a,dcr0,t); a->attr = (pointer)t; + + sprintf(buf,"t%s",name); makevar(buf,&s); - if ( NEXT(ALG) ) { - tv = (V)NEXT(ALG)->v->priv; - for ( vl = CO; NEXT(NEXT(vl)); vl = NEXT(vl) ); - nvl = NEXT(vl); NEXT(vl) = 0; - for ( vl = CO; NEXT(vl) && (NEXT(vl)->v != tv); vl = NEXT(vl) ); - mvl = NEXT(vl); NEXT(vl) = nvl; NEXT(nvl) = mvl; - } - - a->priv = (pointer)VR(s); VR(s)->priv = (pointer)a; - } - MKV(a,x); MKAlg(x,*r); + if ( NEXT(ALG) ) { + tv = (V)NEXT(ALG)->v->priv; + for ( vl = CO; NEXT(NEXT(vl)); vl = NEXT(vl) ); + nvl = NEXT(vl); NEXT(vl) = 0; + for ( vl = CO; NEXT(vl) && (NEXT(vl)->v != tv); vl = NEXT(vl) ); + mvl = NEXT(vl); NEXT(vl) = nvl; NEXT(nvl) = mvl; + } + + a->priv = (pointer)VR(s); VR(s)->priv = (pointer)a; + } + MKV(a,x); MKAlg(x,*r); } int cmpalgp(p,defp) P p,defp; { - DCP dc,dcd; - P t; + DCP dc,dcd; + P t; - for ( dc = DC(p), dcd = DC(defp); dc && dcd; - dc = NEXT(dc), dcd = NEXT(dcd) ) { - if ( cmpq(DEG(dc),DEG(dcd)) ) - break; - t = NID((Num)COEF(dc)) == N_A ? (P)BDY((Alg)COEF(dc)) : COEF(dc); - if ( compp(ALG,t,COEF(dcd)) ) - break; - } - if ( dc || dcd ) - return 1; - else - return 0; + for ( dc = DC(p), dcd = DC(defp); dc && dcd; + dc = NEXT(dc), dcd = NEXT(dcd) ) { + if ( cmpq(DEG(dc),DEG(dcd)) ) + break; + t = NID((Num)COEF(dc)) == N_A ? (P)BDY((Alg)COEF(dc)) : COEF(dc); + if ( compp(ALG,t,COEF(dcd)) ) + break; + } + if ( dc || dcd ) + return 1; + else + return 0; } void Pdefpoly(arg,rp) NODE arg; P *rp; { - asir_assert(ARG0(arg),O_N,"defpoly"); - algptop((P)VR((P)BDY((Alg)ARG0(arg)))->attr,rp); + asir_assert(ARG0(arg),O_N,"defpoly"); + algptop((P)VR((P)BDY((Alg)ARG0(arg)))->attr,rp); } void Pmainalg(arg,r) NODE arg; Alg *r; { - Num c; - V v; - P b; + Num c; + V v; + P b; - c = (Num)(ARG0(arg)); - if ( NID(c) <= N_R ) - *r = 0; - else { - v = VR((P)BDY((Alg)c)); MKV(v,b); MKAlg(b,*r); - } + c = (Num)(ARG0(arg)); + if ( NID(c) <= N_R ) + *r = 0; + else { + v = VR((P)BDY((Alg)c)); MKV(v,b); MKAlg(b,*r); + } } void Palgtorat(arg,rp) NODE arg; Obj *rp; { - asir_assert(ARG0(arg),O_N,"algtorat"); - algtorat((Num)ARG0(arg),rp); + asir_assert(ARG0(arg),O_N,"algtorat"); + algtorat((Num)ARG0(arg),rp); } void Prattoalg(arg,rp) NODE arg; Alg *rp; { - asir_assert(ARG0(arg),O_R,"rattoalg"); - rattoalg((Obj)ARG0(arg),rp); + asir_assert(ARG0(arg),O_R,"rattoalg"); + rattoalg((Obj)ARG0(arg),rp); } void Pgetalg(arg,rp) NODE arg; LIST *rp; { - Obj t; - P p; - VL vl; - Num a; - Alg b; - NODE n0,n; + Obj t; + P p; + VL vl; + Num a; + Alg b; + NODE n0,n; - if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R ) - vl = 0; - else { - t = BDY((Alg)a); - switch ( OID(t) ) { - case O_P: case O_R: - clctvr(ALG,t,&vl); break; - default: - vl = 0; break; - } - } - for ( n0 = 0; vl; vl = NEXT(vl) ) { - NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b; - } - if ( n0 ) - NEXT(n) = 0; - MKLIST(*rp,n0); + if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R ) + vl = 0; + else { + t = BDY((Alg)a); + switch ( OID(t) ) { + case O_P: case O_R: + clctvr(ALG,t,&vl); break; + default: + vl = 0; break; + } + } + for ( n0 = 0; vl; vl = NEXT(vl) ) { + NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b; + } + if ( n0 ) + NEXT(n) = 0; + MKLIST(*rp,n0); } void Pgetalgtree(arg,rp) NODE arg; LIST *rp; { - Obj t; - P p; - VL vl,vl1,vl2; - Num a; - Alg b; - NODE n0,n; + Obj t; + P p; + VL vl,vl1,vl2; + Num a; + Alg b; + NODE n0,n; #if 0 - if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R ) - vl = 0; - else { - t = BDY((Alg)a); - switch ( OID(t) ) { - case O_P: - clctalg((P)t,&vl); break; - case O_R: - clctalg(NM((R)t),&vl1); - clctalg(DN((R)t),&vl2); - mergev(ALG,vl1,vl2,&vl); break; - default: - vl = 0; break; - } - } + if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R ) + vl = 0; + else { + t = BDY((Alg)a); + switch ( OID(t) ) { + case O_P: + clctalg((P)t,&vl); break; + case O_R: + clctalg(NM((R)t),&vl1); + clctalg(DN((R)t),&vl2); + mergev(ALG,vl1,vl2,&vl); break; + default: + vl = 0; break; + } + } #else - get_algtree((Obj)ARG0(arg),&vl); + get_algtree((Obj)ARG0(arg),&vl); #endif - for ( n0 = 0; vl; vl = NEXT(vl) ) { - NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b; - } - if ( n0 ) - NEXT(n) = 0; - MKLIST(*rp,n0); + for ( n0 = 0; vl; vl = NEXT(vl) ) { + NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b; + } + if ( n0 ) + NEXT(n) = 0; + MKLIST(*rp,n0); } void clctalg(p,vl) P p; VL *vl; { - int n,i; - VL tvl; - VN vn,vn1; - P d; - DCP dc; + int n,i; + VL tvl; + VN vn,vn1; + P d; + DCP dc; - for ( n = 0, tvl = ALG; tvl; tvl = NEXT(tvl), n++ ); - vn = (VN) ALLOCA((n+1)*sizeof(struct oVN)); - for ( i = n-1, tvl = ALG; tvl; tvl = NEXT(tvl), i-- ) { - vn[i].v = tvl->v; - vn[i].n = 0; - } - markv(vn,n,p); - for ( i = n-1; i >= 0; i-- ) { - if ( !vn[i].n ) - continue; - d = (P)vn[i].v->attr; - for ( dc = DC(d); dc; dc = NEXT(dc) ) - markv(vn,i,COEF(dc)); - } - vn1 = (VN) ALLOCA((n+1)*sizeof(struct oVN)); - for ( i = 0; i < n; i++ ) { - vn1[i].v = vn[n-1-i].v; vn1[i].n = vn[n-1-i].n; - } - vntovl(vn1,n,vl); + for ( n = 0, tvl = ALG; tvl; tvl = NEXT(tvl), n++ ); + vn = (VN) ALLOCA((n+1)*sizeof(struct oVN)); + for ( i = n-1, tvl = ALG; tvl; tvl = NEXT(tvl), i-- ) { + vn[i].v = tvl->v; + vn[i].n = 0; + } + markv(vn,n,p); + for ( i = n-1; i >= 0; i-- ) { + if ( !vn[i].n ) + continue; + d = (P)vn[i].v->attr; + for ( dc = DC(d); dc; dc = NEXT(dc) ) + markv(vn,i,COEF(dc)); + } + vn1 = (VN) ALLOCA((n+1)*sizeof(struct oVN)); + for ( i = 0; i < n; i++ ) { + vn1[i].v = vn[n-1-i].v; vn1[i].n = vn[n-1-i].n; + } + vntovl(vn1,n,vl); } void Palg(arg,rp) NODE arg; Alg *rp; { - Q a; - VL vl; - P x; - int n; + Q a; + VL vl; + P x; + int n; - a = (Q)ARG0(arg); - if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) ) - *rp = 0; - else { - n = ACNT-QTOS(a)-1; - for ( vl = ALG; vl && n; vl = NEXT(vl), n-- ); - if ( vl ) { - MKV(vl->v,x); MKAlg(x,*rp); - } else - *rp = 0; - } + a = (Q)ARG0(arg); + if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) ) + *rp = 0; + else { + n = ACNT-QTOS(a)-1; + for ( vl = ALG; vl && n; vl = NEXT(vl), n-- ); + if ( vl ) { + MKV(vl->v,x); MKAlg(x,*rp); + } else + *rp = 0; + } } void Palgv(arg,rp) NODE arg; Obj *rp; { - Q a; - VL vl; - P x; - int n; - Alg b; + Q a; + VL vl; + P x; + int n; + Alg b; - a = (Q)ARG0(arg); - if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) ) - *rp = 0; - else { - n = ACNT-QTOS(a)-1; - for ( vl = ALG; vl && n; vl = NEXT(vl), n-- ); - if ( vl ) { - MKV(vl->v,x); MKAlg(x,b); algtorat((Num)b,rp); - } else - *rp = 0; - } + a = (Q)ARG0(arg); + if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) ) + *rp = 0; + else { + n = ACNT-QTOS(a)-1; + for ( vl = ALG; vl && n; vl = NEXT(vl), n-- ); + if ( vl ) { + MKV(vl->v,x); MKAlg(x,b); algtorat((Num)b,rp); + } else + *rp = 0; + } } void algptop(p,r) P p,*r; { - DCP dc,dcr,dcr0; + DCP dc,dcr,dcr0; - if ( NUM(p) ) - *r = (P)p; - else { - for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { - NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); - algptop(COEF(dc),&COEF(dcr)); - } - NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r); - } + if ( NUM(p) ) + *r = (P)p; + else { + for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { + NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); + algptop(COEF(dc),&COEF(dcr)); + } + NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r); + } } void algtorat(n,r) Num n; Obj *r; { - Obj obj; - P nm,dn; + Obj obj; + P nm,dn; - if ( !n || NID(n) <= N_R ) - *r = (Obj)n; - else { - obj = BDY((Alg)n); - if ( ID(obj) <= O_P ) - algptop((P)obj,(P *)r); - else { - algptop(NM((R)obj),&nm); algptop(DN((R)obj),&dn); - divr(CO,(Obj)nm,(Obj)dn,r); - } - } + if ( !n || NID(n) <= N_R ) + *r = (Obj)n; + else { + obj = BDY((Alg)n); + if ( ID(obj) <= O_P ) + algptop((P)obj,(P *)r); + else { + algptop(NM((R)obj),&nm); algptop(DN((R)obj),&dn); + divr(CO,(Obj)nm,(Obj)dn,r); + } + } } void rattoalg(obj,n) Obj obj; Alg *n; { - P nm,dn; - Obj t; + P nm,dn; + Obj t; - if ( !obj || ID(obj) == O_N ) - *n = (Alg)obj; - else if ( ID(obj) == O_P ) { - ptoalgp((P)obj,(P *)&t); MKAlg(t,*n); - } else { - ptoalgp(NM((R)obj),&nm); ptoalgp(DN((R)obj),&dn); - divr(ALG,(Obj)nm,(Obj)dn,&t); MKAlg(t,*n); - } + if ( !obj || ID(obj) == O_N ) + *n = (Alg)obj; + else if ( ID(obj) == O_P ) { + ptoalgp((P)obj,(P *)&t); MKAlg(t,*n); + } else { + ptoalgp(NM((R)obj),&nm); ptoalgp(DN((R)obj),&dn); + divr(ALG,(Obj)nm,(Obj)dn,&t); MKAlg(t,*n); + } } void ptoalgp(p,r) P p,*r; { - DCP dc,dcr,dcr0; + DCP dc,dcr,dcr0; - if ( NUM(p) ) - *r = (P)p; - else { - for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { - NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); - ptoalgp(COEF(dc),&COEF(dcr)); - } - NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r); - } + if ( NUM(p) ) + *r = (P)p; + else { + for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { + NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); + ptoalgp(COEF(dc),&COEF(dcr)); + } + NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r); + } } void Pinvalg_chrem(NODE arg,LIST *r) { - NODE n; + NODE n; - inva_chrem((P)ARG0(arg),(P)ARG1(arg),&n); - MKLIST(*r,n); + inva_chrem((P)ARG0(arg),(P)ARG1(arg),&n); + MKLIST(*r,n); } void invalg_le(Alg a,LIST *r); void Pinvalg_le(NODE arg,LIST *r) { - invalg_le((Alg)ARG0(arg),r); + invalg_le((Alg)ARG0(arg),r); } typedef struct oMono_nf { - DP mono; - DP nf; - Q dn; + DP mono; + DP nf; + Q dn; } *Mono_nf; void invalg_le(Alg a,LIST *r) { - Alg inv; - MAT mobj,sol; - int *rinfo,*cinfo; - P p,dn,dn1,ap; - VL vl,tvl; - Q c1,c2,c3,cont,c,two,iq,dn0,mul,dnsol; - int i,j,n,len,k; - MP mp,mp0; - DP dp,nm,nm1,m,d,u,u1; - NODE b,b1,hlist,mblist,t,s,rev0,rev,hist; - DP *ps; - struct order_spec *spec; - Mono_nf h,h1; - N nq,nr,nl,ng; - Q **mat,**solmat; - Q *w; - int *wi; + Alg inv; + MAT mobj,sol; + int *rinfo,*cinfo; + P p,dn,dn1,ap; + VL vl,tvl; + Q c1,c2,c3,cont,c,two,iq,dn0,mul,dnsol; + int i,j,n,len,k; + MP mp,mp0; + DP dp,nm,nm1,m,d,u,u1; + NODE b,b1,hlist,mblist,t,s,rev0,rev,hist; + DP *ps; + struct order_spec *spec; + Mono_nf h,h1; + N nq,nr,nl,ng; + Q **mat,**solmat; + Q *w; + int *wi; - ap = (P)BDY(a); - asir_assert(ap,O_P,"invalg_le"); + ap = (P)BDY(a); + asir_assert(ap,O_P,"invalg_le"); - /* collecting algebraic numbers */ - clctalg(ap,&vl); + /* collecting algebraic numbers */ + clctalg(ap,&vl); - /* setup */ - ptozp(ap,1,&c,&p); - STOQ(2,two); create_order_spec(0,(Obj)two,&spec); initd(spec); - for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ ); - ps = (DP *)ALLOCA(n*sizeof(DP)); + /* setup */ + ptozp(ap,1,&c,&p); + STOQ(2,two); create_order_spec(0,(Obj)two,&spec); initd(spec); + for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ ); + ps = (DP *)ALLOCA(n*sizeof(DP)); - /* conversion to DP */ - for ( i = 0, tvl = vl; i < n; i++, tvl = NEXT(tvl) ) { - ptod(ALG,vl,tvl->v->attr,&ps[i]); - } - ptod(ALG,vl,p,&dp); - /* index list */ - for ( b = 0, i = 0; i < n; i++ ) { - STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1; - } - /* simplification */ - dp_true_nf(b,dp,ps,1,&nm,&dn); + /* conversion to DP */ + for ( i = 0, tvl = vl; i < n; i++, tvl = NEXT(tvl) ) { + ptod(ALG,vl,tvl->v->attr,&ps[i]); + } + ptod(ALG,vl,p,&dp); + /* index list */ + for ( b = 0, i = 0; i < n; i++ ) { + STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1; + } + /* simplification */ + dp_true_nf(b,dp,ps,1,&nm,&dn); - /* construction of NF table */ + /* construction of NF table */ - /* stdmono: <<0,...,0>> < ... < max */ - for ( hlist = 0, i = 0; i < n; i++ ) { - MKNODE(b1,(pointer)ps[i],hlist); hlist = b1; - } - dp_mbase(hlist,&rev0); - for ( mblist = 0, rev = rev0; rev; rev = NEXT(rev) ) { - MKNODE(b1,BDY(rev),mblist); mblist = b1; - } - dn0 = ONE; - for ( hist = 0, t = mblist; t; t = NEXT(t) ) { - /* searching a predecessor */ - for ( m = (DP)BDY(t), s = hist; s; s = NEXT(s) ) { - h = (Mono_nf)BDY(s); - if ( dp_redble(m,h->mono) ) - break; - } - h1 = (Mono_nf)ALLOCA(sizeof(struct oMono_nf)); - if ( s ) { - dp_subd(m,h->mono,&d); - muld(CO,d,h->nf,&u); - dp_true_nf(b,u,ps,1,&nm1,&dn1); - mulq(h->dn,(Q)dn1,&h1->dn); - } else { - muld(CO,m,nm,&u); - dp_true_nf(b,u,ps,1,&nm1,&dn1); - h1->dn = (Q)dn1; - } - h1->mono = m; - h1->nf = nm1; - MKNODE(b1,(pointer)h1,hist); hist = b1; + /* stdmono: <<0,...,0>> < ... < max */ + for ( hlist = 0, i = 0; i < n; i++ ) { + MKNODE(b1,(pointer)ps[i],hlist); hlist = b1; + } + dp_mbase(hlist,&rev0); + for ( mblist = 0, rev = rev0; rev; rev = NEXT(rev) ) { + MKNODE(b1,BDY(rev),mblist); mblist = b1; + } + dn0 = ONE; + for ( hist = 0, t = mblist; t; t = NEXT(t) ) { + /* searching a predecessor */ + for ( m = (DP)BDY(t), s = hist; s; s = NEXT(s) ) { + h = (Mono_nf)BDY(s); + if ( dp_redble(m,h->mono) ) + break; + } + h1 = (Mono_nf)ALLOCA(sizeof(struct oMono_nf)); + if ( s ) { + dp_subd(m,h->mono,&d); + muld(CO,d,h->nf,&u); + dp_true_nf(b,u,ps,1,&nm1,&dn1); + mulq(h->dn,(Q)dn1,&h1->dn); + } else { + muld(CO,m,nm,&u); + dp_true_nf(b,u,ps,1,&nm1,&dn1); + h1->dn = (Q)dn1; + } + h1->mono = m; + h1->nf = nm1; + MKNODE(b1,(pointer)h1,hist); hist = b1; - /* dn0 = LCM(dn0,h1->dn) */ - gcdn(NM(dn0),NM(h1->dn),&ng); divn(NM(dn0),ng,&nq,&nr); - muln(nq,NM(h1->dn),&nl); NTOQ(nl,1,dn0); - } - /* create a matrix */ - len = length(mblist); - MKMAT(mobj,len,len+1); - mat = (Q **)BDY(mobj); - mat[len-1][len] = dn0; - for ( j = 0, t = hist; j < len; j++, t = NEXT(t) ) { - h = (Mono_nf)BDY(t); - nm1 = h->nf; - divq((Q)dn0,h->dn,&mul); - for ( i = 0, rev = rev0, mp = BDY(nm1); mp && i < len; i++, rev = NEXT(rev) ) - if ( dl_equal(n,BDY((DP)BDY(rev))->dl,mp->dl) ) { - mulq(mul,(Q)mp->c,&mat[i][j]); - mp = NEXT(mp); - } - } + /* dn0 = LCM(dn0,h1->dn) */ + gcdn(NM(dn0),NM(h1->dn),&ng); divn(NM(dn0),ng,&nq,&nr); + muln(nq,NM(h1->dn),&nl); NTOQ(nl,1,dn0); + } + /* create a matrix */ + len = length(mblist); + MKMAT(mobj,len,len+1); + mat = (Q **)BDY(mobj); + mat[len-1][len] = dn0; + for ( j = 0, t = hist; j < len; j++, t = NEXT(t) ) { + h = (Mono_nf)BDY(t); + nm1 = h->nf; + divq((Q)dn0,h->dn,&mul); + for ( i = 0, rev = rev0, mp = BDY(nm1); mp && i < len; i++, rev = NEXT(rev) ) + if ( dl_equal(n,BDY((DP)BDY(rev))->dl,mp->dl) ) { + mulq(mul,(Q)mp->c,&mat[i][j]); + mp = NEXT(mp); + } + } #if 0 - w = (Q *)ALLOCA((len+1)*sizeof(Q)); - wi = (int *)ALLOCA((len+1)*sizeof(int)); - for ( i = 0; i < len; i++ ) { - for ( j = 0, k = 0; j <= len; j++ ) - if ( mat[i][j] ) { - w[k] = mat[i][j]; - wi[k] = j; - k++; - } - removecont_array(w,k); - for ( j = 0; j < k; j++ ) - mat[i][wi[j]] = w[j]; - } + w = (Q *)ALLOCA((len+1)*sizeof(Q)); + wi = (int *)ALLOCA((len+1)*sizeof(int)); + for ( i = 0; i < len; i++ ) { + for ( j = 0, k = 0; j <= len; j++ ) + if ( mat[i][j] ) { + w[k] = mat[i][j]; + wi[k] = j; + k++; + } + removecont_array(w,k); + for ( j = 0; j < k; j++ ) + mat[i][wi[j]] = w[j]; + } #endif - generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo); - solmat = (Q **)BDY(sol); - for ( i = 0, t = rev0, mp0 = 0; i < len; i++, t = NEXT(t) ) - if ( solmat[i][0] ) { - NEXTMP(mp0,mp); - mp->c = (Obj)solmat[i][0]; - mp->dl = BDY((DP)BDY(t))->dl; - } - NEXT(mp) = 0; MKDP(n,mp0,u); - dp_ptozp(u,&u1); - divq((Q)BDY(u)->c,(Q)BDY(u1)->c,&cont); - dtop(ALG,vl,u1,(Obj *)&ap); - MKAlg(ap,inv); - mulq(dnsol,(Q)dn,&c1); - mulq(c1,c,&c2); - divq(c2,cont,&c3); - b = mknode(2,inv,c3); - MKLIST(*r,b); + generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo); + solmat = (Q **)BDY(sol); + for ( i = 0, t = rev0, mp0 = 0; i < len; i++, t = NEXT(t) ) + if ( solmat[i][0] ) { + NEXTMP(mp0,mp); + mp->c = (Obj)solmat[i][0]; + mp->dl = BDY((DP)BDY(t))->dl; + } + NEXT(mp) = 0; MKDP(n,mp0,u); + dp_ptozp(u,&u1); + divq((Q)BDY(u)->c,(Q)BDY(u1)->c,&cont); + dtop(ALG,vl,u1,(Obj *)&ap); + MKAlg(ap,inv); + mulq(dnsol,(Q)dn,&c1); + mulq(c1,c,&c2); + divq(c2,cont,&c3); + b = mknode(2,inv,c3); + MKLIST(*r,b); } void get_algtree(Obj f,VL *r) { - VL vl1,vl2,vl3; - Obj t; - DCP dc; - NODE b; - pointer *a; - pointer **m; - int len,row,col,i,j,l; + VL vl1,vl2,vl3; + Obj t; + DCP dc; + NODE b; + pointer *a; + pointer **m; + int len,row,col,i,j,l; - if ( !f ) *r = 0; - else - switch ( OID(f) ) { - case O_N: - if ( NID((Num)f) != N_A ) *r = 0; - else { - t = BDY((Alg)f); - switch ( OID(t) ) { - case O_P: - clctalg((P)t,r); break; - case O_R: - clctalg(NM((R)t),&vl1); - clctalg(DN((R)t),&vl2); - mergev(ALG,vl1,vl2,r); break; - default: - *r = 0; break; - } - } - break; - case O_P: - vl1 = 0; - for ( dc = DC((P)f); dc; dc = NEXT(dc) ) { - get_algtree((Obj)COEF(dc),&vl2); - mergev(ALG,vl1,vl2,&vl3); - vl1 = vl3; - } - *r = vl1; - break; - case O_R: - get_algtree((Obj)NM((R)f),&vl1); - get_algtree((Obj)DN((R)f),&vl2); - mergev(ALG,vl1,vl2,r); - break; - case O_LIST: - vl1 = 0; - for ( b = BDY((LIST)f); b; b = NEXT(b) ) { - get_algtree((Obj)BDY(b),&vl2); - mergev(ALG,vl1,vl2,&vl3); - vl1 = vl3; - } - *r = vl1; - break; - case O_VECT: - vl1 = 0; - l = ((VECT)f)->len; - a = BDY((VECT)f); - for ( i = 0; i < l; i++ ) { - get_algtree((Obj)a[i],&vl2); - mergev(ALG,vl1,vl2,&vl3); - vl1 = vl3; - } - *r = vl1; - break; - case O_MAT: - vl1 = 0; - row = ((MAT)f)->row; col = ((MAT)f)->col; - m = BDY((MAT)f); - for ( i = 0; i < row; i++ ) - for ( j = 0; j < col; j++ ) { - get_algtree((Obj)m[i][j],&vl2); - mergev(ALG,vl1,vl2,&vl3); - vl1 = vl3; - } - *r = vl1; - break; - default: - *r = 0; - break; - } + if ( !f ) *r = 0; + else + switch ( OID(f) ) { + case O_N: + if ( NID((Num)f) != N_A ) *r = 0; + else { + t = BDY((Alg)f); + switch ( OID(t) ) { + case O_P: + clctalg((P)t,r); break; + case O_R: + clctalg(NM((R)t),&vl1); + clctalg(DN((R)t),&vl2); + mergev(ALG,vl1,vl2,r); break; + default: + *r = 0; break; + } + } + break; + case O_P: + vl1 = 0; + for ( dc = DC((P)f); dc; dc = NEXT(dc) ) { + get_algtree((Obj)COEF(dc),&vl2); + mergev(ALG,vl1,vl2,&vl3); + vl1 = vl3; + } + *r = vl1; + break; + case O_R: + get_algtree((Obj)NM((R)f),&vl1); + get_algtree((Obj)DN((R)f),&vl2); + mergev(ALG,vl1,vl2,r); + break; + case O_LIST: + vl1 = 0; + for ( b = BDY((LIST)f); b; b = NEXT(b) ) { + get_algtree((Obj)BDY(b),&vl2); + mergev(ALG,vl1,vl2,&vl3); + vl1 = vl3; + } + *r = vl1; + break; + case O_VECT: + vl1 = 0; + l = ((VECT)f)->len; + a = BDY((VECT)f); + for ( i = 0; i < l; i++ ) { + get_algtree((Obj)a[i],&vl2); + mergev(ALG,vl1,vl2,&vl3); + vl1 = vl3; + } + *r = vl1; + break; + case O_MAT: + vl1 = 0; + row = ((MAT)f)->row; col = ((MAT)f)->col; + m = BDY((MAT)f); + for ( i = 0; i < row; i++ ) + for ( j = 0; j < col; j++ ) { + get_algtree((Obj)m[i][j],&vl2); + mergev(ALG,vl1,vl2,&vl3); + vl1 = vl3; + } + *r = vl1; + break; + default: + *r = 0; + break; + } } void algobjtorat(Obj f,Obj *r) { - Obj t; - DCP dc,dcr,dcr0; - P p,nm,dn; - R rat; - NODE b,s,s0; - VECT v; - MAT mat; - LIST list; - pointer *a; - pointer **m; - int len,row,col,i,j,l; + Obj t; + DCP dc,dcr,dcr0; + P p,nm,dn; + R rat; + NODE b,s,s0; + VECT v; + MAT mat; + LIST list; + pointer *a; + pointer **m; + int len,row,col,i,j,l; - if ( !f ) *r = 0; - else - switch ( OID(f) ) { - case O_N: - algtorat((Num)f,r); - break; - case O_P: - dcr0 = 0; - for ( dc = DC((P)f); dc; dc = NEXT(dc) ) { - NEXTDC(dcr0,dcr); - algobjtorat((Obj)COEF(dc),&t); - COEF(dcr) = (P)t; - DEG(dcr) = DEG(dc); - } - NEXT(dcr) = 0; MKP(VR((P)f),dcr0,p); *r = (Obj)p; - break; - case O_R: - algobjtorat((Obj)NM((R)f),&t); nm = (P)t; - algobjtorat((Obj)DN((R)f),&t); dn = (P)t; - MKRAT(nm,dn,0,rat); *r = (Obj)rat; - break; - case O_LIST: - s0 = 0; - for ( b = BDY((LIST)f); b; b = NEXT(b) ) { - NEXTNODE(s0,s); - algobjtorat((Obj)BDY(b),&t); - BDY(s) = (pointer)t; - } - NEXT(s) = 0; - MKLIST(list,s0); - *r = (Obj)list; - break; - case O_VECT: - l = ((VECT)f)->len; - a = BDY((VECT)f); - MKVECT(v,l); - for ( i = 0; i < l; i++ ) { - algobjtorat((Obj)a[i],&t); - BDY(v)[i] = (pointer)t; - } - *r = (Obj)v; - break; - case O_MAT: - row = ((MAT)f)->row; col = ((MAT)f)->col; - m = BDY((MAT)f); - MKMAT(mat,row,col); - for ( i = 0; i < row; i++ ) - for ( j = 0; j < col; j++ ) { - algobjtorat((Obj)m[i][j],&t); - BDY(mat)[i][j] = (pointer)t; - } - *r = (Obj)mat; - break; - default: - *r = f; - break; - } + if ( !f ) *r = 0; + else + switch ( OID(f) ) { + case O_N: + algtorat((Num)f,r); + break; + case O_P: + dcr0 = 0; + for ( dc = DC((P)f); dc; dc = NEXT(dc) ) { + NEXTDC(dcr0,dcr); + algobjtorat((Obj)COEF(dc),&t); + COEF(dcr) = (P)t; + DEG(dcr) = DEG(dc); + } + NEXT(dcr) = 0; MKP(VR((P)f),dcr0,p); *r = (Obj)p; + break; + case O_R: + algobjtorat((Obj)NM((R)f),&t); nm = (P)t; + algobjtorat((Obj)DN((R)f),&t); dn = (P)t; + MKRAT(nm,dn,0,rat); *r = (Obj)rat; + break; + case O_LIST: + s0 = 0; + for ( b = BDY((LIST)f); b; b = NEXT(b) ) { + NEXTNODE(s0,s); + algobjtorat((Obj)BDY(b),&t); + BDY(s) = (pointer)t; + } + NEXT(s) = 0; + MKLIST(list,s0); + *r = (Obj)list; + break; + case O_VECT: + l = ((VECT)f)->len; + a = BDY((VECT)f); + MKVECT(v,l); + for ( i = 0; i < l; i++ ) { + algobjtorat((Obj)a[i],&t); + BDY(v)[i] = (pointer)t; + } + *r = (Obj)v; + break; + case O_MAT: + row = ((MAT)f)->row; col = ((MAT)f)->col; + m = BDY((MAT)f); + MKMAT(mat,row,col); + for ( i = 0; i < row; i++ ) + for ( j = 0; j < col; j++ ) { + algobjtorat((Obj)m[i][j],&t); + BDY(mat)[i][j] = (pointer)t; + } + *r = (Obj)mat; + break; + default: + *r = f; + break; + } }