===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/gr,v
retrieving revision 1.10
retrieving revision 1.24
diff -u -p -r1.10 -r1.24
--- OpenXM_contrib2/asir2000/lib/gr 2001/09/06 00:24:07 1.10
+++ OpenXM_contrib2/asir2000/lib/gr 2006/08/09 02:43:38 1.24
@@ -45,8 +45,13 @@
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
*
- * $OpenXM: OpenXM_contrib2/asir2000/lib/gr,v 1.9 2001/09/05 08:09:10 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/lib/gr,v 1.23 2006/07/24 07:31:17 noro Exp $
*/
+
+module gr $
+ /* Empty for now. It will be used in a future. */
+endmodule $
+
extern INIT_COUNT,ITOR_FAIL$
extern REMOTE_MATRIX,REMOTE_NF,REMOTE_VARS$
@@ -126,27 +131,43 @@ def tolex_tl(G0,V,O,W,H)
def tolex(G0,V,O,W)
{
+ Procs = getopt(procs);
+
TM = TE = TNF = 0;
N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
- if ( !ZD )
- error("tolex : ideal is not zero-dimensional!");
- MB = dp_mbase(map(dp_ptod,HM,V));
+ if ( ZD )
+ MB = dp_mbase(map(dp_ptod,HM,V));
+ else
+ MB = 0;
for ( J = 0; ; J++ ) {
M = lprime(J);
if ( !valid_modulus(HM,M) )
continue;
- T0 = time()[0]; GM = tolexm(G0,V,O,W,M); TM += time()[0] - T0;
- dp_ord(2);
- DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
- D = newvect(N); TL = [];
- do
- TL = cons(dp_dtop(dp_vtoe(D),W),TL);
- while ( nextm(D,DL,N) );
- L = npos_check(DL); NPOSV = L[0]; DIM = L[1];
- T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],1)[0];
+ T0 = time()[0];
+ if ( ZD ) {
+ GM = tolexm(G0,V,O,W,M);
+ dp_ord(2);
+ DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
+ D = newvect(N); TL = [];
+ do
+ TL = cons(dp_dtop(dp_vtoe(D),W),TL);
+ while ( nextm(D,DL,N) );
+ } else {
+ GM = dp_gr_mod_main(G0,W,0,M,2);
+ dp_ord(2);
+ for ( T = GM, S = 0; T != []; T = cdr(T) )
+ for ( D = dp_ptod(car(T),V); D; D = dp_rest(D) )
+ S += dp_ht(D);
+ TL = dp_terms(S,V);
+ }
+ TM += time()[0] - T0;
+ T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],ZD)[0];
TNF += time()[0] - T0;
T0 = time()[0];
- R = tolex_main(V,O,NF,GM,M,MB);
+ if ( type(Procs) != -1 )
+ R = tolex_d_main(V,O,NF,GM,M,MB,Procs);
+ else
+ R = tolex_main(V,O,NF,GM,M,MB);
TE += time()[0] - T0;
if ( R ) {
if ( dp_gr_print() )
@@ -316,12 +337,20 @@ def dptov(P,W,MB)
def tolex_main(V,O,NF,GM,M,MB)
{
- DIM = length(MB);
- DV = newvect(DIM);
+ if ( MB ) {
+ PosDim = 0;
+ DIM = length(MB);
+ DV = newvect(DIM);
+ } else
+ PosDim = 1;
for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
S = p_terms(car(T),V,2);
+ if ( PosDim ) {
+ MB = gather_nf_terms(S,NF,V,O);
+ DV = newvect(length(MB));
+ }
dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
- dp_ord(0); NHT = nf_tab_gsl(dp_ptod(LCM*car(S),V),NF);
+ dp_ord(O); NHT = nf_tab_gsl(dp_ptod(LCM*car(S),V),NF);
dptov(NHT[0],DV,MB);
dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
if ( !B )
@@ -338,6 +367,104 @@ def tolex_main(V,O,NF,GM,M,MB)
return SL;
}
+def tolex_d_main(V,O,NF,GM,M,MB,Procs)
+{
+ map(ox_reset,Procs);
+ /* register data in servers */
+ map(ox_cmo_rpc,Procs,"register_data_for_find_base",NF,V,O,MB,M);
+ /* discard return value in stack */
+ map(ox_pop_cmo,Procs);
+ Free = Procs;
+ Busy = [];
+ T = GM;
+ SL = [];
+ while ( T != [] || Busy != [] ){
+ if ( Free == [] || T == [] ) {
+ /* someone is working; wait for data */
+ Ready = ox_select(Busy);
+ Busy = setminus(Busy,Ready);
+ Free = append(Ready,Free);
+ for ( ; Ready != []; Ready = cdr(Ready) )
+ SL = cons(ox_get(car(Ready)),SL);
+ } else {
+ P = car(Free);
+ Free = cdr(Free);
+ Busy = cons(P,Busy);
+ Template = car(T);
+ T = cdr(T);
+ ox_cmo_rpc(P,"find_base",Template);
+ ox_push_cmd(P,262); /* 262 = OX_popCMO */
+ }
+ }
+ return SL;
+}
+
+struct find_base_data { NF,V,O,MB,M,PosDim,DV }$
+extern Find_base$
+
+def register_data_for_find_base(NF,V,O,MB,M)
+{
+ Find_base = newstruct(find_base_data);
+ Find_base->NF = NF;
+ Find_base->V = V;
+ Find_base->O = O;
+ Find_base->M = M;
+ Find_base->MB = MB;
+
+ if ( MB ) {
+ Find_base->PosDim = 0;
+ DIM = length(MB);
+ Find_base->DV = newvect(DIM);
+ } else
+ Find_base->PosDim = 1;
+}
+
+def find_base(S) {
+ NF = Find_base->NF;
+ V = Find_base->V;
+ O = Find_base->O;
+ MB = Find_base->MB;
+ M = Find_base->M;
+ PosDim = Find_base->PosDim;
+ DV = Find_base->DV;
+
+ S = p_terms(S,V,2);
+ if ( PosDim ) {
+ MB = gather_nf_terms(S,NF,V,O);
+ DV = newvect(length(MB));
+ }
+ dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
+ dp_ord(O); NHT = nf_tab_gsl(dp_ptod(car(S),V),NF);
+ dptov(NHT[0],DV,MB);
+ dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
+ if ( !B )
+ return 0;
+ Len = length(S);
+ for ( U = B[1]*car(S), I = 1; I < Len; I++ )
+ U += B[0][I-1]*S[I];
+ R = ptozp(U);
+ return R;
+}
+
+/*
+ * NF = [Pairs,DN]
+ * Pairs = [[NF1,T1],[NF2,T2],...]
+ */
+
+def gather_nf_terms(S,NF,V,O)
+{
+ R = 0;
+ for ( T = S; T != []; T = cdr(T) ) {
+ DT = dp_ptod(car(T),V);
+ for ( U = NF[0]; U != []; U = cdr(U) )
+ if ( car(U)[1] == DT ) {
+ R += tpoly(dp_terms(car(U)[0],V));
+ break;
+ }
+ }
+ return map(dp_ptod,p_terms(R,V,O),V);
+}
+
def reduce_dn(L)
{
NM = L[0]; DN = L[1]; V = vars(NM);
@@ -352,6 +479,9 @@ def minipoly(G0,V,O,P,V0)
if ( !zero_dim(hmlist(G0,V,O),V,O) )
error("tolex : ideal is not zero-dimensional!");
+ Pin = P;
+ P = ptozp(P);
+ CP = sdiv(P,Pin);
G1 = cons(V0-P,G0);
O1 = [[0,1],[O,length(V)]];
V1 = cons(V0,V);
@@ -372,7 +502,7 @@ def minipoly(G0,V,O,P,V0)
TL = cons(V0^J,TL);
NF = gennf(G1,TL,V1,O1,V0,1)[0];
R = tolex_main(V1,O1,NF,[MP],M,MB);
- return R[0];
+ return ptozp(subst(R[0],V0,CP*V0));
}
}
@@ -380,6 +510,15 @@ def minipoly(G0,V,O,P,V0)
def gennf(G,TL,V,O,V0,FLAG)
{
+ F = dp_gr_flags();
+ for ( T = F; T != []; T = cdr(T) ) {
+ Key = car(T); T = cdr(T);
+ if ( Key == "Demand" ) {
+ Dir = car(T); break;
+ }
+ }
+ if ( Dir )
+ return gennf_demand(G,TL,V,O,V0,FLAG,Dir);
N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
@@ -431,6 +570,78 @@ def gennf(G,TL,V,O,V0,FLAG)
return [[map(adj_dn,H,LCM),LCM],PS,GI];
}
+def gennf_demand(G,TL,V,O,V0,FLAG,Dir)
+{
+ N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
+ NTL = length(TL);
+ for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
+ PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
+ }
+ for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
+ DTL = cons(dp_ptod(car(TL),V),DTL);
+ for ( I = Len - 1, GI = []; I >= 0; I-- )
+ GI = cons(I,GI);
+
+ USE_TAB = (FLAG != 0);
+ if ( USE_TAB ) {
+ T0 = time()[0];
+ MB = dp_mbase(HL); DIM = length(MB);
+ U = dp_ptod(V0,V);
+ UTAB = newvect(DIM);
+ for ( I = 0; I < DIM; I++ ) {
+ UTAB[I] = [MB[I],remove_cont(dp_true_nf(GI,U*MB[I],PS,1))];
+ if ( dp_gr_print() )
+ print(".",2);
+ }
+ if ( dp_gr_print() )
+ print("");
+ TTAB = time()[0]-T0;
+ }
+
+ T0 = time()[0];
+ for ( LCM = 1, Index = 0, H = []; DTL != []; Index++ ) {
+ if ( dp_gr_print() )
+ print(".",2);
+ T = car(DTL); DTL = cdr(DTL);
+ if ( L = search_redble(T,H) ) {
+ L = nf_load(Dir,L[0]);
+ DD = dp_subd(T,L[1]);
+ if ( USE_TAB && (DD == U) ) {
+ NF = nf_tab(L[0],UTAB);
+ NF = [NF[0],dp_hc(L[1])*NF[1]*T];
+ } else
+ NF = nf(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
+ } else
+ NF = nf(GI,T,T,PS);
+ NF = remove_cont(NF);
+ nf_save(NF,Dir,Index);
+ H = cons([Index,NF[1]],H);
+ LCM = ilcm(LCM,dp_hc(NF[1]));
+ }
+ TNF = time()[0]-T0;
+ if ( dp_gr_print() )
+ print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
+
+ for ( I = 0; I < NTL; I++ ) {
+ NF = nf_load(Dir,I);
+ NF = adj_dn(NF,LCM);
+ nf_save(NF,Dir,I);
+ }
+ for ( H = [], I = NTL-1; I >= 0; I-- )
+ H = cons(nf_load(Dir,I),H);
+ return [[H,LCM],PS,GI];
+}
+
+def nf_load(Dir,I)
+{
+ return bload(Dir+"/nf"+rtostr(I));
+}
+
+def nf_save(NF,Dir,I)
+{
+ bsave(NF,Dir+"/nf"+rtostr(I));
+}
+
def adj_dn(P,D)
{
return [(idiv(D,dp_hc(P[1])))*P[0],dp_ht(P[1])];
@@ -462,6 +673,8 @@ def vtop(S,L,GSL)
}
}
+/* broken */
+
def leq_nf(TL,NF,LHS,V)
{
TLen = length(NF);
@@ -918,6 +1131,17 @@ def p_true_nf(P,B,V,O) {
return [dp_dtop(L[0],V),L[1]];
}
+def p_nf_mod(P,B,V,O,Mod) {
+ setmod(Mod);
+ dp_ord(O); DP = dp_mod(dp_ptod(P,V),Mod,[]);
+ N = length(B); DB = newvect(N);
+ for ( I = N-1, IL = []; I >= 0; I-- ) {
+ DB[I] = dp_mod(dp_ptod(B[I],V),Mod,[]);
+ IL = cons(I,IL);
+ }
+ return dp_dtop(dp_nf_mod(IL,DP,DB,1,Mod),V);
+}
+
def p_terms(D,V,O)
{
dp_ord(O);
@@ -939,9 +1163,15 @@ def gb_comp(A,B)
LB = length(B);
if ( LA != LB )
return 0;
- A1 = qsort(newvect(LA,A));
- B1 = qsort(newvect(LB,B));
+ A = newvect(LA,A);
+ B = newvect(LB,B);
for ( I = 0; I < LA; I++ )
+ A[I] *= headsgn(A[I]);
+ for ( I = 0; I < LB; I++ )
+ B[I] *= headsgn(B[I]);
+ A1 = qsort(A);
+ B1 = qsort(B);
+ for ( I = 0; I < LA; I++ )
if ( A1[I] != B1[I] && A1[I] != -B1[I] )
break;
return I == LA ? 1 : 0;
@@ -1330,13 +1560,39 @@ def dgr(G,V,O)
Win = "nonhomo";
Lose = P1;
} else {
- Win = "nhomo";
+ Win = "homo";
Lose = P0;
}
ox_reset(Lose);
return [Win,R];
}
+/* competitive Gbase computation : F4 vs. Bucbberger */
+/* P : process list */
+
+def dgrf4mod(G,V,M,O)
+{
+ P = getopt(proc);
+ if ( type(P) == -1 )
+ return dp_f4_mod_main(G,V,M,O);
+ P0 = P[0]; P1 = P[1]; P = [P0,P1];
+ map(ox_reset,P);
+ ox_cmo_rpc(P0,"dp_f4_mod_main",G,V,M,O);
+ ox_cmo_rpc(P1,"dp_gr_mod_main",G,V,0,M,O);
+ map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
+ F = ox_select(P);
+ R = ox_get(F[0]);
+ if ( F[0] == P0 ) {
+ Win = "F4";
+ Lose = P1;
+ } else {
+ Win = "Buchberger";
+ Lose = P0;
+ }
+ ox_reset(Lose);
+ return [Win,R];
+}
+
/* functions for rpc */
def register_matrix(M)
@@ -1413,6 +1669,136 @@ def check_trace(NF,NFIndex,HL)
}
/*
+ * Trace = [Input,[[j1,[[c,i,m,d],...]],[j2,[[...],...]],...]]
+ * if c != 0
+ * g = 0
+ * g = (c*g + m*gi)/d
+ * ...
+ * finally fj = g
+ */
+
+def show_trace(Trace,V)
+{
+ Input = Trace[0];
+ for ( I = 0, T = Input; T != []; T = cdr(T), I++ ) {
+ print("F"+rtostr(I)+"=",0);
+ print(dp_dtop(car(T),V));
+ }
+ Trace = cdr(Trace);
+ for ( T = Trace; T != []; T = cdr(T) ) {
+ HL = car(T);
+ J = car(HL); HL = HL[1];
+ L = length(HL);
+ print("F"+rtostr(J)+"=",0);
+ for ( I = 0; I < L; I++ ) print("(",0);
+ for ( First = 1, S = HL; S != []; S = cdr(S) ) {
+ H = car(S);
+
+ Coeff = H[0];
+ Index = H[1];
+ Monomial = H[2];
+ Denominator = H[3];
+ if ( First ) {
+ if ( Monomial != 1 ) {
+ print("(",0);
+ print(type(Monomial)==9?dp_dtop(Monomial,V):Monomial,0);
+ print(")*",0);
+ }
+ print("F"+rtostr(Index)+")",0);
+ } else {
+ if ( Coeff != 1 ) {
+ print("*(",0); print(Coeff,0); print(")",0);
+ }
+ print("+",0);
+ if ( Monomial != 1 ) {
+ print("(",0);
+ print(type(Monomial)==9?dp_dtop(Monomial,V):Monomial,0);
+ print(")*",0);
+ }
+ print("F"+rtostr(Index)+")",0);
+ if ( Denominator != 1 ) {
+ print("/",0); print(Denominator,0);
+ }
+ }
+ if ( First ) First = 0;
+ }
+ print("");
+ }
+}
+
+def generating_relation(Trace,V)
+{
+ Trace = cdr(Trace);
+ Tab = [];
+ for ( T = Trace; T != []; T = cdr(T) ) {
+ HL = car(T);
+ J = car(HL); HL = HL[1];
+ L = length(HL);
+ LHS = strtov("f"+rtostr(J));
+ Dn = 1;
+ for ( First = 1, S = HL; S != []; S = cdr(S) ) {
+ H = car(S);
+
+ Coeff = H[0];
+ Index = H[1];
+ Monomial = type(H[2])==9?dp_dtop(H[2],V):H[2];
+ Denominator = H[3];
+ F = strtov("f"+rtostr(Index));
+ for ( Z = Tab; Z != []; Z = cdr(Z) )
+ if ( Z[0][0] == F ) break;
+ if ( Z != [] ) Value = Z[0][1];
+ else Value = [F,1];
+ if ( First ) {
+ RHS = Monomial*Value[0];
+ Dn *= Value[1];
+ } else {
+ RHS = RHS*Coeff*Value[1]+Dn*Value[0]*Monomial;
+ Dn = Value[1]*Dn*Denominator;
+ }
+ VVVV = tttttttt;
+ P = ptozp(Dn*VVVV+RHS);
+ RHS = coef(P,0,VVVV);
+ Dn = coef(P,1,VVVV);
+ if ( First ) First = 0;
+ }
+ Tab = cons([LHS,[RHS,Dn]],Tab);
+ }
+ return Tab;
+}
+
+def generating_relation_mod(Trace,V,M)
+{
+ Trace = cdr(Trace);
+ Tab = [];
+ for ( T = Trace; T != []; T = cdr(T) ) {
+ HL = car(T);
+ J = car(HL); HL = HL[1];
+ L = length(HL);
+ LHS = strtov("f"+rtostr(J));
+ Dn = 1;
+ for ( First = 1, S = HL; S != []; S = cdr(S) ) {
+ H = car(S);
+
+ Coeff = H[0];
+ Index = H[1];
+ Monomial = type(H[2])==9?dp_dtop(H[2],V):H[2];
+ F = strtov("f"+rtostr(Index));
+ for ( Z = Tab; Z != []; Z = cdr(Z) )
+ if ( Z[0][0] == F ) break;
+ if ( Z != [] ) Value = Z[0][1];
+ else Value = F;
+ if ( First ) {
+ RHS = (Monomial*Value)%M;
+ } else {
+ RHS = ((RHS*Coeff+Value*Monomial)*inv(H[3],M))%M;
+ }
+ if ( First ) First = 0;
+ }
+ Tab = cons([LHS,RHS],Tab);
+ }
+ return Tab;
+}
+/*
* realloc NFArray so that it can hold * an element as NFArray[Ind].
*/
@@ -1516,5 +1902,74 @@ def compute_coef_by_trace(I,Tr,Coef)
CI = CT;
}
Coef[I] = CI;
+}
+
+extern Gbcheck_DP,Gbcheck_IL$
+
+def register_data_for_gbcheck(DPL)
+{
+ for ( IL = [], I = length(DPL)-1; I >= 0; I-- )
+ IL = cons(I,IL);
+ Gbcheck_DP = newvect(length(DPL),DPL);
+ Gbcheck_IL = IL;
+}
+
+def sp_nf_for_gbcheck(Pair)
+{
+ SP = dp_sp(Gbcheck_DP[Pair[0]],Gbcheck_DP[Pair[1]]);
+ return dp_nf(Gbcheck_IL,SP,Gbcheck_DP,1);
+}
+
+def gbcheck(B,V,O)
+{
+ dp_ord(O);
+ D = map(dp_ptod,B,V);
+ L = dp_gr_checklist(D,length(V));
+ DP = L[0]; Plist = L[1];
+ for ( IL = [], I = size(DP)[0]-1; I >= 0; I-- )
+ IL = cons(I,IL);
+ Procs = getopt(proc);
+ if ( type(Procs) == 4 ) {
+ map(ox_reset,Procs);
+ /* register DP in servers */
+ map(ox_cmo_rpc,Procs,"register_data_for_gbcheck",vtol(DP));
+ /* discard return value in stack */
+ map(ox_pop_cmo,Procs);
+ Free = Procs;
+ Busy = [];
+ T = Plist;
+ while ( T != [] || Busy != [] ){
+ if ( Free == [] || T == [] ) {
+ /* someone is working; wait for data */
+ Ready = ox_select(Busy);
+ Busy = setminus(Busy,Ready);
+ Free = append(Ready,Free);
+ for ( ; Ready != []; Ready = cdr(Ready) ) {
+ if ( ox_get(car(Ready)) ) {
+ map(ox_reset,Procs);
+ return 0;
+ }
+ }
+ } else {
+ P = car(Free);
+ Free = cdr(Free);
+ Busy = cons(P,Busy);
+ Pair = car(T);
+ T = cdr(T);
+ ox_cmo_rpc(P,"sp_nf_for_gbcheck",Pair);
+ ox_push_cmd(P,262); /* 262 = OX_popCMO */
+ }
+ }
+ map(ox_reset,Procs);
+ return 1;
+ } else {
+ for ( T = Plist; T != []; T = cdr(T) ) {
+ Pair = T[0];
+ SP = dp_sp(DP[Pair[0]],DP[Pair[1]]);
+ if ( dp_nf(IL,SP,DP,1) )
+ return 0;
+ }
+ return 1;
+ }
}
end$