diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
index 7a821b8..c7c7745 100644
--- a/.github/workflows/pkgdown.yaml
+++ b/.github/workflows/pkgdown.yaml
@@ -10,7 +10,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v7
 
       - uses: r-lib/actions/setup-pandoc@v2
 
diff --git a/.github/workflows/rcmdcheck.yaml b/.github/workflows/rcmdcheck.yaml
index 6bb6ec9..932cd3d 100644
--- a/.github/workflows/rcmdcheck.yaml
+++ b/.github/workflows/rcmdcheck.yaml
@@ -28,7 +28,7 @@ jobs:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
 
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v7
 
     - uses: r-lib/actions/setup-r@v2
       with:
diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index 9f6c700..8a0dd92 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -15,7 +15,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v7
 
       - uses: r-lib/actions/setup-r@v2
 
@@ -34,7 +34,7 @@ jobs:
           covr::to_cobertura(cov)
         shell: Rscript {0}
 
-      - uses: codecov/codecov-action@v5
+      - uses: codecov/codecov-action@v7
         with:
           # Fail if error if not on PR, or if on PR and token is given
           fail_ci_if_error: ${{ github.event_name != 'pull_request' || secrets.CODECOV_TOKEN }}
diff --git a/NEWS.md b/NEWS.md
index c477895..7490e25 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,10 @@
 - Accessibility functions now treat `NA` opportunity, demand and supply
   values, and `travel_matrix` ids missing from `land_use_data`, as 0.
   Previously, they could make results `NA`.
+- Integer opportunity, demand, supply and population columns are now
+  converted to doubles, so sums no longer overflow. As a result,
+  `cumulative_cutoff()` and `gravity()` (with `decay_binary()`) return doubles
+  instead of integers when opportunities are integers.
 - `concentration_index(type = "corrected")` now uses the bounds of the whole
   accessibility distribution, not of each `group_by` group, to calculate the
   correction factor.
@@ -40,6 +44,8 @@
     `fill_missing_ids = TRUE` and `group_by` leave combinations missing.
   - `constraint = "total"` no longer adds an index to the user's
     `land_use_data`.
+- `cost_to_closest()` no longer returns `Inf` due to integer overflow with
+  large opportunity counts.
 - `cumulative_cutoff()`: cutoffs with floating point noise (e.g. from
   `seq(0, 0.5, by = 0.1)`) could be zeroed when filling missing ids. Travel
   cost columns with non-syntactic names now work.
@@ -77,9 +83,7 @@
   strings. With `constraint = "doubly"`, it is about 2 to 4 times faster and
   uses much less memory.
 - `cost_to_closest()` is about 1.3 to 2.3 times faster.
-- `cumulative_cutoff()` and `cumulative_interval()` return accessibility
-  above the largest integer as doubles instead of erroring.
-  `cumulative_interval()` always sorts rows by interval, id and `group_by`.
+- `cumulative_interval()` always sorts rows by interval, id and `group_by`.
 - `decay_exponential()`, `decay_logistic()`, `decay_power()` and
   `decay_stepped()` are about 2 to 200 times faster.
 - `decay_linear()` requires `cutoff` greater than 0.
@@ -89,9 +93,7 @@
   `poor_below_threshold`, and warns when `poor_below_threshold = FALSE` and
   `accessibility_data` has `Inf` values (e.g. from `cost_to_closest()`).
 - `floating_catchment_area()` and `gravity()` are faster with several decay
-  values (about 4 times with three). `gravity()` with integer opportunities and
-  `decay_binary()` no longer fails when accessibility exceeds the largest
-  integer.
+  values (about 4 times with three).
 - `spatial_availability()` is about 2 times faster.
 
 ## Notes
diff --git a/R/concentration_index.R b/R/concentration_index.R
index d7e4361..94718f0 100644
--- a/R/concentration_index.R
+++ b/R/concentration_index.R
@@ -106,11 +106,15 @@ concentration_index <- function(
 
   warn_extra_access_cols(accessibility_data, opportunity, group_by)
 
-  # integer opportunities and population counts may overflow when multiplied,
-  # so we convert them to doubles first
+  # integer opportunities may overflow when multiplied by the population, so
+  # they are converted to doubles first (merge_by_reference() already converts
+  # the population)
 
-  .numeric_cols <- c(opportunity, population)
-  data[, (.numeric_cols) := lapply(.SD, as.numeric), .SDcols = .numeric_cols]
+  data.table::set(
+    data,
+    j = opportunity,
+    value = as.numeric(data[[opportunity]])
+  )
 
   .opp_colname <- opportunity
   .pop_colname <- population
diff --git a/R/constrained_accessibility.R b/R/constrained_accessibility.R
index da90c19..62ed461 100644
--- a/R/constrained_accessibility.R
+++ b/R/constrained_accessibility.R
@@ -158,7 +158,9 @@
 #' to match; otherwise, the model will not converge.
 #'  - `active` must be `NULL`. Since supply must match demand, their units are
 #'  the same and there is no distinction between 'active' and 'passive' notions.
-#' - Only accepts `detailed_results = TRUE`.
+#' - Should be used with `detailed_results = TRUE`. With
+#'   `detailed_results = FALSE`, it only returns the demand and supply
+#'   marginals, with a warning.
 #'
 #' ### Interpretation:
 #' - Results include OD-level flows (`flow`) along with balancing factors (`A_i`,
@@ -171,12 +173,40 @@
 #' "What is the count of A->1 flows given A->1's region-relative travel impedance,
 #'  demand and supply?"
 #'
+#' @return A data frame of the same class as `travel_matrix` (`data.table`,
+#'   `data.frame` or `tibble`). Its columns depend on `constraint`, `active`
+#'   and `detailed_results`:
+#'
+#' - `constraint = "total"` or `"singly"`, `detailed_results = FALSE`: one
+#'   row per origin (`from_id`, when `active = TRUE`) or destination (`to_id`,
+#'   when `active = FALSE`), plus the `group_by` columns. Accessibility is in a
+#'   column named `supply` (active) or `demand` (passive).
+#' - `constraint = "total"`, `detailed_results = TRUE`: one row per
+#'   origin-destination pair, with the columns `from_id`, `to_id`, the weighted
+#'   opportunities (`weighted_supply` or `weighted_demand`), the share of the
+#'   regional total allocated to the pair (`kappa_total` or `hatkappa_total`),
+#'   the allocated flow (`supply` or `demand`) and the scaling constant
+#'   (`K_total` or `hatK_total`).
+#' - `constraint = "singly"`, `detailed_results = TRUE`: one row per
+#'   origin-destination pair, with the columns `from_id`, `to_id`, the
+#'   allocation share (`kappa_singly` or `hatkappa_singly`), the allocated flow
+#'   (`supply` or `demand`) and the balancing factor (`B_j` or `A_i`).
+#' - `constraint = "doubly"`, `detailed_results = TRUE`: one row per
+#'   origin-destination pair, with the columns `from_id`, `to_id`, `flow`, the
+#'   balancing factors `A_i` and `B_j`, `kappa_doubly` and the final balancing
+#'   `error`.
+#' - `constraint = "doubly"`, `detailed_results = FALSE`: the columns `id`
+#'   and `flow`, holding the demand of each origin followed by the supply of
+#'   each destination (so an id can appear twice).
+#'
+#' Detailed results don't include the `group_by` columns.
+#'
 #' @references
 #' \insertAllCited{}
-
+#'
 #' @family Constrained accessibility
 #'
-#' @examples
+#' @examplesIf identical(tolower(Sys.getenv("NOT_CRAN")), "true")
 #' # Load demo data shipped with the package
 #' data_dir <- system.file("extdata", package = "accessibility")
 #' travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds"))
diff --git a/R/cumulative_cutoff.R b/R/cumulative_cutoff.R
index 3a6c8e2..f32a65f 100644
--- a/R/cumulative_cutoff.R
+++ b/R/cumulative_cutoff.R
@@ -230,8 +230,7 @@ cumulative_by_cutoff <- function(data,
     binned_sums,
     grid_groups,
     groups,
-    n_cutoffs = length(cutoffs),
-    integer_opportunity = is.integer(data[[opportunity]])
+    n_cutoffs = length(cutoffs)
   )
 
   if (!fill_missing_ids) access <- access[.n > 0]
@@ -284,11 +283,9 @@ cumulative_by_filtering <- function(data,
           cum_opps <- cum_opps[cum_opps[, .I[1], by = .od_group]$V1]
         }
 
-        # opportunities are summed as doubles, so integer sums don't overflow
-
         cum_opps <- cum_opps[
           ,
-          .(access = sum(as.numeric(get(.opportunity_colname)))),
+          .(access = sum(get(.opportunity_colname))),
           by = eval(groups, envir = env)
         ]
 
@@ -303,11 +300,6 @@ cumulative_by_filtering <- function(data,
 
   access <- data.table::rbindlist(access)
 
-  is_integer_opportunity <- is.integer(data[[opportunity]])
-  if (is_integer_opportunity && all(access$access <= .Machine$integer.max)) {
-    data.table::set(access, j = "access", value = as.integer(access$access))
-  }
-
   if (fill_missing_ids) {
     unique_values <- lapply(groups, function(x) unique(travel_matrix[[x]]))
     names(unique_values) <- groups
diff --git a/R/cumulative_interval.R b/R/cumulative_interval.R
index dd3daa6..0a8bb99 100644
--- a/R/cumulative_interval.R
+++ b/R/cumulative_interval.R
@@ -104,12 +104,17 @@ cumulative_interval <- function(travel_matrix,
 
   # each interval is summarised over its own cutoffs only, which always include
   # the interval's upper limit (even if it's not a multiple of the increment).
+  # with fractional increments, the last value of the sequence may be the upper
+  # limit with floating point noise (e.g. 0.8999999999999999 instead of 0.9),
+  # so values that close to the limit are dropped before adding the limit.
   # accessibility is calculated once for the union of all cutoffs
 
   interval_cutoffs <- lapply(
     interval,
     function(pair) {
-      unique(c(seq.int(pair[1], pair[2], by = interval_increment), pair[2]))
+      steps <- seq.int(pair[1], pair[2], by = interval_increment)
+      steps <- steps[pair[2] - steps > interval_increment * 1e-6]
+      c(steps, pair[2])
     }
   )
   cutoffs <- unique(unlist(interval_cutoffs))
@@ -162,8 +167,7 @@ cumulative_interval <- function(travel_matrix,
     binned_sums,
     grid_groups,
     groups,
-    n_cutoffs = length(cutoffs),
-    integer_opportunity = is.integer(data[[opportunity]])
+    n_cutoffs = length(cutoffs)
   )
   data.table::set(access_list, j = "cutoffs", value = cutoffs[access_list$.bin])
 
diff --git a/R/gini_index.R b/R/gini_index.R
index adcd3cc..b3de88c 100644
--- a/R/gini_index.R
+++ b/R/gini_index.R
@@ -71,11 +71,15 @@ gini_index <- function(
 
   warn_extra_access_cols(accessibility_data, opportunity, group_by)
 
-  # integer opportunities and population counts may overflow when multiplied,
-  # so we convert them to doubles first
+  # integer opportunities may overflow when multiplied by the population, so
+  # they are converted to doubles first (merge_by_reference() already converts
+  # the population)
 
-  .numeric_cols <- c(opportunity, population)
-  data[, (.numeric_cols) := lapply(.SD, as.numeric), .SDcols = .numeric_cols]
+  data.table::set(
+    data,
+    j = opportunity,
+    value = as.numeric(data[[opportunity]])
+  )
 
   .opp_colname <- opportunity
   .pop_colname <- population
@@ -88,9 +92,7 @@ gini_index <- function(
     by = .groups
   ]
   data[,
-    .cum_total_access := cumsum(
-      as.numeric(get(.pop_colname) * get(.opp_colname))
-    ) /
+    .cum_total_access := cumsum(get(.pop_colname) * get(.opp_colname)) /
       sum(get(.pop_colname) * get(.opp_colname)),
     by = .groups
   ]
diff --git a/R/merge_by_reference.R b/R/merge_by_reference.R
index 437dfda..ee5d591 100644
--- a/R/merge_by_reference.R
+++ b/R/merge_by_reference.R
@@ -21,13 +21,20 @@ merge_by_reference <- function(data,
     match(left_ids, right_ids)
   }
 
-  values <- land_use_data[[opportunity]][match_index]
+  # integer counts (opportunities, demand, supply, population) are converted to
+  # doubles, so that sums and products in the measures don't overflow. they are
+  # converted before matching, because land_use_data is usually much shorter
+  # than the travel matrix
+
+  values <- land_use_data[[opportunity]]
+  if (is.integer(values)) values <- as.numeric(values)
+  values <- values[match_index]
 
   # in accessibility measures, NA values and ids missing from land_use_data
   # count as 0 opportunities (or 0 demand/supply). the inequality and poverty
   # functions keep the NAs
 
-  if (na_as_zero && anyNA(values)) values[is.na(values)] <- 0L
+  if (na_as_zero && anyNA(values)) values[is.na(values)] <- 0
 
   data.table::set(data, j = opportunity, value = values)
 }
diff --git a/R/palma_ratio.R b/R/palma_ratio.R
index b084d40..79bdc2a 100644
--- a/R/palma_ratio.R
+++ b/R/palma_ratio.R
@@ -86,11 +86,13 @@ palma_ratio <- function(
   warn_extra_access_cols(accessibility_data, opportunity, group_by)
 
   # normwt = TRUE treats population as relative weights, so the quantiles don't
-  # depend on its scale (e.g. counts vs. shares)
+  # depend on its scale (e.g. counts vs. shares). integer weights are converted
+  # to doubles, because Hmisc multiplies them by the number of observations,
+  # which may overflow
 
   income_quantiles <- Hmisc::wtd.quantile(
     sociodemographic_data[[income]],
-    weights = sociodemographic_data[[population]],
+    weights = as.numeric(sociodemographic_data[[population]]),
     probs = c(0.4, 0.9),
     normwt = TRUE
   )
diff --git a/R/sum_by_cost_bin.R b/R/sum_by_cost_bin.R
index 8216bd0..963c0e9 100644
--- a/R/sum_by_cost_bin.R
+++ b/R/sum_by_cost_bin.R
@@ -18,9 +18,10 @@ sum_by_cost_bin <- function(data, groups, opportunity, travel_cost, cutoffs) {
   names(binned) <- groups
   binned$.bin <- bin[keep]
 
-  # opportunities are summed as doubles, so integer sums don't overflow
+  # the opportunities are doubles (merge_by_reference() converts integers), so
+  # their sums don't overflow
 
-  binned$.opp <- as.numeric(data[[opportunity]][keep])
+  binned$.opp <- data[[opportunity]][keep]
   binned <- data.table::setDT(binned)
 
   binned_sums <- binned[
@@ -37,12 +38,9 @@ sum_by_cost_bin <- function(data, groups, opportunity, travel_cost, cutoffs) {
 cumulative_sums_grid <- function(binned_sums,
                                  grid_groups,
                                  groups,
-                                 n_cutoffs,
-                                 integer_opportunity) {
+                                 n_cutoffs) {
   # one row per group (from grid_groups) and cutoff. bins that no row falls
-  # into get 0 opportunities and 0 rows. the cumulative sums are doubles, so
-  # they don't overflow, and are converted back to integers when the
-  # opportunities are integers and the sums fit
+  # into get 0 opportunities and 0 rows
 
   grid <- grid_groups[rep(seq_len(nrow(grid_groups)), each = n_cutoffs)]
   data.table::set(
@@ -69,9 +67,5 @@ cumulative_sums_grid <- function(binned_sums,
   grid[, `:=`(access = cumsum(.opp), .n = cumsum(.n)), by = c(groups)]
   grid[, .opp := NULL]
 
-  if (integer_opportunity && all(grid$access <= .Machine$integer.max)) {
-    data.table::set(grid, j = "access", value = as.integer(grid$access))
-  }
-
   return(grid[])
 }
diff --git a/R/theil_t.R b/R/theil_t.R
index 1c4f4f0..cac7f2e 100644
--- a/R/theil_t.R
+++ b/R/theil_t.R
@@ -120,11 +120,15 @@ theil_t <- function(
 
   warn_extra_access_cols(accessibility_data, opportunity, group_by)
 
-  # integer opportunities and population counts may overflow when multiplied,
-  # so we convert them to doubles first
+  # integer opportunities may overflow when multiplied by the population, so
+  # they are converted to doubles first (merge_by_reference() already converts
+  # the population)
 
-  .numeric_cols <- c(opportunity, population)
-  data[, (.numeric_cols) := lapply(.SD, as.numeric), .SDcols = .numeric_cols]
+  data.table::set(
+    data,
+    j = opportunity,
+    value = as.numeric(data[[opportunity]])
+  )
 
   if (is.null(socioeconomic_groups)) {
     result <- theil_without_groups(data, opportunity, population, group_by)
@@ -212,9 +216,7 @@ theil_with_groups <- function(
 
   summarized_data <- data[,
     .(
-      group_tot_access = as.numeric(
-        sum(get(..opportunity) * get(..population))
-      ),
+      group_tot_access = sum(get(..opportunity) * get(..population)),
       group_tot_pop = sum(get(..population)),
       group_theil_t = calc_theil_t(
         get(..opportunity),
diff --git a/man/constrained_accessibility.Rd b/man/constrained_accessibility.Rd
index af2ddb6..1574224 100644
--- a/man/constrained_accessibility.Rd
+++ b/man/constrained_accessibility.Rd
@@ -94,6 +94,36 @@ performance penalty.}
 
 \item{detailed_results}{Logical. Whether to return detailed OD-level results.}
 }
+\value{
+A data frame of the same class as \code{travel_matrix} (\code{data.table},
+\code{data.frame} or \code{tibble}). Its columns depend on \code{constraint}, \code{active}
+and \code{detailed_results}:
+\itemize{
+\item \code{constraint = "total"} or \code{"singly"}, \code{detailed_results = FALSE}: one
+row per origin (\code{from_id}, when \code{active = TRUE}) or destination (\code{to_id},
+when \code{active = FALSE}), plus the \code{group_by} columns. Accessibility is in a
+column named \code{supply} (active) or \code{demand} (passive).
+\item \code{constraint = "total"}, \code{detailed_results = TRUE}: one row per
+origin-destination pair, with the columns \code{from_id}, \code{to_id}, the weighted
+opportunities (\code{weighted_supply} or \code{weighted_demand}), the share of the
+regional total allocated to the pair (\code{kappa_total} or \code{hatkappa_total}),
+the allocated flow (\code{supply} or \code{demand}) and the scaling constant
+(\code{K_total} or \code{hatK_total}).
+\item \code{constraint = "singly"}, \code{detailed_results = TRUE}: one row per
+origin-destination pair, with the columns \code{from_id}, \code{to_id}, the
+allocation share (\code{kappa_singly} or \code{hatkappa_singly}), the allocated flow
+(\code{supply} or \code{demand}) and the balancing factor (\code{B_j} or \code{A_i}).
+\item \code{constraint = "doubly"}, \code{detailed_results = TRUE}: one row per
+origin-destination pair, with the columns \code{from_id}, \code{to_id}, \code{flow}, the
+balancing factors \code{A_i} and \code{B_j}, \code{kappa_doubly} and the final balancing
+\code{error}.
+\item \code{constraint = "doubly"}, \code{detailed_results = FALSE}: the columns \code{id}
+and \code{flow}, holding the demand of each origin followed by the supply of
+each destination (so an id can appear twice).
+}
+
+Detailed results don't include the \code{group_by} columns.
+}
 \description{
 Calculates accessibility using constraints, as proposed in
 \insertCite{soukhov2025family;textual}{accessibility}. Accessibility is
@@ -245,7 +275,9 @@ travel  impedance.
 to match; otherwise, the model will not converge.
 \item \code{active} must be \code{NULL}. Since supply must match demand, their units are
 the same and there is no distinction between 'active' and 'passive' notions.
-\item Only accepts \code{detailed_results = TRUE}.
+\item Should be used with \code{detailed_results = TRUE}. With
+\code{detailed_results = FALSE}, it only returns the demand and supply
+marginals, with a warning.
 }
 }
 
@@ -270,6 +302,7 @@ demand and supply?"
 }
 
 \examples{
+\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf}
 # Load demo data shipped with the package
 data_dir <- system.file("extdata", package = "accessibility")
 travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds"))
@@ -341,7 +374,7 @@ constrained_accessibility(
   supply          = "jobs",
   detailed_results = TRUE
 )
-
+\dontshow{\}) # examplesIf}
 }
 \references{
 \insertAllCited{}
diff --git a/tests/testthat/test-cost_to_closest.R b/tests/testthat/test-cost_to_closest.R
index bf3e355..2f87ccf 100644
--- a/tests/testthat/test-cost_to_closest.R
+++ b/tests/testthat/test-cost_to_closest.R
@@ -331,3 +331,27 @@ test_that("works even if travel_matrix and land_use has specific colnames", {
 
   land_use_data[, opportunity := NULL]
 })
+
+test_that("large integer opportunity counts don't overflow", {
+  # the cumulative opportunities at costs 1, 2 and 3 are 1e9, 2e9 and 3e9, so
+  # 2.5e9 opportunities are first reached at cost 3
+
+  travel_matrix <- data.table::data.table(
+    from_id = "a",
+    to_id = c("a", "b", "c"),
+    travel_time = c(1, 2, 3)
+  )
+  land_use_data <- data.table::data.table(
+    id = c("a", "b", "c"),
+    jobs = rep(1000000000L, 3)
+  )
+
+  result <- cost_to_closest(
+    travel_matrix,
+    land_use_data,
+    opportunity = "jobs",
+    travel_cost = "travel_time",
+    n = 2.5e9
+  )
+  expect_identical(result$travel_time, 3)
+})
diff --git a/tests/testthat/test-cumulative_cutoff.R b/tests/testthat/test-cumulative_cutoff.R
index dc86e6f..b757263 100644
--- a/tests/testthat/test-cumulative_cutoff.R
+++ b/tests/testthat/test-cumulative_cutoff.R
@@ -118,25 +118,25 @@ test_that("result has correct structure", {
   expect_true(ncol(result) == 3)
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   result <- tester(opportunity = "schools")
   expect_true(ncol(result) == 3)
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
-  expect_is(result$schools, "integer")
+  expect_is(result$schools, "numeric")
 
   result <- tester(cutoff = c(15, 30))
   expect_true(ncol(result) == 4)
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
   expect_is(result$travel_time, "numeric")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   suppressWarnings(result <- tester(group_by = character(0)))
   expect_true(ncol(result) == 2)
   expect_is(result$id, "character")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   result <- tester(
     data.table::data.table(
@@ -150,7 +150,7 @@ test_that("result has correct structure", {
   expect_true(nrow(result) == 0)
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   result <- tester(
     data.table::data.table(
@@ -166,7 +166,7 @@ test_that("result has correct structure", {
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
   expect_is(result$travel_time, "numeric")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   result <- tester_with_cost()
   expect_true(ncol(result) == 5)
@@ -174,14 +174,14 @@ test_that("result has correct structure", {
   expect_is(result$mode, "character")
   expect_is(result$travel_time, "numeric")
   expect_is(result$monetary_cost, "numeric")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   suppressWarnings(result <- tester_with_cost(group_by = character()))
   expect_true(ncol(result) == 4)
   expect_is(result$id, "character")
   expect_is(result$travel_time, "numeric")
   expect_is(result$monetary_cost, "numeric")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   result <- tester_with_cost(small_frontier[0])
   expect_true(ncol(result) == 5)
@@ -189,7 +189,7 @@ test_that("result has correct structure", {
   expect_is(result$mode, "character")
   expect_is(result$travel_time, "numeric")
   expect_is(result$monetary_cost, "numeric")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 })
 
 test_that("input data sets remain unchanged", {
@@ -221,7 +221,7 @@ test_that("active and passive accessibility is correctly calculated", {
   expected_result <- data.table::data.table(
     id = rep(selected_ids, each = 2),
     mode = rep(c("transit", "transit2"), 5),
-    jobs = rep(as.integer(c(0, 82, 408, 408, 109)), each = 2)
+    jobs = rep(as.numeric(c(0, 82, 408, 408, 109)), each = 2)
   )
   expect_identical(result, expected_result)
 
@@ -234,7 +234,7 @@ test_that("active and passive accessibility is correctly calculated", {
   expected_result <- data.table::data.table(
     id = rep(selected_ids, each = 2),
     mode = rep(c("transit", "transit2"), 5),
-    population = rep(as.integer(c(0, 5404, 4552, 2363, 4552)), each = 2)
+    population = rep(as.numeric(c(0, 5404, 4552, 2363, 4552)), each = 2)
   )
   expect_identical(result, expected_result)
 
@@ -252,7 +252,7 @@ test_that("active and passive accessibility is correctly calculated", {
     mode = rep(rep(c("transit", "transit2"), each = 2), times = 5),
     travel_time = 120,
     monetary_cost = rep(c(5, 15), 10),
-    jobs = rep(as.integer(c(499, 599, 499, 599, rep(599, 16))))
+    jobs = rep(as.numeric(c(499, 599, 499, 599, rep(599, 16))))
   )
   expect_identical(result, expected_result)
 
@@ -267,7 +267,7 @@ test_that("active and passive accessibility is correctly calculated", {
     mode = rep(rep(c("transit", "transit2"), each = 2), times = 5),
     travel_time = 120,
     monetary_cost = rep(c(5, 15), 10),
-    population = rep(as.integer(c(3435, 5404, 3435, 5404, rep(5404, 16))))
+    population = rep(as.numeric(c(3435, 5404, 3435, 5404, rep(5404, 16))))
   )
   expect_identical(result, expected_result)
 })
@@ -289,7 +289,7 @@ test_that("fill_missing_ids arg works correctly", {
     data.table::data.table(
       id = rep(c("89a88cdb57bffff", "89a88cdb597ffff"), each = 2),
       mode = rep(c("transit", "transit2"), times = 2),
-      jobs = rep(as.integer(c(0, 36567)), each = 2)
+      jobs = rep(as.numeric(c(0, 36567)), each = 2)
     )
   )
 
@@ -304,7 +304,7 @@ test_that("fill_missing_ids arg works correctly", {
     data.table::data.table(
       id = rep("89a88cdb597ffff", each = 2),
       mode = c("transit", "transit2"),
-      jobs = rep(36567L, each = 2)
+      jobs = rep(36567, each = 2)
     )
   )
 
@@ -316,7 +316,7 @@ test_that("fill_missing_ids arg works correctly", {
     id = rep(c("89a88cdb57bffff", "89a88cdb597ffff"), each = 4),
     mode = rep(rep(c("transit", "transit2"), each = 2), times = 2),
     travel_time = rep(c(15, 50), times = 4),
-    jobs = as.integer(c(rep(c(0, 187799), 2), rep(c(3008, 257648), 2)))
+    jobs = as.numeric(c(rep(c(0, 187799), 2), rep(c(3008, 257648), 2)))
   )
 
   expect_identical(result, expected_result)
@@ -342,7 +342,7 @@ test_that("fill_missing_ids arg works correctly", {
     mode = c("transit", "transit2"),
     travel_time = 10,
     monetary_cost = 10,
-    jobs = as.integer(c(323, 0))
+    jobs = as.numeric(c(323, 0))
   )
   data.table::setkey(result, NULL)
   expect_identical(result, expected_result)
diff --git a/tests/testthat/test-cumulative_interval.R b/tests/testthat/test-cumulative_interval.R
index 44f5896..9cbe729 100644
--- a/tests/testthat/test-cumulative_interval.R
+++ b/tests/testthat/test-cumulative_interval.R
@@ -485,6 +485,38 @@ test_that("float cutoffs are all used in the summary", {
   expect_identical(result$jobs, 5)
 })
 
+test_that("the upper limit is counted once when the sequence has noise", {
+  # seq(0, 0.9, 0.3) ends at 0.8999999999999999, not 0.9. the cutoffs are 0,
+  # 0.3, 0.6 and 0.9, and a reaches 1, 1, 1 and 11 jobs, so the mean is 3.5
+
+  toy_matrix <- data.table::data.table(
+    from_id = "a",
+    to_id = c("a", "b"),
+    travel_time = c(0, 0.8)
+  )
+  toy_land_use <- data.table::data.table(id = c("a", "b"), jobs = c(1, 10))
+
+  result <- tester(
+    toy_matrix,
+    toy_land_use,
+    interval = c(0, 0.9),
+    interval_increment = 0.3,
+    summary_function = length,
+    group_by = character(0)
+  )
+  expect_identical(result$jobs, 4)
+
+  result <- tester(
+    toy_matrix,
+    toy_land_use,
+    interval = c(0, 0.9),
+    interval_increment = 0.3,
+    summary_function = mean,
+    group_by = character(0)
+  )
+  expect_identical(result$jobs, 3.5)
+})
+
 test_that("works when nothing is reached within the interval", {
   toy_matrix <- data.table::data.table(
     from_id = c("a", "b"),
@@ -535,3 +567,28 @@ test_that("the upper limit of the interval is always considered", {
   )
   expect_equal(result$jobs, c(32.5, 27.5, 32.5))
 })
+
+test_that("large integer sums don't overflow", {
+  # a reaches 10 jobs within 10 minutes and 10 + the largest integer within 20
+
+  travel_matrix <- data.table::data.table(
+    from_id = "a",
+    to_id = c("a", "b"),
+    travel_time = c(10, 20)
+  )
+  land_use_data <- data.table::data.table(
+    id = c("a", "b"),
+    jobs = c(10L, .Machine$integer.max)
+  )
+
+  result <- cumulative_interval(
+    travel_matrix,
+    land_use_data,
+    opportunity = "jobs",
+    travel_cost = "travel_time",
+    interval = c(10, 20),
+    interval_increment = 10,
+    summary_function = max
+  )
+  expect_identical(result$jobs, .Machine$integer.max + 10)
+})
diff --git a/tests/testthat/test-gravity.R b/tests/testthat/test-gravity.R
index 13aba05..8df7490 100644
--- a/tests/testthat/test-gravity.R
+++ b/tests/testthat/test-gravity.R
@@ -100,18 +100,18 @@ test_that("result has correct structure", {
   expect_true(ncol(result) == 3)
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   result <- tester(opportunity = "schools")
   expect_true(ncol(result) == 3)
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
-  expect_is(result$schools, "integer")
+  expect_is(result$schools, "numeric")
 
   suppressWarnings(result <- tester(group_by = character(0)))
   expect_true(ncol(result) == 2)
   expect_is(result$id, "character")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 
   result <- tester(
     data.table::data.table(
@@ -125,7 +125,7 @@ test_that("result has correct structure", {
   expect_true(nrow(result) == 0)
   expect_is(result$id, "character")
   expect_is(result$mode, "character")
-  expect_is(result$jobs, "integer")
+  expect_is(result$jobs, "numeric")
 })
 
 test_that("input data sets remain unchanged", {
@@ -206,7 +206,7 @@ test_that("accepts custom decay function", {
   ]
   expected_result <- expected_result[
     ,
-    .(jobs = sum(jobs)),
+    .(jobs = sum(as.numeric(jobs))),
     by = .(id = from_id, mode)
   ]
 
@@ -232,7 +232,7 @@ test_that("fill_missing_ids arg works correctly", {
     data.table::data.table(
       id = rep(c("89a88cdb57bffff", "89a88cdb597ffff"), each = 2),
       mode = rep(c("transit", "transit2"), times = 2),
-      jobs = c(22239L, 0L, 36567L, 36567L)
+      jobs = c(22239, 0, 36567, 36567)
     )
   )
 
@@ -246,7 +246,7 @@ test_that("fill_missing_ids arg works correctly", {
     data.table::data.table(
       id = c("89a88cdb57bffff", "89a88cdb597ffff", "89a88cdb597ffff"),
       mode = c("transit", "transit", "transit2"),
-      jobs = c(22239L, 36567L, 36567L)
+      jobs = c(22239, 36567, 36567)
     )
   )
 })
@@ -337,3 +337,26 @@ test_that("raises an error if travel_matrix contains negative travel costs", {
   na_matrix[1, travel_time := NA]
   expect_error(tester(travel_matrix = na_matrix), NA)
 })
+
+test_that("large integer opportunity counts don't overflow or warn", {
+  travel_matrix <- data.table::data.table(
+    from_id = "a",
+    to_id = c("a", "b"),
+    travel_time = c(10, 20)
+  )
+  land_use_data <- data.table::data.table(
+    id = c("a", "b"),
+    jobs = c(1L, .Machine$integer.max)
+  )
+
+  expect_silent(
+    result <- gravity(
+      travel_matrix,
+      land_use_data,
+      opportunity = "jobs",
+      travel_cost = "travel_time",
+      decay_function = decay_binary(30)
+    )
+  )
+  expect_identical(result$jobs, .Machine$integer.max + 1)
+})
diff --git a/tests/testthat/test-merge_by_reference.R b/tests/testthat/test-merge_by_reference.R
index 3c0b3b1..fa88de8 100644
--- a/tests/testthat/test-merge_by_reference.R
+++ b/tests/testthat/test-merge_by_reference.R
@@ -30,7 +30,7 @@ test_that("matches numeric and factor ids", {
   land_use_data <- data.table::data.table(id = c(1, 2), jobs = c(5L, 7L))
 
   tester(data, land_use_data)
-  expect_identical(data$jobs, c(7L, 5L, NA))
+  expect_identical(data$jobs, c(7, 5, NA))
 
   data <- data.table::data.table(to_id = factor(c("b", "a")))
   land_use_data <- data.frame(id = c("a", "b"), jobs = c(1, 2))
@@ -48,10 +48,28 @@ test_that("na_as_zero turns NA values and missing ids into 0", {
   original_land_use <- data.table::copy(land_use_data)
 
   merge_by_reference(data, land_use_data, "jobs", "to_id", na_as_zero = TRUE)
-  expect_identical(data$jobs, c(10L, 0L, 5L, 0L))
+  expect_identical(data$jobs, c(10, 0, 5, 0))
 
   merge_by_reference(data, land_use_data, "jobs", "to_id")
-  expect_identical(data$jobs, c(10L, NA, 5L, NA))
+  expect_identical(data$jobs, c(10, NA, 5, NA))
 
   expect_identical(land_use_data, original_land_use)
 })
+
+test_that("converts integer columns to doubles and keeps other types", {
+  data <- data.table::data.table(to_id = c("a", "b"))
+  land_use_data <- data.table::data.table(
+    id = c("a", "b"),
+    jobs = c(1L, .Machine$integer.max),
+    group = c("low", "high")
+  )
+
+  merge_by_reference(data, land_use_data, "jobs", "to_id")
+  expect_identical(data$jobs, c(1, .Machine$integer.max))
+  expect_identical(sum(data$jobs), .Machine$integer.max + 1)
+
+  merge_by_reference(data, land_use_data, "group", "to_id")
+  expect_identical(data$group, c("low", "high"))
+
+  expect_is(land_use_data$jobs, "integer")
+})
diff --git a/tests/testthat/test-palma_ratio.R b/tests/testthat/test-palma_ratio.R
index dcc56bb..8a98387 100644
--- a/tests/testthat/test-palma_ratio.R
+++ b/tests/testthat/test-palma_ratio.R
@@ -337,3 +337,26 @@ test_that("only returns the group_by combinations present in the data", {
   expect_setequal(paste(result$mode, result$period), c("car am", "bus pm"))
   expect_null(data.table::key(result))
 })
+
+test_that("large integer populations don't overflow", {
+  # with equal populations, the 40% and 90% income quantiles are 4.6 and 9.1,
+  # so ids 1 to 4 are the poorest and id 10 is the wealthiest. the palma ratio
+  # is 10 / mean(1:4) = 4
+
+  accessibility_data <- data.table::data.table(
+    id = as.character(1:10),
+    jobs = as.numeric(1:10)
+  )
+  sociodemographic_data <- data.table::data.table(
+    id = as.character(1:10),
+    population = rep(1000000000L, 10),
+    income_per_capita = as.numeric(1:10)
+  )
+
+  result <- tester(
+    accessibility_data,
+    sociodemographic_data,
+    group_by = character(0)
+  )
+  expect_equal(result$palma_ratio, 4)
+})