diff --git a/DESCRIPTION b/DESCRIPTION
index e465302..bc0b166 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: errors
 Type: Package
 Title: Uncertainty Propagation for R Vectors
-Version: 0.4.5
+Version: 0.4.5.1
 Authors@R: c(
     person("Iñaki", "Ucar", email="iucar@fedoraproject.org",
       role=c("aut", "cph", "cre"), comment=c(ORCID="0000-0001-6403-5550")),
diff --git a/NEWS.md b/NEWS.md
index 9848db3..197c87a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,9 @@
+# errors devel
+
+- Fix PDG rounding (`digits="pdg"`) for uncertainties whose three highest
+  order digits are 950-999, which are now rounded up to 1000 and keep two
+  significant digits (@davidchall #70 addressing #65).
+
 # errors 0.4.5
 
 - Fix `print()` method, which now correctly returns the object invisibly (#64).
diff --git a/R/print.R b/R/print.R
index e01f24a..0cb7de5 100644
--- a/R/print.R
+++ b/R/print.R
@@ -42,13 +42,21 @@ format.errors = function(x,
   stopifnot(notation %in% c("parenthesis", "plus-minus"))
 
   if (is.null(digits)) digits <- getOption("errors.digits", 1)
-  digits <- if (digits == "pdg") digits_pdg(.e(x)) else rep(digits, length(x))
+  ee <- .e(x)
+  if (identical(digits, "pdg")) {
+    hod <- pdg_hod(ee)
+    digits <- digits_pdg(hod)
+    # PDG: 950-999 round up to 1000 and keep two significant digits
+    e <- ifelse(hod >= 950, 10^(get_exponent(ee) + 1), signif(ee, digits))
+  } else {
+    digits <- rep(digits, length(x))
+    e <- signif(ee, digits)
+  }
 
   scipen <- getOption("scipen", 0)
   prepend <- rep("", length(x))
   append <- rep("", length(x))
 
-  e <- signif(.e(x), digits)
   nulle <- e == 0 & !is.na(e)
   eexp <- get_exponent(e)
   xexp <- ifelse(.v(x) == 0, eexp + 1, get_exponent(x))
diff --git a/R/utils.R b/R/utils.R
index 3e74f64..4615441 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -16,15 +16,15 @@ zstd <- function(x, op) {
 
 get_exponent <- function(x) ifelse(.v(x), floor(log10(abs(.v(x)))), 0)
 
-digits_pdg <- function(x) {
+pdg_hod <- function(x) {
   # extract 3 highest order digits
   x <- ifelse(is.finite(x), x, 0)
   x_sci <- formatC(abs(x), digits=2, format="e", decimal.mark=".")
-  x_hod <- as.integer(gsub("(\\.|e.*)", "", x_sci))
-
-  ifelse(x_hod < 355, 2, ifelse(x_hod < 950, 1, 0))
+  as.integer(gsub("(\\.|e.*)", "", x_sci))
 }
 
+digits_pdg <- function(hod) ifelse(hod < 355, 2, ifelse(hod < 950, 1, 2))
+
 propagate <- function(xx, x, y, dx, dy, method=getOption("errors.propagation", "taylor-first-order")) {
   # if y not defined, use a vector of NAs
   if (length(y) == 1 && is.na(y))
diff --git a/tests/testthat/test-print.R b/tests/testthat/test-print.R
index 6281e28..446e23a 100644
--- a/tests/testthat/test-print.R
+++ b/tests/testthat/test-print.R
@@ -47,7 +47,31 @@ test_that("error formatting works properly", {
 
   x <- set_errors(rep(0.827, 3), c(0.119, 0.367, 0.962))
   expect_equal(format(x, notation="plus-minus", digits="pdg"), sapply(list(
-    c("0.83", "0.12"), c("0.8", "0.4"), c("1", "1")),
+    c("0.83", "0.12"), c("0.8", "0.4"), c("0.8", "1.0")),
+    paste, collapse=paste("", .pm, "")))
+
+  # PDG 950-999 rounding up to 1000: boundary, magnitudes, notations
+  x <- set_errors(c(1, 1, 0.827, 827, 1e-5), c(0.95, 0.9495, 0.962, 96.2, 9.62e-6))
+  expect_equal(format(x, notation="plus-minus", digits="pdg"), c(sapply(list(
+    c("1.0", "1.0"), c("1.0", "1.0"), c("0.8", "1.0"), c("830", "100")),
+    paste, collapse=paste("", .pm, "")), paste0("(1.0 ", .pm, " 1.0)e-5")))
+  expect_equal(format(x, notation="parenthesis", digits="pdg"),
+               c("1.0(10)", "1.0(10)", "0.8(10)", "830(100)", "1.0(10)e-5"))
+  expect_equal(format(x, notation="parenthesis", digits="pdg", scientific=TRUE),
+               c("1.0(10)e0", "1.0(10)e0", "8(10)e-1", "8.3(10)e2", "1.0(10)e-5"))
+
+  # PDG band edges (three highest order digits are rounded, not truncated)
+  x <- set_errors(rep(0.827, 7), c(0.100, 0.354, 0.3549, 0.355, 0.949, 0.999, 0.9996))
+  expect_equal(format(x, notation="plus-minus", digits="pdg"), sapply(list(
+    c("0.83", "0.10"), c("0.83", "0.35"), c("0.8", "0.4"), c("0.8", "0.4"),
+    c("0.8", "0.9"), c("0.8", "1.0"), c("0.8", "1.0")),
+    paste, collapse=paste("", .pm, "")))
+
+  # PDG with zero/negative values and zero/non-finite errors
+  x <- set_errors(c(0, 0, -0.827, 0.827, 1, 1), c(0.962, 0.119, 0.962, 0, NA, Inf))
+  expect_equal(format(x, notation="plus-minus", digits="pdg"), sapply(list(
+    c("0.0", "1.0"), c("0.00", "0.12"), c("-0.8", "1.0"), c("0.827", "0"),
+    c("1", "NA"), c("1", "Inf")),
     paste, collapse=paste("", .pm, "")))
 
   x <- set_errors(10, 1)