complexity is again O ( N $^{4}$) due to the preprocessing phase. F A fast algorithm for computing the HHG statis- tic Here we describe a fast algorithm for computing the univariate original (distribution- dependent) HHG test statistic. If D$_{X}$ ( x$_{1}$, x$_{2}$ ) and D$_{Y}$ ( y$_{1}$, y$_{2}$ ) are distance metrics in the variables tested for indepen- dence (for which, again, we have a paired sample with N i.i.d samples), the HHG test requires computing N ( N − 1) different 2 × 2 contingency tables according to the following partitions of the distance-distance plane. For every “origin” sample i , and for every “radius” sample j , the remaining N − 2 samples are classified according to whether their X and Y distances from i are both smaller than the X and Y distances from j to i , or if only X , only Y , or neither are smaller than the respective j to i distances. In the general case, generating contingency tables for all pairs can be done in N $^{2}$log N , as described in (Heller et al., 2013). In the univariate case, an O ( N $^{2}$) algorithm proceeds as follows. Instead of working in the distance-distance plane, the algorithm is specified in terms of the ( x, y ), i.e., the sample plane. It is sufficient to consider the discrete grid expanded from unique X samples and unique Y samples actually observed (these can be identified in O ( N log N )). The double cumulative sum over this N × N grid is computed as in Section 3.1 of the main text, with the only difference being that after A is initialized to all zeros, it is updated sequentially with A ( r$_{i}$, s$_{i}$ ) = A ( r$_{i}$, s$_{i}$ ) + 1 for every sample of ranks in x and ranks in y , ( r$_{i}$, s$_{i}$ ) , i = 1 , 2 , . . . , N , to account for possible ties. 64