// LValue : computation of L( M, n )
// with M a t-motif and n an integer > max w(M)


// addition chains precomputed for exponents up to 40

adch := 
[
	[],
	[ 1 ],
	[ 1, 1 ],
	[ 1, 2 ],
	[ 1, 1, 2 ],
	[ 1, 1, 3 ],
	[ 1, 1, 1, 3 ],
	[ 1, 2, 3 ],
	[ 1, 1, 3, 3 ],
	[ 1, 1, 2, 4 ],
	[ 1, 1, 1, 3, 4 ],
	[ 1, 1, 3, 4 ],
	[ 1, 1, 2, 3, 4 ],
	[ 1, 1, 2, 2, 5 ],
	[ 1, 1, 2, 4, 4 ],
	[ 1, 2, 3, 4 ],
	[ 1, 2, 3, 4, 1 ],
	[ 1, 1, 3, 3, 5 ],
	[ 1, 1, 1, 4, 5, 3 ],
	[ 1, 2, 3, 2, 5 ],
	[ 1, 1, 1, 3, 5, 5 ],
	[ 1, 1, 1, 3, 4, 6 ],
	[ 1, 2, 1, 3, 4, 5 ],
	[ 1, 2, 3, 4, 4 ],
	[ 1, 1, 2, 4, 4, 5 ],
	[ 1, 1, 2, 3, 4, 6 ],
	[ 1, 1, 3, 3, 5, 5 ],
	[ 1, 1, 2, 2, 5, 6 ],
	[ 1, 1, 1, 3, 4, 5, 6 ],
	[ 1, 1, 3, 4, 5, 4 ],
	[ 1, 1, 2, 2, 4, 5, 6 ],
	[ 1, 2, 3, 4, 5 ],
	[ 1, 2, 3, 4, 5, 1 ],
	[ 1, 2, 3, 4, 5, 2 ],
	[ 1, 1, 2, 4, 4, 6, 4 ],
	[ 1, 1, 3, 3, 5, 6 ],
	[ 1, 1, 2, 3, 5, 4, 6 ],
	[ 1, 1, 2, 2, 4, 5, 7 ],
	[ 1, 1, 2, 4, 3, 6, 6 ],
	[ 1, 2, 2, 3, 5, 6 ]
];


// FrobeniusTwist( P, q, n ):
//	 takes a square matrix P of polynomials over a finite field
//	 returns the q^n Frobenius twist of P

function FrobeniusTwist( P, q, n )
	if n eq 0 then return P; end if;
	r := NumberOfRows( P );
	qn := q^n;
	for i in [1..r] do for j in [1..r] do
		cf := Coefficients( P[i,j] );
		P[i,j] := [ c^qn : c in cf ]; // TIME CRITICAL
	end for; end for;
	return P;
end function;


// FrobeniusProduct( P, q, n ):
//   takes a square matrix P over a finite field
//   returns FrobeniusTwist( P, q, n-1 ) * ... FrobeniusTwist( P, q, 0 )
// computed using minimal addition chains

function FrobeniusProduct( P, q, n ) 
	ac := adch[ n ];
	l := # ac + 1;
	Q := []; d := [];
	Q[1] := P; d[1] := 1;
	for i in [ 2 .. l ] do
		a := ac[i-1];
		L := FrobeniusTwist( Q[a], q, d[i-1] );
		Q[i] := L * Q[i-1];  // TIME CRITICAL
		d[i] := d[a] + d[i-1];
	end for;
	return Q[l];
end function;


// PlacesOfBadReduction( M, d )
//   returns the ``minimal primitive'' elements x of GF(q^d) 
//   for which th=x defines a pole of some coefficient of M

function PlacesOfBadReduction( M, d )
	Kt := BaseRing( Parent( M ) );
	k := BaseRing( BaseRing( Kt ) );
	q := # k;
	r := NumberOfRows( M );
	B := { GF( q^d ) | };
	
	// produces big set B of all bad values of th
	for i in [1..r] do for j in [1..r] do
		for cf in Coefficients( M[i][j] ) do
			poles := Roots( Denominator( cf ), GF( q^d ) );
			poles := { x[1] : x in poles };
			B := B join poles;
		end for;
	end for; end for;

	// filter out the nonminimal and imprimitive ones
	Bad := { GF( q^d ) | };
	for g in B do
		a := g;
		for i in [2..d] do
			a := a^q;
			if g ge a then continue g; end if; 
		end for;
		Bad := Bad join { g };
	end for;
	return Bad;
end function;


// ErrorEstimate( M, n, deg : SampleSize := 100, UpperBound := 1000 )
//   returns e such that 
//   L ( M, n ) = L( M, n, deg ) + O( z^e ) 
//   (statistically, no proof)

function ErrorEstimate( M, n, deg : SampleSize := 100, UpperBound := 1000 )
	Kt := BaseRing( Parent( M ) );
	k := BaseRing( BaseRing( Kt ) );
	q := # k;
	r := NumberOfRows( M );

	RR< z > := LaurentSeriesRing( k, UpperBound );
	tt := z^-1;

	errorest := UpperBound;
	for d in [ deg + 1 .. deg + 3 ] do
		kd := GF( q^d );
		kdx< x > := PolynomialRing( kd );
		P := ZeroMatrix( kdx, r, r );
		cf := { RR | };
		valuations := { UpperBound };
		Bad := PlacesOfBadReduction( M, d );
		for sample in [1..SampleSize] do
			gen := Random( kd );
			a := gen;
			for j in [ 2..d ] do
				a := a^q;
				if gen eq a then continue sample; end if;
			end for;
			if gen in Bad then continue sample; end if;
			f := MinimalPolynomial( gen, k );

			// build matrix of frobenius at f
			for i in [1..r] do for j in [1..r] do
				c := Coefficients( M[i,j] );
				e := [ kd!Evaluate( y, gen ) : y in c ];
				P[i,j] := kdx ! e;
			end for; end for;
			P := FrobeniusProduct( P, q, d );

			// char pol of frob at f
			cpt := Kt ! CharacteristicPolynomial( P ); 

			// compute value of inverse Euler factor at f in n
			ff := Evaluate( f, tt )^(-n);
			Lf := 0;
			for i in [0..r] do
				b := Evaluate( Coefficient( cpt, r-i ), tt );
				Lf := Lf + ff^i * b;
			end for;

			if not IsWeaklyZero( Lf - 1 ) then
				val := Valuation( Lf - 1 );
				valuations := valuations join { val };
				lead := LeadingCoefficient( Lf - 1 );
				cf := cf join { lead * z^val };
			end if;
		end for;

		minimalvaluation := Minimum( valuations );

		alleq := #cf eq 1;
		tot := NumberOfPrimePolynomials( q, d );
		divbyq := IsDivisibleBy( tot, q );
		if alleq and divbyq then
			minimalvaluation := minimalvaluation + 1;
		end if;

		errorest := Minimum( errorest, minimalvaluation );

	end for;
	return errorest;
end function;

				
// PartialLValue( M, n, d : Precision := 20 )
//   M a t-motive, n an integer > maximal weight of M, d positive integer
//   returns product of Euler factors of primes of deg = d
//   occuring in L( M, n )

function PartialLValue( M, n, d : Precision := 20 )
	Kt := BaseRing( Parent( M ) );
	k := BaseRing( BaseRing( Kt ) );
	q := # k;
	r := NumberOfRows( M );

	RR< z > := LaurentSeriesRing( k, Precision );
	tt := z^-1;
	L := RR!1;

	kd := GF( q^d );
	kdx< x > := PolynomialRing( kd );
	P := ZeroMatrix( kdx, r, r );

	Bad := PlacesOfBadReduction( M, d );
	
	for gen in kd do

		// test if gen is minimal root of some irr deg d pol f
		a := gen; 
		for i in [2..d] do
			a := a^q;
			if gen ge a then continue gen; end if; // T.C.
		end for; 
		if gen in Bad then continue gen; end if;
		f := MinimalPolynomial( gen, k );	
		
		// build matrix of frobenius at f
		for i in [1..r] do for j in [1..r] do
			c := Coefficients( M[i,j] );
			e := [ kd!Evaluate( y, gen ) : y in c ]; // T.C.
			P[i,j] := kdx ! e;
		end for; end for; 
		P := FrobeniusProduct( P, q, d ); 

		// characteristic polynomial of frobenius at f
		cpt := Kt ! CharacteristicPolynomial( P ); 

		// compute value of inverse Euler factor at f in n
		ff := Evaluate( f, tt )^(-n); 
		Lf := 0;
		for i in [0..r] do
			b := Evaluate( Coefficient( cpt, r-i ), tt );
			Lf := Lf + ff^i * b;
		end for;

		// multiply inverse Euler factor onto the big L
		L := L * Lf;

	end for;

	// invert L to get the final result
	return L^-1;
end function;


// LValue( M, n, d : Precision := 20 )
//  returns approximation to L( M, n ) computed
//  by multiplying the Euler factors at all
//  primes of degree at most d

function LValue( M, n, d : Precision := 30 )
	k := BaseRing( BaseRing( BaseRing( Parent( M ) ) ) );
	RR< z > := PowerSeriesRing( k, Precision );
	L := RR!1;
	for i in [1..d] do
		Li := PartialLValue( M, n, i : Precision := Precision );
		L := L * RR ! Li;
	end for;
	return L;
end function;


// main benchmark example
function benchmark()
	K< th > := RationalFunctionField( GF( 2 ) );
	Kt< t > := PolynomialRing( K );
	M := Matrix( Kt, 2, 2, [ 1, th-t, 1, 0 ] );

	t := 0;
	n := 0;
	while true do
		n := n + 1;
		ct := Cputime();
		PartialLValue( M, 1, 21 : Precision := 20 );
		tn :=  Cputime(ct); t := t + tn;
		print tn, t/n, GetMemoryUsage();
	end while;
	return 0;
end function;