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