# $Source: /u/maple/research/lib/factor/src/RCS/knapsack,v $
#
# Knapsack factorization algorithm. This algorithm is called
# from `factor/combfact` (= Zassenhaus method, searching all
# possible combinations of all modular factors) when it turns
# out that there are many modular factors and when combinations
# of 1, 2, or 3 modular factors don't give any true factors.
#
# Mark van Hoeij, July 2000.

macro(
	compute_answer=`factor/knapsack/comp_ans`,
	cut_away=`factor/knapsack/cut_away`,
	Tr_i=`factor/knapsack/Tr_i`
):

# Input:
# f: the polynomial to be factored
# x: variable
# p: prime number
# a: accuracy of modular factors (they're determined mod p^a)
# B: rootbound(f,x)
# aN: leading term of f
# mf: number of modular factors, and the modular factors themselves
# mf[0]=n mf[1..n] are the modular factors
#
# Ouput: a factorization of f into irreducible factors.
#
`factor/knapsack`:=proc(mf,Q,f,x,B,aN,n)
	local N,p,a,MF,lp,b,BB,k,total,i,max_i,tot_ratio,M,C,inp,L,
	ith,inc,nL,a_b,E,m,LL,ans,q,inC,old_a,f_s,not_rem,j,v;

	userinfo(2,factor,`knapsack factoring started at time`,time());

	N:=degree(f,x);
	q:=Q;
	p,a:=op(ifactors(q)[2][1]);
	# modular factors:
	MF:=[seq(modp1('ConvertOut'(mf[i],x),q),i=1..n)];

	# Calculate the bounds b[i] for Tr_i and calculate
	# the total of all a-b[i]
	Digits:=10;
	lp:=evalf(log(p));
	b[0]:=evalf(log(2*N)/lp);
	BB:=max(2,abs(aN*B));
	k:=evalf(log(BB)/lp);
	total:=0;
	for i to iquo(N,2) while b[i-1] < a-2 do
		b[i]:=ceil(b[0] + i*k);
		total:=total+max(0,a-b[i])
	od;
	max_i:=i-2;

	# This the ratio "amount of info we can input into lattice"
	# divided by "estimate of what we need"
	tot_ratio:=evalf(9*total*log(p)/n^2);
	# If tot_ratio=2 or so that should be enough but just to be
	# on the safe side:
	if tot_ratio < 4 or max_i < 4 then
		# Do additional lifting. Doesn't happen too often.
		if n<18 + 2*(nargs-6)^2 then
			# no more lifting, just return to Zassenhaus method.
			return FAIL
		fi;
		M[0]:=n;
		for i to n do
			M[i]:=modp1('ConvertIn'(MF[i] mod p,x),p)
		od;
		# double a, which increases tot_ratio
		`factor/lift`(f,M,x,p,'i','j','k',`lift to`,q);
		# and try again
		return procname(M,q^2,args[3..-1],x)
	fi;

	# Multiply the 0-1 vectors by this C
	C:=max(1,round(2*sqrt(n)));
	# Changing this number affects the lengths of the input
	# vectors.
	inp:=0.11;

	# L is the "search space".
	# The 0-1 vector of a factor g in Z[x] of f is a vector with
	# 0's and 1's, with a 1 at entry i precisely when modular
	# factor mf[i] is used in g.
	# The search space is a vector space that contains all of
	# these 0-1 vectors. The following L is a basis of that
	# search space. We'll try to reduce the number of elements
	# of L by lattice reductions. Then a Gram-Schmidt computation
	# is used to decide which vectors can be thrown away.
	L:=[seq([seq(`if`(i=j,C,0),j=1..n)],i=1..n)];

	# Some initialization.
	ith:=0;
	nL:=n;

	# Initially use "inc" rows at a time (could increase) when
	# it turns out to be necessary. The lengths of these vectors
	# is such that if we used "Inc" of them we would have
	# determinant around C^(..) * exp(inp*n^2), but since we
	# start with fewer Tr_i's than Inc, we'll start with a
	# smaller determinant. Larger Inc -> smaller vectors.
	if n<45 then
		inc:=2;
		inC:=2.5
	elif n<65 then
		inc:=2;
		inC:=3
	elif n<100 then
		inc:=2;
		inC:=4
	elif n<150 then
		inc:=3;
		inC:=8
	else
		inc:=4;
		inC:=12
	fi;
do
	a_b:=ceil((2*log(nops(L[1]))+inp/max(inC,inc)*max(nL,n/2)^2)/lp);
	if inp=0 then a_b:=a-b[ith+inc] fi;
	if not (ith+inc <= max_i and a-b[ith+inc] >= a_b) then
		inc:=inc+1;
		if ith+inc > max_i then
			inc:=1;
			ith:=0;
			if inp=0 then
				# give up
				break
			fi;
			inp:=0
		fi;
		next
	elif a > 4*a_b + b[inc+ith] and not assigned(old_a) then
		old_a:=a;
		a:=4*a_b + b[inc+ith];
		q:=p^a
	elif a <= 2*a_b + b[inc+ith] and assigned(old_a) then
		a:=min(old_a,a+2*a_b+1);
		q:=p^a
	fi;
	E:=p^(a_b);  # determines the size of the input vectors.
	m:=p^(a-a_b);
	nL:=nops(L);
	# Put in "inc" new vectors, using Tr_i for i from ith+1..ith+inc
	LL:=[seq(subsop(  nops(L[1])+i  =E,[0$(nops(L[1])+inc)]),i=1..inc),
	 seq([op(v),seq(mods( add(v[j]/C*round(
	aN^i*Tr_i(MF[j],x,i,q)
	 /m),j=1..n) ,E),i=ith+1..ith+inc)],v=L)];

	# Bound for the size of the 0-1 vector, which has been multiplied
	# by C and some extra columns have been added:
	M:=C^2*n + nops(L[1]-n)*(n/2)^2;
	userinfo(3,factor,`reducing lattice for traces numbers`
	 ,seq(i,i=ith+1..ith+inc));
	L:=lattice(LL,integer);

	# Found a Short vector?
	f_s:=evalb(min(seq(add(j^2,j=i),i=L)) <= 2*M);
	if f_s then
		# Apply Gram-Schmidt to figure out how many rows in L can
		# be thrown away, and determine how many "long" rows could
		# not be thrown away
		L, not_rem:=cut_away(L,n,M);
	    if nops(L) <= n/4 then
		# rational factors consisting of 3 modular factors are computed
		# by Zassenhaus, so we can have at most n/4 rational factors.
		if nops(L)=1 then
			userinfo(2,factor,`polynomial found irreducible`);
			ans:=f
		else
			ans:=compute_answer(L,n,f,x,mf,Q,aN)
		fi;
		if ans<>FAIL then
			# Clear remember tables to free some memory
			forget(Tr_i);
			return ans
		fi
	    fi
	fi;
	if (not f_s or not_rem>0) and 2*a_b < a-b[ith+inc] then
	    E,m:=p^(a_b), p^(a-2*a_b);
	    # Put in more info from the same Tr_i's because there are
	    # still long vectors remaining
	    LL:=[seq([op(v),seq(mods(
		add(v[j]/C*round(aN^i*Tr_i(MF[j],x,i,q)/m),j=1..n)
	     ,E),i=ith+1..ith+inc)],v=L),
	     seq(subsop(  nops(L[1])+i  =E,[0$(nops(L[1])+inc)]),i=1..inc)];
	    M:=C^2*n + nops(L[1]-n)*(n/2)^2;
	    if max(seq(add(j^2,j=LL[i][-inc..-1]),i=1..nops(L)-inc)) > 2*M then
		# The extra info was not superfluous, so do something with it:
		userinfo(3,factor,`reducing lattice for traces numbers`
		 ,seq(i,i=ith+1..ith+inc));
		L:=lattice(LL,integer);
		L, not_rem:=cut_away(L,n,M)
	    fi
	fi;
	if nops(L) > nL then
		# Our reductions were no good, so we have to use longer
		# vectors.
		inp:=inp*1.2
	fi;
	ith:=ith+inc;
od;
if n>18 + 2*(nargs-6)^2 then
	M:='M';
	M[0]:=n;
	for i to n do
		M[i]:=modp1('ConvertIn'(MF[i] mod p,x),p)
	od;
	# Do more Hensel lifting:
	`factor/lift`(f,M,x,p,'i','j','k',`lift to`,Q);
	# and try again:
	procname(M,Q^2,args[3..-1],x)
else
	# We could implement another strategy because it shouldn't be
	# hopeless given that tot_ratio is definitely high enough, but
	userinfo(2,factor,`giving up, back to Zassenhaus`);
	FAIL
fi
end:


compute_answer:=proc(L,n,f,x,mf,q,aN)
	local FA,S,F,R,i,t,j,Q,v;
	# R:=linalg['rref']([seq(v[1..n],v=L)]);
	# R:=[seq(convert(linalg['col'](R,i),list),i=1..n)];
	R:=LinearAlgebra:-RowSpace(convert([seq(v[1..n],v=L)],Matrix));
	R:=LinearAlgebra:-Transpose(Matrix(map(convert,R,list)));
	R:=[seq(convert(LinearAlgebra:-Row(R,i),list),i=1..n)];
	FA:={op(R)};
	if nops(FA)<>nops(L) or {op(map(op,FA))} <> {0,1} or
	  map(convert,FA,`+`)<>{1} then
		# Matrix R was not a matrix with precisely one 1
		# in each column, rest 0.
		return FAIL
	fi;
	S:=sort([seq({seq(`if`(R[i]=j,i,NULL),i=1..n)},j=FA)],
	 (a,b) -> evalb( nops(a)<nops(b) ) );
	F:=f;
	R:=1;
	for i in S do
		t:=modp1('Constant'(aN,x),q);
		for j in i do
			t := modp1('Multiply'(t,mf[j]),q)
		od;
		t := mods(modp1('ConvertOut'(t,x),q),q);
		divide(t,icontent(t),'t');
		if not divide(F,t,'Q') then
			return FAIL
		fi;
		R:=R*t;
		F:=eval(Q)
	od;
	if has(F,x) then
		error "should not happen"
	fi;
	F*R
end:

# Throw away everything from the lattice that can't contribute
# to the space W of M-short vectors
cut_away:=proc(L,n,M)
	local R,d,i,j;
	userinfo(3,factor,`Computing Gram-Schmidt`);
	UseHardwareFloats:=false;
	Digits:=ceil(nops(L[1])*log(2)/log(10)+max(seq(op(map(length,i)),i=L)));
	R:=LinearAlgebra:-GramSchmidt(map(convert,evalf(L),Vector)):
	d:=[seq(LinearAlgebra:-DotProduct(i,i),i=R)]:
	i:=nops(d);
	# Take 1.2*M, gives some safety margin
	while d[i] > 1.2*M do
		if i=1 then
			# There must be M-short vectors
			error "this should not happen"
		else
			i:=i-1
		fi
	od;
	userinfo(3,factor,i,`of the`,nops(L),
	 `vectors remaining after GramSchmidt`);
	# Return these i vectors and determine how many of them are "long"
	L[1..i], add(`if`(d[j] > 1.2*M,1,0),j=1..i)
end:


# Compute the i'th trace of f mod m, so compute:
#       evala(Trace( RootOf(f,x)^i )) mod m
# We could of course just issue that command but we want
# to reduce mod m in intermediate computations for
# efficiency reasons
Tr_i:=proc(f,x::name,i::nonnegint,m::posint)
	local d,k;
	options remember;
	d:=degree(f,x);
	if lcoeff(f,x)<>1 then
		error "f should be monic"
	elif i=0 then
		d
	else
		# Newton identities
		-add(procname(f,x,k,m)*coeff(f,x,d-i+k),k=1..i-1)
		-i*coeff(f,x,d-i) mod m
	fi
end:

#savelib('`factor/knapsack`', 'cut_away','compute_answer','Tr_i'):