Annotation of OpenXM_contrib2/asir2018/lib/de, Revision 1.1
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 dp_monic_mod,monic_gb;
! 7: localf membership_test;
! 8: localf dp_chrem,intdptoratdp,intdpltoratdpl;
! 9: localf comp_by_ht,dp_gr_mod,gr_chrem;
! 10: localf construct_sqfrbasis;
! 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 */
! 57: Dp = dalgtodp(Ret);
! 58: return dp_dtop(Dp[0],V)/Dp[1];
! 59: } else {
! 60: /* Ret = GB(Id:F) */
! 61: /* compute GB(Id+<f>) */
! 62: Gquo = append(map(ptozp,map(dp_dtop,Ret,V)),Id);
! 63: /* inter-reduction */
! 64: Gquo = nd_gr_postproc(Gquo,V,0,2,0);
! 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:
! 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: if ( type(F) == 1 )
! 81: return [];
! 82: B = sort_lex_dec(B,V);
! 83: V1 = cons(X,V);
! 84: F = nd_nf(F,reverse(B),cons(X,V),2,0);
! 85: D = deg(F,X);
! 86: H = coef(F,D,X);
! 87: 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(subst(car(T)[0],zzz,X),B1),V1,1,1,2);
! 107: R1 = split_lexgb(G,V1);
! 108: R = append(R1,R);
! 109: }
! 110: return R;
! 111: }
! 112: }
! 113:
! 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:
! 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:
! 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;
! 161: M = idiv(isqrt(2*Mod),2);
! 162: R1 = intdpltoratdpl(G,Mod,M);
! 163: if ( R1 ) {
! 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) )
! 167: break;
! 168: else {
! 169: R = R1; Found = 1;
! 170: }
! 171: }
! 172: }
! 173: return GB;
! 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:
! 187: def intdpltoratdpl(G,Mod,M)
! 188: {
! 189: for ( R = []; G != []; G = cdr(G) ) {
! 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);
! 224: ModM1 = Mod*M1;
! 225: for ( ; G != []; G = cdr(G), GM = cdr(GM) ) {
! 226: E = car(G); EM = car(GM);
! 227: E1 = E+(EM-E)*ModM1;
! 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$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>