/* $OpenXM: OpenXM/src/ox_ntl/ntl.rr,v 1.4 2008/09/19 10:55:40 iwane Exp $ */

module ntl;
localf factor$
localf lll$
localf ex_data$
localf ex_data_tmp$
localf mat2list$
localf list2mat$


/* static variables */

/* extern variables */

#if 1
localf test_factor$
localf test_lll$

/*&usage begin: ntl.test_factor(PID, POLY)
 compare on ox_NTL and on asir.

example: 
[1027] PID=ox_launch(ox_ntl);
[1028] F=ntl.ex_data(4);
x^16-136*x^14+6476*x^12-141912*x^10+1513334*x^8-7453176*x^6+13950764*x^4-5596840*x^2+46225
[1029] F = F * subst(F, x, x + 1)$
[1030] ntl.factor(PID, F);
[[1,1],[x^16+16*x^15-16*x^14-1344*x^13-4080*x^12+32576*x^11+157376*x^10-255232*x^9-2062624*x^8-249088*x^7+10702080*x^6+9126912*x^5-18643712*x^4-24167424*x^3+2712576*x^2+10653696*x+2324736,1],[x^16-136*x^14+6476*x^12-141912*x^10+1513334*x^8-7453176*x^6+13950764*x^4-5596840*x^2+46225,1]]
[1031] ntl.test(PID, F);
[CPU,0.121539,0.001354,GC,0.0222,0]

end: */
def test_factor(PID, F)
{
	T0 = time();
	fctr(F);
	T1 = time();
	ntl.factor(PID, F);
	T2 = time();

	return (["CPU", T1[0]-T0[0],T2[0]-T1[0],"GC",T1[1]-T0[1],T2[1]-T1[1]]);
}

def test_lll(PID, F)
{
	T0 = time();
	pari(lllint, F);
	T1 = time();
	ntl.lll(PID, F);
	T2 = time();

	return (["CPU", T1[0]-T0[0],T2[0]-T1[0],"GC",T1[1]-T0[1],T2[1]-T1[1]]);
}

#endif

/*&usage begin: ntl.factor(PID, POLY)

  Factorize polynomial {POLY} over the rationals.

{RETURN}
	list
{POLY}
	univariate polynomial with rational coefficients 

description: 
  Factorizes polynomial {POLY} over the rationals.
  The result is represented by a list, whose elements are a pair represented as

[[num,1],[factor1,multiplicity1],[factor2,multiplicity2],...].

Products of all factor^multiplicity and num is equal to {POLY}.

The number {num} is determined so that ({POLY}/{num}) is an integral polynomial
and its content (GCD of all coefficients) is 1.

author: H.Iwane <iwane@math.sci.kobe-u.ac.jp> 

example: 
[1282] F=(x^5-1)*(x^3+1)^2;
x^11+2*x^8-x^6+x^5-2*x^3-1
[1283] ntl.factor(PID, F);
[[1,1],[x^4+x^3+x^2+x+1,1],[x-1,1],[x+1,2],[x^2-x+1,2]]
[1284] fctr(F);
[[1,1],[x-1,1],[x^4+x^3+x^2+x+1,1],[x+1,2],[x^2-x+1,2]]

ref: fctr
end:
*/
def factor(PID, POLY)
{
	local F, C, LIST, STR, RET, RET_NTL, VAR, I;
	local TYPE;

	/* 入力チェック */
	TYPE = type(POLY);
	if (TYPE == 0 || TYPE == 1) {
		return ([[POLY,1]]);
	} else if (TYPE != 2) {
		error("ntl.factor: invalid argument");
	}


	LIST = vars(POLY);
	if (length(LIST) >= 2) { /* 一変数多項式のみ */
		error("ntl.factor: invalid argument");
	}


	/* NTL で 有理係数多項式は不可 */
	F = ptozp(POLY);

	C = sdiv(POLY, F);

	ox_cmo_rpc(PID, "fctr", F);

	RET_NTL = ox_pop_cmo(PID);	

	/* ERROR Check */
	if (type(RET_NTL) != 4 || length(RET_NTL) < 2) {
		return (RET_NTL);
	}

	RET = cons([RET_NTL[0] * C, 1], RET_NTL[1]);

	return (RET);
}

/*&usage begin: ex_data_tmp(F, N)
  Generate sample irreducible polynomial

example:
[1032] F=t^3-p;       
-p+t^3
[1033] ntl.ex_data_tmp(F,3);
x^27-90*x^24+1089*x^21-62130*x^18+105507*x^15-16537410*x^12-30081453*x^9-1886601330*x^6+73062900*x^3-6859000
[1034] ntl.ex_data_tmp(F,4);
-x^81+459*x^78-76896*x^75+7538094*x^72-347721147*x^69+3240161703*x^66+1032617170332*x^63-37499673798288*x^60+784360767442050*x^57+150576308695750650*x^54+771023617441694964*x^51+67248913649472410184*x^48+13913995714637027898294*x^45+270221527865987051874714*x^42-8828542741395296724347412*x^39-154971101776040822743879716*x^36+4343529580943017469231383983*x^33+4648027555241630173815780123*x^30-1072436585643253024332438894564*x^27+16237394255218510503554781142602*x^24-134104542851048701593527668875195*x^21+727430949790032393675174790142991*x^18-2727255031466780569416130788693624*x^15+7102683996190585423335589883738868*x^12-12524463688445776069953452180105904*x^9+14119870779369458313271232460210576*x^6-8591747198480108022372636571451136*x^3+2346749360904699138972190279765184

ref: ex_data

end: */
def ex_data_tmp(F, N)
{
	FP = subst(F, p, prime(0));
	GP = subst(FP, t, x);

	for (I = 1; I < N; I++) {
		FP = subst(F, p, prime(I));
		GP = res(t, subst(GP, x, x-t), FP);
	}

	return (GP);
}

/*&usage begin: ntl.ex_data(N)
  Generate sample irreducible polynomial

example: 
[1028] ntl.ex_data(3);
x^8-40*x^6+352*x^4-960*x^2+576
[1029] ntl.ex_data(4);
x^16-136*x^14+6476*x^12-141912*x^10+1513334*x^8-7453176*x^6+13950764*x^4-5596840*x^2+46225

ref:
ex_data_tmp

end: */
def ex_data(N)
{
	if (type(N) != 1 && type(N) != 10) {
		print("invalid argument");
		return (0);
	}

	return (ex_data_tmp(t^2-p, N));
}


def mat2list(M)
{
	A = size(M);

	ROW=A[0];
	COL=A[1];

	for (I = 0; I < ROW; I++) {
		for (J = 0; J < COL; J++) {
			A = append(A, [M[I][J]]);
		}
	}

	return (A);
}

def list2mat(L)
{
	if (type(L) != 4) {
		return ("Invalid Argument");
	}

	ROW = L[0];
	if (type(ROW) == 10)
		ROW = int32ton(ROW);
	COL = L[1];
	if (type(COL) == 10)
		COL = int32ton(COL);

	A = newmat(ROW, COL);

	C = 2;
	for (I = 0; I < ROW; I++) {
		for (J = 0; J < COL; J++) {
			A[I][J] = L[C];
			C++;
		}
	}

	
	return (A);
}



/*&usage begin: ntl.lll(PID, MAT)
  the basics of LLL reducation.

{M}
	Matrix which element is Integer.

example: 
[1046] def trans(M) {
  RET = newmat(size(M)[1], size(M)[0]);
    for (I = 0; I < size(M)[0]; I++)
      for (J = 0; J < size(M)[1]; J++)
        RET[J][I] = M[I][J];
  return (RET);
}
[1047] def lllpari(A) {
  return (trans(trans(A) * pari(lllint, trans(A))));
}
[1048] M=newmat(3, 3, [[10,0, 10], [-7,3, 30], [7, 3, 20]]);
[ 10 0 10 ]
[ -7 3 30 ]
[ 7 3 20 ]
[1049] ntl.lll(PID, M);
[ 2 -6 0 ]
[ -1 -3 10 ]
[ 11 3 0 ]
[1050] lllpari(M);
[ 2 -6 0 ]
[ -1 -3 10 ]
[ 11 3 0 ]


ref:
pari(lll)

end: */
def lll(PID, M)
{
	/* parameter check */
	TYPE = type(M);
	if (TYPE != 6) {	/* matrix */
		MES = "ntl.lll: invalid argument: " + TYPE;
		error(MES);
	}

	A = mat2list(M);

	ox_cmo_rpc(PID, "lll", A);

	RET_NTL = ox_pop_cmo(PID);	

	/* return value check */
	if (type(RET_NTL) != 4) { /* list */
		error("ntl.lll: error");
	}
	
	R = list2mat(RET_NTL);

	return (R);
}


endmodule;


end$