=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.120 retrieving revision 1.129 diff -u -p -r1.120 -r1.129 --- OpenXM_contrib2/asir2000/engine/nd.c 2004/12/07 15:15:52 1.120 +++ OpenXM_contrib2/asir2000/engine/nd.c 2006/04/17 04:35:20 1.129 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.119 2004/12/06 09:29:34 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.128 2005/08/03 06:10:47 noro Exp $ */ #include "nd.h" @@ -8,6 +8,7 @@ NM _nm_free_list; ND _nd_free_list; ND_pairs _ndp_free_list; +static int nd_ntrans; static int nd_nalg; #if 0 static int ndv_alloc; @@ -46,6 +47,7 @@ NumberField get_numberfield(); UINT *nd_det_compute_bound(NDV **dm,int n,int j); void nd_det_reconstruct(NDV **dm,int n,int j,NDV d); int nd_monic(int m,ND *p); +NDV plain_vect_to_ndv_q(Q *mat,int col,UINT *s0vect); void nd_free_private_storage() { @@ -1283,7 +1285,7 @@ int ndv_check_candidate(NODE input,int obpe,int oadv,E NODE t,s; union oNDC dn; - ndv_setup(0,0,cand); + ndv_setup(0,0,cand,0); n = length(cand); /* membercheck : list is a subset of Id(cand) ? */ @@ -1536,7 +1538,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; @@ -1552,7 +1554,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; @@ -1573,6 +1576,7 @@ void do_diagonalize(int sugar,int m) } else nd_ps[i] = nfv; } + return 1; } /* return value = 0 => input is not a GB */ @@ -1598,8 +1602,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); } @@ -1645,7 +1655,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; @@ -1658,7 +1668,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); @@ -1671,7 +1682,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) ) { @@ -1690,6 +1702,7 @@ void do_diagonalize_trace(int sugar,int m) } else nd_ps_trace[i] = nfv; } + return 1; } static struct oEGT eg_invdalg; @@ -1720,7 +1733,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); } @@ -1820,7 +1840,7 @@ NODE ndv_reduceall(int m,NODE f) (int (*)(const void *,const void *))ndv_compare); for ( t = f, i = 0; t; i++, t = NEXT(t) ) BDY(t) = (pointer)w[i]; #endif - ndv_setup(m,0,f); + ndv_setup(m,0,f,0); for ( i = 0; i < n; ) { g = ndvtond(m,nd_ps[i]); g = nd_separate_head(g,&head); @@ -2160,7 +2180,7 @@ int ndv_newps(int m,NDV a,NDV aq) return nd_psn++; } -void ndv_setup(int mod,int trace,NODE f) +void ndv_setup(int mod,int trace,NODE f,int dont_sort) { int i,j,td,len,max; NODE s,s0,f0; @@ -2174,8 +2194,9 @@ void ndv_setup(int mod,int trace,NODE f) for ( nd_psn = 0, s = f; s; s = NEXT(s) ) if ( BDY(s) ) nd_psn++; w = (NDV *)ALLOCA(nd_psn*sizeof(NDV)); for ( i = 0, s = f; s; s = NEXT(s) ) if ( BDY(s) ) w[i++] = BDY(s); - qsort(w,nd_psn,sizeof(NDV), - (int (*)(const void *,const void *))ndv_compare); + if ( !dont_sort ) + qsort(w,nd_psn,sizeof(NDV), + (int (*)(const void *,const void *))ndv_compare); nd_pslen = 2*nd_psn; nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV)); nd_ps_trace = (NDV *)MALLOC(nd_pslen*sizeof(NDV)); @@ -2214,20 +2235,25 @@ void ndv_setup(int mod,int trace,NODE f) struct order_spec *append_block(struct order_spec *spec, int nv,int nalg,int ord); +extern VECT current_dl_weight_vector_obj; +static VECT prev_weight_vector_obj; + void preprocess_algcoef(VL vv,VL av,struct order_spec *ord,LIST f, struct order_spec **ord1p,LIST *f1p,NODE *alistp) { - NODE alist,t,s,r0,r; + NODE alist,t,s,r0,r,arg; VL tv; P poly; DP d; Alg alpha,dp; DAlg inv,da,hc; MP m; - int i,nvar,nalg; + int i,nvar,nalg,n; NumberField nf; LIST f1,f2; struct order_spec *current_spec; + VECT obj,obj0; + Obj tmp; for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++); for ( nalg = 0, tv = av; tv; tv = NEXT(tv), nalg++); @@ -2271,8 +2297,46 @@ void preprocess_algcoef(VL vv,VL av,struct order_spec MKLIST(f1,t); *alistp = alist; algobjtorat(f1,f1p); + + /* creating a new weight vector */ + prev_weight_vector_obj = obj0 = current_dl_weight_vector_obj; + n = nvar+nalg+1; + MKVECT(obj,n); + if ( obj0 && obj0->len == nvar ) + for ( i = 0; i < nvar; i++ ) BDY(obj)[i] = BDY(obj0)[i]; + else + for ( i = 0; i < nvar; i++ ) BDY(obj)[i] = (pointer)ONE; + for ( i = 0; i < nalg; i++ ) BDY(obj)[i+nvar] = 0; + BDY(obj)[n-1] = (pointer)ONE; + arg = mknode(1,obj); + Pdp_set_weight(arg,&tmp); } +NODE postprocess_algcoef(VL av,NODE alist,NODE r) +{ + NODE s,t,u0,u; + P p; + VL tv; + Obj obj,tmp; + NODE arg; + + u0 = 0; + for ( t = r; t; t = NEXT(t) ) { + p = (P)BDY(t); + for ( tv = av, s = alist; tv; tv = NEXT(tv), s = NEXT(s) ) { + substr(CO,0,(Obj)p,tv->v,(Obj)BDY(s),&obj); p = (P)obj; + } + if ( OID(p) == O_P || (OID(p) == O_N && NID((Num)p) != N_A) ) { + NEXTNODE(u0,u); + BDY(u) = (pointer)p; + } + } + arg = mknode(1,prev_weight_vector_obj); + Pdp_set_weight(arg,&tmp); + + return u0; +} + void nd_gr(LIST f,LIST v,int m,int f4,struct order_spec *ord,LIST *rp) { VL tv,fv,vv,vc,av; @@ -2304,9 +2368,11 @@ 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++ ); + nd_ntrans = nvar; nd_nalg = nalg; /* #i -> t#i */ if ( nalg ) { @@ -2332,7 +2398,7 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe if ( b ) { NEXTNODE(fd0,fd); BDY(fd) = (pointer)b; } } if ( fd0 ) NEXT(fd) = 0; - ndv_setup(m,0,fd0); + ndv_setup(m,0,fd0,0); x = f4?nd_f4(m):nd_gb(m,ishomo,0); nd_demand = 0; x = ndv_reducebase(x); @@ -2340,21 +2406,96 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe for ( r0 = 0, t = x; t; t = NEXT(t) ) { NEXTNODE(r0,r); BDY(r) = ndvtop(m,CO,vv,BDY(t)); - if ( nalg ) { - p = BDY(r); - for ( tv = av, s = alist; tv; tv = NEXT(tv), s = NEXT(s) ) { - substr(CO,0,(Obj)p,tv->v,(Obj)BDY(s),&obj); p = (P)obj; - } - BDY(r) = p; - } } if ( r0 ) NEXT(r) = 0; + if ( nalg ) + r0 = postprocess_algcoef(av,alist,r0); MKLIST(*rp,r0); #if 0 fprintf(asir_out,"ndv_alloc=%d\n",ndv_alloc); #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,0); + 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; @@ -2385,6 +2526,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru 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 ) { @@ -2427,7 +2569,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); @@ -2441,7 +2583,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru while ( 1 ) { if ( Demand ) nd_demand = 1; - ndv_setup(m,1,fd0); + ndv_setup(m,1,fd0,0); cand = nd_gb_trace(m,ishomo || homo); if ( !cand ) { /* failure */ @@ -2483,16 +2625,10 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru /* dp->p */ nd_bpe = cbpe; nd_setup_parameters(nd_nvar,0); - for ( r = cand; r; r = NEXT(r) ) { + for ( r = cand; r; r = NEXT(r) ) BDY(r) = (pointer)ndvtop(0,CO,vv,BDY(r)); - if ( nalg ) { - poly = BDY(r); - for ( tv = av, s = alist; tv; tv = NEXT(tv), s = NEXT(s) ) { - substr(CO,0,(Obj)poly,tv->v,(Obj)BDY(s),&obj); poly = (P)obj; - } - BDY(r) = poly; - } - } + if ( nalg ) + cand = postprocess_algcoef(av,alist,cand); MKLIST(*rp,cand); } @@ -3361,7 +3497,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; @@ -3482,7 +3618,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; @@ -3504,6 +3640,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); @@ -3891,7 +4034,8 @@ void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec if ( m ) ndv_mod(m,(NDV)BDY(in)); NEXT(in) = 0; - ndv_setup(m,0,in0); + /* dont sort */ + ndv_setup(m,0,in0,1); nd_psn--; nd_scale=2; while ( 1 ) { @@ -3939,6 +4083,29 @@ int nd_to_vect_q(UINT *s0,int n,ND d,Q *r) return i; } +Q *nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_pair pair) +{ + NM m; + NMV mr; + UINT *d,*t,*s; + NDV p; + int i,j,len; + Q *r; + + m = pair->mul; + d = DL(m); + p = nd_ps[pair->index]; + len = LEN(p); + r = (Q *)CALLOC(n,sizeof(Q)); + t = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); + for ( i = j = 0, s = s0, mr = BDY(p); j < len; j++, NMV_ADV(mr) ) { + ndl_add(d,DL(mr),t); + for ( ; !ndl_equal(t,s); s += nd_wpd, i++ ); + r[i] = CQ(mr); + } + return r; +} + IndArray nm_ind_pair_to_vect_compress(int mod,UINT *s0,int n,NM_ind_pair pair) { NM m; @@ -4184,6 +4351,8 @@ NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea } } +/* for preprocessed vector */ + NDV vect_to_ndv_q(Q *vect,int spcol,int col,int *rhead,UINT *s0vect) { int j,k,len; @@ -4214,6 +4383,36 @@ NDV vect_to_ndv_q(Q *vect,int spcol,int col,int *rhead } } +/* for plain vector */ + +NDV plain_vect_to_ndv_q(Q *vect,int col,UINT *s0vect) +{ + int j,k,len; + UINT *p; + Q c; + NDV r; + NMV mr0,mr; + + for ( j = 0, len = 0; j < col; j++ ) if ( vect[j] ) len++; + if ( !len ) return 0; + else { + mr0 = (NMV)GC_malloc(nmv_adv*len); +#if 0 + ndv_alloc += nmv_adv*len; +#endif + mr = mr0; + p = s0vect; + for ( j = k = 0; j < col; j++, p += nd_wpd, k++ ) + if ( c = vect[k] ) { + if ( DN(c) ) + error("afo"); + ndl_copy(p,DL(mr)); CQ(mr) = c; NMV_ADV(mr); + } + MKNDV(nd_nvar,mr0,len,r); + return r; + } +} + int nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket) { ND_pairs t; @@ -4333,6 +4532,14 @@ NODE nd_f4(int m) for ( r = nflist; r; r = NEXT(r) ) { nf = (NDV)BDY(r); ndv_removecont(m,nf); + if ( !m && nd_nalg ) { + ND nf1; + + nf1 = ndvtond(m,nf); + nd_monic(0,&nf1); + nd_removecont(m,nf1); + nf = ndtondv(m,nf1); + } nh = ndv_newps(m,nf,0); d = update_pairs(d,g,nh); g = update_base(g,nh); @@ -4447,6 +4654,7 @@ NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s return r0; } +#if 0 NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,UINT *s0vect,int col, NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred) { @@ -4514,7 +4722,65 @@ NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,UINT *s0vec } return r0; } +#else +void printm(Q **mat,int row,int col) +{ + int i,j; + printf("["); + for ( i = 0; i < row; i++ ) { + for ( j = 0; j < col; j++ ) { + printexpr(CO,mat[i][j]); printf(" "); + } + printf("]\n"); + } +} +NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,UINT *s0vect,int col, + NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred) +{ + int row,a; + int i,j,rank; + NODE r0,r; + ND_pairs sp; + ND spol; + Q **mat; + int *colstat; + int *sugar; + + row = nsp+nred; + /* make the matrix */ + mat = (Q **)ALLOCA(row*sizeof(Q *)); + sugar = (int *)ALLOCA(row*sizeof(int)); + for ( row = a = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { + nd_sp(0,0,sp,&spol); + if ( !spol ) continue; + mat[row] = (Q *)MALLOC(col*sizeof(Q)); + nd_to_vect_q(s0vect,col,spol,mat[row]); + sugar[row] = SG(spol); + row++; + } + for ( i = 0; i < nred; i++, row++ ) { + mat[row] = nm_ind_pair_to_vect(0,s0vect,col,rvect[i]); + sugar[row] = rvect[i]->sugar; + } + /* elimination */ + colstat = (int *)ALLOCA(col*sizeof(int)); + rank = nd_gauss_elim_q(mat,sugar,row,col,colstat); + r0 = 0; + for ( i = 0; i < rank; i++ ) { + for ( j = 0; j < col; j++ ) if ( mat[i][j] ) break; + if ( j == col ) error("nd_f4_red_q_main : cannot happen"); + if ( rhead[j] ) continue; + NEXTNODE(r0,r); BDY(r) = + (pointer)plain_vect_to_ndv_q(mat[i],col,s0vect); + SG((NDV)BDY(r)) = sugar[i]; + } + if ( r0 ) NEXT(r) = 0; + printf("\n"); + return r0; +} +#endif + FILE *nd_write,*nd_read; void nd_send_int(int a) { @@ -5036,6 +5302,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; @@ -5044,6 +5312,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); @@ -5051,8 +5321,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) ) { @@ -5071,7 +5378,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; @@ -5128,12 +5435,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) @@ -5241,24 +5553,31 @@ UINT *nd_det_compute_bound(NDV **dm,int n,int j) DL nd_separate_d(UINT *d,UINT *trans) { - int n,ntrans,td,i,e; + int n,td,i,e,j; DL a; - n = nd_nvar; ntrans = n-nd_nalg; ndl_zero(trans); td = 0; - for ( i = 0; i < ntrans; i++ ) { + for ( i = 0; i < nd_ntrans; i++ ) { e = GET_EXP(d,i); PUT_EXP(trans,i,e); td += MUL_WEIGHT(e,i); } + if ( nd_ntrans+nd_nalg < nd_nvar ) { + /* homogenized */ + i = nd_nvar-1; + e = GET_EXP(d,i); + PUT_EXP(trans,i,e); + td += MUL_WEIGHT(e,i); + } TD(trans) = td; if ( nd_blockmask) ndl_weight_mask(trans); NEWDL(a,nd_nalg); td = 0; - for ( ; i < n; i++ ) { - e = GET_EXP(d,i); - a->d[i-ntrans] = e; + for ( i = 0; i < nd_nalg; i++ ) { + j = nd_ntrans+i; + e = GET_EXP(d,j); + a->d[i] = e; td += e; } a->td = td; @@ -5286,7 +5605,6 @@ int nd_monic(int mod,ND *p) if ( !(nf = get_numberfield()) ) error("nd_monic : current_numberfield is not set"); - n = nd_nvar; ntrans = n-nd_nalg; /* Q coef -> DAlg coef */ NEWNM(ma0); ma = ma0; m = BDY(*p); @@ -5340,10 +5658,10 @@ int nd_monic(int mod,ND *p) dl = mp->dl; td = TD(DL(m)); ndl_copy(DL(m),DL(mr)); - for ( i = ntrans; i < n; i++ ) { - e = dl->d[i-ntrans]; - PUT_EXP(DL(mr),i,e); - td += MUL_WEIGHT(e,i); + for ( i = 0; i < nd_nalg; i++ ) { + e = dl->d[i]; + PUT_EXP(DL(mr),i+nd_ntrans,e); + td += MUL_WEIGHT(e,i+nd_ntrans); } TD(DL(mr)) = td; if ( nd_blockmask) ndl_weight_mask(DL(mr));