=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.101 retrieving revision 1.102 diff -u -p -r1.101 -r1.102 --- OpenXM_contrib2/asir2000/engine/nd.c 2004/09/14 10:00:26 1.101 +++ OpenXM_contrib2/asir2000/engine/nd.c 2004/09/15 01:43:33 1.102 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.100 2004/09/14 09:25:48 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.101 2004/09/14 10:00:26 noro Exp $ */ #include "nd.h" @@ -3196,11 +3196,12 @@ ND ndv_mul_nm(int mod,NM m0,NDV p) } } -ND nd_quo(ND p,ND d) +ND nd_quo(int mod,ND p,NDV d) { NM mq0,mq; + NMV tm; Q q; - int i,nv,sg; + int i,nv,sg,c,c1,c2; ND t,r; if ( !p ) return 0; @@ -3208,14 +3209,21 @@ ND nd_quo(ND p,ND d) nv = NV(p); sg = SG(p); mq0 = 0; + tm = (NMV)ALLOCA(nmv_adv); while ( p ) { NEXTNM(mq0,mq); ndl_sub(HDL(p),HDL(d),DL(mq)); - divq(HCQ(p),HCQ(d),&q); - chsgnq(q,&CQ(mq)); - t = nd_mul_nm_trunc(mq,d,HDL(d)); - CQ(mq) = q; - p = nd_add(0,p,t); + ndl_copy(DL(mq),DL(tm)); + if ( mod ) { + c1 = invm(HCM(d),mod); c2 = HCM(p); + DMAR(c1,c2,0,mod,c); CM(mq) = c; + CM(tm) = mod-c; + } else { + divq(HCQ(p),HCQ(d),&CQ(mq)); + chsgnq(CQ(mq),&CQ(tm)); + } + t = ndv_mul_nmv_trunc(mod,tm,d,HDL(d)); + p = nd_add(mod,p,t); } NEXT(mq) = 0; for ( i = 0, mq = mq0; mq; mq = NEXT(mq), i++ ); @@ -4647,18 +4655,18 @@ NDV ndv_load(int index) return d; } -void nd_det(MAT f,P *rp) +void nd_det(int mod,MAT f,P *rp) { VL fv,tv; - int n,i,j,max,e,nvar,sgn,k0,l0,len0,len,k,l; + int n,i,j,max,e,nvar,sgn,k0,l0,len0,len,k,l,a; pointer **m; Q mone; - ND **dm; - ND *t,*mi,*mj; - ND d,s,mij,mjj,m1,m2,u; - NDV dv; + NDV **dm; + NDV *t,*mi,*mj; + NDV d,s,mij,mjj; + ND u; + NMV nmv; PGeoBucket bucket; - NM nm; struct order_spec *ord; create_order_spec(0,0,&ord); @@ -4676,11 +4684,15 @@ void nd_det(MAT f,P *rp) max = MAX(e,max); } nd_setup_parameters(nvar,1024); - dm = (ND **)almat_pointer(n,n); + dm = (NDV **)almat_pointer(n,n); for ( i = 0, max = 0; i < n; i++ ) - for ( j = 0; j < n; j++ ) - dm[i][j] = ptond(CO,fv,m[i][j]); - d = ptond(CO,fv,(P)ONE); + for ( j = 0; j < n; j++ ) { + dm[i][j] = ptondv(CO,fv,m[i][j]); + if ( mod ) ndv_mod(mod,dm[i][j]); + if ( dm[i][j] && !LEN(dm[i][j]) ) dm[i][j] = 0; + } + d = ptondv(CO,fv,(P)ONE); + if ( mod ) ndv_mod(mod,d); chsgnq(ONE,&mone); for ( j = 0, sgn = 1; j < n; j++ ) { for ( i = j; i < n && !dm[i][j]; i++ ); @@ -4688,11 +4700,11 @@ void nd_det(MAT f,P *rp) *rp = 0; return; } - k0 = i; l0 = j; len0 = nd_length(dm[k0][l0]); + k0 = i; l0 = j; len0 = LEN(dm[k0][l0]); for ( k = j; k < n; k++ ) for ( l = j; l < n; l++ ) - if ( dm[k][l] && (len = nd_length(dm[k][l])) < len0 ) { - k0 = k; l0 = l; len0 = len; + if ( dm[k][l] && LEN(dm[k][l]) < len0 ) { + k0 = k; l0 = l; len0 = LEN(dm[k][l]); } if ( k0 != j ) { t = dm[j]; dm[j] = dm[k0]; dm[k0] = t; @@ -4706,36 +4718,45 @@ void nd_det(MAT f,P *rp) } for ( i = j+1, mj = dm[j], mjj = mj[j]; i < n; i++ ) { mi = dm[i]; mij = mi[j]; - nd_mul_c_q(mij,mone); + if ( mod ) + ndv_mul_c(mod,mij,mod-1); + else + ndv_mul_c_q(mij,mone); for ( k = j+1; k < n; k++ ) { bucket = create_pbucket(); if ( mi[k] ) - for ( nm = BDY(mjj); nm; nm = NEXT(nm) ) { - u = nd_mul_nm_trunc(nm,mi[k],DL(BDY(d))); - add_pbucket(0,bucket,u); + nmv = BDY(mjj); len = LEN(mjj); + for ( a = 0; a < len; a++, NMV_ADV(nmv) ) { + u = ndv_mul_nmv_trunc(mod,nmv,mi[k],DL(BDY(d))); + add_pbucket(mod,bucket,u); } if ( mj[k] && mij ) { - for ( nm = BDY(mij); nm; nm = NEXT(nm) ) { - u = nd_mul_nm_trunc(nm,mj[k],DL(BDY(d))); - add_pbucket(0,bucket,u); + nmv = BDY(mij); len = LEN(mij); + for ( a = 0; a < len; a++, NMV_ADV(nmv) ) { + u = ndv_mul_nmv_trunc(mod,nmv,mj[k],DL(BDY(d))); + add_pbucket(mod,bucket,u); } } - s = normalize_pbucket(0,bucket); - mi[k] = nd_quo(s,d); + u = normalize_pbucket(mod,bucket); + u = nd_quo(mod,u,d); + mi[k] = ndtondv(mod,u); } } d = mjj; } if ( sgn < 0 ) - nd_mul_c_q(d,mone); - dv = ndtondv(0,d); - *rp = ndvtop(0,CO,fv,dv); + if ( mod ) + ndv_mul_c(mod,d,mod-1); + else + ndv_mul_c_q(d,mone); + *rp = ndvtop(mod,CO,fv,d); } -ND nd_mul_nm_trunc(NM m0,ND p,UINT *d) +ND ndv_mul_nmv_trunc(int mod,NMV m0,NDV p,UINT *d) { NM mr,mr0; - NM m,tnm; + NM tnm; + NMV m; UINT *d0,*dt,*dm; int c,n,td,i,c1,c2,len; Q q; @@ -4743,19 +4764,30 @@ ND nd_mul_nm_trunc(NM m0,ND p,UINT *d) if ( !p ) return 0; else { - n = NV(p); m = BDY(p); + n = NV(p); m = BDY(p); len = LEN(p); d0 = DL(m0); - len = LEN(p); - mr0 = 0; td = TD(d); - q = CQ(m0); + mr0 = 0; NEWNM(tnm); - for ( ; m;m = NEXT(m) ) { - ndl_add(DL(m),d0,DL(tnm)); - if ( ndl_reducible(DL(tnm),d) ) { - NEXTNM(mr0,mr); - mulq(CQ(m),q,&CQ(mr)); - ndl_add(DL(m),d0,DL(mr)); + if ( mod ) { + c = CM(m0); + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + ndl_add(DL(m),d0,DL(tnm)); + if ( ndl_reducible(DL(tnm),d) ) { + NEXTNM(mr0,mr); + c1 = CM(m); DMAR(c1,c,0,mod,c2); CM(mr) = c2; + ndl_add(DL(m),d0,DL(mr)); + } + } + } else { + q = CQ(m0); + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + ndl_add(DL(m),d0,DL(tnm)); + if ( ndl_reducible(DL(tnm),d) ) { + NEXTNM(mr0,mr); + mulq(CQ(m),q,&CQ(mr)); + ndl_add(DL(m),d0,DL(mr)); + } } } if ( !mr0 )