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