=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.121 retrieving revision 1.127 diff -u -p -r1.121 -r1.127 --- OpenXM_contrib2/asir2000/engine/nd.c 2004/12/09 08:56:43 1.121 +++ OpenXM_contrib2/asir2000/engine/nd.c 2005/08/03 05:01:01 1.127 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.120 2004/12/07 15:15:52 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.126 2005/02/09 14:30:47 noro Exp $ */ #include "nd.h" @@ -1537,7 +1537,7 @@ ND normalize_pbucket(int mod,PGeoBucket g) return r; } -void do_diagonalize(int sugar,int m) +int do_diagonalize(int sugar,int m) { int i,nh,stat; NODE r,g,t; @@ -1553,7 +1553,8 @@ void do_diagonalize(int sugar,int m) nfv = nd_ps[i]; s = ndvtond(m,nfv); s = nd_separate_head(s,&head); - nd_nf(m,s,nd_ps,1,&dn,&nf); + stat = nd_nf(m,s,nd_ps,1,&dn,&nf); + if ( !stat ) return 0; if ( !m ) { NTOQ(NM(dn.z),SGN(dn.z),num); mulq(HCQ(head),num,&q); HCQ(head) = q; @@ -1574,6 +1575,7 @@ void do_diagonalize(int sugar,int m) } else nd_ps[i] = nfv; } + return 1; } /* return value = 0 => input is not a GB */ @@ -1599,8 +1601,14 @@ NODE nd_gb(int m,int ishomo,int checkonly) again: l = nd_minp(d,&d); if ( SG(l) != sugar ) { - if ( ishomo ) do_diagonalize(sugar,m); - + if ( ishomo ) { + stat = do_diagonalize(sugar,m); + if ( !stat ) { + NEXT(l) = d; d = l; + d = nd_reconstruct(0,d); + goto again; + } + } sugar = SG(l); if ( DP_Print ) fprintf(asir_out,"%d",sugar); } @@ -1646,7 +1654,7 @@ again: return g; } -void do_diagonalize_trace(int sugar,int m) +int do_diagonalize_trace(int sugar,int m) { int i,nh,stat; NODE r,g,t; @@ -1659,7 +1667,8 @@ void do_diagonalize_trace(int sugar,int m) /* for nd_ps */ s = ndvtond(m,nd_ps[i]); s = nd_separate_head(s,&head); - nd_nf_pbucket(m,s,nd_ps,1,&nf); + stat = nd_nf_pbucket(m,s,nd_ps,1,&nf); + if ( !stat ) return 0; nf = nd_add(m,head,nf); ndv_free(nd_ps[i]); nd_ps[i] = ndtondv(m,nf); @@ -1672,7 +1681,8 @@ void do_diagonalize_trace(int sugar,int m) nfv = nd_ps_trace[i]; s = ndvtond(0,nfv); s = nd_separate_head(s,&head); - nd_nf(0,s,nd_ps_trace,1,&dn,&nf); + stat = nd_nf(0,s,nd_ps_trace,1,&dn,&nf); + if ( !stat ) return 0; NTOQ(NM(dn.z),SGN(dn.z),num); mulq(HCQ(head),num,&q); HCQ(head) = q; if ( DN(dn.z) ) { @@ -1691,6 +1701,7 @@ void do_diagonalize_trace(int sugar,int m) } else nd_ps_trace[i] = nfv; } + return 1; } static struct oEGT eg_invdalg; @@ -1721,7 +1732,14 @@ NODE nd_gb_trace(int m,int ishomo) again: l = nd_minp(d,&d); if ( SG(l) != sugar ) { - if ( ishomo ) do_diagonalize_trace(sugar,m); + if ( ishomo ) { + stat = do_diagonalize_trace(sugar,m); + if ( !stat ) { + NEXT(l) = d; d = l; + d = nd_reconstruct(1,d); + goto again; + } + } sugar = SG(l); if ( DP_Print ) fprintf(asir_out,"%d",sugar); } @@ -2348,6 +2366,7 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe break; } nd_nalg = 0; + av = 0; if ( !m ) { get_algtree((Obj)f,&av); for ( nalg = 0, tv = av; tv; tv = NEXT(tv), nalg++ ); @@ -2395,6 +2414,86 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe #endif } +void nd_gr_postproc(LIST f,LIST v,int m,struct order_spec *ord,int do_check,LIST *rp) +{ + VL tv,fv,vv,vc,av; + NODE fd,fd0,r,r0,t,x,s,xx,alist; + int e,max,nvar,i; + NDV b; + int ishomo,nalg; + Alg alpha,dp; + P p; + LIST f1,f2; + Obj obj; + NumberField nf; + struct order_spec *ord1; + + get_vars((Obj)f,&fv); pltovl(v,&vv); + for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); + switch ( ord->id ) { + case 1: + if ( ord->nv != nvar ) + error("nd_check : invalid order specification"); + break; + default: + break; + } + nd_nalg = 0; + av = 0; + if ( !m ) { + get_algtree((Obj)f,&av); + for ( nalg = 0, tv = av; tv; tv = NEXT(tv), nalg++ ); + nd_ntrans = nvar; + nd_nalg = nalg; + /* #i -> t#i */ + if ( nalg ) { + preprocess_algcoef(vv,av,ord,f,&ord1,&f1,&alist); + ord = ord1; + f = f1; + } + nvar += nalg; + } + nd_init_ord(ord); + for ( t = BDY(f), max = 0; t; t = NEXT(t) ) + for ( tv = vv; tv; tv = NEXT(tv) ) { + e = getdeg(tv->v,(P)BDY(t)); + max = MAX(e,max); + } + nd_setup_parameters(nvar,max); + ishomo = 1; + for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { + b = (pointer)ptondv(CO,vv,(P)BDY(t)); + if ( ishomo ) + ishomo = ishomo && ndv_ishomo(b); + if ( m ) ndv_mod(m,b); + if ( b ) { NEXTNODE(fd0,fd); BDY(fd) = (pointer)b; } + } + if ( fd0 ) NEXT(fd) = 0; + ndv_setup(m,0,fd0); + for ( x = 0, i = 0; i < nd_psn; i++ ) + x = update_base(x,i); + if ( do_check ) { + x = nd_gb(m,ishomo,1); + if ( !x ) { + *rp = 0; + return; + } + } else { + for ( t = x; t; t = NEXT(t) ) + BDY(t) = (pointer)nd_ps[(int)BDY(t)]; + } + x = ndv_reducebase(x); + x = ndv_reduceall(m,x); + for ( r0 = 0, t = x; t; t = NEXT(t) ) { + NEXTNODE(r0,r); + BDY(r) = ndvtop(m,CO,vv,BDY(t)); + } + if ( r0 ) NEXT(r) = 0; + if ( nalg ) + r0 = postprocess_algcoef(av,alist,r0); + MKLIST(*rp,r0); +} + void nd_gr_trace(LIST f,LIST v,int trace,int homo,struct order_spec *ord,LIST *rp) { VL tv,fv,vv,vc,av; @@ -2468,7 +2567,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru if ( in0 ) NEXT(in) = 0; if ( fd0 ) NEXT(fd) = 0; if ( !ishomo && homo ) { - for ( t = in0, wmax = 0; t; t = NEXT(t) ) { + for ( t = in0, wmax = max; t; t = NEXT(t) ) { c = (NDV)BDY(t); len = LEN(c); for ( a = BDY(c), i = 0; i < len; i++, NMV_ADV(a) ) wmax = MAX(TD(DL(a)),wmax); @@ -3396,7 +3495,7 @@ ND nd_quo(int mod,PGeoBucket bucket,NDV d) ND p,t,r; N tnm; - if ( !p ) return 0; + if ( bucket->m < 0 ) return 0; else { nv = NV(d); mq0 = 0; @@ -3517,7 +3616,7 @@ ND nd_dup(ND p) void ndv_mod(int mod,NDV p) { NMV t,d; - int r; + int r,s,u; int i,len,dlen; Obj gfs; @@ -3539,6 +3638,13 @@ void ndv_mod(int mod,NDV p) if ( r ) { if ( SGN(CQ(t)) < 0 ) r = mod-r; + if ( DN(CQ(t)) ) { + s = rem(DN(CQ(t)),mod); + if ( !s ) + error("ndv_mod : division by 0"); + s = invm(s,mod); + DMAR(r,s,0,mod,u); r = u; + } CM(d) = r; ndl_copy(DL(t),DL(d)); NMV_ADV(d); @@ -5071,6 +5177,8 @@ void nd_det(int mod,MAT f,P *rp) int n,i,j,max,e,nvar,sgn,k0,l0,len0,len,k,l,a; pointer **m; Q mone; + P **w; + P mp,r; NDV **dm; NDV *t,*mi,*mj; NDV d,s,mij,mjj; @@ -5079,6 +5187,8 @@ void nd_det(int mod,MAT f,P *rp) UINT *bound; PGeoBucket bucket; struct order_spec *ord; + Q dq,dt,ds; + N gn,qn,dn0,nm,dn; create_order_spec(0,0,&ord); nd_init_ord(ord); @@ -5086,8 +5196,45 @@ void nd_det(int mod,MAT f,P *rp) if ( f->row != f->col ) error("nd_det : non-square matrix"); n = f->row; - for ( nvar = 0, tv = fv; tv; tv = NEXT(tv), nvar++ ); m = f->body; + for ( nvar = 0, tv = fv; tv; tv = NEXT(tv), nvar++ ); + + if ( !nvar ) { + if ( !mod ) + detp(CO,(P **)m,n,rp); + else { + w = (P **)almat_pointer(n,n); + for ( i = 0; i < n; i++ ) + for ( j = 0; j < n; j++ ) + ptomp(mod,(P)m[i][j],&w[i][j]); + detmp(CO,mod,w,n,&mp); + mptop(mp,rp); + } + return; + } + + if ( !mod ) { + w = (P **)almat_pointer(n,n); + dq = ONE; + for ( i = 0; i < n; i++ ) { + dn0 = ONEN; + for ( j = 0; j < n; j++ ) { + if ( !m[i][j] ) continue; + lgp(m[i][j],&nm,&dn); + gcdn(dn0,dn,&gn); divsn(dn0,gn,&qn); muln(qn,dn,&dn0); + } + if ( !UNIN(dn0) ) { + NTOQ(dn0,1,ds); + for ( j = 0; j < n; j++ ) + mulp(CO,(P)m[i][j],(P)ds,&w[i][j]); + mulq(dq,ds,&dt); dq = dt; + } else + for ( j = 0; j < n; j++ ) + w[i][j] = (P)m[i][j]; + } + m = (pointer **)w; + } + for ( i = 0, max = 0; i < n; i++ ) for ( j = 0; j < n; j++ ) for ( tv = fv; tv; tv = NEXT(tv) ) { @@ -5106,7 +5253,7 @@ void nd_det(int mod,MAT f,P *rp) if ( mod ) ndv_mod(mod,d); chsgnq(ONE,&mone); for ( j = 0, sgn = 1; j < n; j++ ) { - if ( DP_Print ) fprintf(stderr,"j=%d\n",j); + if ( DP_Print ) fprintf(stderr,".",j); for ( i = j; i < n && !dm[i][j]; i++ ); if ( i == n ) { *rp = 0; @@ -5163,12 +5310,17 @@ void nd_det(int mod,MAT f,P *rp) } d = mjj; } + if ( DP_Print ) fprintf(stderr,"\n",k); if ( sgn < 0 ) if ( mod ) ndv_mul_c(mod,d,mod-1); else ndv_mul_c_q(d,mone); - *rp = ndvtop(mod,CO,fv,d); + r = ndvtop(mod,CO,fv,d); + if ( !mod && !UNIQ(dq) ) + divsp(CO,r,(P)dq,rp); + else + *rp = r; } ND ndv_mul_nmv_trunc(int mod,NMV m0,NDV p,UINT *d)