diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index aab5bd68..c462c35a 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -15,7 +15,7 @@ permissions:
 
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - NaN-safe mode ${{ matrix.nansafe && 'enabled' || 'disabled' }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - unsafe mode ${{ matrix.unsafe && 'enabled' || 'disabled' }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -29,7 +29,7 @@ jobs:
           - ubuntu-latest
           - windows-latest
           - macOS-latest
-        nansafe:
+        unsafe:
           - true
           - false
     steps:
@@ -38,10 +38,10 @@ jobs:
         with:
           version: ${{ matrix.version }}
       - uses: julia-actions/cache@v3
-      - name: Set NaN-safe mode preference
+      - name: Set unsafe mode preference
         run: |
           open(joinpath(pwd(), "Project.toml"), "a") do io
-            println(io, "\n[preferences.ForwardDiff]\nnansafe_mode = ${{ matrix.nansafe }}\n")
+            println(io, "\n[preferences.ForwardDiff]\nunsafe_mode = ${{ matrix.unsafe }}\n")
           end
         shell: julia --color=yes {0}
       - uses: julia-actions/julia-buildpkg@v1
diff --git a/Project.toml b/Project.toml
index 396b3f7f..05c4830f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "1.4.7"
+version = "1.5.0"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"
diff --git a/docs/src/user/advanced.md b/docs/src/user/advanced.md
index 69468751..6e36a16c 100644
--- a/docs/src/user/advanced.md
+++ b/docs/src/user/advanced.md
@@ -117,41 +117,36 @@ Note that it is usually best to pick a chunk size which divides evenly into the
 dimension. Otherwise, ForwardDiff has to construct and utilize an extra "remainder" chunk to
 complete the calculation.
 
-## Fixing NaN/Inf Issues
+## NaN-Safe Mode
 
-ForwardDiff's default behavior is to return `NaN` for undefined derivatives (or otherwise
-mirror the behavior of the function in `Base`, if it would return an error). This is
-usually the correct thing to do, but in some cases can erroneously "poison" values which
-aren't sensitive to the input and thus cause ForwardDiff to incorrectly return `NaN` or
-`Inf` derivatives. For example:
+A dual number whose perturbation component is zero is unperturbed: it must not propagate
+derivative information, regardless of the derivative of the function it is passed to.
+Naively propagating derivatives by multiplying them with the perturbation component
+violates this property whenever the derivative is not finite, since `0.0 * Inf` and
+`0.0 * NaN` are `NaN`. ForwardDiff therefore checks whether each perturbation component is
+zero before propagating derivative information ("`NaN`-safe mode"):
 
 ```julia-repl
-# the dual number's perturbation component is zero, so this
-# variable should not propagate derivative information
 julia> log(ForwardDiff.Dual{:tag}(0.0, 0.0))
-Dual{:tag}(-Inf,NaN) # oops, this NaN should be 0.0
-```
-
-Here, ForwardDiff computes the derivative of `log(0.0)` as `NaN` and then propagates
-this derivative by multiplying it by the perturbation component. Usually, ForwardDiff can
-rely on the identity `x * 0.0 == 0.0` to prevent the derivatives from propagating when
-the perturbation component is `0.0`. However, this identity doesn't hold if `isnan(y)`
-or `isinf(y)`, in which case a `NaN` derivative will be propagated instead.
+Dual{:tag}(-Inf,0.0)
 
-It is possible to fix this behavior by checking that the perturbation component is zero
-before attempting to propagate derivative information, but this check can noticeably
-decrease performance (~5%-10% on our benchmarks).
+julia> ForwardDiff.gradient(v -> sqrt(v[1]) + v[2], [0.0, 1.0])
+2-element Vector{Float64}:
+ Inf
+  1.0
+```
 
-In order to preserve performance in the majority of use cases, ForwardDiff disables this
-check by default. If your code is affected by this `NaN` behavior, you can enable
-ForwardDiff's `NaN`-safe mode by using the
+Similar to `@fastmath`, you can trade correctness for performance by disabling `NaN`-safe
+mode: derivatives may then be `NaN` even though they are well-defined (see the example
+below). Unsafe mode can make code faster, in particular code that is dominated by cheap
+arithmetic operations. To enable it, use the
 [Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl) API to set
-the `nansafe_mode` preference to true, for example via:
+the `unsafe_mode` preference to true, for example via:
 
 ```julia-repl
 julia> using ForwardDiff, Preferences
 
-julia> set_preferences!(ForwardDiff, "nansafe_mode" => true)
+julia> set_preferences!(ForwardDiff, "unsafe_mode" => true)
 ```
 
 Note that Julia has to be restarted and ForwardDiff has to be reloaded after changing
@@ -162,16 +157,27 @@ Alternatively, you can set the preference before loading ForwardDiff, for exampl
 ```julia-repl
 julia> using Preferences, UUIDs
 
-julia> set_preferences!(UUID("f6369f11-7733-5829-9624-2563aa707210"), "nansafe_mode" => true)
+julia> set_preferences!(UUID("f6369f11-7733-5829-9624-2563aa707210"), "unsafe_mode" => true)
 
 julia> using ForwardDiff
 
 julia> log(ForwardDiff.Dual{:tag}(0.0, 0.0))
-Dual{:tag}(-Inf,0.0)
+Dual{:tag}(-Inf,NaN)
+
+julia> ForwardDiff.gradient(v -> sqrt(v[1]) + v[2], [0.0, 1.0])
+2-element Vector{Float64}:
+  Inf
+ NaN
 ```
 
+!!! compat "ForwardDiff 1.5"
+    Before ForwardDiff 1.5.0, `NaN`-safe mode was disabled by default and could be enabled by
+    setting the `nansafe_mode` preference to true. The `nansafe_mode` preference is
+    deprecated: setting it emits a warning, and `nansafe_mode = false` disables `NaN`-safe
+    mode only if the `unsafe_mode` preference is not set.
+
 In the future, we plan on allowing users and downstream library authors to dynamically
-enable [`NaN`-safe mode via the `AbstractConfig`
+toggle [`NaN`-safe mode via the `AbstractConfig`
 API](https://github.com/JuliaDiff/ForwardDiff.jl/issues/181).
 
 ## Hessian of a vector-valued function
diff --git a/src/ForwardDiff.jl b/src/ForwardDiff.jl
index b16b986b..58f2b381 100644
--- a/src/ForwardDiff.jl
+++ b/src/ForwardDiff.jl
@@ -24,4 +24,18 @@ include("hessian.jl")
 
 export DiffResults
 
+function __init__()
+    if NANSAFE_MODE_PREFERENCE !== nothing
+        Base.depwarn(
+            "The `nansafe_mode` preference of ForwardDiff is deprecated. " *
+            "NaN-safe mode is enabled by default since ForwardDiff 1.5.0, and the preference to disable it is now called `unsafe_mode`. " *
+            "For now, `nansafe_mode = false` still disables NaN-safe mode unless `unsafe_mode` is set. " *
+            "Please remove the `nansafe_mode` preference with `Preferences.delete_preferences!(ForwardDiff, \"nansafe_mode\"; force = true)` and restart Julia.",
+            :__init__;
+            force = true,
+        )
+    end
+    return nothing
+end
+
 end # module
diff --git a/src/dual.jl b/src/dual.jl
index 48d5d0e3..ed8f133f 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -416,6 +416,10 @@ end
 
 Base.iszero(x::Dual) = iszero(value(x)) && iszero(partials(x))  # shortcut, equivalent to x == zero(x)
 
+@inline _isfinite_all(x::Dual) = _isfinite_all(x.value) & all(_isfinite_all, x.partials.values)
+@inline _is_strong_divisor(x::Dual) = _is_strong_divisor(x.value) & all(_isfinite_all, x.partials.values)
+@inline _is_strong_zero_type(::Type{<:Dual{<:Any,V}}) where {V} = _is_strong_zero_type(V)
+
 for pred in [:isequal, :(==)]
     @eval begin
         @define_binary_dual_op(
@@ -619,90 +623,37 @@ end
     calc_hypot(x, y, z, Tz),
 )
 
-# fma #
-#-----#
-
-@generated function calc_fma_xyz(x::Dual{T,<:Any,N},
-                                 y::Dual{T,<:Any,N},
-                                 z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_fma_partial(value(x), partials(y, $i), _fma_partial(value(y), partials(x, $i), partials(z, $i)))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = fma(value(x), value(y), value(z))
-        return Dual{T}(v, $ex)
-    end
-end
-
-@inline function calc_fma_xy(x::Dual{T}, y::Dual{T}, z::Real) where T
-    vx, vy = value(x), value(y)
-    result = fma(vx, vy, z)
-    return Dual{T}(result, _mul_partials(partials(x), partials(y), vy, vx))
-end
+# fma/muladd #
+#------------#
 
-@generated function calc_fma_xz(x::Dual{T,<:Any,N},
-                                y::Real,
-                                z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_fma_partial(y, partials(x, $i), partials(z, $i))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = fma(value(x), y, value(z))
-        Dual{T}(v, $ex)
-    end
+# `f` is `fma` or `muladd`
+@inline function calc_muladd_xyz(f::F, x::Dual{T}, y::Dual{T}, z::Dual{T}) where {F,T}
+    vx, vy = x.value, y.value
+    return Dual{T}(f(vx, vy, z.value), _muladd_partials(f, y.partials, x.partials, z.partials, vx, vy))
 end
 
-@define_ternary_dual_op(
-    Base.fma,
-    calc_fma_xyz(x, y, z),                         # xyz_body
-    calc_fma_xy(x, y, z),                          # xy_body
-    calc_fma_xz(x, y, z),                          # xz_body
-    Base.fma(y, x, z),                             # yz_body
-    Dual{Tx}(fma(value(x), y, z), partials(x) * y), # x_body
-    Base.fma(y, x, z),                              # y_body
-    Dual{Tz}(fma(x, y, value(z)), partials(z))      # z_body
-)
-
-# muladd #
-#--------#
-
-@generated function calc_muladd_xyz(x::Dual{T,<:Any,N},
-                                    y::Dual{T,<:Any,N},
-                                    z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_muladd_partial(value(x), partials(y, $i), _muladd_partial(value(y), partials(x, $i), partials(z, $i)))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = muladd(value(x), value(y), value(z))
-        return Dual{T}(v, $ex)
-    end
+@inline function calc_muladd_xy(f::F, x::Dual{T}, y::Dual{T}, z::Real) where {F,T}
+    vx, vy = x.value, y.value
+    return Dual{T}(f(vx, vy, z), _mul_partials(x.partials, y.partials, vy, vx))
 end
 
-@inline function calc_muladd_xy(x::Dual{T}, y::Dual{T}, z::Real) where T
-    vx, vy = value(x), value(y)
-    result = muladd(vx, vy, z)
-    return Dual{T}(result, _mul_partials(partials(x), partials(y), vy, vx))
+@inline function calc_muladd_xz(f::F, x::Dual{T}, y::Real, z::Dual{T}) where {F,T}
+    return Dual{T}(f(x.value, y, z.value), _muladd_partials(f, x.partials, z.partials, y))
 end
 
-@generated function calc_muladd_xz(x::Dual{T,<:Any,N},
-                                   y::Real,
-                                   z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_muladd_partial(y, partials(x, $i), partials(z, $i))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = muladd(value(x), y, value(z))
-        Dual{T}(v, $ex)
-    end
+for f in (:fma, :muladd)
+    @eval @define_ternary_dual_op(
+        Base.$f,
+        calc_muladd_xyz($f, x, y, z),                  # xyz_body
+        calc_muladd_xy($f, x, y, z),                   # xy_body
+        calc_muladd_xz($f, x, y, z),                   # xz_body
+        Base.$f(y, x, z),                              # yz_body
+        Dual{Tx}($f(x.value, y, z), x.partials * y),   # x_body
+        Base.$f(y, x, z),                              # y_body
+        Dual{Tz}($f(x, y, z.value), z.partials)        # z_body
+    )
 end
 
-@define_ternary_dual_op(
-    Base.muladd,
-    calc_muladd_xyz(x, y, z),                         # xyz_body
-    calc_muladd_xy(x, y, z),                          # xy_body
-    calc_muladd_xz(x, y, z),                          # xz_body
-    Base.muladd(y, x, z),                             # yz_body
-    Dual{Tx}(muladd(value(x), y, z), partials(x) * y), # x_body
-    Base.muladd(y, x, z),                             # y_body
-    Dual{Tz}(muladd(x, y, value(z)), partials(z))      # z_body
-)
-
 # sin/cos #
 #--------#
 
diff --git a/src/partials.jl b/src/partials.jl
index e6bf3447..4d76bd6c 100644
--- a/src/partials.jl
+++ b/src/partials.jl
@@ -84,16 +84,44 @@ Base.convert(::Type{Partials{N,V}}, partials::Partials{N,V}) where {N,V} = parti
 @inline Base.:-(partials::Partials) = Partials(minus_tuple(partials.values))
 @inline Base.:*(x::Real, partials::Partials) = partials*x
 
-@inline function Base.:*(partials::Partials, x::Real)
-    return Partials(scale_tuple(partials.values, x))
+@inline function Base.:*(partials::Partials{N,V}, x::Real) where {N,V}
+    if _mul_is_strong(V, x)
+        return Partials(scale_tuple(*, partials.values, x))
+    else
+        return Partials(scale_tuple(_strong_mul, partials.values, x))
+    end
+end
+
+@inline function Base.:/(partials::Partials{N,V}, x::Real) where {N,V}
+    if _div_is_strong(V, x)
+        return Partials(scale_tuple(/, partials.values, x))
+    else
+        return Partials(scale_tuple(_strong_div, partials.values, x))
+    end
 end
 
-@inline function Base.:/(partials::Partials, x::Real)
-    return Partials(div_tuple_by_scalar(partials.values, x))
+@inline function _mul_partials(a::Partials{N,A}, b::Partials{N,B}, x_a, x_b) where {N,A,B}
+    if _mul_is_strong(A, x_a) & _mul_is_strong(B, x_b)
+        return Partials(mul_tuples(*, a.values, b.values, x_a, x_b))
+    else
+        return Partials(mul_tuples(_strong_mul, a.values, b.values, x_a, x_b))
+    end
 end
 
-@inline function _mul_partials(a::Partials{N}, b::Partials{N}, x_a, x_b) where N
-    return Partials(mul_tuples(a.values, b.values, x_a, x_b))
+# `x_a * a + c` and `x_a * a + x_b * b + c` with `f = fma` or `f = muladd`
+@inline function _muladd_partials(f::F, a::Partials{N,A}, c::Partials{N}, x_a) where {F,N,A}
+    if _mul_is_strong(A, x_a)
+        return Partials(muladd_tuples(f, a.values, c.values, x_a))
+    else
+        return Partials(muladd_tuples((x, p, acc) -> _strong_muladd(f, x, p, acc), a.values, c.values, x_a))
+    end
+end
+@inline function _muladd_partials(f::F, a::Partials{N,A}, b::Partials{N,B}, c::Partials{N}, x_a, x_b) where {F,N,A,B}
+    if _mul_is_strong(A, x_a) & _mul_is_strong(B, x_b)
+        return Partials(muladd_tuples(f, a.values, b.values, c.values, x_a, x_b))
+    else
+        return Partials(muladd_tuples((x, p, acc) -> _strong_muladd(f, x, p, acc), a.values, b.values, c.values, x_a, x_b))
+    end
 end
 
 # Divide by `bval` instead of multiplying by `inv(bval)` or dividing by `bval^2`, which may overflow or underflow
@@ -104,35 +132,43 @@ end
 # NaN/Inf-safe methods #
 #----------------------#
 
-if NANSAFE_MODE_ENABLED
-    # A dual number with a zero partial is just an unperturbed non-dual number
-    # Hence when propagated the resulting dual number is unperturbed as well,
-    # ie., its partial is zero as well, regardless of the primal value
-    # However, standard floating point multiplication/division would return `NaN`
-    # if the primal is not-finite/zero
-    @inline function _mul_partial(partial::Real, x::Real)
-        y = partial * x
-        return iszero(partial) ? zero(y) : y
-    end
-    @inline function _div_partial(partial::Real, x::Real)
-        y = partial / x
-        return iszero(partial) ? zero(y) : y
-    end
-    @inline function _fma_partial(x::Real, partial::Real, acc::Real)
-        y = fma(x, partial, acc)
-        return iszero(partial) ? oftype(y, acc) : y
+# Zero perturbations in input yield zero perturbations in output
+@inline function _strong_mul(partial::Real, x::Real)
+    y = partial * x
+    return iszero(partial) ? zero(y) : y
+end
+@inline function _strong_div(partial::Real, x::Real)
+    y = partial / x
+    return iszero(partial) ? zero(y) : y
+end
+@inline function _strong_muladd(f::F, x::Real, partial::Real, acc::Real) where {F}
+    y = f(x, partial, acc)
+    return iszero(partial) ? oftype(y, acc) : y
+end
+
+# Whether plain `*`/`/` already treat zero partials of type `V` as strong zeros
+if !UNSAFE_MODE_ENABLED
+    @inline function _mul_is_strong(::Type{V}, x::Real) where {V}
+        return _is_strong_zero_type(V) && _is_strong_zero_type(typeof(x)) && _isfinite_all(x)
     end
-    @inline function _muladd_partial(x::Real, partial::Real, acc::Real)
-        y = muladd(x, partial, acc)
-        return iszero(partial) ? oftype(y, acc) : y
+    @inline function _div_is_strong(::Type{V}, x::Real) where {V}
+        return _is_strong_zero_type(V) && _is_strong_zero_type(typeof(x)) && _is_strong_divisor(x)
     end
 else
-    @inline _mul_partial(partial::Real, x::Real) = partial * x
-    @inline _div_partial(partial::Real, x::Real) = partial / x
-    @inline _fma_partial(x::Real, partial::Real, acc::Real) = fma(x, partial, acc)
-    @inline _muladd_partial(x::Real, partial::Real, acc::Real) = muladd(x, partial, acc)
+    @inline _mul_is_strong(::Type, ::Real) = true
+    @inline _div_is_strong(::Type, ::Real) = true
 end
 
+# Types whose zeros stay zero when multiplied by finite values of such types
+@inline _is_strong_zero_type(::Type) = false
+@inline _is_strong_zero_type(::Type{<:Union{Base.IEEEFloat,Base.BitInteger,Bool,BigInt,BigFloat,Rational{<:Union{Base.BitInteger,Bool,BigInt}},AbstractIrrational}}) = true
+
+# Whether `x` and, for dual numbers, all of its partials are finite
+@inline _isfinite_all(x::Real) = isfinite(x)
+
+# Whether zero divided by `x` is zero
+@inline _is_strong_divisor(x::Real) = !(isnan(x) | iszero(x))
+
 # edge cases where N == 0 #
 #-------------------------#
 
@@ -154,6 +190,9 @@ end
 @inline _mul_partials(a::Partials{0,A}, b::Partials{N,B}, afactor, bfactor) where {N,A,B} = bfactor * b
 @inline _mul_partials(a::Partials{N,A}, b::Partials{0,B}, afactor, bfactor) where {N,A,B} = afactor * a
 
+@inline _muladd_partials(f::F, a::Partials{0,A}, c::Partials{0,C}, x_a) where {F,A,C} = Partials{0,promote_type(A,C)}(tuple())
+@inline _muladd_partials(f::F, a::Partials{0,A}, b::Partials{0,B}, c::Partials{0,C}, x_a, x_b) where {F,A,B,C} = Partials{0,promote_type(A,B,C)}(tuple())
+
 ##################################
 # Generated Functions on NTuples #
 ##################################
@@ -210,12 +249,8 @@ end
     return tupexpr(i -> :(rand(V)), N)
 end
 
-@generated function scale_tuple(tup::NTuple{N}, x) where N
-    return tupexpr(i -> :(_mul_partial(tup[$i], x)), N)
-end
-
-@generated function div_tuple_by_scalar(tup::NTuple{N}, x) where N
-    return tupexpr(i -> :(_div_partial(tup[$i], x)), N)
+@generated function scale_tuple(f::F, tup::NTuple{N}, x) where {F,N}
+    return tupexpr(i -> :(f(tup[$i], x)), N)
 end
 
 @generated function add_tuples(a::NTuple{N}, b::NTuple{N})  where N
@@ -230,8 +265,16 @@ end
     return tupexpr(i -> :(-tup[$i]), N)
 end
 
-@generated function mul_tuples(a::NTuple{N}, b::NTuple{N}, afactor, bfactor) where N
-    return tupexpr(i -> :(_mul_partial(a[$i], afactor) + _mul_partial(b[$i], bfactor)), N)
+@generated function mul_tuples(f::F, a::NTuple{N}, b::NTuple{N}, afactor, bfactor) where {F,N}
+    return tupexpr(i -> :(f(a[$i], afactor) + f(b[$i], bfactor)), N)
+end
+
+@generated function muladd_tuples(f::F, a::NTuple{N}, c::NTuple{N}, x_a) where {F,N}
+    return tupexpr(i -> :(f(x_a, a[$i], c[$i])), N)
+end
+
+@generated function muladd_tuples(f::F, a::NTuple{N}, b::NTuple{N}, c::NTuple{N}, x_a, x_b) where {F,N}
+    return tupexpr(i -> :(f(x_a, a[$i], f(x_b, b[$i], c[$i]))), N)
 end
 
 ###################
diff --git a/src/prelude.jl b/src/prelude.jl
index 9e037afa..a64e31a8 100644
--- a/src/prelude.jl
+++ b/src/prelude.jl
@@ -1,4 +1,5 @@
-const NANSAFE_MODE_ENABLED = @load_preference("nansafe_mode", false)
+const NANSAFE_MODE_PREFERENCE = @load_preference("nansafe_mode", nothing) # deprecated
+const UNSAFE_MODE_ENABLED = @load_preference("unsafe_mode", NANSAFE_MODE_PREFERENCE === false)
 const DEFAULT_CHUNK_THRESHOLD = @load_preference("default_chunk_threshold", 12)
 
 const AMBIGUOUS_TYPES = (AbstractFloat, Irrational, Integer, Rational, Real, RoundingMode)
diff --git a/test/DerivativeTest.jl b/test/DerivativeTest.jl
index ad171524..6770fe9e 100644
--- a/test/DerivativeTest.jl
+++ b/test/DerivativeTest.jl
@@ -116,7 +116,7 @@ end
 
 @testset "NaN-safe mode" begin
     x = ForwardDiff.derivative(log ∘ zero, 1.0)
-    if ForwardDiff.NANSAFE_MODE_ENABLED
+    if !ForwardDiff.UNSAFE_MODE_ENABLED
         @test iszero(x)
     else
         @test isnan(x)
diff --git a/test/GradientTest.jl b/test/GradientTest.jl
index bf121239..a665fe1f 100644
--- a/test/GradientTest.jl
+++ b/test/GradientTest.jl
@@ -149,7 +149,7 @@ end
 end
 
 @testset "exponential function at base zero" begin
-    if ForwardDiff.NANSAFE_MODE_ENABLED
+    if !ForwardDiff.UNSAFE_MODE_ENABLED
         @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0, -0.5]), [-Inf, -Inf])
         @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0,  0.0]), [NaN, -Inf])
         @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0,  0.5]), [Inf, 0.0])
@@ -214,7 +214,7 @@ end
 end
 
 @testset "gradient for exponential with NaNMath" begin
-    if ForwardDiff.NANSAFE_MODE_ENABLED
+    if !ForwardDiff.UNSAFE_MODE_ENABLED
         @test isequal(ForwardDiff.gradient(x -> NaNMath.pow(x[1],x[2]), [NaN, 1.0]), [1.0, NaN])
     else
         @test isequal(ForwardDiff.gradient(x -> NaNMath.pow(x[1],x[2]), [NaN, 1.0]), [NaN, NaN])
@@ -222,7 +222,7 @@ end
     @test ForwardDiff.gradient(x -> NaNMath.pow(x[1], x[2]), [1.0, 1.0]) == [1.0, 0.0]
     @test isnan(ForwardDiff.gradient((x) -> NaNMath.pow(x[1], x[2]), [-1.0, 0.5])[1])
 
-    if ForwardDiff.NANSAFE_MODE_ENABLED
+    if !ForwardDiff.UNSAFE_MODE_ENABLED
         @test isequal(ForwardDiff.gradient(x -> x[1]^x[2], [NaN, 1.0]), [1.0, NaN])
     else
         @test isequal(ForwardDiff.gradient(x -> x[1]^x[2], [NaN, 1.0]), [NaN, NaN])
@@ -309,7 +309,7 @@ end
     y2 = ForwardDiff.gradient(f, x, ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{1}()))
     y3 = ForwardDiff.gradient(f, x, ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{2}()))
     for y in (y1, y2, y3)
-        if ForwardDiff.NANSAFE_MODE_ENABLED
+        if !ForwardDiff.UNSAFE_MODE_ENABLED
             @test y == [0.0, Inf]
         else
             @test isequal(y, [NaN, Inf])
@@ -323,7 +323,7 @@ end
     b2 = ForwardDiff.gradient(g, a, ForwardDiff.GradientConfig(g, a, ForwardDiff.Chunk{1}()))
     b3 = ForwardDiff.gradient(g, a, ForwardDiff.GradientConfig(g, a, ForwardDiff.Chunk{2}()))
     for b in (b1, b2, b3)
-        if ForwardDiff.NANSAFE_MODE_ENABLED
+        if !ForwardDiff.UNSAFE_MODE_ENABLED
             @test b == [Inf, -Inf]
         else
             @test isequal(b, [NaN, NaN])
diff --git a/test/PartialsTest.jl b/test/PartialsTest.jl
index bed8b48e..10e4a546 100644
--- a/test/PartialsTest.jl
+++ b/test/PartialsTest.jl
@@ -143,8 +143,8 @@ samerng() = MersenneTwister(1)
         # x / x has zero derivative
         @test iszero(ForwardDiff._div_partials(PARTIALS, PARTIALS, X, X))
 
-        if ForwardDiff.NANSAFE_MODE_ENABLED
-            for f in ((p -> NaN * p), (p -> Inf * p), (p -> -Inf * p), (p -> p / 0), (p -> p / NaN), (p -> p / Inf), (p -> p / -Inf))
+        if !ForwardDiff.UNSAFE_MODE_ENABLED
+            for f in ((p -> NaN * p), (p -> Inf * p), (p -> -Inf * p), (p -> big(Inf) * p), (p -> NaN32 * p), (p -> p / 0), (p -> p / NaN), (p -> p / Inf), (p -> p / -Inf), (p -> p / big(0.0)))
                 # Only zero partials
                 @test iszero(@inferred(f(ALLZERO)))
 
@@ -162,6 +162,40 @@ samerng() = MersenneTwister(1)
     end
 end
 
+# A number type that claims to be finite but yields `NaN` when multiplied by zero
+struct FiniteNaN <: Real end
+Base.isfinite(::FiniteNaN) = true
+Base.:*(x::Float64, ::FiniteNaN) = x * NaN
+Base.:*(::FiniteNaN, x::Float64) = NaN * x
+Base.fma(x::FiniteNaN, y::Float64, z::Float64) = x * y + z
+Base.muladd(x::FiniteNaN, y::Float64, z::Float64) = x * y + z
+
+@testset "strong zeros" begin
+    z = ForwardDiff.UNSAFE_MODE_ENABLED ? NaN : 0.0
+    @test isequal(Partials((0.0, 1.0)) * FiniteNaN(), Partials((z, NaN)))
+    @test isequal(Partials((0.0, 1.0)) * ForwardDiff.Dual(1.0, Inf), Partials((ForwardDiff.Dual(0.0, z), ForwardDiff.Dual(1.0, Inf))))
+    @test isequal(Partials((0.0, 1.0)) / ForwardDiff.Dual(1.0, Inf), Partials((ForwardDiff.Dual(0.0, z), ForwardDiff.Dual(1.0, -Inf))))
+    @test isequal(Partials((0.0, 1.0)) / ForwardDiff.Dual(0.0, 1.0), Partials((ForwardDiff.Dual(z, z), ForwardDiff.Dual(Inf, -Inf))))
+    @test Partials((0.0, 1.0)) / ForwardDiff.Dual(Inf, 1.0) == Partials((ForwardDiff.Dual(0.0, 0.0), ForwardDiff.Dual(0.0, 0.0)))
+    x = ForwardDiff.Dual{:a}(ForwardDiff.Dual{:b}(2.0, Inf), -3.0)
+    @test isequal(Partials((0.0, 1.0)) / x, Partials((ForwardDiff.Dual{:a}(ForwardDiff.Dual{:b}(0.0, z), ForwardDiff.Dual{:b}(0.0, z)), 1.0 / x)))
+    for x in (ForwardDiff.Dual{:a}(ForwardDiff.Dual{:b}(1e-200, 1.0), 1.0), ForwardDiff.Dual{:a}(ForwardDiff.Dual{:b}(1e-310, 1.0), 1.0))
+        @test iszero((Partials((0.0, 1.0)) / x)[1])
+    end
+    for x in (ForwardDiff.Dual{:b}(1e-200, 1.0), ForwardDiff.Dual{:b}(1e-310, 1.0))
+        @test iszero((Partials((ForwardDiff.Dual{:b}(0.0, 0.0), ForwardDiff.Dual{:b}(1.0, 0.0))) / x)[1])
+    end
+    @test isequal(ForwardDiff._mul_partials(Partials((0.0, 1.0)), Partials((1.0, 0.0)), ForwardDiff.Dual(1.0, Inf), 2.0), Partials((ForwardDiff.Dual(2.0, z), ForwardDiff.Dual(1.0, Inf))))
+    @test Partials((0.0, 1.0)) * π === Partials((0.0, Float64(π)))
+    @test Partials((0.0, 1.0)) / π === Partials((0.0, 1 / π))
+    a, b, c = Partials((0.0, 1.0)), Partials((1.0, 0.0)), Partials((1.0, 2.0))
+    for f in (fma, muladd), x in (2.0, Inf, NaN, ForwardDiff.Dual(1.0, Inf), FiniteNaN())
+        @test isequal(ForwardDiff._muladd_partials(f, a, c, x), x * a + c)
+        @test isequal(ForwardDiff._muladd_partials(f, a, b, c, x, 2.0), x * a + 2.0 * b + c)
+        @test isequal(ForwardDiff._muladd_partials(f, b, a, c, 2.0, x), 2.0 * b + x * a + c)
+    end
+end
+
 io = IOBuffer()
 show(io, MIME("text/plain"), Partials((1, 2, 3)))
 str = String(take!(io))
diff --git a/test/PreferencesTest.jl b/test/PreferencesTest.jl
new file mode 100644
index 00000000..b062aa31
--- /dev/null
+++ b/test/PreferencesTest.jl
@@ -0,0 +1,41 @@
+module PreferencesTest
+
+using Test
+using ForwardDiff
+
+const DEPRECATION_MESSAGE = "The `nansafe_mode` preference of ForwardDiff is deprecated. " *
+    "NaN-safe mode is enabled by default since ForwardDiff 1.5.0, and the preference to disable it is now called `unsafe_mode`. " *
+    "For now, `nansafe_mode = false` still disables NaN-safe mode unless `unsafe_mode` is set. " *
+    "Please remove the `nansafe_mode` preference with `Preferences.delete_preferences!(ForwardDiff, \"nansafe_mode\"; force = true)` and restart Julia."
+
+@testset "deprecated `nansafe_mode` preference" begin
+    # ForwardDiff is an extra so that it is loaded from the test environment
+    # `__clear__` discards preferences inherited from the test environment
+    # Julia 1.10 hashes only the values of preferences for precompile caches, so we set
+    # `default_chunk_threshold` (to its default) to avoid reusing the cache of the test process
+    for (preferences, unsafe_mode_enabled) in (
+        ("nansafe_mode = false\n__clear__ = [\"unsafe_mode\"]", true),
+        ("nansafe_mode = false\nunsafe_mode = false", false),
+        ("nansafe_mode = true\n__clear__ = [\"unsafe_mode\"]", false),
+    )
+        mktempdir() do dir
+            write(joinpath(dir, "Project.toml"), """
+                [extras]
+                ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+
+                [preferences.ForwardDiff]
+                default_chunk_threshold = 12
+                $preferences
+                """)
+            cmd = `$(Base.julia_cmd()) --startup-file=no --depwarn=no --project=$dir -e 'using ForwardDiff; print(ForwardDiff.UNSAFE_MODE_ENABLED)'`
+            load_path = join(["@", Base.load_path()...], Sys.iswindows() ? ';' : ':')
+            out = IOBuffer()
+            err = IOBuffer()
+            run(pipeline(addenv(cmd, "JULIA_LOAD_PATH" => load_path); stdout = out, stderr = err))
+            @test String(take!(out)) == string(unsafe_mode_enabled)
+            @test occursin(DEPRECATION_MESSAGE, String(take!(err)))
+        end
+    end
+end
+
+end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index e39f5e46..f86b0972 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,6 @@
 using ForwardDiff, Test, Random
 
-@info "Testing ForwardDiff with NaN-safe mode $(ForwardDiff.NANSAFE_MODE_ENABLED ? "enabled" : "disabled")"
+@info "Testing ForwardDiff with unsafe mode $(ForwardDiff.UNSAFE_MODE_ENABLED ? "enabled" : "disabled")"
 
 SEED = trunc(Int, time())
 println("##### Random.seed!($SEED), on VERSION == $VERSION")
@@ -58,6 +58,11 @@ Random.seed!(SEED)
         t = @elapsed include("AllocationsTest.jl")
         println("##### done (took $t seconds).")
     end
+    @testset "Preferences" begin
+        println("##### Testing preferences...")
+        t = @elapsed include("PreferencesTest.jl")
+        println("##### done (took $t seconds).")
+    end
     @testset "QA" begin
         println("##### QA testing...")
         t = @elapsed include("QATest.jl")