![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to de CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib|
Return to de CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.1 noro 1: module de;
2:
3: localf split_lexgb;
4: localf sort_lex_dec,sort_lex_inc;
5: localf inverse_or_split, linear_dim;
1.4 ! noro 6: localf dp_monic_mod,monic_gb;
1.3 noro 7: localf membership_test;
1.1 noro 8: localf dp_chrem,intdptoratdp,intdpltoratdpl;
9: localf comp_by_ht,dp_gr_mod,gr_chrem;
1.4 ! noro 10: localf construct_sqfrbasis;
1.1 noro 11:
12: /*
13: * G : a 0-dim lex gb, reduced
14: * if there exists a non-monic g = g(x')x^k+...
15: * then g(x') is a zero divisor and Id(G) splits
16: */
17:
18: def split_lexgb(G,V)
19: {
20: G = map(ptozp,G);
21: G = sort_lex_dec(G,V);
22: TG = G;
23: for ( ; TG != []; TG = cdr(TG) ) {
24: F = car(TG);
25: for ( TV = V; !deg(F,car(TV)); TV = cdr(TV) );
26: VF = car(TV);
27: DF = deg(F,VF);
28: CF = coef(F,DF,VF);
29: if ( type(CF) == 2 ) {
30: L = inverse_or_split(V,G,CF);
31: R = map(split_lexgb,L,V);
32: return append(R[0],R[1]);
33: }
34: }
35: return [G];
36: }
37:
38: /*
39: * V : a variable list
40: * Id : a 0-dim radical ideal represented by a lex basis
41: * F : a poly
42: * if F is a unit of Q[V]/Id, then 1/F is returned
43: * else F is a zero divisor and Id = (Id:F)\cap(Id+<F>)
44: * [GB(Id:F),GB(Id+<F>)] is returned.
45: */
46:
47: def inverse_or_split(V,Id,F)
48: {
49: Id = map(ptozp,Id);
50: N = length(V);
51: dp_ord(2);
52: set_field(Id,V,2);
53: DF = dptodalg(dp_ptod(F,V));
54: Ret = inv_or_split_dalg(DF);
55: if ( type(Ret) == 1 ) {
56: /* Ret = 1/F */
1.4 ! noro 57: Dp = dalgtodp(Ret);
! 58: return dp_dtop(Dp[0],V)/Dp[1];
1.1 noro 59: } else {
60: /* Ret = GB(Id:F) */
61: /* compute GB(Id+<f>) */
62: Gquo = append(map(ptozp,map(dp_dtop,Ret,V)),Id);
1.3 noro 63: /* inter-reduction */
64: Gquo = nd_gr_postproc(Gquo,V,0,2,0);
1.1 noro 65: DTotal = linear_dim(Id,V,2);
66: Dquo = linear_dim(Gquo,V,2);
67: Drem = DTotal-Dquo;
68: B = cons(F,Id);
69: Grem = gr_chrem(B,V,2,Drem);
70: return [map(ptozp,Gquo),map(ptozp,Grem)];
71: }
72: }
73:
1.4 ! noro 74: /* add F(X,V) to Id(B) */
! 75: /* returns a list of splitted ideals */
! 76: /* B should be a triangular basis */
! 77:
! 78: def construct_sqfrbasis(F,X,B,V)
! 79: {
! 80: B = sort_lex_dec(B,V);
! 81: V1 = cons(X,V);
! 82: F = nd_nf(F,reverse(B),cons(X,V),2,0);
! 83: D = deg(F,X);
! 84: H = coef(F,D,X);
! 85: if ( type(H) == 1 )
! 86: return [];
! 87: else if ( type(H) == 2 ) {
! 88: Ret = inverse_or_split(V,B,H);
! 89: if ( type(Ret) == 4 ) {
! 90: /* H != 0 on Id_nz, H = 0 on Id_z */
! 91: B0=construct_sqfrbasis(F,X,Ret[0],V);
! 92: B1=construct_sqfrbasis(F,X,Ret[1],V);
! 93: return append(B0,B1);
! 94: } else
! 95: F = nd_nf(F*Ret,reverse(B),cons(X,V),2,0);
! 96: }
! 97: B1 = cons(F,B);
! 98: /* F is monic */
! 99: M = minipoly(B1,V1,2,X,zzz);
! 100: S = sqfr(M); S = cdr(S);
! 101: if ( length(S) == 1 && car(S)[1] == 1 )
! 102: return [cons(F,B)];
! 103: else {
! 104: R = [];
! 105: for ( T = S; T != []; T = cdr(T) ) {
! 106: G = nd_gr_trace(cons(car(T),B1),V1,1,1,2);
! 107: R1 = split_lexgb(G,V1);
! 108: R = append(R1,R);
! 109: }
! 110: return R;
! 111: }
! 112: }
! 113:
1.1 noro 114: def sort_lex_dec(B,V)
115: {
116: dp_ord(2);
117: B = map(dp_ptod,B,V);
118: B = vtol(qsort(ltov(B),comp_by_ht));
119: B = map(dp_dtop,B,V);
120: return reverse(B);
121: }
122:
123: def sort_lex_inc(B,V)
124: {
125: dp_ord(2);
126: B = map(dp_ptod,B,V);
127: B = vtol(qsort(ltov(B),comp_by_ht));
128: B = map(dp_dtop,B,V);
129: return B;
130: }
131:
132: def linear_dim(G,V,Ord)
133: {
134: dp_ord(Ord);
135: MB = dp_mbase(map(dp_ptod,G,V));
136: return length(MB);
137: }
138:
1.3 noro 139: def membership_test(B,G,V,O)
140: {
141: B = map(ptozp,B);
142: G = map(ptozp,G);
143: for ( T = B; T != []; T = cdr(T) )
144: if ( nd_nf(car(T),G,V,O,0) ) return 0;
145: return 1;
146: }
147:
1.1 noro 148: def gr_chrem(B,V,O,Dim)
149: {
150: B = map(ptozp,B);
151: G = []; HS = []; Mod = 1;
152: for ( I = 0; ; I++ ) {
153: P = lprime(I);
154: GM = nd_gr(B,V,P,O);
155: if ( linear_dim(GM,V,O) != Dim ) continue;
156: L = monic_gb(GM,V,O,P); GM = L[0]; HSM = L[1];
157: G = dp_chrem(G,HS,Mod,GM,HSM,P);
158: Mod *= P;
159: if ( G != [] )
160: HS = HSM;
1.4 ! noro 161: M = idiv(isqrt(2*Mod),2);
! 162: R1 = intdpltoratdpl(G,Mod,M);
1.1 noro 163: if ( R1 ) {
1.3 noro 164: if ( Found && R == R1
165: && (GB=nd_gr_postproc(map(dp_dtop,R,V),V,0,O,1))
166: && membership_test(B,GB,V,O) )
1.1 noro 167: break;
168: else {
169: R = R1; Found = 1;
170: }
171: }
172: }
1.3 noro 173: return GB;
1.1 noro 174: }
175:
176: def comp_by_ht(A,B)
177: {
178: HA = dp_ht(A); HB = dp_ht(B);
179: if ( HA > HB )
180: return 1;
181: else if ( HA < HB )
182: return -1;
183: else
184: return 0;
185: }
186:
1.4 ! noro 187: def intdpltoratdpl(G,Mod,M)
1.1 noro 188: {
1.4 ! noro 189: for ( R = []; G != []; G = cdr(G) ) {
1.1 noro 190: T = intdptoratdp(car(G),Mod,M);
191: if ( !T )
192: return 0;
193: else
194: R = cons(T,R);
195: }
196: R = reverse(R);
197: return vtol(qsort(newvect(length(R),R),comp_by_ht));
198: }
199:
200: def intdptoratdp(F,Mod,M)
201: {
202: for ( T = F, N = 0; T; T = dp_rest(T), N++ );
203: C = newvect(N);
204: for ( I = 0, T = F; I < N; T = dp_rest(T), I++ ) {
205: L = inttorat(dp_hc(T),Mod,M);
206: if ( !L )
207: return 0;
208: else
209: C[I] = (L[0]/L[1])*dp_ht(T);
210: }
211: for ( R = 0, I = N-1; I >= 0; I-- )
212: R += C[I];
213: return R;
214: }
215:
216: def dp_chrem(G,HS,Mod,GM,HSM,P)
217: {
218: if ( G == [] )
219: return GM;
220: if ( HS != HSM )
221: return [];
222: R = [];
223: M1 = inv(Mod,P);
1.4 ! noro 224: ModM1 = Mod*M1;
1.1 noro 225: for ( ; G != []; G = cdr(G), GM = cdr(GM) ) {
226: E = car(G); EM = car(GM);
1.4 ! noro 227: E1 = E+(EM-E)*ModM1;
1.1 noro 228: R = cons(E1,R);
229: }
230: return reverse(R);
231: }
232:
233: def monic_gb(G,V,O,P)
234: {
235: dp_ord(O); setmod(P);
236: D = map(dp_ptod,G,V);
237: D = map(dp_monic_mod,D,P);
238: D = vtol(qsort(newvect(length(D),D),comp_by_ht));
239: return [D,map(dp_ht,D)];
240: }
241:
242: def dp_monic_mod(F,P)
243: {
244: FP = dp_mod(F,P,[]);
245: return dp_rat(FP/dp_hc(FP));
246: }
247:
248: endmodule;
249: end$