diff --git a/DESCRIPTION b/DESCRIPTION
index 5dff4ca..10934b6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: natcpp
 Title: Fast C++ Primitives for the 'NeuroAnatomy Toolbox'
-Version: 0.2
+Version: 0.3.0
 Authors@R: 
     person(given = "Gregory",
            family = "Jefferis",
@@ -18,9 +18,10 @@ Description: Fast functions implemented in C++ via 'Rcpp' to support the
 License: GPL (>= 3)
 URL: https://github.com/natverse/natcpp, https://natverse.org/natcpp/
 BugReports: https://github.com/natverse/natcpp/issues
-Imports: 
+Imports:
+    Matrix,
     Rcpp (>= 1.0.6)
-Suggests: 
+Suggests:
     spelling,
     testthat (>= 3.0.0)
 LinkingTo: 
@@ -28,4 +29,4 @@ LinkingTo:
 Config/testthat/edition: 3
 Encoding: UTF-8
 Language: en-GB
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
diff --git a/NAMESPACE b/NAMESPACE
index 504bfff..355a815 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,5 +10,7 @@ export(c_sub2ind)
 export(c_topntail)
 export(c_topntail_list)
 export(c_total_cable)
+export(c_weighted_jaccard_dense)
+export(c_weighted_jaccard_sparse)
 import(Rcpp)
 useDynLib(natcpp, .registration = TRUE)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index c90d42d..6474e92 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -135,3 +135,25 @@ c_EdgeListFromSegList <- function(L) {
     .Call(`_natcpp_c_EdgeListFromSegList`, L)
 }
 
+#' @noRd
+weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE) {
+    .Call(`_natcpp_weighted_jaccard_sparse_fill`, x, pattern, transpose)
+}
+
+#' Dense weighted Jaccard similarity via C++
+#'
+#' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using
+#' an adaptive dense accumulation strategy. For small output matrices, uses a
+#' feature-oriented loop; for large outputs, switches to a column-oriented
+#' loop for better cache performance. Intended for internal use by
+#' \code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}.
+#'
+#' @param x A dgCMatrix (sparse column-compressed matrix)
+#' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare
+#'   rows
+#' @return A dense numeric similarity matrix
+#' @export
+c_weighted_jaccard_dense <- function(x, transpose = FALSE) {
+    .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose)
+}
+
diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R
new file mode 100644
index 0000000..eb4eb35
--- /dev/null
+++ b/R/weighted_jaccard.R
@@ -0,0 +1,51 @@
+#' Sparse weighted Jaccard similarity via C++
+#'
+#' Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a
+#' sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute
+#' min-sums only for column (or row) pairs that share at least one non-zero
+#' feature, then normalises to similarity.
+#'
+#' @param x A dgCMatrix (sparse column-compressed matrix)
+#' @param transpose If \code{FALSE} (default), compare columns; if
+#'   \code{TRUE}, compare rows.
+#' @return A sparse dgCMatrix similarity matrix
+#' @export
+#' @seealso \code{\link{c_weighted_jaccard_dense}} for the dense equivalent
+#' @examples
+#' \dontrun{
+#' library(Matrix)
+#' m <- sparseMatrix(i = c(1,2,1,2,3,3), j = c(1,1,2,2,2,3),
+#'                   x = c(4,2,1,3,3,1), dims = c(3,3))
+#' c_weighted_jaccard_sparse(m)
+#' }
+c_weighted_jaccard_sparse <- function(x, transpose = FALSE) {
+  crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod
+  n <- if (transpose) nrow(x) else ncol(x)
+
+  if (length(x@x) == 0L) {
+    return(Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1,
+                                dims = c(n, n)))
+  }
+
+  # Sparsity pattern from binarised crossprod
+  b <- x
+  b@x <- rep(1, length(b@x))
+  A <- as(crossfun(b), "generalMatrix")
+
+  # Column/row totals for denominator
+  totals <- if (transpose) Matrix::rowSums(x) else Matrix::colSums(x)
+
+  # Fill in min_sums using C++
+  A@x <- weighted_jaccard_sparse_fill(x, A, transpose = transpose)
+
+  # Convert min_sums to similarity: sim = ms / (s_i + s_j - ms)
+  col_idx <- rep(seq_along(diff(A@p)), diff(A@p))
+  row_idx <- A@i + 1L
+  denom <- totals[row_idx] + totals[col_idx] - A@x
+  nonzero <- denom != 0
+  A@x[nonzero] <- A@x[nonzero] / denom[nonzero]
+  A@x[!nonzero] <- 0
+  Matrix::diag(A) <- 1
+
+  A
+}
diff --git a/man/c_weighted_jaccard_dense.Rd b/man/c_weighted_jaccard_dense.Rd
new file mode 100644
index 0000000..d63a502
--- /dev/null
+++ b/man/c_weighted_jaccard_dense.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{c_weighted_jaccard_dense}
+\alias{c_weighted_jaccard_dense}
+\title{Dense weighted Jaccard similarity via C++}
+\usage{
+c_weighted_jaccard_dense(x, transpose = FALSE)
+}
+\arguments{
+\item{x}{A dgCMatrix (sparse column-compressed matrix)}
+
+\item{transpose}{If \code{FALSE}, compare columns; if \code{TRUE}, compare
+rows}
+}
+\value{
+A dense numeric similarity matrix
+}
+\description{
+Compute the full weighted Jaccard similarity matrix for a dgCMatrix using
+an adaptive dense accumulation strategy. For small output matrices, uses a
+feature-oriented loop; for large outputs, switches to a column-oriented
+loop for better cache performance. Intended for internal use by
+\code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}.
+}
diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd
new file mode 100644
index 0000000..7f41e52
--- /dev/null
+++ b/man/c_weighted_jaccard_sparse.Rd
@@ -0,0 +1,34 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/weighted_jaccard.R
+\name{c_weighted_jaccard_sparse}
+\alias{c_weighted_jaccard_sparse}
+\title{Sparse weighted Jaccard similarity via C++}
+\usage{
+c_weighted_jaccard_sparse(x, transpose = FALSE)
+}
+\arguments{
+\item{x}{A dgCMatrix (sparse column-compressed matrix)}
+
+\item{transpose}{If \code{FALSE} (default), compare columns; if
+\code{TRUE}, compare rows.}
+}
+\value{
+A sparse dgCMatrix similarity matrix
+}
+\description{
+Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a
+sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute
+min-sums only for column (or row) pairs that share at least one non-zero
+feature, then normalises to similarity.
+}
+\examples{
+\dontrun{
+library(Matrix)
+m <- sparseMatrix(i = c(1,2,1,2,3,3), j = c(1,1,2,2,2,3),
+                  x = c(4,2,1,3,3,1), dims = c(3,3))
+c_weighted_jaccard_sparse(m)
+}
+}
+\seealso{
+\code{\link{c_weighted_jaccard_dense}} for the dense equivalent
+}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 3415f85..a0e9c71 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -135,6 +135,31 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// weighted_jaccard_sparse_fill
+NumericVector weighted_jaccard_sparse_fill(const S4& x, const S4& pattern, bool transpose);
+RcppExport SEXP _natcpp_weighted_jaccard_sparse_fill(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const S4& >::type x(xSEXP);
+    Rcpp::traits::input_parameter< const S4& >::type pattern(patternSEXP);
+    Rcpp::traits::input_parameter< bool >::type transpose(transposeSEXP);
+    rcpp_result_gen = Rcpp::wrap(weighted_jaccard_sparse_fill(x, pattern, transpose));
+    return rcpp_result_gen;
+END_RCPP
+}
+// c_weighted_jaccard_dense
+NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose);
+RcppExport SEXP _natcpp_c_weighted_jaccard_dense(SEXP xSEXP, SEXP transposeSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const S4& >::type x(xSEXP);
+    Rcpp::traits::input_parameter< bool >::type transpose(transposeSEXP);
+    rcpp_result_gen = Rcpp::wrap(c_weighted_jaccard_dense(x, transpose));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_natcpp_c_seglengths", (DL_FUNC) &_natcpp_c_seglengths, 4},
@@ -147,6 +172,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_natcpp_c_topntail", (DL_FUNC) &_natcpp_c_topntail, 1},
     {"_natcpp_c_topntail_list", (DL_FUNC) &_natcpp_c_topntail_list, 1},
     {"_natcpp_c_EdgeListFromSegList", (DL_FUNC) &_natcpp_c_EdgeListFromSegList, 1},
+    {"_natcpp_weighted_jaccard_sparse_fill", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 3},
+    {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 2},
     {NULL, NULL, 0}
 };
 
diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp
new file mode 100644
index 0000000..297bd32
--- /dev/null
+++ b/src/weighted_jaccard.cpp
@@ -0,0 +1,278 @@
+#include <Rcpp.h>
+#include <algorithm>
+#include <vector>
+
+using namespace Rcpp;
+
+//' @noRd
+// [[Rcpp::export]]
+NumericVector weighted_jaccard_sparse_fill(
+    const S4& x, const S4& pattern, bool transpose = false) {
+
+  // Original matrix slots (dgCMatrix: nr x nc)
+  IntegerVector dims = x.slot("Dim");
+  IntegerVector xi = x.slot("i");
+  IntegerVector xp = x.slot("p");
+  NumericVector xx = x.slot("x");
+
+  const int nr = dims[0];
+  const int nc = dims[1];
+  const int nfeat = transpose ? nc : nr;
+
+  // Pattern matrix slots (dgCMatrix: ncomp x ncomp)
+  IntegerVector Ai = pattern.slot("i");
+  IntegerVector Ap = pattern.slot("p");
+  IntegerVector Adims = pattern.slot("Dim");
+  const int ncomp = Adims[0];
+  const int annz = Ai.size();
+
+  // Output: min_sums values aligned with pattern's (i, p) structure
+  NumericVector out(annz, 0.0);
+
+  // Build CSR for feature dimension: for each feature f, the items (columns
+  // or rows being compared) that are nonzero in that feature, with their values.
+  std::vector<int> feat_count(nfeat, 0);
+  if (!transpose) {
+    for (int col = 0; col < nc; ++col)
+      for (int idx = xp[col]; idx < xp[col + 1]; ++idx)
+        feat_count[xi[idx]]++;
+  } else {
+    for (int col = 0; col < nc; ++col)
+      feat_count[col] = xp[col + 1] - xp[col];
+  }
+
+  std::vector<int> feat_offsets(nfeat + 1, 0);
+  for (int f = 0; f < nfeat; ++f)
+    feat_offsets[f + 1] = feat_offsets[f] + feat_count[f];
+
+  const int total_nnz = feat_offsets[nfeat];
+  std::vector<int> feat_items(total_nnz);   // which item (col/row being compared)
+  std::vector<double> feat_vals(total_nnz);  // value of x at that position
+  std::vector<int> feat_fill(feat_offsets.begin(), feat_offsets.end());
+
+  if (!transpose) {
+    for (int col = 0; col < nc; ++col) {
+      for (int idx = xp[col]; idx < xp[col + 1]; ++idx) {
+        const int row = xi[idx];
+        const int dest = feat_fill[row]++;
+        feat_items[dest] = col;
+        feat_vals[dest] = xx[idx];
+      }
+    }
+  } else {
+    for (int col = 0; col < nc; ++col) {
+      for (int idx = xp[col]; idx < xp[col + 1]; ++idx) {
+        const int dest = feat_fill[col]++;
+        feat_items[dest] = xi[idx];
+        feat_vals[dest] = xx[idx];
+      }
+    }
+  }
+
+  // Also need item->features CSR: for each compared item, which features
+  // it participates in and with what value.
+  std::vector<int> item_count(ncomp, 0);
+  for (int f = 0; f < nfeat; ++f) {
+    for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx)
+      item_count[feat_items[idx]]++;
+  }
+
+  std::vector<int> item_offsets(ncomp + 1, 0);
+  for (int c = 0; c < ncomp; ++c)
+    item_offsets[c + 1] = item_offsets[c] + item_count[c];
+
+  std::vector<int> item_feats(total_nnz);   // which feature
+  std::vector<double> item_vals(total_nnz);  // value
+  std::vector<int> item_fill(item_offsets.begin(), item_offsets.end());
+
+  for (int f = 0; f < nfeat; ++f) {
+    for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) {
+      const int item = feat_items[idx];
+      const int dest = item_fill[item]++;
+      item_feats[dest] = f;
+      item_vals[dest] = feat_vals[idx];
+    }
+  }
+
+  // row_to_pos[row] = position in out[] for entry (row, current_col).
+  // Populated once per output column, then cleared.
+  std::vector<int> row_to_pos(ncomp, -1);
+
+  // Iterate over output columns
+  for (int cb = 0; cb < ncomp; ++cb) {
+    // Populate row_to_pos for pattern column cb
+    for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx)
+      row_to_pos[Ai[idx]] = idx;
+
+    // For each feature f where cb is nonzero
+    for (int fi = item_offsets[cb]; fi < item_offsets[cb + 1]; ++fi) {
+      const int f = item_feats[fi];
+      const double vb = item_vals[fi];
+
+      // For each other item ca also in this feature, accumulate min into (ca, cb)
+      for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) {
+        const int ca = feat_items[idx];
+        const int pos = row_to_pos[ca];
+        if (pos >= 0)
+          out[pos] += std::min(feat_vals[idx], vb);
+      }
+    }
+
+    // Clear row_to_pos
+    for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx)
+      row_to_pos[Ai[idx]] = -1;
+  }
+
+  return out;
+}
+
+//' Dense weighted Jaccard similarity via C++
+//'
+//' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using
+//' an adaptive dense accumulation strategy. For small output matrices, uses a
+//' feature-oriented loop; for large outputs, switches to a column-oriented
+//' loop for better cache performance. Intended for internal use by
+//' \code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}.
+//'
+//' @param x A dgCMatrix (sparse column-compressed matrix)
+//' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare
+//'   rows
+//' @return A dense numeric similarity matrix
+//' @export
+// [[Rcpp::export]]
+NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false) {
+  IntegerVector dims = x.slot("Dim");
+  IntegerVector xi = x.slot("i");
+  IntegerVector xp = x.slot("p");
+  NumericVector xx = x.slot("x");
+
+  const int nr = dims[0];
+  const int nc = dims[1];
+  const int ncomp = transpose ? nr : nc;
+  const int nfeat = transpose ? nc : nr;
+
+  // Totals per compared item
+  std::vector<double> totals(ncomp, 0.0);
+
+  // Build feature CSR: for each feature f, which items are nonzero
+  std::vector<int> feat_count(nfeat, 0);
+  if (!transpose) {
+    for (int col = 0; col < nc; ++col) {
+      for (int idx = xp[col]; idx < xp[col + 1]; ++idx) {
+        totals[col] += xx[idx];
+        feat_count[xi[idx]]++;
+      }
+    }
+  } else {
+    for (int col = 0; col < nc; ++col) {
+      feat_count[col] = xp[col + 1] - xp[col];
+      for (int idx = xp[col]; idx < xp[col + 1]; ++idx) {
+        totals[xi[idx]] += xx[idx];
+      }
+    }
+  }
+
+  std::vector<int> feat_offsets(nfeat + 1, 0);
+  for (int f = 0; f < nfeat; ++f)
+    feat_offsets[f + 1] = feat_offsets[f] + feat_count[f];
+
+  const int total_nnz = feat_offsets[nfeat];
+  std::vector<int> feat_items(total_nnz);
+  std::vector<double> feat_vals(total_nnz);
+  std::vector<int> feat_fill(feat_offsets.begin(), feat_offsets.end());
+
+  if (!transpose) {
+    for (int col = 0; col < nc; ++col) {
+      for (int idx = xp[col]; idx < xp[col + 1]; ++idx) {
+        const int row = xi[idx];
+        const int dest = feat_fill[row]++;
+        feat_items[dest] = col;
+        feat_vals[dest] = xx[idx];
+      }
+    }
+  } else {
+    for (int col = 0; col < nc; ++col) {
+      for (int idx = xp[col]; idx < xp[col + 1]; ++idx) {
+        const int dest = feat_fill[col]++;
+        feat_items[dest] = xi[idx];
+        feat_vals[dest] = xx[idx];
+      }
+    }
+  }
+
+  NumericMatrix out(ncomp, ncomp);
+
+  // Threshold: ncomp^2 * 8 bytes > ~12MB (typical L3 share per core)
+  const bool use_colwise = (static_cast<std::size_t>(ncomp) * ncomp * 8 > 12 * 1024 * 1024);
+
+  if (use_colwise) {
+    // Column-oriented: build item->features CSR, iterate output columns
+    std::vector<int> item_count(ncomp, 0);
+    for (int f = 0; f < nfeat; ++f)
+      for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx)
+        item_count[feat_items[idx]]++;
+
+    std::vector<int> item_offsets(ncomp + 1, 0);
+    for (int c = 0; c < ncomp; ++c)
+      item_offsets[c + 1] = item_offsets[c] + item_count[c];
+
+    std::vector<int> item_feats(total_nnz);
+    std::vector<double> item_vals(total_nnz);
+    std::vector<int> item_fill(item_offsets.begin(), item_offsets.end());
+
+    for (int f = 0; f < nfeat; ++f) {
+      for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) {
+        const int item = feat_items[idx];
+        const int dest = item_fill[item]++;
+        item_feats[dest] = f;
+        item_vals[dest] = feat_vals[idx];
+      }
+    }
+
+    for (int cb = 0; cb < ncomp; ++cb) {
+      double* col = &out[static_cast<std::size_t>(cb) * ncomp];
+      for (int fi = item_offsets[cb]; fi < item_offsets[cb + 1]; ++fi) {
+        const int f = item_feats[fi];
+        const double vb = item_vals[fi];
+        for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) {
+          col[feat_items[idx]] += std::min(feat_vals[idx], vb);
+        }
+      }
+    }
+  } else {
+    // Feature-oriented: iterate features, scatter into output
+    for (int f = 0; f < nfeat; ++f) {
+      const int start = feat_offsets[f];
+      const int end = feat_offsets[f + 1];
+      const int k = end - start;
+
+      for (int a = 0; a < k; ++a) {
+        const int ca = feat_items[start + a];
+        const double va = feat_vals[start + a];
+        out[static_cast<std::size_t>(ca) * ncomp + ca] += va;
+        for (int b = a + 1; b < k; ++b) {
+          const int cb = feat_items[start + b];
+          const double mv = std::min(va, feat_vals[start + b]);
+          out[static_cast<std::size_t>(ca) * ncomp + cb] += mv;
+          out[static_cast<std::size_t>(cb) * ncomp + ca] += mv;
+        }
+      }
+    }
+  }
+
+  // Convert min_sums to similarity
+  for (int cb = 0; cb < ncomp; ++cb) {
+    for (int ca = 0; ca < ncomp; ++ca) {
+      if (ca == cb) {
+        out[static_cast<std::size_t>(cb) * ncomp + ca] = 1.0;
+      } else {
+        const std::size_t pos = static_cast<std::size_t>(cb) * ncomp + ca;
+        const double ms = out[pos];
+        const double denom = totals[ca] + totals[cb] - ms;
+        out[pos] = (denom > 0.0) ? ms / denom : 0.0;
+      }
+    }
+  }
+
+  return out;
+}
diff --git a/tests/testthat/test-weighted-jaccard.R b/tests/testthat/test-weighted-jaccard.R
new file mode 100644
index 0000000..e1c1f11
--- /dev/null
+++ b/tests/testthat/test-weighted-jaccard.R
@@ -0,0 +1,65 @@
+test_that("c_weighted_jaccard_dense computes correct similarity", {
+  m <- Matrix::sparseMatrix(
+    i = c(1L, 2L, 1L, 2L, 3L, 3L),
+    j = c(1L, 1L, 2L, 2L, 2L, 3L),
+    x = c(4, 2, 1, 3, 3, 1),
+    dims = c(3L, 3L)
+  )
+
+  # Naive reference: sim[a,b] = sum(pmin(col_a, col_b)) / sum(pmax(col_a, col_b))
+  dm <- as.matrix(m)
+  n <- ncol(dm)
+  ref <- matrix(0, n, n)
+  for (a in seq_len(n)) {
+    for (b in seq_len(n)) {
+      mins <- sum(pmin(dm[, a], dm[, b]))
+      maxs <- sum(pmax(dm[, a], dm[, b]))
+      ref[a, b] <- if (maxs > 0) mins / maxs else 0
+    }
+  }
+  diag(ref) <- 1
+
+  sim <- c_weighted_jaccard_dense(m, transpose = FALSE)
+  expect_equal(sim, ref, tolerance = 1e-12)
+})
+
+test_that("c_weighted_jaccard_dense transpose works", {
+  m <- Matrix::sparseMatrix(
+    i = c(1L, 2L, 1L, 2L, 3L, 3L),
+    j = c(1L, 1L, 2L, 2L, 2L, 3L),
+    x = c(4, 2, 1, 3, 3, 1),
+    dims = c(3L, 3L)
+  )
+
+  dm <- as.matrix(m)
+  nr <- nrow(dm)
+  ref_t <- matrix(0, nr, nr)
+  for (a in seq_len(nr)) {
+    for (b in seq_len(nr)) {
+      mins <- sum(pmin(dm[a, ], dm[b, ]))
+      maxs <- sum(pmax(dm[a, ], dm[b, ]))
+      ref_t[a, b] <- if (maxs > 0) mins / maxs else 0
+    }
+  }
+  diag(ref_t) <- 1
+
+  sim_t <- c_weighted_jaccard_dense(m, transpose = TRUE)
+  expect_equal(sim_t, ref_t, tolerance = 1e-12)
+})
+
+test_that("c_weighted_jaccard_sparse matches dense", {
+  m <- Matrix::sparseMatrix(
+    i = c(1L, 2L, 1L, 2L, 3L, 3L),
+    j = c(1L, 1L, 2L, 2L, 2L, 3L),
+    x = c(4, 2, 1, 3, 3, 1),
+    dims = c(3L, 3L)
+  )
+
+  sim_dense <- c_weighted_jaccard_dense(m, transpose = FALSE)
+  sim_sparse <- c_weighted_jaccard_sparse(m, transpose = FALSE)
+  expect_equal(as.matrix(sim_sparse), sim_dense, tolerance = 1e-12)
+
+  sim_dense_t <- c_weighted_jaccard_dense(m, transpose = TRUE)
+  sim_sparse_t <- c_weighted_jaccard_sparse(m, transpose = TRUE)
+  expect_equal(as.matrix(sim_sparse_t), sim_dense_t, tolerance = 1e-12)
+})