cmax2 := 2*cmax
cmin2 := 2*cmin
chi := cmin2 -cmax
FOR i := 1 TO cmax DO
step in max-direction
IF chi ≥ 0 THEN
chi := chi - cmax2
step in min-direction
ENDIF
chi := chi + cmin2
ENDDO
The generalisation to three dimensions is achieved by
combining two of these 2d-algorithms with different χ’s
for max-mid and max-min in just one loop over the maxdirection

// Initialize variables
cmax := maximum of (dx, dy, dz)
cmed := median of (dx, dy, dz)
cmin := minimum of (dx, dy, dz)

cmax2 := 2 * cmax
cmed2 := 2 * cmed
cmin2 := 2 * cmin

chi_mid := cmed2 - cmax // Decision variable for mid-direction
chi_min := cmin2 - cmax // Decision variable for min-direction

FOR i := 1 TO cmax DO

step in max-direction

// Handle mid-direction step
IF chi_mid ≥ 0 THEN
chi_mid := chi_mid - cmax2
step in mid-direction
ENDIF
chi_mid := chi_mid + cmed2

// Handle min-direction step
IF chi_min ≥ 0 THEN
chi_min := chi_min - cmax2
step in min-direction
ENDIF
chi_min := chi_min + cmin2

ENDDO

2d bresenham:
_device_ _host_ SU3<floatT> BresenhamLoopXYL(int &x,int &y, int &z, int &t, gSite site, int &mu, int &nu, int &rho)

{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
int cmax,cmin,cmax_dir,cmin_dir;
if (x>=y) {cmax=x; cmin=y; cmax_dir=mu; cmin_dir=nu;}
else {cmax=y; cmin=x; cmax_dir=nu; cmin_dir=mu;}
typedef GIndexer<All, HaloDepth> GInd;
int cmax2=2*cmax;
int cmin2=2*cmin;
int chi=cmin2-cmax; 
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>=0)
{
chi-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmin_dir)); //step in cmin
site=GInd::site_up(site,cmin_dir);
}
chi+=cmin2;
}
for(int iz=0;iz<z;iz++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //step in cmin
site=GInd::site_up(site,3-mu-nu);
}
return temp; 
}

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

{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
int cmax,cmin,cmax_dir,cmin_dir;
if (x>=y) {cmax=x; cmin=y; cmax_dir=mu; cmin_dir=nu;}
else {cmax=y; cmin=x; cmax_dir=nu; cmin_dir=mu;}
typedef GIndexer<All, HaloDepth> GInd;
int cmax2=2*cmax;
int cmin2=2*cmin;
int chi=cmin2-cmax; 
for(int iz=0;iz<z;iz++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //step in cmin
site=GInd::site_up(site,3-mu-nu);
}
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>=0)
{
chi-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmin_dir)); //step in cmin
site=GInd::site_up(site,cmin_dir);
}
chi+=cmin2;
}
return temp; 
}
_device_ _host_ floatT operator()(gSite site) {

typedef GIndexer<All, HaloDepth> GInd;
/// Define a SU(3) matrix
floatT result = 0;
SU3<floatT> temp1,temp2;
int mu,nu;
switch (_ex_dir) {
case 0: mu=1;nu=2;break;
case 1: mu=0;nu=2;break;
case 2: mu=0;nu=1;break;
case 3: for(nu=1;nu<3;nu++)
{
for(mu=0;mu<nu;mu++)
{
temp1=BresenhamLoopXYL(_x,_y,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_x,_y,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
//replacing x<->y
temp1=BresenhamLoopXYL(_y,_x,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_y,_x,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
}
}
return result/6.0;
break;
default: for (nu=1;nu<3;nu++)
{
for(mu=0;mu<nu;mu++)
{
temp1=BresenhamLoopXYL(_x,_y,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_x,_y,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
//replacing x<->y
temp1=BresenhamLoopXYL(_y,_x,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_y,_x,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
}
}
return result/6.0;
break;
}        
temp1=BresenhamLoopXYL(_x,_y,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_x,_y,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));

//replacing x<->y

temp1=BresenhamLoopXYL(_y,_x,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_y,_x,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
return result/2.0;
}

};