File: [local] / OpenXM_contrib2 / asir2000 / lib / sp (download)
Revision 1.15, Fri Jun 23 08:57:47 2006 UTC (18 years, 3 months ago) by noro
Branch: MAIN
CVS Tags: R_1_3_1-2, RELEASE_1_2_3_12, DEB_REL_1_2_3-9 Changes since 1.14: +10 -1
lines
Fixed a bug in ufctrhint2().
|
/*
* 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.15 2006/06/23 08:57:47 noro Exp $
*/
/*
sp : functions related to algebraic number fields
Revision History:
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;
L = flatmf(fctr(P)); X = var(P);
AL = []; ADL = [];
while ( 1 ) {
L = sort_by_deg(L);
for ( T = L, H = []; T != []; H = cons(car(T),H), T = cdr(T) )
if ( deg(car(T),X) > 1 )
break;
if ( T == [] ) {
if ( dp_gr_print() ) {
print(["GCDTIME = ",GCDTIME]);
print(["UFTIME = ",UFTIME]);
print(["RESTIME = ",RESTIME]);
}
return [L,ADL];
} else {
A = newalg(car(T));
R = pdiva(car(T),X-A);
AL = cons(A,AL);
ADL = cons([A,defpoly(A)],ADL);
L = aflist(append(H,append([X-A,R],cdr(T))),AL);
}
}
}
/*
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) ) {
T = af_sp(car(L),AL,1);
DC = append(DC,T);
}
return DC;
}
def sort_by_deg(F)
{
for ( T = F, S = []; T != []; T = cdr(T) )
if ( type(car(T)) != NUM )
S = cons(car(T),S);
N = length(S); W = newvect(N);
for ( I = 0; I < N; I++ )
W[I] = S[I];
V = var(W[0]);
for ( I = 0; I < N; I++ ) {
for ( J = I + 1, J0 = I; J < N; J++ )
if ( deg(W[J0],V) > deg(W[J],V) )
J0 = J;
if ( J0 != I ) {
T = W[I]; W[I] = W[J0]; W[J0] = T;
}
}
if ( ASCENT )
for ( I = N-1, S = []; I >= 0; I-- )
S = cons(W[I],S);
else
for ( I = 0, S = []; I < N; I++ )
S = cons(W[I],S);
return S;
}
def flatmf(L) {
for ( S = []; L != []; L = cdr(L) )
if ( type(F=car(car(L))) != NUM )
S = append(S,[F]);
return S;
}
def af(P,AL)
{
RESTIME=UFTIME=GCDTIME=SQTIME=0;
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);
for ( ; G != []; G = cdr(G) )
L = cons([car(G),M],L);
}
if ( dp_gr_print() )
print(["GCDTIME = ",GCDTIME,"UFTIME = ",UFTIME,"RESTIME = ",RESTIME,"SQTIME=",SQTIME]);
return L;
}
def af_sp(P,AL,HINT)
{
if ( !P || type(P) == NUM )
return [P];
P1 = simpcoef(simpalg(P));
return af_spmain(P1,AL,1,HINT,P1,[]);
}
def af_spmain(P,AL,INIT,HINT,PP,SHIFT)
{
if ( !P || type(P) == NUM )
return [P];
P = simpcoef(simpalg(P));
if ( DEG(P) == 1 )
return [simpalg(P)];
if ( AL == [] ) {
TTT = time()[0];
F = flatmf(ufctrhint_heuristic(P,HINT,PP,SHIFT));
UFTIME+=time()[0]-TTT;
return F;
}
A0 = car(AL); P0 = defpoly(A0);
V = var(P); V0 = var(P0);
P = simpcoef(P);
TTT = time()[0];
N = simpcoef(sp_norm(A0,V,subst(P,V,V-INIT*A0),AL));
RESTIME+=time()[0]-TTT;
TTT = time()[0];
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);
if ( !var(C) )
continue;
if ( D == 1 )
DCT = af_spmain(C,cdr(AL),1,HINT*deg(P0,V0),PP,cons([A0,INIT],SHIFT));
else
DCT = af_spmain(C,cdr(AL),1,1,C,[]);
for ( ; DCT != []; DCT = cdr(DCT) ) {
if ( !var(car(DCT)) )
continue;
if ( length(DCSQ) == 1 && length(DCT) == 1 )
U = simpcoef(G);
else {
S = subst(car(DCT),V,A);
if ( pra(G,S,AL) )
U = cr_gcda(S,G);
else
U = S;
}
if ( var(U) == V ) {
G = pdiva(G,U);
if ( D == 1 )
DCR = cons(simpcoef(U),DCR);
else {
T = af_spmain(U,AL,sp_next(INIT),HINT,PP,SHIFT);
DCR = append(DCR,T);
}
}
}
}
return DCR;
}
def sp_next(I)
{
if ( I > 0 )
return -I;
else
return -I+1;
}
extern USE_RES;
def sp_norm(A,V,P,AL)
{
P = simpcoef(simpalg(P));
if (USE_RES)
return sp_norm_res(A,V,P,AL);
else
return sp_norm_ch(A,V,P,AL);
}
def sp_norm_ch(A,V,P,AL)
{
Len = length(AL);
P0 = defpoly(A); V0 = var(P0);
PR = algptorat(P);
if ( nmono(P0) == 2 )
R = res(V0,PR,P0);
else if ( Len == 1 || Len == 3 )
R = res_ch1(V0,V,PR,P0);
else if ( Len == 2 ) {
P1 = defpoly(AL[1]);
R = norm_ch1(V0,V,PR,P0,P1);
} else
R = res(V0,PR,P0);
return rattoalgp(R,cdr(AL));
}
def sp_norm_res(A,V,P,AL)
{
Len = length(AL);
P0 = defpoly(A); V0 = var(P0);
PR = algptorat(P);
R = res(V0,PR,P0);
return rattoalgp(R,cdr(AL));
}
def simpalg(P) {
if ( !P )
return 0;
else if ( type(P) == NUM )
return ntype(P) <= 1 ? P : simpalgn(P);
else if ( type(P) == POLY )
return simpalgp(P);
else if ( type(P) == RAT )
return simpalg(nm(P))/simpalg(dn(P));
}
def simpalgp(P) {
for ( V = var(P), I = deg(P,V), T = 0; I >= 0; I-- )
if ( C = coef(P,I) )
T += simpalg(C)*V^I;
return T;
}
def simpalgn(A) {
if ( ntype(A) <= 1 )
return A;
else if ( type(R=algtorat(A)) == POLY )
return simpalgb(A);
else
return simpalgb(
invalgp(simpalgb(rattoalg(dn(R))))
*simpalgb(rattoalg(nm(R)))
);
}
def simpalgb(P) {
if ( ntype(P) <= 1 )
return P;
else {
A0 = getalg(P);
Used = [];
while ( A0 != [] ) {
S = algtorat(P);
for ( A = A0; A != []; A = cdr(A) )
S = srem(S,defpoly(car(A)));
P = rattoalg(S);
Used = append(Used,[car(A0)]);
A0 = setminus(getalg(P),Used);
}
return P;
}
}
def setminus(A,B) {
for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
for ( S = B, M = car(T); S != []; S = cdr(S) )
if ( car(S) == M )
break;
if ( S == [] )
R = cons(M,R);
}
return R;
}
def getalgp(P) {
if ( type(P) <= 1 )
return getalg(P);
else {
for ( V = var(P), I = deg(P,V), T = []; I >= 0; I-- )
if ( C = coef(P,I) )
T = union_sort(T,getalgp(C));
return T;
}
}
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;
}
}
/* C = union of A and B; A and B is sorted. C should also be sorted. */
def union_sort(A,B)
{
if ( A == [] )
return B;
else if ( B == [] )
return A;
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)
{
if ( ntype(A) <= 1 )
return 1/A;
P0 = defpoly(mainalg(A)); P = algtorat(A);
V = var(P0); G1 = P0;
G2 = DEG(P)>=DEG(P0)?srem(P,P0):P;
for ( H = 1, X = 1, U1 = 0, U2 = 1; deg(G2,V); ) {
D = DEG(G1)-DEG(G2); T = LCOEF(G2)^(D+1);
L = sqr(G1*T,G2); Q = car(L); R = car(cdr(L));
S = U1*T-U2*Q;
M = H^D; M1 = M*X;
G1 = G2; G2 = sdiv(R,M1);
U1 = U2; U2 = sdiv(S,M1);
X = LCOEF(G1); H = sdiv(X^D*H,M);
}
C = invalgp(rattoalg(srem(P*U2,P0)));
return C*rattoalg(U2);
}
def algptorat(P) {
if ( type(P) <= 1 )
return algtorat(P);
else {
for ( V = var(P), I = deg(P,V), T = 0; I >= 0; I-- )
if ( C = coef(P,I) )
T += algptorat(C)*V^I;
return T;
}
}
def rattoalgp(P,M) {
for ( T = M, S = P; T != []; T = cdr(T) )
S = subst(S,algtorat(FIRST(T)),FIRST(T));
return S;
}
def sortfs(L)
{
#define Factor(a) car(a)
#define Mult(a) car(cdr(a))
if ( type(TT(L)) == NUM )
L = cdr(L);
for ( N = 0, T = L; T != []; T = cdr(T), N++ );
P = newvect(N); P1 = newvect(N);
for ( I = 0, T = L, R = []; T != []; T = cdr(T) )
if ( Mult(car(T)) == 1 ) {
R = cons(car(T),R); N--;
} else {
P[I] = car(T); I++;
}
for ( J = 0, V = var(Factor(P[0])); J < N; J++ ) {
for ( K0 = K = J, D = deg(Factor(P[J]),V); K < N; K++ )
if ( deg(Factor(P[K]),V) < D ) {
K0 = K;
D = deg(Factor(P[K]),V);
}
P1[J] = P[K0];
if ( J != K0 )
P[K0] = P[J];
}
for ( I = N - 1; I >= 0; I-- )
R = cons(P1[I],R);
return R;
}
def pdiva(P1,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));
}
def pqra(P1,P2)
{
if ( type(P2) != POLY )
return [P1,0];
else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) )
return [0,P1];
else {
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))];
}
}
def pra(P1,P2,AL)
{
if ( type(P2) != POLY )
return 0;
else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) )
return P1;
else {
F1 = algptorat(P1); F2 = algptorat(P2); ML = map(defpoly,AL);
B = append(reverse(ML),[F2]);
V0 = var(P1);
V = cons(V0,map(algtorat,AL));
G = srem_by_nf(F1,B,V,2);
return simpalg(rattoalgp(G[0]/G[1],AL));
}
}
def sort_alg(VL)
{
N = length(VL); W = newvect(N,VL);
for ( I = 0; I < N; I++ ) {
for ( M = I, J = I + 1; J < N; J++ )
if ( W[J] > W[M] )
M = J;
if ( I != M ) {
T = W[I]; W[I] = W[M]; W[M] = T;
}
}
for ( I = N-1, L = []; I >= 0; I-- )
L = cons(W[I],L);
return L;
}
def asq(P)
{
P = simpalg(P);
if ( type(P) == NUM )
return [[1,1]];
else if ( getalgp(P) == [] )
return sqfr(P);
else {
V = var(P); N = DEG(P); A = newvect(N+1); B = newvect(N+1);
for ( I = 0, F = P; ;I++ ) {
if ( type(F) == NUM )
break;
F1 = diff(F,V);
GCD = cr_gcda(F,F1);
FLAT = pdiva(F,GCD);
if ( type(GCD) == NUM ) {
A[I] = F; B[I] = 1;
break;
}
for ( J = 1, F = GCD; ; J++ ) {
L = pqra(F,FLAT); Q = L[0]; R = L[1];
if ( R )
break;
else
F = Q;
}
A[I] = FLAT; B[I] = J;
}
for ( I = 0, J = 0, L = []; A[I]; I++ ) {
J += B[I];
if ( A[I+1] )
C = pdiva(A[I],A[I+1]);
else
C = A[I];
L = cons([C,J],L);
}
return L;
}
}
def ufctrhint1(P,HINT)
{
if ( deg(P,var(P)) == 168 ) {
SQ = sqfr(P);
if ( length(SQ) == 2 && SQ[1][1] == 1 )
return [[1,1],[P,1]];
else
return ufctrhint(P,HINT);
} else
return ufctrhint(P,HINT);
}
def simpcoef(P) {
return rattoalgp(ptozp(algptorat(P)),getalgp(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]];
for ( S = 0, L = SHIFT, AL = [], K = 1; L != []; L = cdr(L) ) {
A = car(L)[0]; S += A*car(L)[1]; AL = cons(A,AL);
K *= deg(defpoly(A),algtorat(A));
}
PPP = simpcoef(simpalg(subst(PP,V,V-S)));
for ( T = P-coef(P,D)*V^D, G = D; T; T -= coef(T,DT)*V^DT )
G = igcd(G,DT=deg(T,V));
if ( G == 1 ) {
if ( K*deg(PPP,V) != deg(P,V) )
PPP = cr_gcda(PPP,P);
return ufctrhint2(P,HINT,PPP,AL);
} else {
for ( S = 0, T = P; T; T -= coef(T,DT)*V^DT ) {
DT = deg(T,V);
S += coef(T,DT)*V^(DT/G);
}
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);
T = ufctrhint2(H,HINT,HH,AL);
DC = append(DC,T);
}
return DC;
}
}
def ufctrhint2(P,HINT,PP,AL)
{
if ( deg(P,var(P)) == HINT )
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);
return resfmain(P,L[2],L[0],DL);
}
def res_det(V,P1,P2)
{
D1 = deg(P1,V); D2 = deg(P2,V);
M = newmat(D1+D2,D1+D2);
for ( J = 0; J <= D2; J++ )
M[0][J] = coef(P2,D2-J,V);
for ( I = 1; I < D1; I++ )
for ( J = 0; J <= D2; J++ )
M[I][I+J] = M[0][J];
for ( J = 0; J <= D1; J++ )
M[D1][J] = coef(P1,D1-J,V);
for ( I = 1; I < D2; I++ )
for ( J = 0; J <= D1; J++ )
M[D1+I][I+J] = M[D1][J];
return det(M);
}
def norm_ch1(V0,VM,P,P0,PR) {
D = deg(P,V0); D0 = deg(P0,V0); DM = deg(P,VM); N = DM*D0;
X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
Min = -idiv(N,2);
C = coef(P,D,V0);
for ( I = J = 0; I <= N; J++ ) {
if ( PRINT )
print([J,N]);
T=J+Min;
if ( subst(C,VM,T) ) {
U[I] = srem(res(V0,subst(P,VM,T),P0),PR);
X[I++] = T;
}
}
for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
for ( J = 0, T = U[I]; J < I; J++ )
T = sdiv(T-V[J],X[I]-X[J]);
V[I] = T;
M *= (VM-X[I-1]);
S += T*M;
}
return S;
}
def norm_ch2(V0,VM,P,P0,PR) {
D = deg(P,V0); D0 = deg(P0,V0); DM = deg(P,VM); N = DM*D0;
X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
Min = -idiv(N,2);
C = coef(P,D,V0);
for ( I = J = 0; I <= N; J++ ) {
T=J+Min;
if ( subst(C,VM,T) ) {
U[I] = srem(res_det(V0,subst(P,VM,T),P0),PR);
X[I++] = T;
}
}
for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
for ( J = 0, T = U[I]; J < I; J++ )
T = sdiv(T-V[J],X[I]-X[J]);
V[I] = T;
M *= (VM-X[I-1]);
S += T*M;
}
return S;
}
def res_ch1(V0,VM,P,P0) {
D = deg(P,V0); D0 = deg(P0,V0); N = deg(P,VM)*D0+deg(P0,VM)*D;
X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
Min = -idiv(N,2);
C = coef(P,D,V0); C0 = coef(P0,D0,V0);
for ( I = J = 0; I <= N; J++ ) {
if ( PRINT )
print([J,N]);
T=J+Min;
if ( subst(C,VM,T) && subst(C0,VM,T) ) {
U[I] = res(V0,subst(P,VM,T),subst(P0,VM,T));
X[I++] = T;
}
}
for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
for ( J = 0, T = U[I]; J < I; J++ )
T = sdiv(T-V[J],X[I]-X[J]);
V[I] = T;
M *= (VM-X[I-1]);
S += T*M;
}
return S;
}
def res_ch(V0,VM,P,P0) {
D = deg(P,V0); D0 = deg(P0,V0); N = deg(P,VM)*D0+deg(P0,VM)*D;
X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
Min = -idiv(N,2);
C = coef(P,D,V0); C0 = coef(P0,D0,V0);
for ( I = J = 0; I <= N; J++ ) {
T=J+Min;
if ( subst(C,VM,T) && subst(C0,VM,T) ) {
U[I] = res_det(V0,subst(P,VM,T),subst(P0,VM,T));
X[I++] = T;
}
}
for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
for ( J = 0, T = U[I]; J < I; J++ )
T = sdiv(T-V[J],X[I]-X[J]);
V[I] = T;
M *= (VM-X[I-1]);
S += T*M;
}
return S;
}
def norm_ch2_lag(V,VM,P,P0,PR) {
D0 = deg(P0,V); DM = deg(P,VM); N = DM*D0;
Min = -idiv(N,2);
for ( A = 1, I = 0; I <= N; I++ )
A *= (VM-I-Min);
for ( I = 0, S = 0; I <= N; I++ ) {
R = res_det(V,subst(P,VM,I+Min),P0);
R = srem(R,PR);
T = sdiv(A,VM-I-Min);
S += R*T/subst(T,VM,I+Min);
}
return S;
}
def norm_ch_lag(V,VM,P,P0) {
D0 = deg(P0,V); DM = deg(P,VM); N = DM*D0;
Min = -idiv(N,2);
for ( A = 1, I = 0; I <= N; I++ )
A *= (VM-I-Min);
for ( I = 0, S = 0; I <= N; I++ ) {
R = res_det(V,subst(P,VM,I+Min),P0);
T = sdiv(A,VM-I-Min);
S += R*T/subst(T,VM,I+Min);
}
return S;
}
def cr_gcda(P1,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);
NEXT = length(EXT);
if ( deg(P1,V) < deg(P2,V) ) {
T = P1; P1 = P2; P2 = T;
}
G1 = ptozp(algptorat(P1)); G2 = ptozp(algptorat(P2));
for ( ML = VL = [], T = reverse(EXT); T != []; T = cdr(T) ) {
ML = cons(defpoly(car(T)),ML);
VL = cons(algptorat(car(T)),VL);
}
DL = [coef(G1,deg(G1,V),V),coef(G2,deg(G2,V),V)];
for ( T = EXT; T != []; T = cdr(T) ) {
DL = cons(discr(sp_minipoly(car(T),EXT)),DL);
C = LCOEF(defpoly(car(T)));
if ( C != 1 && C != -1 )
DL = cons(C,DL);
}
TIME = time()[0];
for ( D = deg(P1,V)+1, I = 0; ; I++ ) {
MOD = lprime(I);
for ( J = 0; J < length(DL); J++ )
if ( !(DL[J] % MOD) )
break;
if ( J != length(DL) )
continue;
Ord = 2; NOSUGAR = 1;
T = ag_mod(G1 % MOD,G2 % MOD,ML,VL,MOD);
if ( dp_gr_print() )
print(".");
if ( !T )
continue;
T = (T*inv(coef(T,deg(T,V),V),MOD))%MOD;
if ( deg(T,V) > D )
continue;
else if ( deg(T,V) < D ) {
IMAGE = T; M = MOD; D = deg(T,V);
} else {
L = cr(IMAGE,M,T,MOD); IMAGE = L[0]; M = L[1];
}
F = intptoratp(IMAGE,M,calcb(M));
if ( F != [] ) {
F = ptozp(F);
DIV = rattoalgp(F,EXT);
if ( type(DIV) == 1 )
return 1;
/*
if ( srem_simp(G1,F,V,ML) )
continue;
if ( srem_simp(G2,F,V,ML) )
continue;
*/
if ( srem_by_nf(G1,reverse(cons(F,ML)),cons(V,VL),2)[0] )
continue;
if ( srem_by_nf(G2,reverse(cons(F,ML)),cons(V,VL),2)[0] )
continue;
TIME = time()[0]-TIME;
if ( dp_gr_print() )
print([TIME]);
GCDTIME += TIME;
return DIV;
}
}
}
def srem_simp(F1,F2,V,D)
{
D2 = deg(F2,V); C = coef(F2,D2);
while ( (D1 = deg(F1,V)) >= D2 ) {
F1 -= coef(F1,D1)/C*V^(D1-D2)*F2;
F1 = simp_by_dp(F1,D);
}
return F1;
}
def member(E,L)
{
for ( ; L != []; L = cdr(L) )
if ( E == car(L) )
return 1;
return 0;
}
def discr(P) {
V = var(P);
return res(V,P,diff(P,V));
}
def sp_minipoly(A,EXT)
{
while ( car(EXT) != A )
EXT = cdr(EXT);
for ( M = x-A; EXT != []; EXT = cdr(EXT) )
M = sp_norm(car(EXT),x,M,EXT);
F = sqfr(M);
return F[1][0];
}
def cr(F1,M1,F2,M2)
{
K = inv(M1 % M2,M2);
M3 = M1*M2;
F3 = (F1 + (F2-(F1%M2))*K*M1) % M3;
return [F3,M3];
}
#define ABS(a) ((a)>=0?(a):(-a))
#if 0
def calcb(M) {
setprec(800);
return pari(floor,eval((M/2)^(1/2)));
}
#endif
def calcb(M) {
N = 2*M;
T = sp_sqrt(N);
if ( T^2 <= N && N < (T+1)^2 )
return idiv(T,2);
else
error("afo");
}
def sp_sqrt(A) {
for ( J = 0, T = A; T >= 2^27; J++ ) {
T = idiv(T,2^27)+1;
}
for ( I = 0; T >= 2; I++ ) {
S = idiv(T,2);
if ( T = S+S )
T = S;
else
T = S+1;
}
X = (2^27)^idiv(J,2)*2^idiv(I,2);
while ( 1 ) {
if ( (Y=X^2) < A )
X += X;
else if ( Y == A )
return X;
else
break;
}
while ( 1 )
if ( (Y = X^2) <= A )
return X;
else
X = idiv(A + Y,2*X);
}
def intptoratp(P,M,B) {
if ( type(P) == 1 ) {
L = inttorat(P,M,B);
if ( L == 0 )
return [];
else
return L[0]/L[1];
} else {
V = var(P);
S = 0;
while ( P ) {
D = deg(P,V);
C = coef(P,D,V);
T = intptoratp(C,M,B);
if ( T == [] )
return [];
S += T*V^D;
P -= C*V^D;
}
return S;
}
}
def ltoalg(L) {
F = L[0]; V = reverse(L[1]);
N = length(V)-1;
for ( I = 0, G = F; I < N; I++ ) {
D = car(G);
A = newalg(D); V = var(D);
for ( G = reverse(cdr(G)), T = []; G != []; G = cdr(G) )
T = cons(subst(car(G),V,A),T);
G = T;
}
return G;
}
/*
def ag_mod(F1,F2,D,MOD)
{
if ( length(D) == 1 )
return ag_mod_single(F1,F2,D,MOD);
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
for ( T = reverse(D), S = []; T != []; T = cdr(T) ) {
U = car(T);
S = cons((inv(LCOEF(U),MOD)*U) % MOD,S);
}
D = S;
while ( 1 ) {
F = srem_simp_mod(F1,F2,V,D,MOD);
if ( !F )
return F2;
if ( !deg(F,V) )
return 1;
C = LCOEF(F);
INV = inverse_by_gr_mod(C,D,MOD);
if ( !INV )
return 0;
F = simp_by_dp_mod(F*INV,D,MOD);
F = (inv(LCOEF(F),MOD)*F) % MOD;
F1 = F2; F2 = F;
}
}
*/
def ag_mod(F1,F2,D,VL,MOD)
{
if ( length(D) == 1 )
return ag_mod_single6(F1,F2,D,MOD);
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
for ( T = reverse(D), S = []; T != []; T = cdr(T) ) {
U = car(T);
S = cons((inv(LCOEF(U),MOD)*U) % MOD,S);
}
D = S;
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);
dp_gr_flags(FLAGS);
if ( length(G) == 1 )
return 1;
if ( length(G) == N ) {
for ( T = G; T != []; T = cdr(T) )
if ( member(V,vars(car(T))) )
return car(T);
}
}
}
def srem_simp_mod(F1,F2,V,D,MOD)
{
D2 = deg(F2,V); C = coef(F2,D2);
while ( (D1 = deg(F1,V)) >= D2 ) {
F1 -= coef(F1,D1)/C*V^(D1-D2)*F2;
F1 = simp_by_dp_mod(F1,D,MOD);
}
return F1;
}
def ag_mod_single(F1,F2,D,MOD)
{
TD = TI = TM = 0;
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
FLAGS = dp_gr_flags(); dp_gr_flags(["Reverse",1,"NoSugar",1]);
G = dp_gr_mod_main([F1,F2,D],[V,var(D)],0,MOD,2);
dp_gr_flags(FLAGS);
if ( length(G) == 1 )
return 1;
if ( length(G) != 2 )
return 0;
if ( vars(G[0]) == [var(D)] )
return G[1];
else
return G[0];
}
def ag_mod_single2(F1,F2,D,MOD)
{
TD = TI = TM = 0;
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
while ( 1 ) {
T0 = time()[0];
F = srem((srem(F1,F2) % MOD),D) % MOD;
TD += time()[0] - T0;
if ( !F ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI]);
return F2;
}
if ( !deg(F,V) ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI]);
return 1;
}
C = LCOEF(F);
T0 = time()[0];
INV = inva_mod(D,MOD,C);
TI += time()[0] - T0;
if ( !INV )
return 0;
T0 = time()[0];
F = remc_mod((INV*F) % MOD,D,MOD);
TM += time()[0] - T0;
F1 = F2; F2 = F;
}
}
def ag_mod_single3(F1,F2,D,MOD)
{
TD = TI = TM = 0;
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
while ( 1 ) {
if ( !D2 )
return 1;
while ( D1 >= D2 ) {
F = srem((coef(F2,D2,V)*F1-coef(F1,D1,V)*F2*V^(D1-D2))%MOD,D)%MOD;
F1 = F; D1 = deg(F1,V);
}
if ( !F1 ) {
INV = inva_mod(D,MOD,coef(F2,D2,V));
if ( dp_gr_print() )
print(".");
return srem((INV*F2) % MOD,D)%MOD;
} else {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
}
}
def ag_mod_single4(F1,F2,D,MOD)
{
if ( !F1 )
return F2;
if ( !F2 )
return F1;
TD = TI = TM = TR = 0;
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
while ( 1 ) {
T0 = time()[0]; R = srem(F1,F2); TR += time()[0] - T0;
T0 = time()[0]; F = srem(R % MOD,D) % MOD; TD += time()[0] - T0;
if ( !F ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
return F2;
}
if ( !deg(F,V) ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
return 1;
}
C = LCOEF(F);
T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
if ( !INV )
return 0;
T0 = time()[0]; F = srem((INV*F) % MOD,D) % MOD; TM += time()[0] - T0;
F1 = F2; F2 = F;
}
}
def ag_mod_single5(F1,F2,D,MOD)
{
TD = TI = TM = TR = 0;
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
while ( 1 ) {
T0 = time()[0];
D1 = deg(F1,V); D2 = deg(F2,V); F = F1;
while ( D1 >= D2 ) {
R = (F-coef(F,D1,V)*F2*V^(D1-D2))%MOD;
D1 = deg(R,V); HC = coef(R,D1,V);
F = (R - HC*V^D1) + (srem(HC,D)%MOD)*V^D1;
}
TR += time()[0] - T0;
T0 = time()[0]; F = srem(R % MOD,D) % MOD; TD += time()[0] - T0;
if ( !F ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
return F2;
}
if ( !deg(F,V) ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
return 1;
}
C = LCOEF(F);
T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
if ( !INV )
return 0;
T0 = time()[0]; F = srem((INV*F) % MOD,D) % MOD; TM += time()[0] - T0;
F1 = F2; F2 = F;
}
}
def ag_mod_single6(F1,F2,D,MOD)
{
TD = TI = TM = TR = 0;
V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
if ( D1 < D2 ) {
T = F1; F1 = F2; F2 = T;
T = D1; D1 = D2; D2 = T;
}
F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
while ( 1 ) {
T0 = time()[0];
D1 = deg(F1,V); D2 = deg(F2,V); F = F1;
while ( D1 >= D2 ) {
R = (F-coef(F,D1,V)*F2*V^(D1-D2))%MOD;
D1 = deg(R,V); HC = coef(R,D1,V);
/* F = (R - HC*V^D1) + (srem_mod(HC,D,MOD))*V^D1; */
F = remc_mod(R,D,MOD);
}
TR += time()[0] - T0;
T0 = time()[0]; F = remc_mod(R%MOD,D,MOD); TD += time()[0] - T0;
if ( !F ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
return F2;
}
if ( !deg(F,V) ) {
if ( dp_gr_print() )
print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
return 1;
}
C = LCOEF(F);
T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
if ( !INV )
return 0;
T0 = time()[0]; F = remc_mod((INV*F)%MOD,D,MOD); TM += time()[0] - T0;
F1 = F2; F2 = F;
}
}
def inverse_by_gr_mod(C,D,MOD)
{
Ord = 2;
dp_gr_flags(["NoSugar",1]);
G = dp_gr_mod_main(cons(x*C-1,D),cons(x,vars(D)),0,MOD,Ord);
dp_gr_flags(["NoSugar",0]);
if ( length(G) == 1 )
return 1;
else if ( length(G) == length(D)+1 ) {
for ( T = G; T != []; T = cdr(T) )
if ( member(x,vars(car(G))) )
break;
T = car(G);
if ( type(coef(T,1,x)) != NUM )
return 0;
else
return coef(T,0,x);
} else
return 0;
}
def simp_by_dp(F,D)
{
for ( T = D; T != []; T = cdr(T) )
F = srem(F,car(T));
return F;
}
def simp_by_dp_mod(F,D,MOD)
{
F %= MOD;
for ( T = D; T != []; T = cdr(T) )
F = srem(F,car(T)) % MOD;
return F;
}
def remc_mod(P,D,M)
{
V = var(P);
if ( !V || V == var(D) )
return srem_mod(P,D,M);
for ( I = deg(P,V), S = 0; I >= 0; I-- )
if ( C = coef(P,I,V) )
S += srem_mod(C,D,M)*V^I;
return S;
}
def rem_mod(C,D,M)
{
V = var(D);
D2 = deg(D,V);
while ( (D1 = deg(C,V)) >= D2 ) {
C -= (D*V^(D1-D2)*coef(C,D1,V))%M;
C %= M;
}
return C;
}
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++;
T = nfctr_mod(N,M);
if ( T < Len ) {
Len = T; M0 = M;
}
}
}
S = spm(L,V,M0);
T = resfctr_mod(F,S,M0);
return [T,S,M0];
}
def resfctr_mod(F,L,M)
{
for ( T = L, R = []; T != []; T = cdr(T) ) {
U = car(T); MP = U[0]; W = U[1];
for ( A = W, B = F; A != []; A = cdr(cdr(A)) )
B = sremm(subst(B,A[0],A[1]),MP,M);
C = res(var(MP),B,MP) % M;
R = cons(flatten(cdr(modfctr(C,M))),R);
}
return reverse(R);
}
def flatten(L)
{
for ( T = L, R = []; T != []; T = cdr(T) )
R = cons(car(car(T)),R);
return R;
}
def spm(L,V,M)
{
if ( length(V) == 1 ) {
U = modfctr(car(L),M);
for ( T = cdr(U), R = []; T != []; T = cdr(T) ) {
S = car(T);
R = cons([subst(S[0],var(S[0]),a_),[var(S[0]),a_]],R);
}
return R;
}
L1 = spm(cdr(L),cdr(V),M);
F0 = car(L); V0 = car(V); VR = cdr(V);
for ( T = L1, R = []; T != []; T = cdr(T) ) {
S = car(T);
F1 = subst(F0,S[1]);
U = fctr_mod(F1,V0,S[0],M);
VS = var(S[0]);
for ( W = U; W != []; W = cdr(W) ) {
A = car(car(W));
if ( deg(A,V0) == 1 ) {
A = monic_mod(A,V0,S[0],M);
R = cons([S[0],append([V0,-coef(A,0,V0)],S[1])],R);
} else {
B = pe_mod(A,S[0],M);
MP = B[0]; VMP = var(MP); NV = B[1];
for ( C = S[1], D = []; C != []; C = cdr(cdr(C)) ) {
G = subst(sremm(subst(C[1],VS,NV[1]),MP,M),VMP,VS);
D = append([C[0],G],D);
}
R = cons([subst(MP,VMP,VS),
append([B[2][0],subst(B[2][1],VMP,VS)],D)],R);
}
}
}
return R;
}
def pe_mod(F,G,M)
{
V = var(G); W = car(setminus(vars(F),[V]));
NG = deg(G,V); NF = deg(F,W); N = NG*NF;
X = prim;
while ( 1 ) {
D = mrandompoly(N,X,M);
if ( irred_check(D,M) )
break;
}
L = fctr_mod(G,V,D,M);
for ( T = L; T != []; T = cdr(T) ) {
U = car(car(T));
if ( deg(U,V) == 1 )
break;
}
U = monic_mod(U,V,D,M); RV = -coef(U,0,V);
L = fctr_mod(sremm(subst(F,V,RV),D,M),W,D,M);
for ( T = L; T != []; T = cdr(T) ) {
U = car(car(T));
if ( deg(U,W) == 1 )
break;
}
U = monic_mod(U,W,D,M); RW = -coef(U,0,W);
return [D,[V,RV],[W,RW]];
}
def fctr_mod(F,V,D,M)
{
if ( V != x ) {
F = subst(F,V,x); V0 = V; V = x;
} else
V0 = x;
F = monic_mod(F,V,D,M);
L = sqfr_mod(F,V,D,M);
for ( R = [], T = L; T != []; T = cdr(T) ) {
S = car(T); A = S[0]; E = S[1];
B = ddd_mod(A,V,D,M);
R = append(append_mult(B,E),R);
}
if ( V0 != x ) {
for ( R = reverse(R), T = []; R != []; R = cdr(R) )
T = cons([subst(car(R)[0],x,V0),car(R)[1]],T);
R = T;
}
return R;
}
def append_mult(L,E)
{
for ( T = L, R = []; T != []; T = cdr(T) )
R = cons([car(T),E],R);
return R;
}
def sqfr_mod(F,V,D,M)
{
setmod(M);
F = sremm(F,D,M);
F1 = sremm(diff(F,V),D,M);
F1 = sremm(F1*inva_mod(D,M,LCOEF(F1)),D,M);
if ( F1 ) {
F2 = ag_mod_single4(F,F1,[D],M);
FLAT = sremm(sdivm(F,F2,M,V),D,M);
I = 0; L = [];
while ( deg(FLAT,V) ) {
while ( 1 ) {
QR = sqrm(F,FLAT,M,V);
if ( !sremm(QR[1],D,M) ) {
F = sremm(QR[0],D,M); I++;
} else
break;
}
if ( !deg(F,V) )
FLAT1 = 1;
else
FLAT1 = ag_mod_single4(F,FLAT,[D],M);
G = sremm(sdivm(FLAT,FLAT1,M,V),D,M);
FLAT = FLAT1;
L = cons([G,I],L);
}
}
if ( deg(F,V) ) {
T = sqfr_mod(pthroot_p_mod(F,V,D,M),V,D,M);
for ( R = []; T != []; T = cdr(T) ) {
H = car(T); R = cons([H[0],M*H[1]],R);
}
} else
R = [];
return append(L,R);
}
def pthroot_p_mod(F,V,D,M)
{
for ( T = F, R = 0; T; ) {
D1 = deg(T,V); C = coef(T,D1,V); T -= C*V^D1;
R += pthroot_n_mod(C,D,M)*V^idiv(D1,M);
}
return R;
}
def pthroot_n_mod(C,D,M)
{
pwr_n_mod(C,D,M,deg(D,var(D))-1);
}
def pwr_n_mod(C,D,M,N)
{
if ( N == 0 )
return 1;
else if ( N == 1 )
return C;
else {
QR = iqr(N,2);
T = pwr_n_mod(C,D,M,QR[0]);
S = sremm(T^2,D,M);
if ( QR[1] )
return sremm(S*C,D,M);
else
return S;
}
}
def pwr_p_mod(P,A,V,D,M,N)
{
if ( N == 0 )
return 1;
else if ( N == 1 )
return P;
else {
QR = iqr(N,2);
T = pwr_p_mod(P,A,V,D,M,QR[0]);
S = sremm(sremm(sremm(T^2,D,M),A,M,V),D,M);
if ( QR[1] )
return sremm(sremm(sremm(S*P,D,M),A,M,V),D,M);
else
return S;
}
}
def qmat_mod(F,V,D,M)
{
R = tab_mod(F,V,D,M);
Q = newmat(N,N);
for ( J = 0; J < N; J++ )
for ( I = 0, T = R[J]; I < N; I++ ) {
Q[I][J] = coef(T,I);
}
for ( I = 0; I < N; I++ )
Q[I][I] = (Q[I][I]+(M-1))%M;
return Q;
}
def tab_mod(F,V,D,M)
{
MD = M^deg(D,var(D));
N = deg(F,V);
F = sremm(F*inva_mod(D,M,coef(F,N,V)),D,M);
R = newvect(N); R[0] = 1;
R[1] = pwr_mod(V,F,V,D,M,MD);
for ( I = 2; I < N; I++ )
R[I] = sremm(sremm(R[1]*R[I-1],F,M),D,M);
return R;
}
def ddd_mod(F,V,D,M)
{
if ( deg(F,V) == 1 )
return [F];
TAB = tab_mod(F,V,D,M);
for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
for ( T = 0, K = 0; K <= deg(W,V); K++ )
if ( C = coef(W,K,V) )
T = sremm(T+TAB[K]*C,D,M);
W = T;
GCD = ag_mod_single4(F,monic_mod(W-V,V,D,M),[D],M);
if ( deg(GCD,V) ) {
L = append(berlekamp(GCD,V,I,TAB,D,M),L);
F = sremm(sdivm(F,GCD,M,V),D,M);
W = sremm(sremm(W,F,M,V),D,M);
}
}
if ( deg(F,V) )
return cons(F,L);
else
return L;
}
def monic_mod(F,V,D,M) {
if ( !F || !deg(F,V) )
return F;
return sremm(F*inva_mod(D,M,coef(F,deg(F,V),V)),D,M);
}
def berlekamp(F,V,E,TAB,D,M)
{
N = deg(F,V);
Q = newmat(N,N);
for ( J = 0; J < N; J++ ) {
T = sremm(sremm(TAB[J],F,M,V),D,M);
for ( I = 0; I < N; I++ ) {
Q[I][J] = coef(T,I);
}
}
for ( I = 0; I < N; I++ )
Q[I][I] = (Q[I][I]+(M-1))%M;
L = nullspace(Q,D,M); MT = L[0]; IND = L[1];
NF0 = N/E;
PS = null_to_poly(MT,IND,V,M);
R = newvect(NF0); R[0] = monic_mod(F,V,D,M);
for ( I = 1, NF = 1; NF < NF0 && I < NF0; I++ ) {
PSI = PS[I];
MP = minipoly_mod(PSI,F,V,D,M);
ROOT = find_root(MP,V,D,M); NR = length(ROOT);
for ( J = 0; J < NF; J++ ) {
if ( deg(R[J],V) == E )
continue;
for ( K = 0; K < NR; K++ ) {
GCD = ag_mod_single4(R[J],PSI-ROOT[K],[D],M);
if ( deg(GCD,V) > 0 && deg(GCD,V) < deg(R[J],V) ) {
Q = sremm(sdivm(R[J],GCD,M,V),D,M);
R[J] = Q; R[NF++] = GCD;
}
}
}
}
return vtol(R);
}
def null_to_poly(MT,IND,V,M)
{
N = size(MT)[0];
for ( I = 0, J = 0; I < N; I++ )
if ( IND[I] )
J++;
R = newvect(J);
for ( I = 0, L = 0; I < N; I++ ) {
if ( !IND[I] )
continue;
for ( J = K = 0, T = 0; J < N; J++ )
if ( !IND[J] )
T += MT[K++][I]*V^J;
else if ( J == I )
T += (M-1)*V^I;
R[L++] = T;
}
return R;
}
def minipoly_mod(P,F,V,D,M)
{
L = [[1,1]]; P0 = P1 = 1;
while ( 1 ) {
P0 *= V;
P1 = sremm(sremm(P*P1,F,M,V),D,M);
L1 = lnf_mod(P0,P1,L,V,D,M); NP0 = L1[0]; NP1 = L1[1];
if ( !NP1 )
return NP0;
else
L = lnf_insert([NP0,NP1],L,V);
}
}
def lnf_mod(P0,P1,L,V,D,M)
{
NP0 = P0; NP1 = P1;
for ( T = L; T != []; T = cdr(T) ) {
Q = car(T);
D1 = deg(NP1,V);
if ( D1 == deg(Q[1],V) ) {
C = coef(Q[1],D1,V);
INV = inva_mod(D,M,M-C); H = sremm(coef(NP1,D1,V)*INV,D,M);
NP0 = sremm(NP0+Q[0]*H,D,M);
NP1 = sremm(NP1+Q[1]*H,D,M);
}
}
return [NP0,NP1];
}
def lnf_insert(P,L,V)
{
if ( L == [] )
return [P];
else {
P0 = car(L);
if ( deg(P0[1],V) > deg(P[1],V) )
return cons(P0,lnf_insert(P,cdr(L),V));
else
return cons(P,L);
}
}
def find_root(P,V,D,M)
{
L = c_z(P,V,1,D,M);
for ( T = L, U = []; T != []; T = cdr(T) ) {
S = monic_mod(car(T),V,D,M); U = cons(-coef(S,0,V),U);
}
return U;
}
def c_z(F,V,E,D,M)
{
N = deg(F,V);
if ( N == E )
return [F];
Q = M^deg(D,var(D));
K = idiv(N,E);
L = [F];
while ( 1 ) {
W = mrandomgfpoly(2*E,V,D,M);
if ( M == 2 ) {
W = monic_mod(tr_mod(W,F,V,D,M,N-1),V,D,M);
} else {
/* W = monic_mod(pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2))-1,V,D,M); */
/* T = pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2)); */
T = pwr_mod(W,F,V,D,M,idiv(Q^E-1,2));
W = monic_mod(T-1,V,D,M);
}
if ( !W )
continue;
G = ag_mod_single4(F,W,[D],M);
if ( deg(G,V) && deg(G,V) < N ) {
L1 = c_z(G,V,E,D,M);
L2 = c_z(sremm(sdivm(F,G,M,V),D,M),V,E,D,M);
return append(L1,L2);
}
}
}
def tr_mod(P,F,V,D,M,N)
{
for ( I = 1, S = P, W = P; I <= N; I++ ) {
W = sremm(sremm(W^2,F,M,V),D,M);
S = sremm(S+W,D,M);
}
return S;
}
def mrandomgfpoly(N,V,D,M)
{
W = var(D); ND = deg(D,W);
for ( I = N-2, S = V^(N-1); I >= 0; I-- )
S += randompoly(ND,W,M)*V^I;
return S;
}
def randompoly(N,V,M)
{
for ( I = 0, S = 0; I < N; I++ )
S += (random()%M)*V^I;
return S;
}
def mrandompoly(N,V,M)
{
for ( I = N-1, S = V^N; I >=0; I-- )
S += (random()%M)*V^I;
return S;
}
def srem_by_nf(P,B,V,O) {
dp_ord(O); DP = dp_ptod(P,V);
N = length(B); DB = newvect(N);
for ( I = N-1, IL = []; I >= 0; I-- ) {
DB[I] = dp_ptod(B[I],V);
IL = cons(I,IL);
}
L = dp_true_nf(IL,DP,DB,1);
return [dp_dtop(L[0],V),L[1]];
}
end$