File: [local] / OpenXM_contrib2 / asir2000 / lib / gr (download)
Revision 1.10, Thu Sep 6 00:24:07 2001 UTC (23 years ago) by noro
Branch: MAIN
Changes since 1.9: +10 -5
lines
Fixed a bug in tracetogen().
Changed the argument of tracetogen().
|
/*
* Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
* All rights reserved.
*
* FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
* non-exclusive and royalty-free license to use, copy, modify and
* redistribute, solely for non-commercial and non-profit purposes, the
* computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
* conditions of this Agreement. For the avoidance of doubt, you acquire
* only a limited right to use the SOFTWARE hereunder, and FLL or any
* third party developer retains all rights, including but not limited to
* copyrights, in and to the SOFTWARE.
*
* (1) FLL does not grant you a license in any way for commercial
* purposes. You may use the SOFTWARE only for non-commercial and
* non-profit purposes only, such as academic, research and internal
* business use.
* (2) The SOFTWARE is protected by the Copyright Law of Japan and
* international copyright treaties. If you make copies of the SOFTWARE,
* with or without modification, as permitted hereunder, you shall affix
* to all such copies of the SOFTWARE the above copyright notice.
* (3) An explicit reference to this SOFTWARE and its copyright owner
* shall be made on your publication or presentation in any form of the
* results obtained by use of the SOFTWARE.
* (4) In the event that you modify the SOFTWARE, you shall notify FLL by
* e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
* for such modification or the source code of the modified part of the
* SOFTWARE.
*
* THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
* MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
* EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
* FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
* RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
* MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
* UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
* OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
* DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
* DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
* ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
* FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
* DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
* SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
* OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
* 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.10 2001/09/06 00:24:07 noro Exp $
*/
extern INIT_COUNT,ITOR_FAIL$
extern REMOTE_MATRIX,REMOTE_NF,REMOTE_VARS$
#define MAX(a,b) ((a)>(b)?(a):(b))
#define HigherDim 0
#define ZeroDim 1
#define MiniPoly 2
/* toplevel functions for Groebner basis computation */
def gr(B,V,O)
{
G = dp_gr_main(B,V,0,1,O);
return G;
}
def hgr(B,V,O)
{
G = dp_gr_main(B,V,1,1,O);
return G;
}
def gr_mod(B,V,O,M)
{
G = dp_gr_mod_main(B,V,0,M,O);
return G;
}
def hgr_mod(B,V,O,M)
{
G = dp_gr_mod_main(B,V,1,M,O);
return G;
}
/* toplevel functions for change-of-ordering */
def lex_hensel(B,V,O,W,H)
{
G = dp_gr_main(B,V,H,1,O);
return tolex(G,V,O,W);
}
def lex_hensel_gsl(B,V,O,W,H)
{
G = dp_gr_main(B,V,H,1,O);
return tolex_gsl(G,V,O,W);
}
def gr_minipoly(B,V,O,P,V0,H)
{
G = dp_gr_main(B,V,H,1,O);
return minipoly(G,V,O,P,V0);
}
def lex_tl(B,V,O,W,H)
{
G = dp_gr_main(B,V,H,1,O);
return tolex_tl(G,V,O,W,H);
}
def tolex_tl(G0,V,O,W,H)
{
N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
for ( I = 0; ; I++ ) {
M = lprime(I);
if ( !valid_modulus(HM,M) )
continue;
if ( ZD ) {
if ( G3 = dp_gr_main(G0,W,H,-M,3) )
for ( J = 0; ; J++ )
if ( G2 = dp_gr_main(G3,W,0,-lprime(J),2) )
return G2;
} else if ( G2 = dp_gr_main(G0,W,H,-M,2) )
return G2;
}
}
def tolex(G0,V,O,W)
{
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));
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];
TNF += time()[0] - T0;
T0 = time()[0];
R = tolex_main(V,O,NF,GM,M,MB);
TE += time()[0] - T0;
if ( R ) {
if ( dp_gr_print() )
print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
return R;
}
}
}
def tolex_gsl(G0,V,O,W)
{
TM = TE = TNF = 0;
N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
MB = dp_mbase(map(dp_ptod,HM,V));
if ( !ZD )
error("tolex_gsl : ideal is not zero-dimensional!");
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];
if ( NPOSV >= 0 ) {
V0 = W[NPOSV];
T0 = time()[0]; NFL = gennf(G0,TL,V,O,V0,1);
TNF += time()[0] - T0;
T0 = time()[0];
R = tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB);
TE += time()[0] - T0;
} else {
T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],1)[0];
TNF += time()[0] - T0;
T0 = time()[0];
R = tolex_main(V,O,NF,GM,M,MB);
TE += time()[0] - T0;
}
if ( R ) {
if ( dp_gr_print() )
print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
return R;
}
}
}
def termstomat(NF,TERMS,MB,MOD)
{
DN = NF[1];
NF = NF[0];
N = length(MB);
M = length(TERMS);
MAT = newmat(N,M);
W = newvect(N);
Len = length(NF);
for ( I = 0; I < M; I++ ) {
T = TERMS[I];
for ( K = 0; K < Len; K++ )
if ( T == NF[K][1] )
break;
dptov(NF[K][0],W,MB);
for ( J = 0; J < N; J++ )
MAT[J][I] = W[J];
}
return [henleq_prep(MAT,MOD),DN];
}
def tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB)
{
NF = NFL[0]; PS = NFL[1]; GI = NFL[2];
V0 = W[NPOSV]; N = length(W);
DIM = length(MB);
DV = newvect(DIM);
TERMS = gather_terms(GM,W,M,NPOSV);
Len = length(TERMS);
dp_ord(O); RHS = termstomat(NF,map(dp_ptod,TERMS,V),MB,M);
for ( T = GM; T != []; T = cdr(T) )
if ( vars(car(T)) == [V0] )
break;
dp_ord(0); NHT = nf_tab_gsl(dp_ptod(V0^deg(car(T),V0),V),NF);
dptov(NHT[0],DV,MB);
B = hen_ttob_gsl([DV,NHT[1]],RHS,TERMS,M);
if ( !B )
return 0;
for ( I = 0, U = B[1]*V0^deg(car(T),V0); I < Len; I++ )
U += B[0][I]*TERMS[I];
DN0 = diff(U,V0);
dp_ord(O); DN0NF = nf_tab_gsl(dp_ptod(DN0,V),NF);
SL = [[V0,U,DN0]];
for ( I = N-1, LCM = 1; I >= 0; I-- ) {
if ( I == NPOSV )
continue;
V1 = W[I];
dp_ord(O); L = nf(GI,DN0NF[0]*dp_ptod(-LCM*V1,V),DN0NF[1],PS);
L = remove_cont(L);
dptov(L[0],DV,MB);
dp_ord(O); B = hen_ttob_gsl([DV,L[1]],RHS,TERMS,M);
if ( !B )
return 0;
for ( K = 0, R = 0; K < Len; K++ )
R += B[0][K]*TERMS[K];
LCM *= B[1];
SL = cons(cons(V1,[R,LCM]),SL);
if ( dp_gr_print() )
print(["DN",B[1]]);
}
return SL;
}
def hen_ttob_gsl(LHS,RHS,TERMS,M)
{
LDN = LHS[1]; RDN = RHS[1]; LCM = ilcm(LDN,RDN);
L1 = idiv(LCM,LDN); R1 = idiv(LCM,RDN);
T0 = time()[0];
S = henleq_gsl(RHS[0],LHS[0]*L1,M);
if ( dp_gr_print() )
print(["henleq_gsl",time()[0]-T0]);
N = length(TERMS);
return [S[0],S[1]*R1];
}
def gather_terms(GM,W,M,NPOSV)
{
N = length(W); V0 = W[NPOSV];
for ( T = GM; T != []; T = cdr(T) ) {
if ( vars(car(T)) == [V0] )
break;
}
U = car(T); DU = diff(U,V0);
R = tpoly(cdr(p_terms(U,W,2)));
for ( I = 0; I < N; I++ ) {
if ( I == NPOSV )
continue;
V1 = W[I];
for ( T = GM; T != []; T = cdr(T) )
if ( member(V1,vars(car(T))) )
break;
P = car(T);
R += tpoly(p_terms(srem(DU*coef(P,0,V1),U,M),W,2));
}
return p_terms(R,W,2);
}
def tpoly(L)
{
for ( R = 0; L != []; L = cdr(L) )
R += car(L);
return R;
}
def dptov(P,W,MB)
{
N = size(W)[0];
for ( I = 0; I < N; I++ )
W[I] = 0;
for ( I = 0, S = MB; P; P = dp_rest(P) ) {
HM = dp_hm(P); C = dp_hc(HM); T = dp_ht(HM);
for ( ; T != car(S); S = cdr(S), I++ );
W[I] = C;
I++; S = cdr(S);
}
}
def tolex_main(V,O,NF,GM,M,MB)
{
DIM = length(MB);
DV = newvect(DIM);
for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
S = p_terms(car(T),V,2);
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);
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);
LCM *= B[1];
for ( U = LCM*car(S), I = 1; I < Len; I++ )
U += B[0][I-1]*S[I];
R = ptozp(U);
SL = cons(R,SL);
if ( dp_gr_print() )
print(["DN",B[1]]);
}
return SL;
}
def reduce_dn(L)
{
NM = L[0]; DN = L[1]; V = vars(NM);
T = remove_cont([dp_ptod(NM,V),DN]);
return [dp_dtop(T[0],V),T[1]];
}
/* a function for computation of minimal polynomial */
def minipoly(G0,V,O,P,V0)
{
if ( !zero_dim(hmlist(G0,V,O),V,O) )
error("tolex : ideal is not zero-dimensional!");
G1 = cons(V0-P,G0);
O1 = [[0,1],[O,length(V)]];
V1 = cons(V0,V);
W = append(V,[V0]);
N = length(V1);
dp_ord(O1);
HM = hmlist(G1,V1,O1);
MB = dp_mbase(map(dp_ptod,HM,V1));
dp_ord(O);
for ( J = 0; ; J++ ) {
M = lprime(J);
if ( !valid_modulus(HM,M) )
continue;
MP = minipolym(G0,V,O,P,V0,M);
for ( D = deg(MP,V0), TL = [], J = 0; J <= D; J++ )
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];
}
}
/* subroutines */
def gennf(G,TL,V,O,V0,FLAG)
{
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);
}
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);
T = car(DTL); DTL = cdr(DTL);
H = [nf(GI,T,T,PS)];
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; DTL != []; ) {
if ( dp_gr_print() )
print(".",2);
T = car(DTL); DTL = cdr(DTL);
if ( L = search_redble(T,H) ) {
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);
H = cons(NF,H);
LCM = ilcm(LCM,dp_hc(NF[1]));
}
TNF = time()[0]-T0;
if ( dp_gr_print() )
print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
return [[map(adj_dn,H,LCM),LCM],PS,GI];
}
def adj_dn(P,D)
{
return [(idiv(D,dp_hc(P[1])))*P[0],dp_ht(P[1])];
}
def hen_ttob(T,NF,LHS,V,MOD)
{
if ( length(T) == 1 )
return car(T);
T0 = time()[0]; M = etom(leq_nf(T,NF,LHS,V)); TE = time()[0] - T0;
T0 = time()[0]; U = henleq(M,MOD); TH = time()[0] - T0;
if ( dp_gr_print() ) {
print("(etom="+rtostr(TE)+" hen="+rtostr(TH)+")");
}
return U ? vtop(T,U,LHS) : 0;
}
def vtop(S,L,GSL)
{
U = L[0]; H = L[1];
if ( GSL ) {
for ( A = 0, I = 0; S != []; S = cdr(S), I++ )
A += U[I]*car(S);
return [A,H];
} else {
for ( A = H*car(S), S = cdr(S), I = 0; S != []; S = cdr(S), I++ )
A += U[I]*car(S);
return ptozp(A);
}
}
def leq_nf(TL,NF,LHS,V)
{
TLen = length(NF);
T = newvect(TLen); M = newvect(TLen);
for ( I = 0; I < TLen; I++ ) {
T[I] = dp_ht(NF[I][1]);
M[I] = dp_hc(NF[I][1]);
}
Len = length(TL); INDEX = newvect(Len); COEF = newvect(Len);
for ( L = TL, J = 0; L != []; L = cdr(L), J++ ) {
D = dp_ptod(car(L),V);
for ( I = 0; I < TLen; I++ )
if ( D == T[I] )
break;
INDEX[J] = I; COEF[J] = strtov("u"+rtostr(J));
}
if ( !LHS ) {
COEF[0] = 1; NM = 0; DN = 1;
} else {
NM = LHS[0]; DN = LHS[1];
}
for ( J = 0, S = -NM; J < Len; J++ ) {
DNJ = M[INDEX[J]];
GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD;
S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J];
DN *= CS;
}
for ( D = S, E = []; D; D = dp_rest(D) )
E = cons(dp_hc(D),E);
BOUND = LHS ? 0 : 1;
for ( I = Len - 1, W = []; I >= BOUND; I-- )
W = cons(COEF[I],W);
return [E,W];
}
def nf_tab(F,TAB)
{
for ( NM = 0, DN = 1, I = 0; F; F = dp_rest(F) ) {
T = dp_ht(F);
for ( ; TAB[I][0] != T; I++);
NF = TAB[I][1]; N = NF[0]; D = NF[1];
G = igcd(DN,D); DN1 = idiv(DN,G); D1 = idiv(D,G);
NM = D1*NM + DN1*dp_hc(F)*N; DN *= D1;
}
return [NM,DN];
}
def nf_tab_gsl(A,NF)
{
DN = NF[1];
NF = NF[0];
TLen = length(NF);
for ( R = 0; A; A = dp_rest(A) ) {
HM = dp_hm(A); C = dp_hc(HM); T = dp_ht(HM);
for ( I = 0; I < TLen; I++ )
if ( NF[I][1] == T )
break;
R += C*NF[I][0];
}
return remove_cont([R,DN]);
}
def redble(D1,D2,N)
{
for ( I = 0; I < N; I++ )
if ( D1[I] > D2[I] )
break;
return I == N ? 1 : 0;
}
def tolexm(G,V,O,W,M)
{
N = length(V); Len = length(G);
dp_ord(O); setmod(M); PS = newvect(Len);
for ( I = 0, T = G; T != []; T = cdr(T), I++ )
PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
for ( I = Len-1, HL = []; I >= 0; I-- )
HL = cons(dp_ht(PS[I]),HL);
G2 = tolexm_main(PS,HL,V,W,M,ZeroDim);
L = map(dp_dtop,G2,V);
return L;
}
def tolexm_main(PS,HL,V,W,M,FLAG)
{
N = length(W); D = newvect(N); Len = size(PS)[0];
for ( I = Len - 1, GI = []; I >= 0; I-- )
GI = cons(I,GI);
MB = dp_mbase(HL); DIM = length(MB);
U = dp_mod(dp_ptod(W[N-1],V),M,[]);
UTAB = newvect(DIM);
for ( I = 0; I < DIM; I++ ) {
if ( dp_gr_print() )
print(".",2);
UTAB[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
}
if ( dp_gr_print() )
print("");
T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
H = G = [[T,T]];
DL = []; G2 = [];
TNF = 0;
while ( 1 ) {
if ( dp_gr_print() )
print(".",2);
S = nextm(D,DL,N);
if ( !S )
break;
T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
T0 = time()[0];
if ( L = search_redble(T,H) ) {
DD = dp_mod(dp_subd(T,L[1]),M,[]);
if ( DD == U )
NT = dp_nf_tab_mod(L[0],UTAB,M);
else
NT = dp_nf_mod(GI,L[0]*DD,PS,1,M);
} else
NT = dp_nf_mod(GI,T,PS,1,M);
TNF += time()[0] - T0;
H = cons([NT,T],H);
T0 = time()[0];
L = dp_lnf_mod([NT,T],G,M); N1 = L[0]; N2 = L[1];
TLNF += time()[0] - T0;
if ( !N1 ) {
G2 = cons(N2,G2);
if ( FLAG == MiniPoly )
break;
D1 = newvect(N);
for ( I = 0; I < N; I++ )
D1[I] = D[I];
DL = cons(D1,DL);
} else
G = insert(G,L);
}
if ( dp_gr_print() )
print("tolexm(nfm="+rtostr(TNF)+" lnfm="+rtostr(TLNF)+")");
return G2;
}
def minipolym(G,V,O,P,V0,M)
{
N = length(V); Len = length(G);
dp_ord(O); setmod(M); PS = newvect(Len);
for ( I = 0, T = G; T != []; T = cdr(T), I++ )
PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
for ( I = Len-1, HL = []; I >= 0; I-- )
HL = cons(dp_ht(PS[I]),HL);
for ( I = Len - 1, GI = []; I >= 0; I-- )
GI = cons(I,GI);
MB = dp_mbase(HL); DIM = length(MB); UT = newvect(DIM);
U = dp_mod(dp_ptod(P,V),M,[]);
for ( I = 0; I < DIM; I++ )
UT[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
T = dp_mod(<<0>>,M,[]); TT = dp_mod(dp_ptod(1,V),M,[]);
G = H = [[TT,T]]; TNF = TLNF = 0;
for ( I = 1; ; I++ ) {
T = dp_mod(<<I>>,M,[]);
T0 = time()[0]; NT = dp_nf_tab_mod(H[0][0],UT,M); TNF += time()[0] - T0;
H = cons([NT,T],H);
T0 = time()[0]; L = dp_lnf_mod([NT,T],G,M); TLNF += time()[0] - T0;
if ( !L[0] ) {
if ( dp_gr_print() ) print(["nfm",TNF,"lnfm",TLNF]);
return dp_dtop(L[1],[V0]);
} else
G = insert(G,L);
}
}
def nextm(D,DL,N)
{
for ( I = N-1; I >= 0; ) {
D[I]++;
for ( T = DL; T != []; T = cdr(T) )
if ( car(T) == D )
return 1;
else if ( redble(car(T),D,N) )
break;
if ( T != [] ) {
for ( J = N-1; J >= I; J-- )
D[J] = 0;
I--;
} else
break;
}
if ( I < 0 )
return 0;
else
return 1;
}
def search_redble(T,G)
{
for ( ; G != []; G = cdr(G) )
if ( dp_redble(T,car(G)[1]) )
return car(G);
return 0;
}
def insert(G,A)
{
if ( G == [] )
return [A];
else if ( dp_ht(car(A)) > dp_ht(car(car(G))) )
return cons(A,G);
else
return cons(car(G),insert(cdr(G),A));
}
#if 0
def etom(L) {
E = L[0]; W = L[1];
LE = length(E); LW = length(W);
M = newmat(LE,LW+1);
for(J=0;J<LE;J++) {
for ( T = E[J]; T && (type(T) == 2); )
for ( V = var(T), I = 0; I < LW; I++ )
if ( V == W[I] ) {
M[J][I] = coef(T,1,V);
T = coef(T,0,V);
}
M[J][LW] = T;
}
return M;
}
#endif
def etom(L) {
E = L[0]; W = L[1];
LE = length(E); LW = length(W);
M = newmat(LE,LW+1);
for(J=0;J<LE;J++) {
for ( I = 0, T = E[J]; I < LW; I++ ) {
M[J][I] = coef(T,1,W[I]); T = coef(T,0,W[I]);
}
M[J][LW] = T;
}
return M;
}
def calcb_old(M) {
N = 2*M;
T = gr_sqrt(N);
if ( T^2 <= N && N < (T+1)^2 )
return idiv(T,2);
else
error("afo");
}
def calcb_special(PK,P,K) { /* PK = P^K */
N = 2*PK;
T = sqrt_special(N,2,P,K);
if ( T^2 <= N && N < (T+1)^2 )
return idiv(T,2);
else
error("afo");
}
def sqrt_special(A,C,M,K) { /* A = C*M^K */
L = idiv(K,2); B = M^L;
if ( K % 2 )
C *= M;
D = 2^K; X = idiv((gr_sqrt(C*D^2)+1)*B,D)+1;
while ( 1 )
if ( (Y = X^2) <= A )
return X;
else
X = idiv(A + Y,2*X);
}
def gr_sqrt(A) {
for ( J = 0, T = A; T >= 2^27; J++ ) {
T = idiv(T,2^27)+1;
}
for ( I = 0; T >= 2; I++ ) {
S = idiv(T,2);
if ( T = S+S )
T = S;
else
T = S+1;
}
X = (2^27)^idiv(J,2)*2^idiv(I,2);
while ( 1 ) {
if ( (Y=X^2) < A )
X += X;
else if ( Y == A )
return X;
else
break;
}
while ( 1 )
if ( (Y = X^2) <= A )
return X;
else
X = idiv(A + Y,2*X);
}
#define ABS(a) ((a)>=0?(a):(-a))
def inttorat_asir(C,M,B)
{
if ( M < 0 )
M = -M;
C %= M;
if ( C < 0 )
C += M;
U1 = 0; U2 = M; V1 = 1; V2 = C;
while ( V2 >= B ) {
L = iqr(U2,V2); Q = L[0]; R2 = L[1];
R1 = U1 - Q*V1;
U1 = V1; U2 = V2;
V1 = R1; V2 = R2;
}
if ( ABS(V1) >= B )
return 0;
else
if ( V1 < 0 )
return [-V2,-V1];
else
return [V2,V1];
}
def intvtoratv(V,M,B) {
if ( !B )
B = 1;
N = size(V)[0];
W = newvect(N);
if ( ITOR_FAIL >= 0 ) {
if ( V[ITOR_FAIL] ) {
T = inttorat(V[ITOR_FAIL],M,B);
if ( !T ) {
if ( dp_gr_print() ) {
print("F",2);
}
return 0;
}
}
}
for ( I = 0, DN = 1; I < N; I++ )
if ( V[I] ) {
T = inttorat((V[I]*DN) % M,M,B);
if ( !T ) {
ITOR_FAIL = I;
if ( dp_gr_print() ) {
#if 0
print("intvtoratv : failed at I = ",0); print(ITOR_FAIL);
#endif
print("F",2);
}
return 0;
} else {
for( J = 0; J < I; J++ )
W[J] *= T[1];
W[I] = T[0]; DN *= T[1];
}
}
return [W,DN];
}
def nf(B,G,M,PS)
{
for ( D = 0; G; ) {
for ( U = 0, L = B; L != []; L = cdr(L) ) {
if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
GCD = igcd(dp_hc(G),dp_hc(R));
CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
U = CG*G-dp_subd(G,R)*CR*R;
if ( !U )
return [D,M];
D *= CG; M *= CG;
break;
}
}
if ( U )
G = U;
else {
D += dp_hm(G); G = dp_rest(G);
}
}
return [D,M];
}
def remove_cont(L)
{
if ( type(L[1]) == 1 ) {
T = remove_cont([L[0],L[1]*<<0>>]);
return [T[0],dp_hc(T[1])];
} else if ( !L[0] )
return [0,dp_ptozp(L[1])];
else if ( !L[1] )
return [dp_ptozp(L[0]),0];
else {
A0 = dp_ptozp(L[0]); A1 = dp_ptozp(L[1]);
C0 = idiv(dp_hc(L[0]),dp_hc(A0)); C1 = idiv(dp_hc(L[1]),dp_hc(A1));
GCD = igcd(C0,C1); M0 = idiv(C0,GCD); M1 = idiv(C1,GCD);
return [M0*A0,M1*A1];
}
}
def union(A,B)
{
for ( T = B; T != []; T = cdr(T) )
A = union1(A,car(T));
return A;
}
def union1(A,E)
{
if ( A == [] )
return [E];
else if ( car(A) == E )
return A;
else
return cons(car(A),union1(cdr(A),E));
}
def setminus(A,B) {
for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
for ( S = B, M = car(T); S != []; S = cdr(S) )
if ( car(S) == M )
break;
if ( S == [] )
R = cons(M,R);
}
return R;
}
def member(A,L) {
for ( ; L != []; L = cdr(L) )
if ( A == car(L) )
return 1;
return 0;
}
/* several functions for computation of normal forms etc. */
def p_nf(P,B,V,O) {
dp_ord(O); DP = dp_ptod(P,V);
N = length(B); DB = newvect(N);
for ( I = N-1, IL = []; I >= 0; I-- ) {
DB[I] = dp_ptod(B[I],V);
IL = cons(I,IL);
}
return dp_dtop(dp_nf(IL,DP,DB,1),V);
}
def p_true_nf(P,B,V,O) {
dp_ord(O); DP = dp_ptod(P,V);
N = length(B); DB = newvect(N);
for ( I = N-1, IL = []; I >= 0; I-- ) {
DB[I] = dp_ptod(B[I],V);
IL = cons(I,IL);
}
L = dp_true_nf(IL,DP,DB,1);
return [dp_dtop(L[0],V),L[1]];
}
def p_terms(D,V,O)
{
dp_ord(O);
for ( L = [], T = dp_ptod(D,V); T; T = dp_rest(T) )
L = cons(dp_dtop(dp_ht(T),V),L);
return reverse(L);
}
def dp_terms(D,V)
{
for ( L = [], T = D; T; T = dp_rest(T) )
L = cons(dp_dtop(dp_ht(T),V),L);
return reverse(L);
}
def gb_comp(A,B)
{
LA = length(A);
LB = length(B);
if ( LA != LB )
return 0;
A1 = qsort(newvect(LA,A));
B1 = qsort(newvect(LB,B));
for ( I = 0; I < LA; I++ )
if ( A1[I] != B1[I] && A1[I] != -B1[I] )
break;
return I == LA ? 1 : 0;
}
def zero_dim(G,V,O) {
dp_ord(O);
HL = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
for ( L = []; HL != []; HL = cdr(HL) )
if ( length(vars(car(HL))) == 1 )
L = cons(car(HL),L);
return length(vars(L)) == length(V) ? 1 : 0;
}
def hmlist(G,V,O) {
dp_ord(O);
return map(dp_dtop,map(dp_hm,map(dp_ptod,G,V)),V);
}
def valid_modulus(HL,M) {
V = vars(HL);
for ( T = HL; T != []; T = cdr(T) )
if ( !dp_mod(dp_ptod(car(T),V),M,[]) )
break;
return T == [] ? 1 : 0;
}
def npos_check(DL) {
N = size(car(DL))[0];
if ( length(DL) != N )
return [-1,0];
D = newvect(N);
for ( I = 0; I < N; I++ ) {
for ( J = 0; J < N; J++ )
D[J] = 0;
D[I] = 1;
for ( T = DL; T != []; T = cdr(T) )
if ( D == car(T) )
break;
if ( T != [] )
DL = setminus(DL,[car(T)]);
}
if ( length(DL) != 1 )
return [-1,0];
U = car(DL);
for ( I = 0, J = 0, I0 = -1; I < N; I++ )
if ( U[I] ) {
I0 = I; J++;
}
if ( J != 1 )
return [-1,0];
else
return [I0,U[I0]];
}
def mult_mat(L,TAB,MB)
{
A = L[0]; DN0 = L[1];
for ( NM = 0, DN = 1, I = 0; A; A = dp_rest(A) ) {
H = dp_ht(A);
for ( ; MB[I] != H; I++ );
NM1 = TAB[I][0]; DN1 = TAB[I][1]; I++;
GCD = igcd(DN,DN1); C = DN1/GCD; C1 = DN/GCD;
NM = C*NM + C1*dp_hc(A)*NM1;
DN *= C;
}
Z=remove_cont([NM,DN*DN0]);
return Z;
}
def sepm(MAT)
{
S = size(MAT); N = S[0]; M = S[1]-1;
A = newmat(N,M); B = newvect(N);
for ( I = 0; I < N; I++ )
for ( J = 0, T1 = MAT[I], T2 = A[I]; J < M; J++ )
T2[J] = T1[J];
for ( I = 0; I < N; I++ )
B[I] = MAT[I][M];
return [A,B];
}
def henleq(M,MOD)
{
SIZE = size(M); ROW = SIZE[0]; COL = SIZE[1];
W = newvect(COL);
L = sepm(M); A = L[0]; B = L[1];
COUNT = INIT_COUNT?INIT_COUNT:idiv(max_mag(M),54);
if ( !COUNT )
COUNT = 1;
TINV = TC = TR = TS = TM = TDIV = 0;
T0 = time()[0];
L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
TS += time()[0] - T0;
COL1 = COL - 1;
AA = newmat(COL1,COL1); BB = newvect(COL1);
for ( I = 0; I < COL1; I++ ) {
for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL1; J++ )
T[J] = S[J];
BB[I] = B[INDEX[I]];
}
if ( COL1 != ROW ) {
RESTA = newmat(ROW-COL1,COL1); RESTB = newvect(ROW-COL1);
for ( ; I < ROW; I++ ) {
for ( J = 0, T = RESTA[I-COL1], S = A[INDEX[I]]; J < COL1; J++ )
T[J] = S[J];
RESTB[I-COL1] = B[INDEX[I]];
}
} else
RESTA = RESTB = 0;
MOD2 = idiv(MOD,2);
for ( I = 0, C = BB, X = 0, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
I++, PK *= MOD ) {
if ( COUNT == CCC ) {
CCC = 0;
T0 = time()[0];
ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
TR += time()[0]-T0;
if ( ND ) {
T0 = time()[0];
F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
TM += time()[0]-T0;
if ( zerovector(T) ) {
T0 = time()[0]; T = RESTA*F+LCM*RESTB; TM += time()[0]-T0;
if ( zerovector(T) ) {
#if 0
if ( dp_gr_print() ) print(["init",TS,"pinv",TINV,"c",TC,"div",TDIV,"rat",TR,"mul",TM]);
#endif
if ( dp_gr_print() ) print("end",2);
return [F,LCM];
} else
return 0;
}
} else {
#if 0
if ( dp_gr_print() ) print(I);
#endif
}
} else {
#if 0
if ( dp_gr_print() ) print([I,TINV,TC,TDIV]);
#endif
if ( dp_gr_print() ) print(".",2);
CCC++;
}
T0 = time()[0];
XT = sremainder(INV*sremainder(-C,MOD),MOD);
XT = map(adj_sgn,XT,MOD,MOD2);
TINV += time()[0] - T0;
X += XT*PK;
T0 = time()[0];
C += mul_mat_vect_int(AA,XT);
TC += time()[0] - T0;
T0 = time()[0]; C = map(idiv,C,MOD); TDIV += time()[0] - T0;
}
}
def henleq_prep(A,MOD)
{
SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
AA = newmat(COL,COL);
for ( I = 0; I < COL; I++ )
for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL; J++ )
T[J] = S[J];
if ( COL != ROW ) {
RESTA = newmat(ROW-COL,COL);
for ( ; I < ROW; I++ )
for ( J = 0, T = RESTA[I-COL], S = A[INDEX[I]]; J < COL; J++ )
T[J] = S[J];
} else
RESTA = 0;
return [[A,AA,RESTA],L];
}
def henleq_gsl(L,B,MOD)
{
AL = L[0]; INVL = L[1];
A = AL[0]; AA = AL[1]; RESTA = AL[2];
INV = INVL[0]; INDEX = INVL[1];
SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
BB = newvect(COL);
for ( I = 0; I < COL; I++ )
BB[I] = B[INDEX[I]];
if ( COL != ROW ) {
RESTB = newvect(ROW-COL);
for ( ; I < ROW; I++ )
RESTB[I-COL] = B[INDEX[I]];
} else
RESTB = 0;
COUNT = INIT_COUNT?INIT_COUNT:idiv(MAX(max_mag(A),max_mag_vect(B)),54);
if ( !COUNT )
COUNT = 1;
MOD2 = idiv(MOD,2);
X = newvect(size(AA)[0]);
for ( I = 0, C = BB, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
I++, PK *= MOD ) {
if ( zerovector(C) )
if ( zerovector(RESTA*X+RESTB) ) {
if ( dp_gr_print() ) print("end",0);
return [X,1];
} else
return 0;
else if ( COUNT == CCC ) {
CCC = 0;
ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
if ( ND ) {
F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
if ( zerovector(T) ) {
T = RESTA*F+LCM*RESTB;
if ( zerovector(T) ) {
if ( dp_gr_print() ) print("end",0);
return [F,LCM];
} else
return 0;
}
} else {
}
} else {
if ( dp_gr_print() ) print(".",2);
CCC++;
}
XT = sremainder(INV*sremainder(-C,MOD),MOD);
XT = map(adj_sgn,XT,MOD,MOD2);
X += XT*PK;
C += mul_mat_vect_int(AA,XT);
C = map(idiv,C,MOD);
}
}
def adj_sgn(A,M,M2)
{
return A > M2 ? A-M : A;
}
def zerovector(C)
{
if ( !C )
return 1;
for ( I = size(C)[0]-1; I >= 0 && !C[I]; I-- );
if ( I < 0 )
return 1;
else
return 0;
}
def solvem(INV,COMP,V,MOD)
{
T = COMP*V;
N = size(T)[0];
for ( I = 0; I < N; I++ )
if ( T[I] % MOD )
return 0;
return modvect(INV*V,MOD);
}
def modmat(A,MOD)
{
if ( !A )
return 0;
S = size(A); N = S[0]; M = S[1];
MAT = newmat(N,M);
for ( I = 0, NZ = 0; I < N; I++ )
for ( J = 0, T1 = A[I], T2 = MAT[I]; J < M; J++ ) {
T2[J] = T1[J] % MOD;
NZ = NZ || T2[J];
}
return NZ?MAT:0;
}
def modvect(A,MOD)
{
if ( !A )
return 0;
N = size(A)[0];
VECT = newvect(N);
for ( I = 0, NZ = 0; I < N; I++ ) {
VECT[I] = A[I] % MOD;
NZ = NZ || VECT[I];
}
return NZ?VECT:0;
}
def qrmat(A,MOD)
{
if ( !A )
return [0,0];
S = size(A); N = S[0]; M = S[1];
Q = newmat(N,M); R = newmat(N,M);
for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ )
for ( J = 0, TA = A[I], TQ = Q[I], TR = R[I]; J < M; J++ ) {
L = iqr(TA[J],MOD); TQ[J] = L[0]; TR[J] = L[1];
NZQ = NZQ || TQ[J]; NZR = NZR || TR[J];
}
return [NZQ?Q:0,NZR?R:0];
}
def qrvect(A,MOD)
{
if ( !A )
return [0,0];
N = size(A)[0];
Q = newvect(N); R = newvect(N);
for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ ) {
L = iqr(A[I],MOD); Q[I] = L[0]; R[I] = L[1];
NZQ = NZQ || Q[I]; NZR = NZR || R[I];
}
return [NZQ?Q:0,NZR?R:0];
}
def max_mag(M)
{
R = size(M)[0];
U = 1;
for ( I = 0; I < R; I++ ) {
A = max_mag_vect(M[I]);
U = MAX(A,U);
}
return U;
}
def max_mag_vect(V)
{
R = size(V)[0];
U = 1;
for ( I = 0; I < R; I++ ) {
A = dp_mag(V[I]*<<0>>);
U = MAX(A,U);
}
return U;
}
def gsl_check(B,V,S)
{
N = length(V);
U = S[N-1]; M = U[1]; D = U[2];
W = setminus(V,[var(M)]);
H = uc(); VH = append(W,[H]);
for ( T = B; T != []; T = cdr(T) ) {
A = car(T);
AH = dp_dtop(dp_homo(dp_ptod(A,W)),VH);
for ( I = 0, Z = S; I < N-1; I++, Z = cdr(Z) ) {
L = car(Z); AH = ptozp(subst(AH,L[0],L[1]/L[2]));
}
AH = ptozp(subst(AH,H,D));
R = srem(AH,M);
if ( dp_gr_print() )
if ( !R )
print([A,"ok"]);
else
print([A,"bad"]);
if ( R )
break;
}
return T == [] ? 1 : 0;
}
def vs_dim(G,V,O)
{
HM = hmlist(G,V,O); ZD = zero_dim(HM,V,O);
if ( ZD ) {
MB = dp_mbase(map(dp_ptod,HM,V));
return length(MB);
} else
error("vs_dim : ideal is not zero-dimensional!");
}
def dgr(G,V,O)
{
P = getopt(proc);
if ( type(P) == -1 )
return gr(G,V,O);
P0 = P[0]; P1 = P[1]; P = [P0,P1];
map(ox_reset,P);
ox_cmo_rpc(P0,"dp_gr_main",G,V,0,1,O);
ox_cmo_rpc(P1,"dp_gr_main",G,V,1,1,O);
map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
F = ox_select(P);
R = ox_get(F[0]);
if ( F[0] == P0 ) {
Win = "nonhomo";
Lose = P1;
} else {
Win = "nhomo";
Lose = P0;
}
ox_reset(Lose);
return [Win,R];
}
/* functions for rpc */
def register_matrix(M)
{
REMOTE_MATRIX = M; return 0;
}
def register_nfv(L)
{
REMOTE_NF = L[0]; REMOTE_VARS = L[1]; return 0;
}
def r_ttob(T,M)
{
return hen_ttob(T,REMOTE_NF,0,REMOTE_VARS,M);
}
def r_ttob_gsl(L,M)
{
return cons(L[2],hen_ttob(L[0],REMOTE_NF,L[1],REMOTE_VARS,M));
}
def get_matrix()
{
REMOTE_MATRIX;
}
extern NFArray$
/*
* HL = [[c,i,m,d],...]
* if c != 0
* g = 0
* g = (c*g + m*gi)/d
* ...
* finally compare g with NF
* if g == NF then NFArray[NFIndex] = g
*
* if c = 0 then HL consists of single history [0,i,0,0],
* which means that dehomogenization of NFArray[i] should be
* eqall to NF.
*/
def check_trace(NF,NFIndex,HL)
{
if ( !car(HL)[0] ) {
/* dehomogenization */
DH = dp_dehomo(NFArray[car(HL)[1]]);
if ( NF == DH ) {
realloc_NFArray(NFIndex);
NFArray[NFIndex] = NF;
return 0;
} else
error("check_trace(dehomo)");
}
for ( G = 0, T = HL; T != []; T = cdr(T) ) {
H = car(T);
Coeff = H[0];
Index = H[1];
Monomial = H[2];
Denominator = H[3];
Reducer = NFArray[Index];
G = (Coeff*G+Monomial*Reducer)/Denominator;
}
if ( NF == G ) {
realloc_NFArray(NFIndex);
NFArray[NFIndex] = NF;
return 0;
} else
error("check_trace");
}
/*
* realloc NFArray so that it can hold * an element as NFArray[Ind].
*/
def realloc_NFArray(Ind)
{
if ( Ind == size(NFArray)[0] ) {
New = newvect(Ind + 100);
for ( I = 0; I < Ind; I++ )
New[I] = NFArray[I];
NFArray = New;
}
}
/*
* create NFArray and initialize it by List.
*/
def register_input(List)
{
Len = length(List);
NFArray = newvect(Len+100,List);
}
/*
tracetogen(): preliminary version
dp_gr_main() returns [GB,GBIndex,Trace].
GB : groebner basis
GBIndex : IndexList (corresponding to Trace)
Trace : [InputList,Trace0,Trace1,...]
TraceI : [Index,TraceList]
TraceList : [[Coef,Index,Monomial,Denominator],...]
Poly <- 0
Poly <- (Coef*Poly+Monomial*PolyList[Index])/Denominator
*/
def tracetogen(G)
{
GB = G[0]; GBIndex = G[1]; Trace = G[2];
InputList = Trace[0];
Trace = cdr(Trace);
/* number of initial basis */
Nini = length(InputList);
/* number of generated basis */
Ngen = length(Trace);
N = Nini + Ngen;
/* stores traces */
Tr = vector(N);
/* stores coeffs */
Coef = vector(N);
/* XXX create dp_ptod(1,V) */
HT = dp_ht(InputList[0]);
One = dp_subd(HT,HT);
for ( I = 0; I < Nini; I++ ) {
Tr[I] = [1,I,One,1];
C = vector(Nini);
C[I] = One;
Coef[I] = C;
}
for ( ; I < N; I++ )
Tr[I] = Trace[I-Nini][1];
for ( T = GBIndex; T != []; T = cdr(T) )
compute_coef_by_trace(car(T),Tr,Coef);
return Coef;
}
def compute_coef_by_trace(I,Tr,Coef)
{
if ( Coef[I] )
return;
/* XXX */
Nini = size(Coef[0])[0];
/* initialize coef vector */
CI = vector(Nini);
for ( T = Tr[I]; T != []; T = cdr(T) ) {
/* Trace = [Coef,Index,Monomial,Denominator] */
Trace = car(T);
C = Trace[0];
Ind = Trace[1];
Mon = Trace[2];
Den = Trace[3];
if ( !Coef[Ind] )
compute_coef_by_trace(Ind,Tr,Coef);
/* XXX */
CT = newvect(Nini);
for ( J = 0; J < Nini; J++ )
CT[J] = (C*CI[J]+Mon*Coef[Ind][J])/Den;
CI = CT;
}
Coef[I] = CI;
}
end$