diff --git a/DESCRIPTION b/DESCRIPTION
index f0cd18884..341c26afc 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: DescTools
 Type: Package
 Title: Tools for Descriptive Statistics
-Version: 0.99.52.2
+Version: 0.99.52.3
 Date: 2023-12-10
 Authors@R: c(
   person("Andri", "Signorell", email = "andri@signorell.net", role = c("aut", "cre")),
diff --git a/NEWS b/NEWS
index 87b37c6f7..f5005194a 100644
--- a/NEWS
+++ b/NEWS
@@ -5,10 +5,10 @@ DescTools 0.99.53 (2023-12-30)
 ------------------------------
 
 NEW FUNCTIONS ADDED:
-  * 
+  * New implementation of 2-sample HodgesLehmann estimator.
 
 UPDATED FUNCTIONS:
-  * 
+  * pckg.yaml has a dTri entry now.
 
 BUGFIXES:
   * HodgesLehman() mixed up the one sample and two sample case.
diff --git a/R/RcppExports.R b/R/RcppExports.R
index f63cc566e..d48b040fd 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -57,11 +57,11 @@ check.gompertz <- function(shape, rate) {
     .Call(`_DescTools_check_gompertz`, shape, rate)
 }
 
-hl2qest <- function(x, y) {
-    .Call(`_DescTools_hl2qest`, x, y)
-}
-
 hlqest <- function(x) {
     .Call(`_DescTools_hlqest`, x)
 }
 
+hl2qest <- function(x, y) {
+    .Call(`_DescTools_hl2qest`, x, y)
+}
+
diff --git a/src/DescTools.dll b/src/DescTools.dll
index 0f06eb9b1..2dcb6dfc4 100644
Binary files a/src/DescTools.dll and b/src/DescTools.dll differ
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 8708408c6..9f4ff6412 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -175,26 +175,26 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// hl2qest
-double hl2qest(NumericVector x, NumericVector y);
-RcppExport SEXP _DescTools_hl2qest(SEXP xSEXP, SEXP ySEXP) {
+// hlqest
+double hlqest(NumericVector x);
+RcppExport SEXP _DescTools_hlqest(SEXP xSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP);
-    Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP);
-    rcpp_result_gen = Rcpp::wrap(hl2qest(x, y));
+    rcpp_result_gen = Rcpp::wrap(hlqest(x));
     return rcpp_result_gen;
 END_RCPP
 }
-// hlqest
-double hlqest(NumericVector x);
-RcppExport SEXP _DescTools_hlqest(SEXP xSEXP) {
+// hl2qest
+double hl2qest(NumericVector x, NumericVector y);
+RcppExport SEXP _DescTools_hl2qest(SEXP xSEXP, SEXP ySEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP);
-    rcpp_result_gen = Rcpp::wrap(hlqest(x));
+    Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP);
+    rcpp_result_gen = Rcpp::wrap(hl2qest(x, y));
     return rcpp_result_gen;
 END_RCPP
 }
diff --git a/src/RcppExports.o b/src/RcppExports.o
index b0b58db03..310acab6c 100644
Binary files a/src/RcppExports.o and b/src/RcppExports.o differ
diff --git a/src/hl2qest.cpp b/src/hl2qest.cpp
deleted file mode 100644
index 233f6c8a7..000000000
--- a/src/hl2qest.cpp
+++ /dev/null
@@ -1,175 +0,0 @@
-
-// Port of Monahans algorithm for Hodges-Lehman estimator
-// https://dl.acm.org/doi/10.1145/1271.319414
-// 
-  // by Cyril Flurin Moser 
-
-// 2023-11-29
-
-
-
-#include <Rcpp.h>
-
-#include <iostream>
-#include <vector>
-#include<numeric>
-
-#include <algorithm>
-#include <cstdio>
-#include <cstdlib>
-
-using namespace std;
-using namespace Rcpp;
-
- 
-void *safe_malloc2(size_t n) {
-  void *p = malloc(n);
-  if (p == NULL) {
-    stop("Fatal: failed to allocate enough memory.\n");
-  }
-  return p;
-};
-
-
-
-// [[Rcpp::export]]
-double hl2qest(NumericVector x, NumericVector y) {
-  
-  int m = x.size();
-  int n = y.size();
-  
-  std::sort(x.begin(), x.end(), [](double &a, double &b){ return a<b; });
-  std::sort(y.begin(), y.end(), [](double &a, double &b){ return a<b; });
-  
-  int nn = m * n;
-  int k1 = (nn + 1) / 2;
-  int k2 = (nn + 2) / 2;
-  
-  if (n == 1) {
-    return (x[k1 - 1] + x[k2 - 1]) / 2. - y[0];
-  }
-  if (m == 1) {
-    return x[0] - (y[k1 - 1] + y[k2 - 1]) / 2.;
-  }
-  int *lb = (int *)safe_malloc2(sizeof(int) * n);
-  int *rb = (int *)safe_malloc2(sizeof(int) * n);
-  int *q = (int *)safe_malloc2(sizeof(int) * n);
-  
-  double retval = 0.0;
-  double am, amx, amn;
-  
-  for (int i = 0; i < n; i++) {
-    lb[i] = std::max(1, (m + 1) - (nn + 1 - k1) / (i + 1));
-    rb[i] = std::min(m, k2 / (n - i));
-  }
-  
-  int sm = nn;  // number in set s at step m
-  int l = 0;    // number of pairs less than those in set s at step m
-  int j = 0;
-  int sq = 0;
-  
-  int method = 1;
-  
-  for (;;) {
-    switch (method) {
-      case 1:
-        am = x[((m + 1) / 2) - 1] - y[((n + 1) / 2) - 1];
-        break;
-        case 2:
-          amx = x[0] - y[n - 1];
-          amn = x[m - 1] - y[0];
-          for (int i = 0; i < n; i++) {
-            if (lb[i] <= rb[i]) {  // skip this row if no element in it is in set s on this step
-              int lbi = lb[i] - 1;
-              amn = std::min(amn, x[lbi] - y[i]);  // get the smallest in this row
-              int rbi = rb[i] - 1;
-              amx = std::max(amx, x[rbi] - y[i]);  // get the largest in this row
-            }
-          }
-          am = (amx + amn) / 2.;
-          if ((am <= amn) || (am > amx)) {  // roundoff
-            am = amx;
-          }
-          if ((amn == amx) || (sm == 2)) {
-            retval = am;
-            goto cleanup;
-          }
-          break;
-          default:
-            int k = std::rand() % sm;
-          for (int i = 0; i < n; i++) {
-            j = i;
-            if (k <= (rb[i] - lb[i]))
-              break;
-            k = k - rb[i] + lb[i] - 1;
-          }
-          int mdlrow = (lb[j] + rb[j]) / 2 - 1;
-          am = x[mdlrow] - y[j];
-          break;
-    }
-    // Partition Step
-    j = 0;
-    sq = 0;
-    for (int i = 1; i <= n; i++) {
-      q[i - 1] = j;
-      while (j < m) {
-        if ((x[j] - y[i - 1]) >= am)
-          break;
-        j++;
-      }
-      q[i - 1] = j;
-      sq += j;
-    }
-    // start branching
-    if (sq == l) {
-      method = 2;
-    } else {
-      if ((sq == (k2 - 1)) || (sq == k1))
-        break;
-      if (sq > k1) {
-        for (int i = 0; i < n; i++) {
-          rb[i] = q[i];
-        }
-      } else {
-        for (int i = 0; i < n; i++) {
-          lb[i] = q[i] + 1;
-        }
-      }
-      l = 0;
-      sm = 0;
-      for (int i = 0; i < n; i++) {
-        l += lb[i] - 1;
-        sm += rb[i] - lb[i] + 1;
-      }
-      method = 3;
-    }
-    if (sm == 2)
-      method = 2;
-  }
-  
-  amn = x[m - 1] - y[0];
-  amx = x[0] - y[n - 1];
-  for (int i = 0; i < n; i++) {
-    int iq = q[i];
-    if (iq > 0)
-      amx = std::max(amx, x[iq - 1] - y[i]);
-    if (iq < m)
-      amn = std::min(amn, x[iq] - y[i]);
-  }
-  retval = (amn + amx) / 2.;
-  if (k1 < k2) {
-    goto cleanup;
-  } else if (sq == k1) {
-    retval = amx;
-  } else if (sq == k1 - 1) {
-    retval = amn;
-  }
-  
-  cleanup:
-    free(lb);
-    free(rb);
-    free(q);
-  
-  return retval;
-}
-;
diff --git a/src/hl2qest.o b/src/hl2qest.o
deleted file mode 100644
index 16c59a227..000000000
Binary files a/src/hl2qest.o and /dev/null differ
diff --git a/src/hlqest.cpp b/src/hlqest.cpp
index f2962043f..3144776ce 100644
--- a/src/hlqest.cpp
+++ b/src/hlqest.cpp
@@ -30,28 +30,6 @@ void *safe_malloc(size_t n) {
   return p;
 };
 
-// double rng(int ixx) {
-//   int a = 16807;
-//   int b15 = 32768;
-//   int b16 = 65536;
-//   int p = 2147483647;
-//   static int ix = 0;
-//   if (ix == 0) {
-//     ix = ixx;
-//   }
-//   int xhi = ix / b16;
-//   int xalo = (ix - xhi * b16) * a;
-//   int leftlo = xalo / b16;
-//   int fhi = xhi * a + leftlo;
-//   int k = fhi / b15;
-//   ix = (((xalo - leftlo * b16) - p) + (fhi - k * b15) * b16) + k;
-//   if (ix < 0) {
-//     ix = ix + p;
-//   }
-//   return ((float)ix) * 4.656612875E-10;
-// }
-// ;
-
 
 int rng(int max) //range : [0, max]
 {
@@ -60,20 +38,11 @@ int rng(int max) //range : [0, max]
 }  
 
 
-// NumericVector stl_sort(NumericVector x) {
-//   NumericVector y = clone(x);
-//   std::sort(y.begin(), y.end());
-//   return y;
-// }
-
-
-
 // [[Rcpp::export]]
 double hlqest(NumericVector x) {
   
   int n = x.size();
   
-  // x = stl_sort(x);
   std::sort(x.begin(), x.end(), [](double &a, double &b){ return a<b; });
   
   if (n <= 2) {  // special cases for n=1 and n=2
@@ -258,3 +227,150 @@ double hlqest(NumericVector x) {
 }
 ;
 
+
+
+// [[Rcpp::export]]
+double hl2qest(NumericVector x, NumericVector y) {
+  
+  int m = x.size();
+  int n = y.size();
+  
+  std::sort(x.begin(), x.end(), [](double &a, double &b){ return a<b; });
+  std::sort(y.begin(), y.end(), [](double &a, double &b){ return a<b; });
+  
+  int nn = m * n;
+  int k1 = (nn + 1) / 2;
+  int k2 = (nn + 2) / 2;
+  
+  if (n == 1) {
+    return (x[k1 - 1] + x[k2 - 1]) / 2. - y[0];
+  }
+  if (m == 1) {
+    return x[0] - (y[k1 - 1] + y[k2 - 1]) / 2.;
+  }
+  int *lb = (int *)safe_malloc(sizeof(int) * n);
+  int *rb = (int *)safe_malloc(sizeof(int) * n);
+  int *q = (int *)safe_malloc(sizeof(int) * n);
+  
+  double retval = 0.0;
+  double am, amx, amn;
+  
+  for (int i = 0; i < n; i++) {
+    lb[i] = std::max(1, (m + 1) - (nn + 1 - k1) / (i + 1));
+    rb[i] = std::min(m, k2 / (n - i));
+  }
+  
+  int sm = nn;  // number in set s at step m
+  int l = 0;    // number of pairs less than those in set s at step m
+  int j = 0;
+  int sq = 0;
+  
+  int method = 1;
+  
+  for (;;) {
+    switch (method) {
+    case 1:
+      am = x[((m + 1) / 2) - 1] - y[((n + 1) / 2) - 1];
+      break;
+    case 2:
+      amx = x[0] - y[n - 1];
+      amn = x[m - 1] - y[0];
+      for (int i = 0; i < n; i++) {
+        if (lb[i] <= rb[i]) {  // skip this row if no element in it is in set s on this step
+          int lbi = lb[i] - 1;
+          amn = std::min(amn, x[lbi] - y[i]);  // get the smallest in this row
+          int rbi = rb[i] - 1;
+          amx = std::max(amx, x[rbi] - y[i]);  // get the largest in this row
+        }
+      }
+      am = (amx + amn) / 2.;
+      if ((am <= amn) || (am > amx)) {  // roundoff
+        am = amx;
+      }
+      if ((amn == amx) || (sm == 2)) {
+        retval = am;
+        goto cleanup;
+      }
+      break;
+    default:
+//      int k = std::rand() % sm;
+      int k = rng(sm);  // k is a random int from 0 to sm
+
+    for (int i = 0; i < n; i++) {
+      j = i;
+      if (k <= (rb[i] - lb[i]))
+        break;
+      k = k - rb[i] + lb[i] - 1;
+    }
+    int mdlrow = (lb[j] + rb[j]) / 2 - 1;
+    am = x[mdlrow] - y[j];
+    break;
+    }
+    // Partition Step
+    j = 0;
+    sq = 0;
+    for (int i = 1; i <= n; i++) {
+      q[i - 1] = j;
+      while (j < m) {
+        if ((x[j] - y[i - 1]) >= am)
+          break;
+        j++;
+      }
+      q[i - 1] = j;
+      sq += j;
+    }
+    // start branching
+    if (sq == l) {
+      method = 2;
+    } else {
+      if ((sq == (k2 - 1)) || (sq == k1))
+        break;
+      if (sq > k1) {
+        for (int i = 0; i < n; i++) {
+          rb[i] = q[i];
+        }
+      } else {
+        for (int i = 0; i < n; i++) {
+          lb[i] = q[i] + 1;
+        }
+      }
+      l = 0;
+      sm = 0;
+      for (int i = 0; i < n; i++) {
+        l += lb[i] - 1;
+        sm += rb[i] - lb[i] + 1;
+      }
+      method = 3;
+    }
+    if (sm == 2)
+      method = 2;
+  }
+  
+  amn = x[m - 1] - y[0];
+  amx = x[0] - y[n - 1];
+  for (int i = 0; i < n; i++) {
+    int iq = q[i];
+    if (iq > 0)
+      amx = std::max(amx, x[iq - 1] - y[i]);
+    if (iq < m)
+      amn = std::min(amn, x[iq] - y[i]);
+  }
+  retval = (amn + amx) / 2.;
+  if (k1 < k2) {
+    goto cleanup;
+  } else if (sq == k1) {
+    retval = amx;
+  } else if (sq == k1 - 1) {
+    retval = amn;
+  }
+  
+  cleanup:
+    free(lb);
+    free(rb);
+    free(q);
+  
+  return retval;
+}
+;
+
+
diff --git a/src/hlqest.o b/src/hlqest.o
index 79203f416..fd16c81a1 100644
Binary files a/src/hlqest.o and b/src/hlqest.o differ