_device_ _host_ SU3<floatT> BresenhamLoopXYU(int &x,int &y, int &z,int &t, gSite site, int &mu, int &nu, int &rho, int &t_dir)

{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
int cmax,cmin,cmid,cmax_dir,cmin_dir, cmid_dir;
if (x>=y & x>=z & y>=z) {cmax=x; cmid=y; cmin=z; cmax_dir=mu; cmid_dir=nu; cmin_dir=rho;}
else if (x>=y & x>=z & z>=y) {cmax=x; cmid=z; cmin=y; cmax_dir=mu; cmid_dir=rho; cmin_dir=nu;}
else if (y>=x & y>=z & x>=z) {cmax=y; cmid=x; cmin=z; cmax_dir=nu; cmid_dir=mu; cmin_dir=rho;}
else if (y>=x & y>=z & z>=x) {cmax=y; cmid=z; cmin=x; cmax_dir=nu; cmid_dir=rho; cmin_dir=mu;}
else if (z>=x & z>=y & x>=y) {cmax=z; cmid=x; cmin=y; cmax_dir=rho; cmid_dir=mu; cmin_dir=nu;}
else if (z>=x & z>=y & y>=x) {cmax=z; cmid=y; cmin=x; cmax_dir=rho; cmid_dir=nu; cmin_dir=mu;}
typedef GIndexer<All, HaloDepth> GInd;
int cmax2=2*cmax; //extension of 2d bresenham algorithm https://arxiv.org/pdf/hep-lat/0005018
int cmin2=2*cmin;
int cmid2=2*cmid;
int chi_min=cmin2-cmax; 
int chi_mid=cmid2-cmax; 
for(int it=0;it<t;it++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,t_dir)); //step in time dir
site=GInd::site_up(site,t_dir);
}
for (int ix=0;ix<cmax;ix++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmax_dir));//step in cmax
site=GInd::site_up(site,cmax_dir);

if(chi_mid>=0)
{
chi_mid-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmid_dir)); //step in cmid
site=GInd::site_up(site,cmid_dir);
}
chi_mid+=cmid2;

if(chi_min>=0)
{
chi_min-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmin_dir)); //step in cmin
site=GInd::site_up(site,cmin_dir);
}
chi_min+=cmin2;

}

return temp; 
}