# v1.3	pass in u,n (now local in findrel_v2.1)
# v1.2	assignment for fix when in `if 2*ordg1 >= ordg0 + ordg2' now adds `or ordg1 = 0'
# v1.1	_c used as constant (can be used within a proc but not assigned outside)

with(LinearAlgebra, CharacteristicPolynomial):
with(padic, ordp):
read "/Users/heba/Desktop/Implementations/Giles_Levy_implementation/code/multrecs_v1.3.txt":

ReplaceConstant := proc(f,n,p)
	local A, constant;
	A := Normal(f) mod p;
	A := Normal( (Rem(numer(A),n^p-n-constant,n) mod p)/(Rem(denom(A),n^p-n-constant,n) mod p) ) mod p;
	if has(A,n) then
		error "n should have disappeared"
	fi;
	subs(constant=_c, A)
end:


curvp := proc(f, p, u, n)       # input a rec eqn in u, p
local g, g0, g1, g2, ordg0, ordg1, ordg2, fix, i, k, m, mm, CP, pcurv, a0, a1, invar;
# option trace;
g := primpart(f, {u(n), u(n+1), u(n+2)});
g0, g1, g2 := seq(coeff(g, u(n+i)), i=0..2);
ordg0, ordg1, ordg2 := ordp(content(g0, n), p), ordp(content(g1, n), p), ordp(content(g2,n), p);
if (2*ordg1 >= ordg0 + ordg2 or ordg1=0) and type(ordg2 - ordg0, even) then
	fix := multrecs([seq(coeff(g, u(n+1-k))/coeff(g, u(n+2)), k=0..1)], [-p^((ordg2-ordg0)/2)], u, n);
else
	return []
fi;

try
	g := fix/coeff(fix, u(n+2)) mod p;
catch:
	return []
end try;

g0, g1 := seq(coeff(g, u(n+i)), i=0..1);

m[0] := Matrix([[0,1],[-g0, -g1]]);
seq( assign(m[i], eval(m[0], n=n+i)), i=1..p-1);
mm := m[1] . m[0];
for i from 2 to p-1 do mm := m[i] . mm od;

try    
	CP := CharacteristicPolynomial(mm, lambda) mod p;
catch:
	return []
end try;

a0, a1 := seq(ReplaceConstant(coeff(CP, lambda, i),n,p), i=0..1);

try
	if a0 <> 0 then
		g := Normal(a1^2/a0) mod p;
		[(Expand(numer(g)) mod p) / (Expand(denom(g)) mod p)]
	else
		[]
	fi;
catch:
	return []
end try;

end: