![[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);
! 60: Gquo = nd_gr(Gquo,V,0,2);
! 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$