=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/sp,v retrieving revision 1.1 retrieving revision 1.18 diff -u -p -r1.1 -r1.18 --- OpenXM_contrib2/asir2000/lib/sp 1999/12/03 07:39:11 1.1 +++ OpenXM_contrib2/asir2000/lib/sp 2018/04/06 07:40:44 1.18 @@ -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.17 2016/02/04 04:17:21 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 SpOrd$ +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,11 @@ def flatmf(L) { def af(P,AL) { RESTIME=UFTIME=GCDTIME=SQTIME=0; - S = reverse(asq(P,AL)); + V = var(P); + LC = coef(P,deg(P,V),V); + if ( ntype(LC) != 1 ) + P = simpalg(1/LC*P); + 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); @@ -100,6 +231,10 @@ def af_sp(P,AL,HINT) { if ( !P || type(P) == NUM ) return [P]; + V = var(P); + LC = coef(P,deg(P,V),V); + if ( ntype(LC) != 1 ) + P = simpalg(1/LC*P); P1 = simpcoef(simpalg(P)); return af_spmain(P1,AL,1,HINT,P1,[]); } @@ -124,7 +259,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 +277,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 +405,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 +508,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 +520,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 +559,7 @@ def sort_alg(VL) return L; } -def asq(P,AL) +def asq(P) { P = simpalg(P); if ( type(P) == NUM ) @@ -423,7 +572,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 +616,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 +643,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 +653,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); } @@ -503,6 +667,15 @@ def ufctrhint2(P,HINT,PP,AL) return [[P,1]]; if ( AL == [] ) return ufctrhint(P,HINT); + + /* if P != norm(PP) then call the generic ufctrhint() */ + for ( T = AL, E = 1; T != []; T = cdr(T) ) { + D = defpoly(car(T)); E *= deg(D,var(D)); + } + if ( E*deg(PP,var(PP)) != deg(P,var(P)) ) + return ufctrhint(P,HINT); + + /* P = norm(PP) */ L = resfctr(algptorat(PP),map(defpoly,AL),map(algtorat,AL),P); for ( T = reverse(L[1]), DL = []; T != []; T = cdr(T) ) DL = cons(deg(car(car(T)),a_),DL); @@ -645,27 +818,18 @@ 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) ) + if ( !P1 ) + return P2; + if ( !P2 ) + return P1; + if ( !var(P1) || !var(P2) ) return 1; - AL = union(getalgp(P1),getalgp(P2)); - if ( AL == [] ) + 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; @@ -690,7 +854,7 @@ def cr_gcda(P1,P2,EXT) break; if ( J != length(DL) ) continue; - Ord = 2; NOSUGAR = 1; + SpOrd = 2; NOSUGAR = 1; T = ag_mod(G1 % MOD,G2 % MOD,ML,VL,MOD); if ( dp_gr_print() ) print("."); @@ -747,18 +911,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)); @@ -916,7 +1068,7 @@ def ag_mod(F1,F2,D,VL,MOD) VL = cons(V,VL); B = append([F1,F2],D); N = length(VL); while ( 1 ) { FLAGS = dp_gr_flags(); dp_gr_flags(["Reverse",1,"NoSugar",1]); - G = dp_gr_mod_main(B,VL,0,MOD,Ord); + G = dp_gr_mod_main(B,VL,0,MOD,SpOrd); dp_gr_flags(FLAGS); if ( length(G) == 1 ) return 1; @@ -1150,9 +1302,9 @@ def ag_mod_single6(F1,F2,D,MOD) def inverse_by_gr_mod(C,D,MOD) { - Ord = 2; + SpOrd = 2; dp_gr_flags(["NoSugar",1]); - G = dp_gr_mod_main(cons(x*C-1,D),cons(x,vars(D)),0,MOD,Ord); + G = dp_gr_mod_main(cons(x*C-1,D),cons(x,vars(D)),0,MOD,SpOrd); dp_gr_flags(["NoSugar",0]); if ( length(G) == 1 ) return 1; @@ -1211,18 +1363,27 @@ def resfctr(F,L,V,N) N = ptozp(N); V0 = var(N); DN = diff(N,V0); - for ( I = 0, J = 2, Len = deg(N,V0)+1; I < 5; J++ ) { - M = prime(J); - G = gcd(N,DN,M); - if ( !deg(G,V0) ) { - I++; - T = nfctr_mod(N,M); - if ( T < Len ) { - Len = T; M0 = M; - } - } - } - S = spm(L,V,M0); + LC = coef(N,deg(N,V0),V0); + LCD = coef(DN,deg(DN,V0),V0); + J = 2; + while ( 1 ) { + Len = deg(N,V0)+1; + for ( I = 0; I < 5; J++ ) { + M = prime(J); + if ( !(LC%M) || !(LCD%M)) + continue; + G = gcd(N,DN,M); + if ( !deg(G,V0) ) { + I++; + T = nfctr_mod(N,M); + if ( T < Len ) { + Len = T; M0 = M; + } + } + } + S = spm(L,V,M0); + if ( S ) break; + } T = resfctr_mod(F,S,M0); return [T,S,M0]; } @@ -1236,7 +1397,7 @@ def resfctr_mod(F,L,M) C = res(var(MP),B,MP) % M; R = cons(flatten(cdr(modfctr(C,M))),R); } - return R; + return reverse(R); } def flatten(L) @@ -1252,11 +1413,13 @@ def spm(L,V,M) U = modfctr(car(L),M); for ( T = cdr(U), R = []; T != []; T = cdr(T) ) { S = car(T); + if ( S[1] > 1 ) return 0; R = cons([subst(S[0],var(S[0]),a_),[var(S[0]),a_]],R); } return R; } L1 = spm(cdr(L),cdr(V),M); + if ( !L1 ) return 0; F0 = car(L); V0 = car(V); VR = cdr(V); for ( T = L1, R = []; T != []; T = cdr(T) ) { S = car(T); @@ -1264,6 +1427,7 @@ def spm(L,V,M) U = fctr_mod(F1,V0,S[0],M); VS = var(S[0]); for ( W = U; W != []; W = cdr(W) ) { + if ( car(W)[1] > 1 ) return 0; A = car(car(W)); if ( deg(A,V0) == 1 ) { A = monic_mod(A,V0,S[0],M);