diff --git a/DESCRIPTION b/DESCRIPTION
index 826b3ba..1d79c31 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: natcpp
 Title: Fast C++ Primitives for the 'NeuroAnatomy Toolbox'
-Version: 0.3.1
+Version: 0.3.1.9000
 Authors@R: 
     person(given = "Gregory",
            family = "Jefferis",
@@ -31,4 +31,4 @@ LinkingTo:
 Config/testthat/edition: 3
 Encoding: UTF-8
 Language: en-GB
-Config/roxygen2/version: 8.0.0
+Config/roxygen2/version: 8.1.0
diff --git a/NEWS.md b/NEWS.md
index 0e675df..d9cf5f5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,11 @@
+# natcpp 0.3.1.9000
+
+* the `threads` argument of `c_weighted_jaccard_sparse()` and
+  `c_weighted_jaccard_dense()` now defaults to `NULL`, applying a package-wide
+  thread policy (respecting `getOption("Ncpus")` and `OMP_THREAD_LIMIT`, else a
+  conservative 2) rather than a hard-coded 4. Pass `threads = 0` for all cores,
+  or an integer to override. No new package dependency.
+
 # natcpp 0.3.1
 
 * fix package loading on CRAN's Debian clang r-devel check by conditionally
diff --git a/R/RcppExports.R b/R/RcppExports.R
index c113355..c800624 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -140,27 +140,7 @@ weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_
     .Call(`_natcpp_weighted_jaccard_sparse_fill`, x, pattern, transpose, display_progress, threads)
 }
 
-#' Dense weighted Jaccard similarity via C++
-#'
-#' Compute the full weighted Jaccard similarity matrix for a
-#' \link[Matrix:dgCMatrix-class]{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.
-#'
-#' @param x A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix)
-#' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare
-#'   rows
-#' @param threads Number of threads (default 4). Set to 0 for all cores.
-#' @param triangle If \code{TRUE}, return only the lower triangle as a flat
-#'   numeric vector in \code{\link{dist}} layout. If \code{FALSE} (default),
-#'   return a full square matrix.
-#' @param distance If \code{TRUE}, return distance (\code{1 - similarity})
-#'   instead of similarity. Default \code{FALSE}.
-#' @return A dense numeric similarity matrix, or a numeric vector in
-#'   \code{\link{dist}} layout when \code{triangle = TRUE}.
-#' @export
-c_weighted_jaccard_dense <- function(x, transpose = FALSE, threads = 4L, triangle = FALSE, distance = FALSE) {
+weighted_jaccard_dense_impl <- function(x, transpose = FALSE, threads = 4L, triangle = FALSE, distance = FALSE) {
     .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose, threads, triangle, distance)
 }
 
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000..3d310a9
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,30 @@
+#' Resolve the number of threads for natcpp's parallel functions
+#'
+#' Internal helper implementing the package-wide default thread policy. A
+#' function that accepts a \code{threads} argument should default it to
+#' \code{NULL} and pass it through here.
+#'
+#' The policy follows the advice given by Dirk Eddelbuettel in
+#' \url{https://github.com/Rdatatable/data.table/issues/5658}: an explicit
+#' user value always wins; otherwise fall back to \code{getOption("Ncpus")}
+#' (which CRAN sets for its check farm), then the \code{OMP_THREAD_LIMIT}
+#' environment variable, and finally a conservative default of 2. Deliberately
+#' avoids grabbing every available core by default, and needs no dependency on
+#' \pkg{parallel}. To use all cores, pass \code{threads = 0} (the C++ back ends
+#' interpret 0 as \code{hardware_concurrency()}).
+#'
+#' @param threads User-supplied thread count, or \code{NULL} to apply the
+#'   default policy. A supplied value (including 0) is returned unchanged.
+#' @return An integer thread count.
+#' @noRd
+natcpp_threads <- function(threads = NULL) {
+  if (!is.null(threads)) return(as.integer(threads)[1])
+  n <- getOption("Ncpus")
+  if (is.null(n)) {
+    tl <- Sys.getenv("OMP_THREAD_LIMIT")
+    if (nzchar(tl)) n <- tl
+  }
+  n <- suppressWarnings(as.integer(n))[1]
+  if (length(n) == 0L || is.na(n) || n < 1L) n <- 2L
+  n
+}
diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R
index 870a61e..9f1e778 100644
--- a/R/weighted_jaccard.R
+++ b/R/weighted_jaccard.R
@@ -12,8 +12,10 @@
 #'   \code{TRUE}, compare rows.
 #' @param display_progress Whether to show a text progress bar (default
 #'   \code{TRUE}).
-#' @param threads Number of threads for parallel computation (default 4).
-#'   Set to 0 to use all available cores.
+#' @param threads Number of threads for parallel computation. The default
+#'   \code{NULL} applies the package thread policy (respecting
+#'   \code{getOption("Ncpus")} and the \code{OMP_THREAD_LIMIT} environment
+#'   variable, else 2). Set to 0 to use all available cores.
 #' @param triangle If \code{TRUE}, return a symmetric \code{dsCMatrix}
 #'   (upper triangle only). If \code{FALSE} (default), return a general
 #'   \link[Matrix:dgCMatrix-class]{dgCMatrix}.
@@ -33,8 +35,9 @@
 #' c_weighted_jaccard_sparse(m)
 #' }
 c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = TRUE,
-                                      threads = 4L, triangle = FALSE,
+                                      threads = NULL, triangle = FALSE,
                                       distance = FALSE) {
+  threads <- natcpp_threads(threads)
   if (distance)
     warning("distance=TRUE with sparse output produces a mostly-dense matrix; ",
             "consider using c_weighted_jaccard_dense() with triangle=TRUE instead.")
@@ -85,3 +88,43 @@ c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = T
 
   A
 }
+
+#' Dense weighted Jaccard similarity via C++
+#'
+#' Compute the full weighted Jaccard similarity matrix for a
+#' \link[Matrix:dgCMatrix-class]{dgCMatrix}, returning a dense matrix (or a
+#' \code{\link{dist}}-layout vector).
+#'
+#' @details Uses an adaptive dense accumulation strategy: for small output
+#'   matrices a feature-oriented loop, switching to a column-oriented loop for
+#'   larger outputs for better cache performance.
+#'
+#' @param x A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix)
+#' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare
+#'   rows
+#' @param threads Number of threads for parallel computation. The default
+#'   \code{NULL} applies the package thread policy (respecting
+#'   \code{getOption("Ncpus")} and the \code{OMP_THREAD_LIMIT} environment
+#'   variable, else 2). Set to 0 to use all available cores.
+#' @param triangle If \code{TRUE}, return only the lower triangle as a flat
+#'   numeric vector in \code{\link{dist}} layout. If \code{FALSE} (default),
+#'   return a full square matrix.
+#' @param distance If \code{TRUE}, return distance (\code{1 - similarity})
+#'   instead of similarity. Default \code{FALSE}.
+#' @return A dense numeric similarity matrix, or a numeric vector in
+#'   \code{\link{dist}} layout when \code{triangle = TRUE}.
+#' @export
+#' @seealso \code{\link{c_weighted_jaccard_sparse}} for the sparse 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_dense(m)
+#' }
+c_weighted_jaccard_dense <- function(x, transpose = FALSE, threads = NULL,
+                                     triangle = FALSE, distance = FALSE) {
+  threads <- natcpp_threads(threads)
+  weighted_jaccard_dense_impl(x, transpose = transpose, threads = threads,
+                              triangle = triangle, distance = distance)
+}
diff --git a/man/c_weighted_jaccard_dense.Rd b/man/c_weighted_jaccard_dense.Rd
index e1b1f9c..ca2121d 100644
--- a/man/c_weighted_jaccard_dense.Rd
+++ b/man/c_weighted_jaccard_dense.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
+% Please edit documentation in R/weighted_jaccard.R
 \name{c_weighted_jaccard_dense}
 \alias{c_weighted_jaccard_dense}
 \title{Dense weighted Jaccard similarity via C++}
@@ -7,7 +7,7 @@
 c_weighted_jaccard_dense(
   x,
   transpose = FALSE,
-  threads = 4L,
+  threads = NULL,
   triangle = FALSE,
   distance = FALSE
 )
@@ -18,7 +18,10 @@ c_weighted_jaccard_dense(
 \item{transpose}{If \code{FALSE}, compare columns; if \code{TRUE}, compare
 rows}
 
-\item{threads}{Number of threads (default 4). Set to 0 for all cores.}
+\item{threads}{Number of threads for parallel computation. The default
+\code{NULL} applies the package thread policy (respecting
+\code{getOption("Ncpus")} and the \code{OMP_THREAD_LIMIT} environment
+variable, else 2). Set to 0 to use all available cores.}
 
 \item{triangle}{If \code{TRUE}, return only the lower triangle as a flat
 numeric vector in \code{\link{dist}} layout. If \code{FALSE} (default),
@@ -33,8 +36,22 @@ A dense numeric similarity matrix, or a numeric vector in
 }
 \description{
 Compute the full weighted Jaccard similarity matrix for a
-\link[Matrix:dgCMatrix-class]{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.
+\link[Matrix:dgCMatrix-class]{dgCMatrix}, returning a dense matrix (or a
+\code{\link{dist}}-layout vector).
+}
+\details{
+Uses an adaptive dense accumulation strategy: for small output
+  matrices a feature-oriented loop, switching to a column-oriented loop for
+  larger outputs for better cache performance.
+}
+\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_dense(m)
+}
+}
+\seealso{
+\code{\link{c_weighted_jaccard_sparse}} for the sparse equivalent
 }
diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd
index cadadda..3f17993 100644
--- a/man/c_weighted_jaccard_sparse.Rd
+++ b/man/c_weighted_jaccard_sparse.Rd
@@ -8,7 +8,7 @@ c_weighted_jaccard_sparse(
   x,
   transpose = FALSE,
   display_progress = TRUE,
-  threads = 4L,
+  threads = NULL,
   triangle = FALSE,
   distance = FALSE
 )
@@ -22,8 +22,10 @@ c_weighted_jaccard_sparse(
 \item{display_progress}{Whether to show a text progress bar (default
 \code{TRUE}).}
 
-\item{threads}{Number of threads for parallel computation (default 4).
-Set to 0 to use all available cores.}
+\item{threads}{Number of threads for parallel computation. The default
+\code{NULL} applies the package thread policy (respecting
+\code{getOption("Ncpus")} and the \code{OMP_THREAD_LIMIT} environment
+variable, else 2). Set to 0 to use all available cores.}
 
 \item{triangle}{If \code{TRUE}, return a symmetric \code{dsCMatrix}
 (upper triangle only). If \code{FALSE} (default), return a general
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 7eb45d1..e6ff2f2 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1,8 +1,8 @@
 // Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-#include <Rcpp.h>
 #include <RcppThread.h>
+#include <Rcpp.h>
 
 using namespace Rcpp;
 
diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp
index 396b542..fe25706 100644
--- a/src/weighted_jaccard.cpp
+++ b/src/weighted_jaccard.cpp
@@ -268,27 +268,10 @@ static FeatureCSR build_feature_csr(const IntegerVector& xi, const IntegerVector
   return csr;
 }
 
-//' Dense weighted Jaccard similarity via C++
-//'
-//' Compute the full weighted Jaccard similarity matrix for a
-//' \link[Matrix:dgCMatrix-class]{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.
-//'
-//' @param x A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix)
-//' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare
-//'   rows
-//' @param threads Number of threads (default 4). Set to 0 for all cores.
-//' @param triangle If \code{TRUE}, return only the lower triangle as a flat
-//'   numeric vector in \code{\link{dist}} layout. If \code{FALSE} (default),
-//'   return a full square matrix.
-//' @param distance If \code{TRUE}, return distance (\code{1 - similarity})
-//'   instead of similarity. Default \code{FALSE}.
-//' @return A dense numeric similarity matrix, or a numeric vector in
-//'   \code{\link{dist}} layout when \code{triangle = TRUE}.
-//' @export
-// [[Rcpp::export]]
+// Internal back end for the dense weighted Jaccard similarity. The exported,
+// documented R interface c_weighted_jaccard_dense() (R/weighted_jaccard.R)
+// resolves the thread count and calls this.
+// [[Rcpp::export(weighted_jaccard_dense_impl)]]
 SEXP c_weighted_jaccard_dense(const S4& x, bool transpose = false,
                               int threads = 4, bool triangle = false,
                               bool distance = false) {
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
new file mode 100644
index 0000000..b441c07
--- /dev/null
+++ b/tests/testthat/test-utils.R
@@ -0,0 +1,26 @@
+test_that("natcpp_threads honours explicit values", {
+  expect_identical(natcpp_threads(1L), 1L)
+  expect_identical(natcpp_threads(8), 8L)
+  # 0 (all cores) is passed through unchanged, not treated as "unset"
+  expect_identical(natcpp_threads(0L), 0L)
+})
+
+test_that("natcpp_threads applies the default policy", {
+  old_opt <- options(Ncpus = 3L)
+  on.exit(options(old_opt), add = TRUE)
+  expect_identical(natcpp_threads(NULL), 3L)
+
+  options(Ncpus = NULL)
+  old_env <- Sys.getenv("OMP_THREAD_LIMIT", unset = NA)
+  on.exit(
+    if (is.na(old_env)) Sys.unsetenv("OMP_THREAD_LIMIT")
+    else Sys.setenv(OMP_THREAD_LIMIT = old_env),
+    add = TRUE
+  )
+
+  Sys.setenv(OMP_THREAD_LIMIT = "5")
+  expect_identical(natcpp_threads(NULL), 5L)
+
+  Sys.unsetenv("OMP_THREAD_LIMIT")
+  expect_identical(natcpp_threads(NULL), 2L)
+})