=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/sp,v retrieving revision 1.5 retrieving revision 1.18 diff -u -p -r1.5 -r1.18 --- OpenXM_contrib2/asir2000/lib/sp 2000/04/20 02:20:16 1.5 +++ OpenXM_contrib2/asir2000/lib/sp 2018/04/06 07:40:44 1.18 @@ -1,18 +1,82 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/lib/sp,v 1.1.1.1 1999/12/03 07:39:11 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: - 2000/04/20 noro fixed bugs around gathering algebraic numbers + 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; @@ -40,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) ) { @@ -85,6 +211,10 @@ def flatmf(L) { def af(P,AL) { RESTIME=UFTIME=GCDTIME=SQTIME=0; + 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]; @@ -101,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,[]); } @@ -482,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]]; @@ -518,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); @@ -662,8 +820,13 @@ def norm_ch_lag(V,VM,P,P0) { 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; + V = var(P1); EXT = union_sort(getalgtreep(P1),getalgtreep(P2)); if ( EXT == [] ) return gcd(P1,P2); @@ -691,7 +854,7 @@ def cr_gcda(P1,P2) 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("."); @@ -905,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; @@ -1139,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; @@ -1200,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]; } @@ -1225,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) @@ -1241,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); @@ -1253,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);