/*
 * 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/asir2018/lib/sp,v 1.1 2018/09/19 05:45:08 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 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;
	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;
    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);
		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];
    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,[]);
}

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(factpol,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;
		SpOrd = 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,SpOrd);
		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)
{
	SpOrd = 2;
	dp_gr_flags(["NoSugar",1]);
	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;
	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);
  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];
}

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);
      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);
		F1 = subst(F0,S[1]);
		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);
				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$