=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/sp,v retrieving revision 1.1 retrieving revision 1.12 diff -u -p -r1.1 -r1.12 --- OpenXM_contrib2/asir2000/lib/sp 1999/12/03 07:39:11 1.1 +++ OpenXM_contrib2/asir2000/lib/sp 2004/01/07 08:23:11 1.12 @@ -1,17 +1,82 @@ -/* $OpenXM: OpenXM/src/asir99/lib/sp,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */ +/* + * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED + * All rights reserved. + * + * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, + * non-exclusive and royalty-free license to use, copy, modify and + * redistribute, solely for non-commercial and non-profit purposes, the + * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and + * conditions of this Agreement. For the avoidance of doubt, you acquire + * only a limited right to use the SOFTWARE hereunder, and FLL or any + * third party developer retains all rights, including but not limited to + * copyrights, in and to the SOFTWARE. + * + * (1) FLL does not grant you a license in any way for commercial + * purposes. You may use the SOFTWARE only for non-commercial and + * non-profit purposes only, such as academic, research and internal + * business use. + * (2) The SOFTWARE is protected by the Copyright Law of Japan and + * international copyright treaties. If you make copies of the SOFTWARE, + * with or without modification, as permitted hereunder, you shall affix + * to all such copies of the SOFTWARE the above copyright notice. + * (3) An explicit reference to this SOFTWARE and its copyright owner + * shall be made on your publication or presentation in any form of the + * results obtained by use of the SOFTWARE. + * (4) In the event that you modify the SOFTWARE, you shall notify FLL by + * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification + * for such modification or the source code of the modified part of the + * SOFTWARE. + * + * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL + * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND + * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' + * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY + * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. + * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, + * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY + * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL + * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES + * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES + * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY + * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF + * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART + * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY + * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, + * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. + * + * $OpenXM: OpenXM_contrib2/asir2000/lib/sp,v 1.11 2003/02/07 09:25:58 noro Exp $ +*/ /* sp : functions related to algebraic number fields Revision History: - 99/08/24 noro modified for 1999 release version + 2001/10/12 noro if USE_PARI_FACTOR is nonzero, pari factor is called + 2000/03/10 noro fixed several bugs around gathering algebraic numbers + 1999/08/24 noro modified for 1999 release version */ #include "defs.h" extern ASCENT,GCDTIME,UFTIME,RESTIME,SQTIME,PRINT$ extern Ord$ +extern USE_PARI_FACTOR$ +/* gen_sp can handle non-monic poly */ + +def gen_sp(P) +{ + P = ptozp(P); + V = var(P); + D = deg(P,V); + LC = coef(P,D,V); + F = LC^(D-1)*subst(P,V,V/LC); + /* F must be monic */ + L = sp(F); + return cons(map(subst,car(L),V,LC*V),cdr(L)); +} + def sp(P) { RESTIME=UFTIME=GCDTIME=SQTIME=0; @@ -39,6 +104,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 +211,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 +251,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 +269,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 +397,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 +500,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 +512,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 +551,7 @@ def sort_alg(VL) return L; } -def asq(P,AL) +def asq(P) { P = simpalg(P); if ( type(P) == NUM ) @@ -423,7 +564,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; @@ -467,6 +608,21 @@ def simpcoef(P) { } def ufctrhint_heuristic(P,HINT,PP,SHIFT) { + if ( USE_PARI_FACTOR ) + return pari_ufctr(P); + else + return asir_ufctrhint_heuristic(P,HINT,PP,SHIFT); +} + +def pari_ufctr(P) { + F = pari(factor,P); + S = size(F); + for ( I = S[0]-1, R = []; I >= 0; I-- ) + R = cons(vtol(F[I]),R); + return cons([1,1],R); +} + +def asir_ufctrhint_heuristic(P,HINT,PP,SHIFT) { V = var(P); D = deg(P,V); if ( D == HINT ) return [[P,1]]; @@ -479,7 +635,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 +645,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 +801,16 @@ 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 == [] ) + if ( !var(P1) ) + return P2; + else if ( !var(P2) ) + return P1; + V = var(P1); + 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; @@ -747,18 +892,6 @@ def member(E,L) 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) { V = var(P); return res(V,P,diff(P,V)); @@ -1211,8 +1344,12 @@ def resfctr(F,L,V,N) N = ptozp(N); V0 = var(N); DN = diff(N,V0); + LC = coef(N,deg(N,V0),V0); + LCD = coef(DN,deg(DN,V0),V0); for ( I = 0, J = 2, Len = deg(N,V0)+1; I < 5; J++ ) { M = prime(J); + if ( !(LC%M) || !(LCD%M)) + continue; G = gcd(N,DN,M); if ( !deg(G,V0) ) { I++;