CheckT[f1_,f2_,f3_,f4_,y_,b_,tau_][t_]:=
Function[k,Inverse[tau].(y+{k,t}.b)]/@CheckT1[f1,f2,f3,f4,y,b,tau][t]

CheckT1[f1_,f2_,f3_,f4_,y_,b_,tau_][t_]:=
With[{ineq=CheckT2[f1,f2,f3,f4,y,b,tau][t]},
If[
ToString[ineq]=="False",{},
If[
ToString[ineq[[0]]]=="Or",
Function[i,ineq[[i,2]]]/@Range[Length[ineq]],
{ineq[[2]]}
]
]
]

CheckT2[f1_,f2_,f3_,f4_,y_,b_,tau_][t_]:=
Reduce[
(f1[Inverse[tau].(y+b[[1]]+b[[2]] t)]-f1[Inverse[tau].(y+b[[2]] t)])x+f1[Inverse[tau].(y+b[[2]] t)]≤0&&
(f2[Inverse[tau].(y+b[[1]]+b[[2]] t)]-f2[Inverse[tau].(y+b[[2]] t)])x+f2[Inverse[tau].(y+b[[2]] t)]≤0&&
(f3[Inverse[tau].(y+b[[1]]+b[[2]] t)]-f3[Inverse[tau].(y+b[[2]] t)])x+f3[Inverse[tau].(y+b[[2]] t)]≤0&&
(f4[Inverse[tau].(y+b[[1]]+b[[2]] t)]-f4[Inverse[tau].(y+b[[2]] t)])x+f4[Inverse[tau].(y+b[[2]] t)]≤0
(*f2[Inverse[tau].(y+b[[1]] x+b[[2]] t)]≤0&&
f3[Inverse[tau].(y+b[[1]] x+b[[2]] t)]≤0&&
f4[Inverse[tau].(y+b[[1]] x+b[[2]] t)]≤0*),Integers]

FindH[b_]:=Module[{c,bn,m},
c=Orthogonalize[b[[;;-2]]];
bn=b[[-1]];
m=Total[Function[v,Projection[bn,v]]/@c];
Norm[bn-m]
]

Approx2[pt_,b_]:=If[Length[b]==1,
Module[{bn,v,q,nb,nv},
bn=b[[-1]];
v=Projection[pt,bn];
q=pt-v;
nb=Norm[bn];
nv=((bn.v)/(nb^2));
Round[nv]],
With[{k=Round[m/.With[{p=pt,b1=b[[1]],b2=b[[2]]},
Solve[(m^2)b2.b2+p.p-2m p.b2-(((m^2)(((b1.b2))^2)-2m(b1.b2)(p.b1)+(((p.b1))^2))/(b1.b1))==0,m]]
[[1]]]},{Approx2[pt-k b[[2]],b[[;;-2]]],k}]
]

fun1[p_,a_,eps_,k_][{x_,y_}]:=x+(Sqrt[2])y-(a^2)(p^k)-a eps (p^k);
fun2[p_,a_,eps_,k_][{x_,y_}]:=-x-(Sqrt[2])y+(a^2)(p^k)-a eps (p^k);
fun3[p_,a_,eps_,k_][{x_,y_}]:=x-(Sqrt[2])y-(p^k);
fun4[p_,a_,eps_,k_][{x_,y_}]:=-x+(Sqrt[2])y-(p^k);

HTLenstra[p_,a_,eps_][k_]:=Module[{f1,f2,f3,f4,v0,v1,v2,v,b,bn,ob1,ob2,ob,x,h,tau,R,y},
f1=fun1[p,a,eps,k];
f2=fun2[p,a,eps,k];
f3=fun3[p,a,eps,k];
f4=fun4[p,a,eps,k];
v0=((p^k)/2){1+(a^2)+a eps,(1/(Sqrt[2]))(-1+(a^2)+a eps)};
v1=((p^k)/2){1+(a^2)-a eps,(1/(Sqrt[2]))(-1+(a^2)-a eps)};
v2=((p^k)/2){-1+(a^2)-a eps,(1/(Sqrt[2]))(1+(a^2)-a eps)};
v={v1,v2};
tau=Inverse[Transpose[{v1-v0,v2-v0}]];
x=tau.((v0+v1+v2)/3);
R=(Sqrt[2]);
ob1=tau.{1,0};
ob2=tau.{0,1};
ob={ob1,ob2};
b=SortBy[(1/(10^10))LatticeReduce[Floor[(10^10)ob]],Norm];
bn=b[[2]];
h=FindH[b];
y=Approx2[x,b].b;
Round/@Join@@(CheckT[f1,f2,f3,f4,y,b,tau]/@Range[-Floor[R/h],Floor[R/h]])
];HTLenstra[2,1/3,1/10^10][17]