=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/sp,v retrieving revision 1.1.1.1 retrieving revision 1.6 diff -u -p -r1.1.1.1 -r1.6 --- OpenXM_contrib2/asir2000/lib/sp 1999/12/03 07:39:11 1.1.1.1 +++ OpenXM_contrib2/asir2000/lib/sp 2000/04/20 02:30:35 1.6 @@ -1,10 +1,11 @@ -/* $OpenXM: OpenXM/src/asir99/lib/sp,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */ +/* $OpenXM$ */ /* sp : functions related to algebraic number fields Revision History: - 99/08/24 noro modified for 1999 release version + 2000/03/10 noro fixed several bugs around gathering algebraic numbers + 1999/08/24 noro modified for 1999 release version */ #include "defs.h" @@ -39,6 +40,68 @@ def sp(P) } } +/* + Input: + F=F(x,a1,...,an) + DL = [[an,dn(an,...,a1)],...,[a2,d2(a2,a1)],[a1,d1(a1)]] + 'ai' denotes a root of di(t). + Output: + irreducible factorization of F over Q(a1,...,an) + [[F1(x,a1,...,an),e1],...,[Fk(x,a1,...,an),ek]] + 'ej' denotes the multiplicity of Fj. +*/ + +def af_noalg(F,DL) +{ + DL = reverse(DL); + N = length(DL); + Tab = newvect(N); + /* Tab = [[a1,r1],...]; ri is a root of di(t,r(i-1),...,r1). */ + AL = []; + for ( I = 0; I < N; I++ ) { + T = DL[I]; + for ( J = 0, DP = T[1]; J < I; J++ ) + DP = subst(DP,Tab[J][0],Tab[J][1]); + B = newalg(DP); + Tab[I] = [T[0],B]; + F = subst(F,T[0],B); + AL = cons(B,AL); + } + FL = af(F,AL); + for ( T = FL, R = []; T != []; T = cdr(T) ) + R = cons([conv_noalg(T[0][0],Tab),T[0][1]],R); + return reverse(R); +} + +/* + Input: + F=F(x) univariate polynomial over the rationals + Output: + [FL,DL] + DL = [[an,dn(an,...,a1)],...,[a2,d2(a2,a1)],[a1,d1(a1)]] + 'ai' denotes a root of di(t). + FL = [F1,F2,...] + irreducible factors of F over Q(a1,...,an) +*/ + +def sp_noalg(F) +{ + L = sp(F); + FL = map(algptorat,L[0]); + for ( T = L[1], DL = []; T != []; T = cdr(T) ) + DL = cons([algtorat(T[0][0]),T[0][1]],DL); + return [FL,reverse(DL)]; +} + +def conv_noalg(F,Tab) +{ + N = size(Tab)[0]; + F = algptorat(F); + for ( I = N-1; I >= 0; I-- ) + F = subst(F,algtorat(Tab[I][1]),Tab[I][0]); + return F; +} + def aflist(L,AL) { for ( DC = []; L != []; L = cdr(L) ) { @@ -84,7 +147,7 @@ def flatmf(L) { def af(P,AL) { RESTIME=UFTIME=GCDTIME=SQTIME=0; - S = reverse(asq(P,AL)); + S = reverse(asq(P)); for ( L = []; S != []; S = cdr(S) ) { FM = car(S); F = FM[0]; M = FM[1]; G = af_sp(F,AL,1); @@ -124,7 +187,7 @@ def af_spmain(P,AL,INIT,HINT,PP,SHIFT) N = simpcoef(sp_norm(A0,V,subst(P,V,V-INIT*A0),AL)); RESTIME+=time()[0]-TTT; TTT = time()[0]; - DCSQ = sortfs(asq(N,AL)); + DCSQ = sortfs(asq(N)); SQTIME+=time()[0]-TTT; for ( G = P, A = V+INIT*A0, DCR = []; DCSQ != []; DCSQ = cdr(DCSQ) ) { C = TT(DCSQ); D = TS(DCSQ); @@ -142,7 +205,7 @@ def af_spmain(P,AL,INIT,HINT,PP,SHIFT) else { S = subst(car(DCT),V,A); if ( pra(G,S,AL) ) - U = cr_gcda(S,G,AL); + U = cr_gcda(S,G); else U = S; } @@ -270,26 +333,40 @@ def getalgp(P) { else { for ( V = var(P), I = deg(P,V), T = []; I >= 0; I-- ) if ( C = coef(P,I) ) - T = union(T,getalgp(C)); + T = union_sort(T,getalgp(C)); return T; } } -def union(A,B) -{ - for ( T = B; T != []; T = cdr(T) ) - A = union1(A,car(T)); - return A; +def getalgtreep(P) { + if ( type(P) <= 1 ) + return getalgtree(P); + else { + for ( V = var(P), I = deg(P,V), T = []; I >= 0; I-- ) + if ( C = coef(P,I) ) + T = union_sort(T,getalgtreep(C)); + return T; + } } -def union1(A,E) +/* C = union of A and B; A and B is sorted. C should also be sorted. */ + +def union_sort(A,B) { if ( A == [] ) - return [E]; - else if ( car(A) == E ) + return B; + else if ( B == [] ) return A; - else - return cons(car(A),union1(cdr(A),E)); + else { + A0 = car(A); + B0 = car(B); + if ( A0 == B0 ) + return cons(A0,union_sort(cdr(A),cdr(B))); + else if ( A0 > B0 ) + return cons(A0,union_sort(cdr(A),B)); + else + return cons(B0,union_sort(A,cdr(B))); + } } def invalgp(A) @@ -359,7 +436,7 @@ def sortfs(L) def pdiva(P1,P2) { - A = union(getalgp(P1),getalgp(P2)); + A = union_sort(getalgp(P1),getalgp(P2)); P1 = algptorat(P1); P2 = algptorat(P2); return simpalg(rattoalgp(sdiv(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2),A)); } @@ -371,7 +448,7 @@ def pqra(P1,P2) else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) ) return [0,P1]; else { - A = union(getalgp(P1),getalgp(P2)); + A = union_sort(getalgp(P1),getalgp(P2)); P1 = algptorat(P1); P2 = algptorat(P2); L = sqr(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2); return [simpalg(rattoalgp(L[0],A)),simpalg(rattoalgp(L[1],A))]; @@ -410,7 +487,7 @@ def sort_alg(VL) return L; } -def asq(P,AL) +def asq(P) { P = simpalg(P); if ( type(P) == NUM ) @@ -423,7 +500,7 @@ def asq(P,AL) if ( type(F) == NUM ) break; F1 = diff(F,V); - GCD = cr_gcda(F,F1,AL); + GCD = cr_gcda(F,F1); FLAT = pdiva(F,GCD); if ( type(GCD) == NUM ) { A[I] = F; B[I] = 1; @@ -479,7 +556,7 @@ def ufctrhint_heuristic(P,HINT,PP,SHIFT) { G = igcd(G,DT=deg(T,V)); if ( G == 1 ) { if ( K*deg(PPP,V) != deg(P,V) ) - PPP = cr_gcda(PPP,P,AL); + PPP = cr_gcda(PPP,P); return ufctrhint2(P,HINT,PPP,AL); } else { for ( S = 0, T = P; T; T -= coef(T,DT)*V^DT ) { @@ -489,7 +566,7 @@ def ufctrhint_heuristic(P,HINT,PP,SHIFT) { L = fctr(S); for ( DC = [car(L)], L = cdr(L); L != []; L = cdr(L) ) { H = subst(car(car(L)),V,V^G); - HH = cr_gcda(PPP,H,AL); + HH = cr_gcda(PPP,H); T = ufctrhint2(H,HINT,HH,AL); DC = append(DC,T); } @@ -645,27 +722,13 @@ def norm_ch_lag(V,VM,P,P0) { return S; } -def cr_gcda(P1,P2,EXT) +def cr_gcda(P1,P2) { if ( !(V = var(P1)) || !var(P2) ) return 1; - AL = union(getalgp(P1),getalgp(P2)); - if ( AL == [] ) + EXT = union_sort(getalgtreep(P1),getalgtreep(P2)); + if ( EXT == [] ) return gcd(P1,P2); - T = newvect(length(EXT)); - for ( TAL = AL; TAL != []; TAL = cdr(TAL) ) { - A = getalg(car(TAL)); - for ( TA = A; TA != []; TA = cdr(TA) ) { - B = car(TA); - for ( TEXT = EXT, I = 0; TEXT != []; TEXT = cdr(TEXT), I++ ) - if ( car(TEXT) == B ) - T[I] = B; - } - } - for ( I = length(EXT)-1, S = []; I >= 0; I-- ) - if ( T[I] ) - S = cons(T[I],S); - EXT = S; NEXT = length(EXT); if ( deg(P1,V) < deg(P2,V) ) { T = P1; P1 = P2; P2 = T; @@ -745,18 +808,6 @@ def member(E,L) if ( E == car(L) ) return 1; return 0; -} - -def getallalg(A) -{ - T = cdr(vars(defpoly(A))); - if ( T == [] ) - return [A]; - else { - for ( S = [A]; T != []; T = cdr(T) ) - S = union(S,getallalg(rattoalg(car(T)))); - return S; - } } def discr(P) {