/*
 * 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/nf,v 1.1 2018/09/19 05:45:08 noro Exp $
*/
extern IDIVV_REM,IDIVV_HIST;

def nf_mod(B,G,MOD,DIR,PS)
{
	setmod(MOD);
	for ( D = 0, H = []; G; ) {
		for ( U = 0, L = B; L != []; L = cdr(L) ) {
			if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
				R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
				CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
				U = G-CR*R;
				print(".",2);
				if ( !U )
					return D;
				break;
			}
		}
		if ( U )
			G = U;
		else {
			D += dp_hm(G); G = dp_rest(G); print(!D,2);
		}
	}
	return D;
}

def nf_mod_ext(B,G,MOD,DIR,PS)
{
	setmod(MOD);
	for ( D = 0, H = []; G; ) {
		for ( U = 0, L = B; L != []; L = cdr(L) ) {
			if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
				R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
				CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
				U = G-CR*R;
				H = cons([CR,strtov("t"+rtostr(car(L)))],H);
				print([length(H)]);
				if ( !U )
					return [D,reverse(H)];
				break;
			}
		}
		if ( U )
			G = U;
		else {
			D += dp_hm(G); G = dp_rest(G);
		}
	}
	return [D,reverse(H)];
}

def nf(B,G,M,PS)
{
	for ( D = 0, H = []; G; ) {
		for ( U = 0, L = B; L != []; L = cdr(L) ) {
			if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
				GCD = igcd(dp_hc(G),dp_hc(R));
				CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
				U = CG*G-dp_subd(G,R)*CR*R;
				H = cons(car(L),H);
				if ( !U )
					return [D,M,reverse(H)];
				D *= CG; M *= CG;
				break;
			}
		}
		if ( U )
			G = U;
		else {
			D += dp_hm(G); G = dp_rest(G);
		}
	}
	return [D,M,reverse(H)];
}
extern M_NM,M_DN;

def nf_demand(B,G,DIR,PSH)
{
	T1 = T2 = T3 = T4 = 0;
	Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
	for ( D = 0, H = []; G; ) {
		for ( U = 0, L = B; L != []; L = cdr(L) ) {
			if ( dp_redble(G,PSH[car(L)]) > 0 ) {
				R = bload(DIR+rtostr(car(L)));
				H = cons(car(L),H);
				print([length(H),dp_mag(dp_hm(G))]);
				GCD = igcd(dp_hc(G),dp_hc(R));
				CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
				T0 = time()[0]; A1 = CG*G; T1 += time()[0]-T0;
				T0 = time()[0]; A2 = dp_subd(G,R)*CR*R; T2 += time()[0]-T0;
				T0 = time()[0]; U = A1-A2; T3 += time()[0]-T0;
				if ( !U )
					return [D,reverse(H),T1,T2,T3,T4];
				D *= CG;
				break;
			}
		}
		if ( U ) {
			G = U;
			if ( D ) {
				if ( dp_mag(dp_hm(D)) > Mag ) {
					T0 = time()[0];
					print("ptozp2");
					T = dp_ptozp2(D,G); D = T[0]; G = T[1];
					T4 += time()[0]-T0;
					Mag = idiv(dp_mag(dp_hm(D))*M_NM,M_DN);
				}
			} else {
				if ( dp_mag(dp_hm(G)) > Mag ) {
					T0 = time()[0];
					print("ptozp");
					G = dp_ptozp(G);
					T4 += time()[0]-T0;
					Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
				}
			}
		} else {
			D += dp_hm(G); G = dp_rest(G);
		}
	}
	return [D,reverse(H),T1,T2,T3,T4];
}

def nf_demand_d(B,G,DIR,NDIR,PDIR,PSH,PROC)
{
	INDEX = 0;
	T1 = T2 = 0;
/*	dp_gr_flags(["Dist",PROC]); */
	if ( PROC != [] ) {
		PROC = newvect(length(PROC),PROC);
		NPROC = size(PROC)[0];
	} else
		NPROC = 0;
	Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
	Kara = dp_set_kara()*27; /* XXX */
	D = [0,0]; R = [1,G];
	print("");
	for ( H = []; R[1]; ) {
		for ( U = 0, L = B; L != []; L = cdr(L) ) {
			if ( dp_redble(R[1],PSH[car(L)]) > 0 ) {
				Red = bload(DIR+rtostr(car(L)));
				H = cons(car(L),H);
/*				D0=dp_mag(D[0]*<<0>>); D1=dp_mag(dp_hm(D[1]));
				R0=dp_mag(R[0]*<<0>>); R1=dp_mag(dp_hm(R[1]));
				print([car(L),length(H),[D0,D1,dp_nt(D[1])],[R0,R1,dp_nt(R[1])]],2); */

				GCD = igcd(dp_hc(R[1]),dp_hc(Red));
				CR = idiv(dp_hc(Red),GCD);
				CRed = idiv(dp_hc(R[1]),GCD);

				T0 = time()[3];
				if ( (PROC != []) && (dp_mag(dp_hm(Red)) > Kara) ) {
					print("d",2);
					rpc(PROC[0],"dp_imul_index",CRed,car(L));
					U = CR*R[1] - dp_subd(R[1],Red)*rpcrecv(PROC[0]);
				} else {
					print("l",2);
					U = CR*R[1] - dp_subd(R[1],Red)*CRed*Red;
				}
				TT = time()[3]-T0; T1 += TT; /* print(TT); */
				if ( !U )
					return [D,reverse(H),T1,T2];
				break;
			}
		}
		if ( U ) {
			if ( (HMAG=dp_mag(dp_hm(U))) > Mag ) {
				T0 = time()[3];
				if ( HMAG > Kara ) {
					print("D",2);
					T = dp_ptozp_d(U,PROC,NPROC);
				} else {
					print("L",2);
					T = dp_ptozp(U);
				}
				TT = time()[3]-T0; T2 += TT; /* print(TT); */
				Cont = idiv(dp_hc(U),dp_hc(T));
				D0 = kmul(CR,D[0]);
				R0 = kmul(Cont,R[0]);
				S = igcd(D0,R0);
				D = [idiv(D0,S),D[1]];
				R = [idiv(R0,S),T];
				Mag = idiv(dp_mag(dp_hm(R[1]))*M_NM,M_DN);
			} else {
				D = [kmul(CR,D[0]),D[1]];
				R = [R[0],U];
			}
		} else {
			C = kmul(dp_hc(R[1]),R[0]);
			T = igcd(D[0],C);
			D = [T,idiv(D[0],T)*D[1]+idiv(C,T)*dp_ht(R[1])];
			R = [R[0],dp_rest(R[1])];
		}
	}
	return [D[1],reverse(H),T1,T2];
}

extern PTOZP,DPCV,NFSDIR;

def imulv(A,B)
{
	return A*B;
}

def dp_imul(P,N)
{
	return N*P;
}

def mul(A,B)
{
	return A*B;
}

def dp_imul_index(C,INDEX)
{
	T = C*bload(NFSDIR+rtostr(INDEX));
	return T;
}

def reg_dp_dtov(P)
{
	PTOZP = P;
	DPCV = dp_dtov(PTOZP);
}

def reg_dp_iqr(G)
{
	N = size(DPCV)[0];
	for ( R = [], I = 0; I < N; I++ ) {
		DPCV[I] = T = irem(DPCV[I],G);
		if ( T )
			R = cons(T,R);
	}
	return R;
}

def reg_dp_cont(P)
{
	PTOZP = P;	
	C = dp_dtov(P);
	return igcd(C);
}

def reg_dp_idiv(GCD)
{
	return dp_idiv(PTOZP,GCD);
}

def dp_cont_d(G,PROC,NPROC)
{
	C = dp_sep(G,NPROC);
	N = size(C)[0];
	for ( I = 0; I < N; I++ )
		rpc(PROC[I],"reg_dp_cont",C[I]);
	R = map(rpcrecv,PROC);
	GCD = igcd(R);
	return GCD;
}

/*
def iqrv(C,D)
{
	return map(iqr,C,D);
}
*/

#if 1
def dp_ptozp_d(G,PROC,NPROC)
{
	C = dp_dtov(G); L = size(C)[0];
	T0 = time()[3];
	D0 = igcd_estimate(C);
	TT = time()[3]-T0; TG1 += TT;
	CS = sepvect(C,NPROC+1);
	N = size(CS)[0]; N1 = N-1;
	QR = newvect(N); Q = newvect(L); R = newvect(L);
	T0 = time()[3];
	for ( I = 0; I < N1; I++ )
		rpc0(PROC[I],"iqrv",CS[I],D0);
	QR[N1] = iqrv(CS[N1],D0);
	for ( I = 0; I < N1; I++ )
		QR[I] = rpcrecv(PROC[I]);
	TT = time()[3]-T0; TD += TT;
	for ( I = J = 0; I < N; I++ ) {
		T = QR[I]; M = size(T)[0];
		for ( K = 0; K < M; K++, J++ ) {
			Q[J] = T[K][0]; R[J] = T[K][1];
		}
	}
	T0 = time()[3];
	D1 = igcd(R);
	GCD = igcd(D0,D1);
	A = idiv(D0,GCD);
	TT = time()[3]-T0; TG2 += TT;
	T0 = time()[3];
	for ( I = 0; I < L; I++ )
		Q[I] = kmul(A,Q[I])+idiv(R[I],GCD);
	TT = time()[3]-T0; TM += TT;
	print([TG1,TD,TG2,TM,dp_mag(D0*<<0>>),dp_mag(D1*<<0>>),dp_mag(GCD*<<0>>)],2);
	return dp_vtod(Q,G);
}
#endif

#if 0
def dp_ptozp_d(G,PROC,NPROC)
{
	C = dp_sep(G,NPROC+1);
	N = size(C)[0]; N1 = N-1;
	R = newvect(N);
	for ( I = 0; I < N1; I++ )
		rpc(PROC[I],"reg_igcdv",dp_dtov(C[I]));
	R[N1] = reg_igcdv(dp_dtov(C[N1]));
	for ( I = 0; I < N1; I++ )
		R[I] = rpcrecv(PROC[I]);
	GCD = igcd(R);
	for ( I = 0; I < N1; I++ )
		rpc(PROC[I],"idivv_final",GCD);
	S = dp_vtod(idivv_final(GCD),C[N1]);
	for ( I = 0; I < N1; I++ )
		S += dp_vtod(rpcrecv(PROC[I]),C[I]);
	return S;
}
#endif

#if 0
def dp_ptozp_d(G,PROC,NPROC)
{
	C = dp_sep(G,NPROC+1);
	N = size(C)[0]; N1 = N-1;
	R = newvect(N);
	for ( I = 0; I < N1; I++ )
		rpc(PROC[I],"reg_dp_cont",C[I]);
	R[N1] = igcd(dp_dtov(C[N1]));
	for ( I = 0; I < N1; I++ )
		R[I] = rpcrecv(PROC[I]);
	GCD = igcd(R);
	for ( I = 0; I < N1; I++ )
		rpc(PROC[I],"reg_dp_idiv",GCD);
	S = dp_idiv(C[N1],GCD);
	for ( I = 0; I < N1; I++ )
		S += rpcrecv(PROC[I]);
	return S;
}
#endif

#if 0
def dp_ptozp_d(G,PROC,NPROC)
{
	C = dp_sep(G,NPROC);
	N = size(C)[0]; T = newvect(N);
	for ( I = 0; I < N; I++ )
		T[I] = PROC[I];
	PROC = T;
	for ( I = 0; I < N; I++ )
		rpc(PROC[I],"reg_dp_dtov",C[I]);
	A = dp_dtov(G); A = isort(A); L = size(A)[0];
	map(rpcrecv,PROC);
	while ( 1 ) {
		GCD = igcd(A[0],A[1]);
		for ( I = 2; I < L; I++ ) {
			GT = igcd(GCD,A[I]);
			if ( GCD == GT )
				break;
			else
				GCD = GT;
		}
		if ( I == L )
			break;
		else {
			map(rpc,PROC,"reg_dp_iqr",GCD);
			R = map(rpcrecv,PROC);
			for ( J = 0, U = []; J < N; J++ )
				U = append(R[J],U);
			L = length(U);
			if ( L == 0 )
				break;
			else 
				U = cons(GCD,U);
			A = isort(newvect(L+1,U));
			print(["length=",L,"minmag=",dp_mag(A[0]*<<0>>)]);
		}
	}
	for ( I = 0; I < N; I++ )
		rpc(PROC[I],"reg_dp_idiv",GCD);
	R = map(rpcrecv,PROC);
	for ( I = 0, S = 0; I < N; I++ )
		S += R[I];
	return S;
}
#endif

def dp_ptozp2_d(D,G,PROC,NPROC)
{
	T = dp_ptozp_d(D+G,PROC,NPROC);
	for ( S = 0; D; D = dp_rest(D), T = dp_rest(T) )
		S += dp_hm(T);
	return [S,T];
}

def genindex(N)
{
	R = [];
	for ( I = 0; I < N; I++ )
		R = cons(I,R);
	return reverse(R);
}

def nftest_ext_mod(N1,N2,N,MOD,DIR,PSH)
{
	S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
	R = nf_mod_ext(genindex(N-1),S,MOD,DIR,PSH);
	return R;
}

def nftest_mod(N1,N2,N,MOD,DIR,PSH)
{
	S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
	R = nf_mod(genindex(N-1),S,MOD,DIR,PSH);
	return dp_mdtod(R);
}

def nftest(N1,N2,N,DIR,PSH)
{
	S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
	R = nf_demand(genindex(N-1),S,DIR,PSH);
	return R;
}

def nftest_d(N1,N2,N,DIR,NDIR,PDIR,PSH,PROC)
{
	S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
	R = nf_demand_d(genindex(N-1),S,DIR,NDIR,PDIR,PSH,PROC);
	return R;
}

def abs(A)
{
	return A >= 0 ? A : -A;
}

def sort(A)
{
	N = size(A)[0];
	for ( I = 0; I < N; I++ ) {
		for ( M = I, J = I + 1; J < N; J++ )
			if ( A[J] < A[M] )
				M = J;
		if ( I != M ) {
			T = A[I]; A[I] = A[M]; A[M] = T;
		}
	}
	return A;
}

def igcdv(C)
{
	C = sort(C); N = size(C)[0];
	G = igcd(C[0],C[1]);
	for ( I = 2; I < N; I++ ) {
		T = time();
		G = igcd(G,C[I]);
	/*	print([I,dp_mag(G*<<0>>),time()[0]-T[0]]); */
	}
	return G;
}

def igcdv2(C)
{
	C = sort(C); N = size(C)[0];
	G = igcd(C[0],C[1]);
	for ( I = 2; I < N; I++ ) {
		T = time();
		C[I] %= G;
		GT = igcd(G,C[I]);
		if ( G == GT ) {
			for ( J = I+1, NZ=0; J < N; J++ ) {
				C[J] %= G;
				if ( C[J] )
					NZ = 1;
			}
			if ( !NZ )
				break;
		} else
			G = GT;
		print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
	}
	return G;
}

def igcdv_d(C,P,NP)
{
	C = isort(C); N = size(C)[0];
	G = igcd(C[0],C[1]);
	for ( I = 2; I < N; I++ ) {
		T = time();
		C[I] %= G;
		GT = igcd(G,C[I]);
		if ( G == GT ) {
			for ( J = I+1, NZ=0; J < N; J++ ) {
				C[J] %= G;
				if ( C[J] )
					NZ = 1;
			}
			if ( !NZ )
				break;
		} else
			G = GT;
		print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
	}
	return G;
}

def dup(A)
{
	N = size(A)[0]; V = newvect(N);
	for ( I = 0; I < N; I++ )
		V[I] = A[I];
	return V;
}

def idivv_init(A)
{
	IDIVV_REM = dup(A);
}

def igcd_cofactor(D)
{
	T = map(iqr,IDIVV_REM,D);
	N = size(T)[0];	
	Q = newvect(N); R = newvect(N);
	for ( I = 0; I < N; I++ ) {
		Q[I] = T[I][0]; R[I] = T[I][1];
	}
	D1 = igcdv(dup(R));
	if ( !D1 ) {
		IDIVV_HIST = [[Q,D]];
		return D;
	} else {
		Q1 = map(idiv,R,D1);
		IDIVV_HIST = [[Q,D],[Q1,D1]];
		return D1;
	}
}

def idivv_step(D)
{
	T = map(iqr,IDIVV_REM,D);
	N = size(T)[0];	
	Q = newvect(N); R = newvect(N);
	for ( I = 0, NZ = 0; I < N; I++ ) {
		Q[I] = T[I][0]; R[I] = T[I][1];
		if ( R[I] )
			NZ = 1;
	}
	IDIVV_REM = R;
	IDIVV_HIST = cons([Q,D],IDIVV_HIST);
	return NZ;
}

def igcdv3(C)
{
	idivv_init(C);
	C = isort(C); N = size(C)[0];
	D = igcd(C[0],C[1]);
	for ( I = 2; I < N; I++ ) {
		T = time();
		C[I] %= G;
		GT = igcd(G,C[I]);
		if ( G == GT ) {
			NZ = idivv_step(G);
			if ( !NZ )
				break;
		} else
			G = GT;
	}
	CF = idivv_final(G);
	return [G,CF];
}

def igcdv4(C)
{
	idivv_init(C);
	C = isort(C); N = size(C)[0];
	D = igcd_estimate(C);
	G = igcd_cofactor(D);
	G = igcd(G,D);
	CF = idivv_final(G);	
	return [G,CF];
}

def igcd_estimate(C)
{
	C = isort(C); N = size(C)[0];
	M = idiv(N,2);
	for ( I = 0, S = 0; I < M; I++ )
		S += C[I]>0?C[I]:-C[I];
	for ( T = 0; I < N; I++ )
		T += C[I]>0?C[I]:-C[I];
	if ( !S && !T )
		return igcdv(C);
	else
		return igcd(S,T);
}

def reg_igcdv(C)
{
	D = igcd_estimate(C);
	T = map(iqr,C,D); N = size(T)[0];	
	Q = newvect(N); R = newvect(N);
	for ( I = 0; I < N; I++ ) {
		Q[I] = T[I][0]; R[I] = T[I][1];
	}
	D1 = igcdv(dup(R));
	if ( !D1 ) {
		IDIVV_HIST = [[Q,D]];
		return D;
	} else {
		Q1 = map(idiv,R,D1);
		IDIVV_HIST = [[Q,D],[Q1,D1]];
		return igcd(D,D1);
	}
}

def idivv_final(D)
{
	for ( T = IDIVV_HIST, S = 0; T != []; T = cdr(T) ) {
		H = car(T); S += map(kmul,H[0],idiv(H[1],D));
	}
	return S;
}

def dp_vtod(C,P)
{
	for ( R = 0, I = 0; P; P = dp_rest(P), I++ )
		R += C[I]*dp_ht(P);
	return R;
}

def dp_nt(P)
{
	for ( I = 0; P; P = dp_rest(P), I++ );
	return I;
}
end$