diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 00000000..700707ce
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"
diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml
new file mode 100644
index 00000000..4d0004e8
--- /dev/null
+++ b/.github/workflows/Invalidations.yml
@@ -0,0 +1,40 @@
+name: Invalidations
+
+on:
+  pull_request:
+
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: always.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  evaluate:
+    # Only run on PRs to the default branch.
+    # In the PR trigger above branches can be specified only explicitly whereas this check should work for master, main, or any other default branch
+    if: github.base_ref == github.event.repository.default_branch
+    runs-on: ubuntu-latest
+    steps:
+    - uses: julia-actions/setup-julia@v1
+      with:
+        version: '1'
+    - uses: actions/checkout@v3
+    - uses: julia-actions/julia-buildpkg@v1
+    - uses: julia-actions/julia-invalidations@v1
+      id: invs_pr
+
+    - uses: actions/checkout@v3
+      with:
+        ref: ${{ github.event.repository.default_branch }}
+    - uses: julia-actions/julia-buildpkg@v1
+    - uses: julia-actions/julia-invalidations@v1
+      id: invs_default
+    
+    - name: Report invalidation counts
+      run: |
+        echo "Invalidations on default branch: ${{ steps.invs_default.outputs.total }} (${{ steps.invs_default.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
+        echo "This branch: ${{ steps.invs_pr.outputs.total }} (${{ steps.invs_pr.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
+    - name: Check if the PR does increase number of invalidations
+      if: steps.invs_pr.outputs.total > steps.invs_default.outputs.total
+      run: exit 1
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index f83be7ce..ef9fcadd 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -7,6 +7,7 @@ on:
     branches:
       - master
     tags: '*'
+  workflow_dispatch:
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
@@ -16,7 +17,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
           - 'nightly'
         os:
@@ -24,12 +25,12 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v3
         env:
           cache-name: cache-artifacts
         with:
@@ -42,14 +43,14 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v3
         with:
           file: lcov.info
   docs:
     name: Documentation
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: '1'
diff --git a/Project.toml b/Project.toml
index ece577e5..c55b8442 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.24"
+version = "0.11-DEV"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"
@@ -12,28 +12,37 @@ NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
+[weakdeps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[extensions]
+ForwardDiffStaticArraysExt = "StaticArrays"
+
 [compat]
-Calculus = "0.2, 0.3, 0.4, 0.5"
+Calculus = "0.5"
 CommonSubexpressions = "0.3"
-DiffResults = "0.0.1, 0.0.2, 0.0.3, 0.0.4, 1.0.1"
-DiffRules = "1.4.0"
-DiffTests = "0.0.1, 0.1"
+DiffResults = "1.1"
+DiffRules = "1.4"
+DiffTests = "0.1"
 LogExpFunctions = "0.3"
-NaNMath = "0.2.2, 0.3"
+NaNMath = "1"
 Preferences = "1"
-SpecialFunctions = "0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 1.0, 2"
-StaticArrays = "0.8.3, 0.9, 0.10, 0.11, 0.12, 1.0"
-julia = "1"
+SIMD = "3"
+SpecialFunctions = "1, 2"
+StaticArrays = "1.5"
+julia = "1.6"
 
 [extras]
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 DiffTests = "de460e47-3fe3-5279-bb4a-814414816d5d"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Calculus", "DiffTests", "SparseArrays", "Test", "InteractiveUtils"]
+test = ["Calculus", "DiffTests", "SparseArrays", "StaticArrays", "Test", "InteractiveUtils"]
diff --git a/README.md b/README.md
index 683d5f1f..c9e5154b 100644
--- a/README.md
+++ b/README.md
@@ -8,43 +8,62 @@
 
 ForwardDiff implements methods to take **derivatives**, **gradients**, **Jacobians**, **Hessians**, and higher-order derivatives of native Julia functions (or any callable object, really) using **forward mode automatic differentiation (AD)**.
 
-While performance can vary depending on the functions you evaluate, the algorithms implemented by ForwardDiff **generally outperform non-AD algorithms in both speed and accuracy.**
+While performance can vary depending on the functions you evaluate, the algorithms implemented by ForwardDiff generally outperform non-AD algorithms (such as finite-differencing) in both speed and accuracy.
 
 Here's a simple example showing the package in action:
 
 ```julia
 julia> using ForwardDiff
 
-julia> f(x::Vector) = sum(sin, x) + prod(tan, x) * sum(sqrt, x);
+julia> f(x::Vector) = sin(x[1]) + prod(x[2:end]);  # returns a scalar
 
-julia> x = rand(5) # small size for example's sake
-5-element Array{Float64,1}:
- 0.986403
- 0.140913
- 0.294963
- 0.837125
- 0.650451
+julia> x = vcat(pi/4, 2:4)
+4-element Vector{Float64}:
+ 0.7853981633974483
+ 2.0
+ 3.0
+ 4.0
 
-julia> g = x -> ForwardDiff.gradient(f, x); # g = ∇f
-
-julia> g(x)
-5-element Array{Float64,1}:
- 1.01358
- 2.50014
- 1.72574
- 1.10139
- 1.2445
+julia> ForwardDiff.gradient(f, x)
+4-element Vector{Float64}:
+  0.7071067811865476
+ 12.0
+  8.0
+  6.0
 
 julia> ForwardDiff.hessian(f, x)
-5x5 Array{Float64,2}:
- 0.585111  3.48083  1.7706    0.994057  1.03257
- 3.48083   1.06079  5.79299   3.25245   3.37871
- 1.7706    5.79299  0.423981  1.65416   1.71818
- 0.994057  3.25245  1.65416   0.251396  0.964566
- 1.03257   3.37871  1.71818   0.964566  0.140689
- ```
-
- Trying to switch to the latest version of ForwardDiff? See our [upgrade guide](http://www.juliadiff.org/ForwardDiff.jl/stable/user/upgrade/) for details regarding user-facing changes between releases.
+4×4 Matrix{Float64}:
+ -0.707107  0.0  0.0  0.0
+  0.0       0.0  4.0  3.0
+  0.0       4.0  0.0  2.0
+  0.0       3.0  2.0  0.0
+```
+
+Functions like `f` which map a vector to a scalar are the best case for reverse-mode automatic differentiation,
+but ForwardDiff may still be a good choice if `x` is not too large, as it is much simpler.
+The best case for forward-mode differentiation is a function which maps a scalar to a vector, like this `g`:
+
+```julia
+julia> g(y::Real) = [sin(y), cos(y), tan(y)];  # returns a vector
+
+julia> ForwardDiff.derivative(g, pi/4)
+3-element Vector{Float64}:
+  0.7071067811865476
+ -0.7071067811865475
+  1.9999999999999998
+
+julia> ForwardDiff.jacobian(x) do x  # anonymous function, returns a length-2 vector
+         [sin(x[1]), prod(x[2:end])]
+       end
+2×4 Matrix{Float64}:
+ 0.707107   0.0  0.0  0.0
+ 0.0       12.0  8.0  6.0
+```
+
+See [ForwardDiff's documentation](https://juliadiff.org/ForwardDiff.jl/stable) for full details on how to use this package.
+ForwardDiff relies on [DiffRules](https://github.com/JuliaDiff/DiffRules.jl) for the derivatives of many simple function such as `sin`.
+
+See the [JuliaDiff web page](https://juliadiff.org) for other automatic differentiation packages.
 
 ## Publications
 
diff --git a/docs/src/dev/how_it_works.md b/docs/src/dev/how_it_works.md
index 171492a5..9b68d6ca 100644
--- a/docs/src/dev/how_it_works.md
+++ b/docs/src/dev/how_it_works.md
@@ -39,7 +39,7 @@ central feature that allows ForwardDiff to take derivatives.
 
 In order to implement the above property, elementary numerical functions on a `Dual`
 number are overloaded to evaluate both the original function, *and* evaluate the derivative
-of the function, propogating the derivative via multiplication. For example, `Base.sin`
+of the function, propagating the derivative via multiplication. For example, `Base.sin`
 can be overloaded on `Dual` like so:
 
 ```julia
diff --git a/docs/src/user/advanced.md b/docs/src/user/advanced.md
index a304ff9e..9e2ca8ee 100644
--- a/docs/src/user/advanced.md
+++ b/docs/src/user/advanced.md
@@ -72,7 +72,7 @@ julia> @time gradient!(out, rosenbrock, x, cfg10);
   0.282529 seconds (4 allocations: 160 bytes)
 ```
 
-If you do not explicity provide a chunk size, ForwardDiff will try to guess one for you
+If you do not explicitly provide a chunk size, ForwardDiff will try to guess one for you
 based on your input vector:
 
 ```julia
@@ -84,6 +84,12 @@ julia> @time ForwardDiff.gradient!(out, rosenbrock, x, cfg);
   0.281853 seconds (4 allocations: 160 bytes)
 ```
 
+The underlying heuristic will compute a suitable chunk size smaller or equal
+to the `ForwardDiff.DEFAULT_CHUNK_THRESHOLD` constant. As of ForwardDiff
+v0.10.32 and Julia 1.6, this constant can be configured at load time via
+[Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl) by setting the
+`default_chunk_threshold` value.
+
 If your input dimension is constant across calls, you should explicitly select a chunk size
 rather than relying on ForwardDiff's heuristic. There are two reasons for this. The first is
 that ForwardDiff's heuristic depends only on the input dimension, whereas in reality the
@@ -148,8 +154,24 @@ julia> using ForwardDiff, Preferences
 julia> set_preferences!(ForwardDiff, "nansafe_mode" => true)
 ```
 
+Note that Julia has to be restarted and ForwardDiff has to be reloaded after changing
+this preference.
+
+Alternatively, you can set the preference before loading ForwardDiff, for example via:
+
+```julia
+julia> using Preferences, UUIDs
+
+julia> set_preferences!(UUID("f6369f11-7733-5829-9624-2563aa707210"), "nansafe_mode" => true)
+
+julia> using ForwardDiff
+
+julia> log(ForwardDiff.Dual{:tag}(0.0, 0.0))
+Dual{:tag}(-Inf,0.0)
+```
+
 In the future, we plan on allowing users and downstream library authors to dynamically
-enable [NaN`-safe mode via the `AbstractConfig`
+enable [`NaN`-safe mode via the `AbstractConfig`
 API](https://github.com/JuliaDiff/ForwardDiff.jl/issues/181).
 
 ## Hessian of a vector-valued function
@@ -226,8 +248,8 @@ want to disable this checking.
 
 1. (preferred) Provide an extra `Val{false}()` argument to the differentiation function, e.g.
    ```julia
-   cfg = ForwarDiff.GradientConfig(g, x)
-   ForwarDiff.gradient(f, x, cfg, Val{false}())
+   cfg = ForwardDiff.GradientConfig(g, x)
+   ForwardDiff.gradient(f, x, cfg, Val{false}())
    ```
    If using as part of a library, the tag can be checked manually via
    ```julia
diff --git a/docs/src/user/limitations.md b/docs/src/user/limitations.md
index dc4365e2..1392ae43 100644
--- a/docs/src/user/limitations.md
+++ b/docs/src/user/limitations.md
@@ -12,3 +12,6 @@ function being differentiated):
 - **The target function must be written generically enough to accept numbers of type `T<:Real` as input (or arrays of these numbers).** The function doesn't require a specific type signature, as long as the type signature is generic enough to avoid breaking this rule. This also means that any storage assigned used within the function must be generic as well (see [this comment](https://github.com/JuliaDiff/ForwardDiff.jl/issues/136#issuecomment-237941790) for an example).
 
 - **The types of array inputs must be subtypes of** `AbstractArray` **.** Non-`AbstractArray` array-like types are not officially supported.
+
+ForwardDiff is not natively compatible with rules defined by the [ChainRules.jl](https://github.com/JuliaDiff/ChainRules.jl) ecosystem.
+You can use [ForwardDiffChainRules.jl](https://github.com/ThummeTo/ForwardDiffChainRules.jl) to bridge this gap.
diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl
new file mode 100644
index 00000000..52914cd6
--- /dev/null
+++ b/ext/ForwardDiffStaticArraysExt.jl
@@ -0,0 +1,138 @@
+module ForwardDiffStaticArraysExt
+
+using ForwardDiff, StaticArrays
+using ForwardDiff.LinearAlgebra
+using ForwardDiff.DiffResults
+using ForwardDiff: Dual, partials, GradientConfig, JacobianConfig, HessianConfig, Tag, Chunk,
+                   gradient, hessian, jacobian, gradient!, hessian!, jacobian!,
+                   extract_gradient!, extract_jacobian!, extract_value!,
+                   vector_mode_gradient, vector_mode_gradient!,
+                   vector_mode_jacobian, vector_mode_jacobian!, valtype, value, _lyap_div!
+using DiffResults: DiffResult, ImmutableDiffResult, MutableDiffResult
+
+@generated function dualize(::Type{T}, x::StaticArray) where T
+    N = length(x)
+    dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i}())) for i in 1:N]...)
+    V = StaticArrays.similar_type(x, Dual{T,eltype(x),N})
+    return quote
+        chunk = Chunk{$N}()
+        $(Expr(:meta, :inline))
+        return $V($(dx))
+    end
+end
+
+@inline static_dual_eval(::Type{T}, f, x::StaticArray) where T = f(dualize(T, x))
+
+function LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix}) where {Tg,T<:Real,N}
+    λ,Q = eigen(Symmetric(value.(parent(A))))
+    parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N)
+    Dual{Tg}.(λ, tuple.(parts...))
+end
+
+function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix}) where {Tg,T<:Real,N}
+    λ = eigvals(A)
+    _,Q = eigen(Symmetric(value.(parent(A))))
+    parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N)
+    Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
+end
+
+# Gradient
+@inline ForwardDiff.gradient(f, x::StaticArray)                      = vector_mode_gradient(f, x)
+@inline ForwardDiff.gradient(f, x::StaticArray, cfg::GradientConfig) = gradient(f, x)
+@inline ForwardDiff.gradient(f, x::StaticArray, cfg::GradientConfig, ::Val) = gradient(f, x)
+
+@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray) = vector_mode_gradient!(result, f, x)
+@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::GradientConfig) = gradient!(result, f, x)
+@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::GradientConfig, ::Val) = gradient!(result, f, x)
+
+@generated function extract_gradient(::Type{T}, y::Real, x::S) where {T,S<:StaticArray}
+    result = Expr(:tuple, [:(partials(T, y, $i)) for i in 1:length(x)]...)
+    return quote
+        $(Expr(:meta, :inline))
+        V = StaticArrays.similar_type(S, valtype($y))
+        return V($result)
+    end
+end
+
+@inline function ForwardDiff.vector_mode_gradient(f, x::StaticArray)
+    T = typeof(Tag(f, eltype(x)))
+    return extract_gradient(T, static_dual_eval(T, f, x), x)
+end
+
+@inline function ForwardDiff.vector_mode_gradient!(result, f, x::StaticArray)
+    T = typeof(Tag(f, eltype(x)))
+    return extract_gradient!(T, result, static_dual_eval(T, f, x))
+end
+
+# Jacobian
+@inline ForwardDiff.jacobian(f, x::StaticArray) = vector_mode_jacobian(f, x)
+@inline ForwardDiff.jacobian(f, x::StaticArray, cfg::JacobianConfig) = jacobian(f, x)
+@inline ForwardDiff.jacobian(f, x::StaticArray, cfg::JacobianConfig, ::Val) = jacobian(f, x)
+
+@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray) = vector_mode_jacobian!(result, f, x)
+@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::JacobianConfig) = jacobian!(result, f, x)
+@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::JacobianConfig, ::Val) = jacobian!(result, f, x)
+
+@generated function extract_jacobian(::Type{T}, ydual::StaticArray, x::S) where {T,S<:StaticArray}
+    M, N = length(ydual), length(x)
+    result = Expr(:tuple, [:(partials(T, ydual[$i], $j)) for i in 1:M, j in 1:N]...)
+    return quote
+        $(Expr(:meta, :inline))
+        V = StaticArrays.similar_type(S, valtype(eltype($ydual)), Size($M, $N))
+        return V($result)
+    end
+end
+
+@inline function ForwardDiff.vector_mode_jacobian(f, x::StaticArray)
+    T = typeof(Tag(f, eltype(x)))
+    return extract_jacobian(T, static_dual_eval(T, f, x), x)
+end
+
+function extract_jacobian(::Type{T}, ydual::AbstractArray, x::StaticArray) where T
+    result = similar(ydual, valtype(eltype(ydual)), length(ydual), length(x))
+    return extract_jacobian!(T, result, ydual, length(x))
+end
+
+@inline function ForwardDiff.vector_mode_jacobian!(result, f, x::StaticArray)
+    T = typeof(Tag(f, eltype(x)))
+    ydual = static_dual_eval(T, f, x)
+    result = extract_jacobian!(T, result, ydual, length(x))
+    result = extract_value!(T, result, ydual)
+    return result
+end
+
+@inline function ForwardDiff.vector_mode_jacobian!(result::ImmutableDiffResult, f, x::StaticArray)
+    T = typeof(Tag(f, eltype(x)))
+    ydual = static_dual_eval(T, f, x)
+    result = DiffResults.jacobian!(result, extract_jacobian(T, ydual, x))
+    result = DiffResults.value!(d -> value(T,d), result, ydual)
+    return result
+end
+
+# Hessian
+ForwardDiff.hessian(f, x::StaticArray) = jacobian(y -> gradient(f, y), x)
+ForwardDiff.hessian(f, x::StaticArray, cfg::HessianConfig) = hessian(f, x)
+ForwardDiff.hessian(f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian(f, x)
+
+ForwardDiff.hessian!(result::AbstractArray, f, x::StaticArray) = jacobian!(result, y -> gradient(f, y), x)
+
+ForwardDiff.hessian!(result::MutableDiffResult, f, x::StaticArray) = hessian!(result, f, x, HessianConfig(f, result, x))
+
+ForwardDiff.hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig) = hessian!(result, f, x)
+ForwardDiff.hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian!(result, f, x)
+
+function ForwardDiff.hessian!(result::ImmutableDiffResult, f, x::StaticArray)
+    T = typeof(Tag(f, eltype(x)))
+    d1 = dualize(T, x)
+    d2 = dualize(T, d1)
+    fd2 = f(d2)
+    val = value(T,value(T,fd2))
+    grad = extract_gradient(T,value(T,fd2), x)
+    hess = extract_jacobian(T,partials(T,fd2), x)
+    result = DiffResults.hessian!(result, hess)
+    result = DiffResults.gradient!(result, grad)
+    result = DiffResults.value!(result, val)
+    return result
+end
+
+end
diff --git a/src/ForwardDiff.jl b/src/ForwardDiff.jl
index 93d3b246..9679c1ab 100644
--- a/src/ForwardDiff.jl
+++ b/src/ForwardDiff.jl
@@ -1,20 +1,26 @@
 module ForwardDiff
 
 using DiffRules, DiffResults
-using DiffResults: DiffResult, MutableDiffResult, ImmutableDiffResult
-using StaticArrays
-if VERSION >= v"1.6"
-    using Preferences
-end
+using DiffResults: DiffResult, MutableDiffResult
+using Preferences
 using Random
 using LinearAlgebra
+import SIMD: Vec
 
+using Base: require_one_based_indexing
 import Printf
 import NaNMath
 import SpecialFunctions
 import LogExpFunctions
 import CommonSubexpressions
 
+const SIMDFloat = Union{Float64, Float32}
+const SIMDInt = Union{
+                       Int128, Int64, Int32, Int16, Int8,
+                       UInt128, UInt64, UInt32, UInt16, UInt8,
+                     }
+const SIMDType = Union{SIMDFloat, SIMDInt}
+
 include("prelude.jl")
 include("partials.jl")
 include("dual.jl")
@@ -25,6 +31,10 @@ include("gradient.jl")
 include("jacobian.jl")
 include("hessian.jl")
 
+if !isdefined(Base, :get_extension)
+    include("../ext/ForwardDiffStaticArraysExt.jl")
+end
+
 export DiffResults
 
 end # module
diff --git a/src/apiutils.jl b/src/apiutils.jl
index 971c368c..5a2a1df9 100644
--- a/src/apiutils.jl
+++ b/src/apiutils.jl
@@ -18,19 +18,6 @@ end
 # vector mode function evaluation #
 ###################################
 
-@generated function dualize(::Type{T}, x::StaticArray) where T
-    N = length(x)
-    dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i}())) for i in 1:N]...)
-    V = StaticArrays.similar_type(x, Dual{T,eltype(x),N})
-    return quote
-        chunk = Chunk{$N}()
-        $(Expr(:meta, :inline))
-        return $V($(dx))
-    end
-end
-
-@inline static_dual_eval(::Type{T}, f, x::StaticArray) where T = f(dualize(T, x))
-
 function vector_mode_dual_eval!(f::F, cfg::Union{JacobianConfig,GradientConfig}, x) where {F}
     xdual = cfg.duals
     seed!(xdual, x, cfg.seeds)
diff --git a/src/derivative.jl b/src/derivative.jl
index a60d8d0e..b39e2a48 100644
--- a/src/derivative.jl
+++ b/src/derivative.jl
@@ -22,8 +22,9 @@ stored in `y`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-@inline function derivative(f!, y::AbstractArray, x::Real,
-                            cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK}
+@inline function derivative(f!::F, y::AbstractArray, x::Real,
+                            cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F, T, CHK}
+    require_one_based_indexing(y)
     CHK && checktag(T, f!, x)
     ydual = cfg.duals
     seed!(ydual, y)
@@ -42,6 +43,7 @@ This method assumes that `isa(f(x), Union{Real,AbstractArray})`.
 """
 @inline function derivative!(result::Union{AbstractArray,DiffResult},
                              f::F, x::R) where {F,R<:Real}
+    result isa DiffResult || require_one_based_indexing(result)
     T = typeof(Tag(f, R))
     ydual = f(Dual{T}(x, one(x)))
     result = extract_value!(T, result, ydual)
@@ -58,8 +60,9 @@ called as `f!(y, x)` where the result is stored in `y`.
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
 @inline function derivative!(result::Union{AbstractArray,DiffResult},
-                             f!, y::AbstractArray, x::Real,
-                             cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK}
+                             f!::F, y::AbstractArray, x::Real,
+                             cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F, T, CHK}
+    result isa DiffResult ? require_one_based_indexing(y) : require_one_based_indexing(result, y)
     CHK && checktag(T, f!, x)
     ydual = cfg.duals
     seed!(ydual, y)
@@ -70,6 +73,7 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba
 end
 
 derivative(f, x::AbstractArray) = throw(DimensionMismatch("derivative(f, x) expects that x is a real number. Perhaps you meant gradient(f, x)?"))
+derivative(f, x::Complex) = throw(DimensionMismatch("derivative(f, x) expects that x is a real number (does not support Wirtinger derivatives). Separate real and imaginary parts of the input."))
 
 #####################
 # result extraction #
@@ -78,9 +82,13 @@ derivative(f, x::AbstractArray) = throw(DimensionMismatch("derivative(f, x) expe
 # non-mutating #
 #--------------#
 
-@inline extract_derivative(::Type{T}, y::Dual) where {T}          = partials(T, y, 1)
 @inline extract_derivative(::Type{T}, y::Real) where {T}          = zero(y)
+@inline extract_derivative(::Type{T}, y::Complex) where {T}       = zero(y)
+@inline extract_derivative(::Type{T}, y::Dual) where {T}          = partials(T, y, 1)
 @inline extract_derivative(::Type{T}, y::AbstractArray) where {T} = map(d -> extract_derivative(T,d), y)
+@inline function extract_derivative(::Type{T}, y::Complex{TD}) where {T, TD <: Dual}
+    complex(partials(T, real(y), 1), partials(T, imag(y), 1))
+end
 
 # mutating #
 #----------#
diff --git a/src/dual.jl b/src/dual.jl
index 03330c00..9dacf2b7 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -195,6 +195,38 @@ macro define_ternary_dual_op(f, xyz_body, xy_body, xz_body, yz_body, x_body, y_b
     return esc(defs)
 end
 
+# Support complex-valued functions such as `hankelh1`
+@inline function dual_definition_retval(::Val{T}, val::Real, deriv::Real, partial::Partials) where {T}
+    return Dual{T}(val, deriv * partial)
+end
+@inline function dual_definition_retval(::Val{T}, val::Real, deriv1::Real, partial1::Partials, deriv2::Real, partial2::Partials) where {T}
+    return Dual{T}(val, _mul_partials(partial1, partial2, deriv1, deriv2))
+end
+@inline function dual_definition_retval(::Val{T}, val::Complex, deriv::Union{Real,Complex}, partial::Partials) where {T}
+    reval, imval = reim(val)
+    if deriv isa Real
+        p = deriv * partial
+        return Complex(Dual{T}(reval, p), Dual{T}(imval, zero(p)))
+    else
+        rederiv, imderiv = reim(deriv)
+        return Complex(Dual{T}(reval, rederiv * partial), Dual{T}(imval, imderiv * partial))
+    end
+end
+@inline function dual_definition_retval(::Val{T}, val::Complex, deriv1::Union{Real,Complex}, partial1::Partials, deriv2::Union{Real,Complex}, partial2::Partials) where {T}
+    reval, imval = reim(val)
+    if deriv1 isa Real && deriv2 isa Real
+        p = _mul_partials(partial1, partial2, deriv1, deriv2)
+        return Complex(Dual{T}(reval, p), Dual{T}(imval, zero(p)))
+    else
+        rederiv1, imderiv1 = reim(deriv1)
+        rederiv2, imderiv2 = reim(deriv2)
+        return Complex(
+            Dual{T}(reval, _mul_partials(partial1, partial2, rederiv1, rederiv2)),
+            Dual{T}(imval, _mul_partials(partial1, partial2, imderiv1, imderiv2)),
+        )
+    end
+end
+
 function unary_dual_definition(M, f)
     FD = ForwardDiff
     Mf = M == :Base ? f : :($M.$f)
@@ -206,7 +238,7 @@ function unary_dual_definition(M, f)
         @inline function $M.$f(d::$FD.Dual{T}) where T
             x = $FD.value(d)
             $work
-            return $FD.Dual{T}(val, deriv * $FD.partials(d))
+            return $FD.dual_definition_retval(Val{T}(), val, deriv, $FD.partials(d))
         end
     end
 end
@@ -236,17 +268,17 @@ function binary_dual_definition(M, f)
             begin
                 vx, vy = $FD.value(x), $FD.value(y)
                 $xy_work
-                return $FD.Dual{Txy}(val, $FD._mul_partials($FD.partials(x), $FD.partials(y), dvx, dvy))
+                return $FD.dual_definition_retval(Val{Txy}(), val, dvx, $FD.partials(x), dvy, $FD.partials(y))
             end,
             begin
                 vx = $FD.value(x)
                 $x_work
-                return $FD.Dual{Tx}(val, dvx * $FD.partials(x))
+                return $FD.dual_definition_retval(Val{Tx}(), val, dvx, $FD.partials(x))
             end,
             begin
                 vy = $FD.value(y)
                 $y_work
-                return $FD.Dual{Ty}(val, dvy * $FD.partials(y))
+                return $FD.dual_definition_retval(Val{Ty}(), val, dvy, $FD.partials(y))
             end
         )
     end
@@ -262,6 +294,18 @@ Base.copy(d::Dual) = d
 Base.eps(d::Dual) = eps(value(d))
 Base.eps(::Type{D}) where {D<:Dual} = eps(valtype(D))
 
+# The `base` keyword was added in Julia 1.8:
+# https://github.com/JuliaLang/julia/pull/42428
+if VERSION < v"1.8.0-DEV.725"
+    Base.precision(d::Dual) = precision(value(d))
+    Base.precision(::Type{D}) where {D<:Dual} = precision(valtype(D))
+else
+    Base.precision(d::Dual; base::Integer=2) = precision(value(d); base=base)
+    function Base.precision(::Type{D}; base::Integer=2) where {D<:Dual}
+        precision(valtype(D); base=base)
+    end
+end
+
 function Base.nextfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N}
     ForwardDiff.Dual{T}(nextfloat(d.value), d.partials)
 end
@@ -284,13 +328,14 @@ Base.trunc(d::Dual) = trunc(value(d))
 Base.round(::Type{R}, d::Dual) where {R<:Real} = round(R, value(d))
 Base.round(d::Dual) = round(value(d))
 
-if VERSION ≥ v"1.4"
-    Base.div(x::Dual, y::Dual, r::RoundingMode) = div(value(x), value(y), r)
-else
-    Base.div(x::Dual, y::Dual) = div(value(x), value(y))
-end
+Base.fld(x::Dual, y::Dual) = fld(value(x), value(y))
+
+Base.cld(x::Dual, y::Dual) = cld(value(x), value(y))
+
+Base.exponent(x::Dual) = exponent(value(x))
+
+Base.div(x::Dual, y::Dual, r::RoundingMode) = div(value(x), value(y), r)
 
-Base.hash(d::Dual) = hash(value(d))
 Base.hash(d::Dual, hsh::UInt) = hash(value(d), hsh)
 
 function Base.read(io::IO, ::Type{Dual{T,V,N}}) where {T,V,N}
@@ -336,17 +381,40 @@ for pred in UNARY_PREDICATES
     @eval Base.$(pred)(d::Dual) = $(pred)(value(d))
 end
 
-for pred in BINARY_PREDICATES
+# Before PR#481 this loop ran over this list:
+# BINARY_PREDICATES = Symbol[:isequal, :isless, :<, :>, :(==), :(!=), :(<=), :(>=)]
+# Not a minimal set, as Base defines some in terms of others.
+for pred in [:isless, :<, :>, :(<=), :(>=)]
     @eval begin
         @define_binary_dual_op(
             Base.$(pred),
             $(pred)(value(x), value(y)),
             $(pred)(value(x), y),
-            $(pred)(x, value(y))
+            $(pred)(x, value(y)),
         )
     end
 end
 
+Base.iszero(x::Dual) = iszero(value(x)) && iszero(partials(x))  # shortcut, equivalent to x == zero(x)
+
+for pred in [:isequal, :(==)]
+    @eval begin
+        @define_binary_dual_op(
+            Base.$(pred),
+            $(pred)(value(x), value(y)) && $(pred)(partials(x), partials(y)),
+            $(pred)(value(x), y)        && iszero(partials(x)),
+            $(pred)(x, value(y))        && iszero(partials(y)),
+        )
+    end
+end
+
+@define_binary_dual_op(
+    Base.:(!=),
+    (!=)(value(x), value(y)) || (!=)(partials(x), partials(y)),
+    (!=)(value(x), y)        || !iszero(partials(x)),
+    (!=)(x, value(y))        || !iszero(partials(y)),
+)
+
 ########################
 # Promotion/Conversion #
 ########################
@@ -380,9 +448,9 @@ for R in (Irrational, Real, BigFloat, Bool)
     end
 end
 
-Base.convert(::Type{Dual{T,V,N}}, d::Dual{T}) where {T,V,N} = Dual{T}(convert(V, value(d)), convert(Partials{N,V}, partials(d)))
-Base.convert(::Type{Dual{T,V,N}}, x) where {T,V,N} = Dual{T}(convert(V, x), zero(Partials{N,V}))
-Base.convert(::Type{Dual{T,V,N}}, x::Number) where {T,V,N} = Dual{T}(convert(V, x), zero(Partials{N,V}))
+@inline Base.convert(::Type{Dual{T,V,N}}, d::Dual{T}) where {T,V,N} = Dual{T}(V(value(d)), convert(Partials{N,V}, partials(d)))
+@inline Base.convert(::Type{Dual{T,V,N}}, x) where {T,V,N} = Dual{T}(V(x), zero(Partials{N,V}))
+@inline Base.convert(::Type{Dual{T,V,N}}, x::Number) where {T,V,N} = Dual{T}(V(x), zero(Partials{N,V}))
 Base.convert(::Type{D}, d::D) where {D<:Dual} = d
 
 Base.float(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T,float(V),N}
@@ -486,7 +554,7 @@ for f in (:(Base.:^), :(NaNMath.pow))
             begin
                 v = value(x)
                 expv = ($f)(v, y)
-                if y == zero(y)
+                if y == zero(y) || iszero(partials(x))
                     new_partials = zero(partials(x))
                 else
                     new_partials = partials(x) * y * ($f)(v, y - 1)
@@ -541,6 +609,16 @@ end
 # fma #
 #-----#
 
+@inline function calc_fma_xyz(x::Dual{T,V,N},
+                              y::Dual{T,V,N},
+                              z::Dual{T,V,N}) where {T, V<:SIMDFloat,N}
+    xv, yv, zv = value(x), value(y), value(z)
+    rv = fma(xv, yv, zv)
+    N == 0 && return Dual{T}(rv)
+    xp, yp, zp = Vec(partials(x).values), Vec(partials(y).values), Vec(partials(z).values)
+    parts = Tuple(fma(xv, yp, fma(yv, xp, zp)))
+    Dual{T}(rv, parts)
+end
 @generated function calc_fma_xyz(x::Dual{T,<:Any,N},
                                  y::Dual{T,<:Any,N},
                                  z::Dual{T,<:Any,N}) where {T,N}
@@ -583,6 +661,16 @@ end
 # muladd #
 #--------#
 
+@inline function calc_muladd_xyz(x::Dual{T,V,N},
+                                 y::Dual{T,V,N},
+                                 z::Dual{T,V,N}) where {T, V<:SIMDType,N}
+    xv, yv, zv = value(x), value(y), value(z)
+    rv = muladd(xv, yv, zv)
+    N == 0 && return Dual{T}(rv)
+    xp, yp, zp = Vec(partials(x).values), Vec(partials(y).values), Vec(partials(z).values)
+    parts = Tuple(muladd(xv, yp, muladd(yv, xp, zp)))
+    Dual{T}(rv, parts)
+end
 @generated function calc_muladd_xyz(x::Dual{T,<:Any,N},
                                     y::Dual{T,<:Any,N},
                                     z::Dual{T,<:Any,N}) where {T,N}
@@ -624,6 +712,7 @@ end
 
 # sin/cos #
 #--------#
+
 function Base.sin(d::Dual{T}) where T
     s, c = sincos(value(d))
     return Dual{T}(s, c * partials(d))
@@ -642,11 +731,9 @@ end
 # sincospi #
 #----------#
 
-if VERSION >= v"1.6.0-DEV.292"
-    @inline function Base.sincospi(d::Dual{T}) where T
-        sd, cd = sincospi(value(d))
-        return (Dual{T}(sd, cd * π * partials(d)), Dual{T}(cd, -sd * π * partials(d)))
-    end
+@inline function Base.sincospi(d::Dual{T}) where T
+    sd, cd = sincospi(value(d))
+    return (Dual{T}(sd, cd * π * partials(d)), Dual{T}(cd, -sd * π * partials(d)))
 end
 
 # Symmetric eigvals #
@@ -682,7 +769,7 @@ end
 
 function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
     λ = eigvals(A)
-    _,Q = eigen(SymTridiagonal(value.(parent(A).dv),value.(parent(A).ev)))
+    _,Q = eigen(Symmetric(value.(parent(A))))
     parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N)
     Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
 end
@@ -694,6 +781,23 @@ function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Rea
     Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
 end
 
+# Functions in SpecialFunctions which return tuples #
+# Their derivatives are not defined in DiffRules    #
+#---------------------------------------------------#
+
+function SpecialFunctions.logabsgamma(d::Dual{T,<:Real}) where {T}
+    x = value(d)
+    y, s = SpecialFunctions.logabsgamma(x)
+    return (Dual{T}(y, SpecialFunctions.digamma(x) * partials(d)), s)
+end
+
+# Derivatives wrt to first parameter and precision setting are not supported
+function SpecialFunctions.gamma_inc(a::Real, d::Dual{T,<:Real}, ind::Integer) where {T}
+    x = value(d)
+    p, q = SpecialFunctions.gamma_inc(a, x, ind)
+    ∂p = exp(-x) * x^(a - 1) / SpecialFunctions.gamma(a) * partials(d)
+    return (Dual{T}(p, ∂p), Dual{T}(q, -∂p))
+end
 
 ###################
 # Pretty Printing #
@@ -713,6 +817,4 @@ for op in (:(Base.typemin), :(Base.typemax), :(Base.floatmin), :(Base.floatmax))
     end
 end
 
-if VERSION >= v"1.6.0-rc1"
-    Printf.tofloat(d::Dual) = Printf.tofloat(value(d))
-end
+Printf.tofloat(d::Dual) = Printf.tofloat(value(d))
diff --git a/src/gradient.jl b/src/gradient.jl
index f9f173eb..6d4cc791 100644
--- a/src/gradient.jl
+++ b/src/gradient.jl
@@ -13,7 +13,8 @@ This method assumes that `isa(f(x), Real)`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function gradient(f, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK}
+function gradient(f::F, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {F, T, CHK}
+    require_one_based_indexing(x)
     CHK && checktag(T, f, x)
     if chunksize(cfg) == length(x)
         return vector_mode_gradient(f, x, cfg)
@@ -32,6 +33,7 @@ This method assumes that `isa(f(x), Real)`.
 
 """
 function gradient!(result::Union{AbstractArray,DiffResult}, f::F, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK, F}
+    result isa DiffResult ? require_one_based_indexing(x) : require_one_based_indexing(result, x)
     CHK && checktag(T, f, x)
     if chunksize(cfg) == length(x)
         vector_mode_gradient!(result, f, x, cfg)
@@ -41,29 +43,12 @@ function gradient!(result::Union{AbstractArray,DiffResult}, f::F, x::AbstractArr
     return result
 end
 
-@inline gradient(f, x::StaticArray)                      = vector_mode_gradient(f, x)
-@inline gradient(f, x::StaticArray, cfg::GradientConfig) = gradient(f, x)
-@inline gradient(f, x::StaticArray, cfg::GradientConfig, ::Val) = gradient(f, x)
-
-@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray) = vector_mode_gradient!(result, f, x)
-@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::GradientConfig) = gradient!(result, f, x)
-@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::GradientConfig, ::Val) = gradient!(result, f, x)
-
 gradient(f, x::Real) = throw(DimensionMismatch("gradient(f, x) expects that x is an array. Perhaps you meant derivative(f, x)?"))
 
 #####################
 # result extraction #
 #####################
 
-@generated function extract_gradient(::Type{T}, y::Real, x::S) where {T,S<:StaticArray}
-    result = Expr(:tuple, [:(partials(T, y, $i)) for i in 1:length(x)]...)
-    return quote
-        $(Expr(:meta, :inline))
-        V = StaticArrays.similar_type(S, valtype($y))
-        return V($result)
-    end
-end
-
 function extract_gradient!(::Type{T}, result::DiffResult, y::Real) where {T}
     result = DiffResults.value!(result, y)
     grad = DiffResults.gradient(result)
@@ -115,16 +100,6 @@ function vector_mode_gradient!(result, f::F, x, cfg::GradientConfig{T}) where {T
     return result
 end
 
-@inline function vector_mode_gradient(f, x::StaticArray)
-    T = typeof(Tag(f, eltype(x)))
-    return extract_gradient(T, static_dual_eval(T, f, x), x)
-end
-
-@inline function vector_mode_gradient!(result, f, x::StaticArray)
-    T = typeof(Tag(f, eltype(x)))
-    return extract_gradient!(T, result, static_dual_eval(T, f, x))
-end
-
 ##############
 # chunk mode #
 ##############
diff --git a/src/hessian.jl b/src/hessian.jl
index 7de6109d..9c755c9a 100644
--- a/src/hessian.jl
+++ b/src/hessian.jl
@@ -11,7 +11,8 @@ This method assumes that `isa(f(x), Real)`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function hessian(f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK}
+function hessian(f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F, T,CHK}
+    require_one_based_indexing(x)
     CHK && checktag(T, f, x)
     ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}())
     return jacobian(∇f, x, cfg.jacobian_config, Val{false}())
@@ -27,7 +28,8 @@ This method assumes that `isa(f(x), Real)`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK}
+function hessian!(result::AbstractArray, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T,CHK}
+    require_one_based_indexing(result, x)
     CHK && checktag(T, f, x)
     ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}())
     jacobian!(result, ∇f, x, cfg.jacobian_config, Val{false}())
@@ -61,34 +63,9 @@ because `isa(result, DiffResult)`, `cfg` is constructed as `HessianConfig(f, res
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {T,CHK}
+function hessian!(result::DiffResult, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {F,T,CHK}
     CHK && checktag(T, f, x)
     ∇f! = InnerGradientForHess(result, cfg, f)
     jacobian!(DiffResults.hessian(result), ∇f!, DiffResults.gradient(result), x, cfg.jacobian_config, Val{false}())
     return ∇f!.result
 end
-
-hessian(f, x::StaticArray) = jacobian(y -> gradient(f, y), x)
-hessian(f, x::StaticArray, cfg::HessianConfig) = hessian(f, x)
-hessian(f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian(f, x)
-
-hessian!(result::AbstractArray, f, x::StaticArray) = jacobian!(result, y -> gradient(f, y), x)
-
-hessian!(result::MutableDiffResult, f, x::StaticArray) = hessian!(result, f, x, HessianConfig(f, result, x))
-
-hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig) = hessian!(result, f, x)
-hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian!(result, f, x)
-
-function hessian!(result::ImmutableDiffResult, f, x::StaticArray)
-    T = typeof(Tag(f, eltype(x)))
-    d1 = dualize(T, x)
-    d2 = dualize(T, d1)
-    fd2 = f(d2)
-    val = value(T,value(T,fd2))
-    grad = extract_gradient(T,value(T,fd2), x)
-    hess = extract_jacobian(T,partials(T,fd2), x)
-    result = DiffResults.hessian!(result, hess)
-    result = DiffResults.gradient!(result, grad)
-    result = DiffResults.value!(result, val)
-    return result
-end
diff --git a/src/jacobian.jl b/src/jacobian.jl
index bcda61d7..e7ca96f5 100644
--- a/src/jacobian.jl
+++ b/src/jacobian.jl
@@ -15,7 +15,8 @@ This method assumes that `isa(f(x), AbstractArray)`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function jacobian(f, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK}
+function jacobian(f::F, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T,CHK}
+    require_one_based_indexing(x)
     CHK && checktag(T, f, x)
     if chunksize(cfg) == length(x)
         return vector_mode_jacobian(f, x, cfg)
@@ -32,7 +33,8 @@ stored in `y`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function jacobian(f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK}
+function jacobian(f!::F, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F,T, CHK}
+    require_one_based_indexing(y, x)
     CHK && checktag(T, f!, x)
     if chunksize(cfg) == length(x)
         return vector_mode_jacobian(f!, y, x, cfg)
@@ -52,7 +54,8 @@ This method assumes that `isa(f(x), AbstractArray)`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function jacobian!(result::Union{AbstractArray,DiffResult}, f, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK}
+function jacobian!(result::Union{AbstractArray,DiffResult}, f::F, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T, CHK}
+    result isa DiffResult ? require_one_based_indexing(x) : require_one_based_indexing(result, x)
     CHK && checktag(T, f, x)
     if chunksize(cfg) == length(x)
         vector_mode_jacobian!(result, f, x, cfg)
@@ -72,7 +75,8 @@ This method assumes that `isa(f(x), AbstractArray)`.
 
 Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
 """
-function jacobian!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T,CHK}
+function jacobian!(result::Union{AbstractArray,DiffResult}, f!::F, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F,T,CHK}
+    result isa DiffResult ? require_one_based_indexing(y, x) : require_one_based_indexing(result, y, x)
     CHK && checktag(T, f!, x)
     if chunksize(cfg) == length(x)
         vector_mode_jacobian!(result, f!, y, x, cfg)
@@ -82,35 +86,12 @@ function jacobian!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray
     return result
 end
 
-@inline jacobian(f, x::StaticArray) = vector_mode_jacobian(f, x)
-@inline jacobian(f, x::StaticArray, cfg::JacobianConfig) = jacobian(f, x)
-@inline jacobian(f, x::StaticArray, cfg::JacobianConfig, ::Val) = jacobian(f, x)
-
-@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray) = vector_mode_jacobian!(result, f, x)
-@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::JacobianConfig) = jacobian!(result, f, x)
-@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::JacobianConfig, ::Val) = jacobian!(result, f, x)
-
 jacobian(f, x::Real) = throw(DimensionMismatch("jacobian(f, x) expects that x is an array. Perhaps you meant derivative(f, x)?"))
 
 #####################
 # result extraction #
 #####################
 
-@generated function extract_jacobian(::Type{T}, ydual::StaticArray, x::S) where {T,S<:StaticArray}
-    M, N = length(ydual), length(x)
-    result = Expr(:tuple, [:(partials(T, ydual[$i], $j)) for i in 1:M, j in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        V = StaticArrays.similar_type(S, valtype(eltype($ydual)), Size($M, $N))
-        return V($result)
-    end
-end
-
-function extract_jacobian(::Type{T}, ydual::AbstractArray, x::StaticArray) where T
-    result = similar(ydual, valtype(eltype(ydual)), length(ydual), length(x))
-    return extract_jacobian!(T, result, ydual, length(x))
-end
-
 function extract_jacobian!(::Type{T}, result::AbstractArray, ydual::AbstractArray, n) where {T}
     out_reshaped = reshape(result, length(ydual), n)
     ydual_reshaped = vec(ydual)
@@ -180,27 +161,6 @@ function vector_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T}) wher
     return result
 end
 
-@inline function vector_mode_jacobian(f, x::StaticArray)
-    T = typeof(Tag(f, eltype(x)))
-    return extract_jacobian(T, static_dual_eval(T, f, x), x)
-end
-
-@inline function vector_mode_jacobian!(result, f, x::StaticArray)
-    T = typeof(Tag(f, eltype(x)))
-    ydual = static_dual_eval(T, f, x)
-    result = extract_jacobian!(T, result, ydual, length(x))
-    result = extract_value!(T, result, ydual)
-    return result
-end
-
-@inline function vector_mode_jacobian!(result::ImmutableDiffResult, f, x::StaticArray)
-    T = typeof(Tag(f, eltype(x)))
-    ydual = static_dual_eval(T, f, x)
-    result = DiffResults.jacobian!(result, extract_jacobian(T, ydual, x))
-    result = DiffResults.value!(d -> value(T,d), result, ydual)
-    return result
-end
-
 const JACOBIAN_ERROR = DimensionMismatch("jacobian(f, x) expects that f(x) is an array. Perhaps you meant gradient(f, x)?")
 
 # chunk mode #
diff --git a/src/partials.jl b/src/partials.jl
index fce67b0a..08cf0237 100644
--- a/src/partials.jl
+++ b/src/partials.jl
@@ -52,8 +52,7 @@ Base.:(==)(a::Partials{N}, b::Partials{N}) where {N} = a.values == b.values
 
 const PARTIALS_HASH = hash(Partials)
 
-Base.hash(partials::Partials) = hash(partials.values, PARTIALS_HASH)
-Base.hash(partials::Partials, hsh::UInt64) = hash(hash(partials), hsh)
+Base.hash(partials::Partials, hsh::UInt64) = hash(hash(partials.values, PARTIALS_HASH), hsh)
 
 @inline Base.copy(partials::Partials) = partials
 
@@ -141,6 +140,13 @@ 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
 
+const SIMDFloat = Union{Float64, Float32}
+const SIMDInt = Union{
+                       Int128, Int64, Int32, Int16, Int8,
+                       UInt128, UInt64, UInt32, UInt16, UInt8,
+                     }
+const SIMDType = Union{SIMDFloat, SIMDInt}
+
 ##################################
 # Generated Functions on NTuples #
 ##################################
@@ -164,6 +170,7 @@ end
 @inline rand_tuple(::AbstractRNG, ::Type{Tuple{}}) = tuple()
 @inline rand_tuple(::Type{Tuple{}}) = tuple()
 
+iszero_tuple(tup::NTuple{N,V}) where {N, V<:SIMDType} = sum(Vec(tup) != zero(V)) == 0
 @generated function iszero_tuple(tup::NTuple{N,V}) where {N,V}
     ex = Expr(:&&, [:(z == tup[$i]) for i=1:N]...)
     return quote
@@ -197,29 +204,24 @@ end
     return tupexpr(i -> :(rand(V)), N)
 end
 
-@generated function scale_tuple(tup::NTuple{N}, x) where N
-    return tupexpr(i -> :(tup[$i] * x), N)
-end
+const NT{N,T} = NTuple{N,T}
 
-@generated function div_tuple_by_scalar(tup::NTuple{N}, x) where N
-    return tupexpr(i -> :(tup[$i] / x), N)
-end
+# SIMD implementation
+@inline add_tuples(a::NT{N,T}, b::NT{N,T})               where {N, T<:SIMDType}  = Tuple(Vec(a) + Vec(b))
+@inline sub_tuples(a::NT{N,T}, b::NT{N,T})               where {N, T<:SIMDType}  = Tuple(Vec(a) - Vec(b))
+@inline scale_tuple(tup::NT{N,T}, x::T)                  where {N, T<:SIMDType}  = Tuple(Vec(tup) * x)
+@inline div_tuple_by_scalar(tup::NT{N,T}, x::T)          where {N, T<:SIMDFloat} = Tuple(Vec(tup) / x)
+@inline minus_tuple(tup::NT{N,T})                        where {N, T<:SIMDType}  = Tuple(-Vec(tup))
+@inline mul_tuples(a::NT{N,T}, b::NT{N,T}, af::T, bf::T) where {N, T<:SIMDType}  = Tuple(muladd(af, Vec(a), bf * Vec(b)))
 
-@generated function add_tuples(a::NTuple{N}, b::NTuple{N})  where N
-    return tupexpr(i -> :(a[$i] + b[$i]), N)
-end
 
-@generated function sub_tuples(a::NTuple{N}, b::NTuple{N})  where N
-    return tupexpr(i -> :(a[$i] - b[$i]), N)
-end
-
-@generated function minus_tuple(tup::NTuple{N}) where N
-    return tupexpr(i -> :(-tup[$i]), N)
-end
-
-@generated function mul_tuples(a::NTuple{N}, b::NTuple{N}, afactor, bfactor) where N
-    return tupexpr(i -> :((afactor * a[$i]) + (bfactor * b[$i])), N)
-end
+# Fallback implementations
+@generated add_tuples(a::NT{N}, b::NT{N})         where N = tupexpr(i -> :(a[$i] + b[$i]), N)
+@generated sub_tuples(a::NT{N}, b::NT{N})         where N = tupexpr(i -> :(a[$i] - b[$i]), N)
+@generated scale_tuple(tup::NT{N}, x)             where N = tupexpr(i -> :(tup[$i] * x), N)
+@generated div_tuple_by_scalar(tup::NT{N}, x)     where N = tupexpr(i -> :(tup[$i] / x), N)
+@generated minus_tuple(tup::NT{N})                where N = tupexpr(i -> :(-tup[$i]), N)
+@generated mul_tuples(a::NT{N}, b::NT{N}, af, bf) where N = tupexpr(i -> :(muladd(af, a[$i], bf * b[$i])), N)
 
 ###################
 # Pretty Printing #
diff --git a/src/prelude.jl b/src/prelude.jl
index 2c91fbdd..1b372fdc 100644
--- a/src/prelude.jl
+++ b/src/prelude.jl
@@ -1,15 +1,10 @@
-@static if VERSION >= v"1.6"
-    const NANSAFE_MODE_ENABLED = @load_preference("nansafe_mode", false)
-else
-    const NANSAFE_MODE_ENABLED = false
-end
+const NANSAFE_MODE_ENABLED = @load_preference("nansafe_mode", false)
+const DEFAULT_CHUNK_THRESHOLD = @load_preference("default_chunk_threshold", 12)
 
 const AMBIGUOUS_TYPES = (AbstractFloat, Irrational, Integer, Rational, Real, RoundingMode)
 
 const UNARY_PREDICATES = Symbol[:isinf, :isnan, :isfinite, :iseven, :isodd, :isreal, :isinteger]
 
-const BINARY_PREDICATES = Symbol[:isequal, :isless, :<, :>, :(==), :(!=), :(<=), :(>=)]
-
 const DEFAULT_CHUNK_THRESHOLD = 12
 
 struct Chunk{N} end
diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl
index a16a0103..2a0075a7 100644
--- a/test/AllocationsTest.jl
+++ b/test/AllocationsTest.jl
@@ -4,6 +4,8 @@ using ForwardDiff
 
 include(joinpath(dirname(@__FILE__), "utils.jl"))
 
+convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,Float64,8},4},2}, 1.3)
+
 @testset "Test seed! allocations" begin
     x = rand(1000)
     cfg = ForwardDiff.GradientConfig(nothing, x)
@@ -28,6 +30,11 @@ include(joinpath(dirname(@__FILE__), "utils.jl"))
     alloc = @allocated ForwardDiff.seed!(duals, x, index, seed)
     alloc = @allocated ForwardDiff.seed!(duals, x, index, seed)
     @test alloc == 0
+    
+    alloc = @allocated convert_test_574()
+    alloc = @allocated convert_test_574()
+    @test alloc == 0
+  
 end
 
-end
\ No newline at end of file
+end
diff --git a/test/DerivativeTest.jl b/test/DerivativeTest.jl
index f5645a7e..dfdd8ed2 100644
--- a/test/DerivativeTest.jl
+++ b/test/DerivativeTest.jl
@@ -100,4 +100,8 @@ end
     @test_throws DimensionMismatch ForwardDiff.derivative(sum, fill(2pi, 3))
 end
 
+@testset "complex output" begin
+    @test ForwardDiff.derivative(x -> (1+im)*x, 0) == (1+im)
+end
+
 end # module
diff --git a/test/DualTest.jl b/test/DualTest.jl
index 5a4a8350..938cd0e6 100644
--- a/test/DualTest.jl
+++ b/test/DualTest.jl
@@ -12,6 +12,7 @@ using DiffRules
 import Calculus
 
 struct TestTag end
+struct OuterTestTag end
 
 samerng() = MersenneTwister(1)
 
@@ -26,8 +27,10 @@ dual_isapprox(a::Dual{T,T1,T2}, b::Dual{T3,T4,T5}) where {T,T1,T2,T3,T4,T5} = er
 
 ForwardDiff.:≺(::Type{TestTag()}, ::Int) = true
 ForwardDiff.:≺(::Int, ::Type{TestTag()}) = false
+ForwardDiff.:≺(::Type{TestTag}, ::Type{OuterTestTag}) = true
+ForwardDiff.:≺(::Type{OuterTestTag}, ::Type{TestTag}) = false
 
-for N in (0,3), M in (0,4), V in (Int, Float32)
+@testset "Dual{Z,$V,$N} and Dual{Z,Dual{Z,$V,$M},$N}" for N in (0,3), M in (0,4), V in (Int, Float32)
     println("  ...testing Dual{TestTag(),$V,$N} and Dual{TestTag(),Dual{TestTag(),$V,$M},$N}")
 
     PARTIALS = Partials{N,V}(ntuple(n -> intrand(V), N))
@@ -42,6 +45,13 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
     PRIMAL3 = intrand(V)
     FDNUM3 = Dual{TestTag()}(PRIMAL3, PARTIALS3)
 
+    if !allunique([PRIMAL, PRIMAL2, PRIMAL3])
+        @info "testing with non-unique primals" PRIMAL PRIMAL2 PRIMAL3
+    end
+    if N > 0 && !allunique([PARTIALS, PARTIALS2, PARTIALS3])
+        @info "testing with non-unique partials" PARTIALS PARTIALS2 PARTIALS3
+    end
+
     M_PARTIALS = Partials{M,V}(ntuple(m -> intrand(V), M))
     NESTED_PARTIALS = convert(Partials{N,Dual{TestTag(),V,M}}, PARTIALS)
     NESTED_FDNUM = Dual{TestTag()}(Dual{TestTag()}(PRIMAL, M_PARTIALS), NESTED_PARTIALS)
@@ -105,6 +115,17 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
         @test eps(NESTED_FDNUM) === eps(PRIMAL)
         @test eps(typeof(NESTED_FDNUM)) === eps(V)
 
+        @test precision(FDNUM) === precision(PRIMAL)
+        @test precision(typeof(FDNUM)) === precision(V)
+        @test precision(NESTED_FDNUM) === precision(PRIMAL)
+        @test precision(typeof(NESTED_FDNUM)) === precision(V)
+        if VERSION >= v"1.8.0-DEV.725" # https://github.com/JuliaLang/julia/pull/42428
+            @test precision(FDNUM; base=10) === precision(PRIMAL; base=10)
+            @test precision(typeof(FDNUM); base=10) === precision(V; base=10)
+            @test precision(NESTED_FDNUM; base=10) === precision(PRIMAL; base=10)
+            @test precision(typeof(NESTED_FDNUM); base=10) === precision(V; base=10)
+        end
+
         @test floor(Int, FDNUM) === floor(Int, PRIMAL)
         @test floor(Int, FDNUM2) === floor(Int, PRIMAL2)
         @test floor(Int, NESTED_FDNUM) === floor(Int, PRIMAL)
@@ -137,6 +158,18 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
         @test round(FDNUM2) === round(PRIMAL2)
         @test round(NESTED_FDNUM) === round(PRIMAL)
 
+        @test fld(FDNUM, FDNUM2) === fld(PRIMAL, PRIMAL2)
+        @test fld(FDNUM, PRIMAL2) === fld(PRIMAL, PRIMAL2)
+        @test fld(PRIMAL, FDNUM2) === fld(PRIMAL, PRIMAL2)
+
+        @test exponent(FDNUM) === exponent(PRIMAL)
+        @test exponent(FDNUM2) === exponent(PRIMAL2)
+        @test exponent(NESTED_FDNUM) === exponent(PRIMAL)
+
+        @test cld(FDNUM, FDNUM2) === cld(PRIMAL, PRIMAL2)
+        @test cld(FDNUM, PRIMAL2) === cld(PRIMAL, PRIMAL2)
+        @test cld(PRIMAL, FDNUM2) === cld(PRIMAL, PRIMAL2)
+
         @test div(FDNUM, FDNUM2) === div(PRIMAL, PRIMAL2)
         @test div(FDNUM, PRIMAL2) === div(PRIMAL, PRIMAL2)
         @test div(PRIMAL, FDNUM2) === div(PRIMAL, PRIMAL2)
@@ -145,10 +178,8 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
         @test div(NESTED_FDNUM, PRIMAL2) === div(PRIMAL, PRIMAL2)
         @test div(PRIMAL, NESTED_FDNUM2) === div(PRIMAL, PRIMAL2)
 
-        if VERSION ≥ v"1.4"
-            @test div(FDNUM, FDNUM2, RoundUp) === div(PRIMAL, PRIMAL2, RoundUp)
-            @test div(NESTED_FDNUM, NESTED_FDNUM2, RoundUp) === div(PRIMAL, PRIMAL2, RoundUp)
-        end
+        @test div(FDNUM, FDNUM2, RoundUp) === div(PRIMAL, PRIMAL2, RoundUp)
+        @test div(NESTED_FDNUM, NESTED_FDNUM2, RoundUp) === div(PRIMAL, PRIMAL2, RoundUp)
 
         @test Base.rtoldefault(typeof(FDNUM)) ≡ Base.rtoldefault(typeof(PRIMAL))
         @test Dual{TestTag()}(PRIMAL-eps(V), PARTIALS) ≈ FDNUM
@@ -209,15 +240,27 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
     @test ForwardDiff.isconstant(one(NESTED_FDNUM))
     @test ForwardDiff.isconstant(NESTED_FDNUM) == (N == 0)
 
-    @test isequal(FDNUM, Dual{TestTag()}(PRIMAL, PARTIALS2))
-    @test isequal(PRIMAL, PRIMAL2) == isequal(FDNUM, FDNUM2)
-
-    @test isequal(NESTED_FDNUM, Dual{TestTag()}(Dual{TestTag()}(PRIMAL, M_PARTIALS2), NESTED_PARTIALS2))
-    @test isequal(PRIMAL, PRIMAL2) == isequal(NESTED_FDNUM, NESTED_FDNUM2)
-
-    @test FDNUM == Dual{TestTag()}(PRIMAL, PARTIALS2)
-    @test (PRIMAL == PRIMAL2) == (FDNUM == FDNUM2)
-    @test (PRIMAL == PRIMAL2) == (NESTED_FDNUM == NESTED_FDNUM2)
+    # Recall that FDNUM = Dual{TestTag()}(PRIMAL, PARTIALS) has N partials, 
+    # and FDNUM2 has everything with a 2, and all random numbers nonzero.
+    # M is the length of M_PARTIALS, which affects:
+    # NESTED_FDNUM = Dual{TestTag()}(Dual{TestTag()}(PRIMAL, M_PARTIALS), NESTED_PARTIALS)
+
+    @test (FDNUM == Dual{TestTag()}(PRIMAL, PARTIALS2)) == (PARTIALS == PARTIALS2)
+    @test isequal(FDNUM, Dual{TestTag()}(PRIMAL, PARTIALS2)) == (PARTIALS == PARTIALS2)
+    @test isequal(NESTED_FDNUM, Dual{TestTag()}(Dual{TestTag()}(PRIMAL, M_PARTIALS2), NESTED_PARTIALS2)) == ((M_PARTIALS == M_PARTIALS2) && (NESTED_PARTIALS == NESTED_PARTIALS2))
+
+    if PRIMAL == PRIMAL2
+        @test isequal(FDNUM, Dual{TestTag()}(PRIMAL, PARTIALS2)) == (PARTIALS == PARTIALS2)
+        @test isequal(FDNUM, FDNUM2) == (PARTIALS == PARTIALS2)
+        
+        @test (FDNUM == FDNUM2) == (PARTIALS == PARTIALS2)
+        @test (NESTED_FDNUM == NESTED_FDNUM2) == ((M_PARTIALS == M_PARTIALS2) && (NESTED_PARTIALS == NESTED_PARTIALS2))
+    else
+        @test !isequal(FDNUM, FDNUM2)
+        
+        @test FDNUM != FDNUM2
+        @test NESTED_FDNUM != NESTED_FDNUM2
+    end
 
     @test isless(Dual{TestTag()}(1, PARTIALS), Dual{TestTag()}(2, PARTIALS2))
     @test !(isless(Dual{TestTag()}(1, PARTIALS), Dual{TestTag()}(1, PARTIALS2)))
@@ -322,7 +365,7 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
     @test typeof(WIDE_NESTED_FDNUM) === Dual{TestTag(),Dual{TestTag(),WIDE_T,M},N}
 
     @test value(WIDE_FDNUM) == PRIMAL
-    @test value(WIDE_NESTED_FDNUM) == PRIMAL
+    @test (value(WIDE_NESTED_FDNUM) == PRIMAL) == (M == 0)
 
     @test convert(Dual, FDNUM) === FDNUM
     @test convert(Dual, NESTED_FDNUM) === NESTED_FDNUM
@@ -373,6 +416,8 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
     #----------#
 
     if M > 0 && N > 0
+        # Recall that FDNUM = Dual{TestTag()}(PRIMAL, PARTIALS) has N partials,
+        # all random numbers nonzero, and FDNUM2 another draw. M only affects NESTED_FDNUM.
         @test Dual{1}(FDNUM) / Dual{1}(PRIMAL) === Dual{1}(FDNUM / PRIMAL)
         @test Dual{1}(PRIMAL) / Dual{1}(FDNUM) === Dual{1}(PRIMAL / FDNUM)
         @test_broken Dual{1}(FDNUM) / FDNUM2 === Dual{1}(FDNUM / FDNUM2)
@@ -391,6 +436,7 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
 
     # Exponentiation #
     #----------------#
+
     # If V == Int, the LHS terms are Int's. Large inputs cause integer overflow
     # within the generic fallback of `isapprox`, resulting in a DomainError.
     # Promote to Float64 to avoid issues.
@@ -420,8 +466,8 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
     @test abs(NESTED_FDNUM) === NESTED_FDNUM
 
     if V != Int
-        for (M, f, arity) in DiffRules.diffrules(filter_modules = nothing)
-            if f in (:hankelh1, :hankelh1x, :hankelh2, :hankelh2x, :/, :rem2pi)
+        @testset "$f" for (M, f, arity) in DiffRules.diffrules(filter_modules = nothing)
+            if f in (:/, :rem2pi)
                 continue  # Skip these rules
             elseif !(isdefined(@__MODULE__, M) && isdefined(getfield(@__MODULE__, M), f))
                 continue  # Skip rules for methods not defined in the current scope
@@ -438,9 +484,20 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
                 end
                 @eval begin
                     x = rand() + $modifier
-                    dx = $M.$f(Dual{TestTag()}(x, one(x)))
-                    @test value(dx) == $M.$f(x)
-                    @test partials(dx, 1) == $deriv
+                    dx = @inferred $M.$f(Dual{TestTag()}(x, one(x)))
+                    actualval = $M.$f(x)
+                    @assert actualval isa Real || actualval isa Complex
+                    if actualval isa Real
+                        @test dx isa Dual{TestTag()}
+                        @test value(dx) == actualval
+                        @test partials(dx, 1) == $deriv
+                    else
+                        @test dx isa Complex{<:Dual{TestTag()}}
+                        @test value(real(dx)) == real(actualval)
+                        @test value(imag(dx)) == imag(actualval)
+                        @test partials(real(dx), 1) == real($deriv)
+                        @test partials(imag(dx), 1) == imag($deriv)
+                    end
                 end
             elseif arity == 2
                 derivs = DiffRules.diffrule(M, f, :x, :y)
@@ -453,14 +510,35 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
                 end
                 @eval begin
                     x, y = $x, $y
-                    dx = $M.$f(Dual{TestTag()}(x, one(x)), y)
-                    dy = $M.$f(x, Dual{TestTag()}(y, one(y)))
+                    dx = @inferred $M.$f(Dual{TestTag()}(x, one(x)), y)
+                    dy = @inferred $M.$f(x, Dual{TestTag()}(y, one(y)))
                     actualdx = $(derivs[1])
                     actualdy = $(derivs[2])
-                    @test value(dx) == $M.$f(x, y)
-                    @test value(dy) == value(dx)
-                    @test partials(dx, 1) ≈ actualdx nans=true
-                    @test partials(dy, 1) ≈ actualdy nans=true
+                    actualval = $M.$f(x, y)
+                    @assert actualval isa Real || actualval isa Complex
+                    if actualval isa Real
+                        @test dx isa Dual{TestTag()}
+                        @test dy isa Dual{TestTag()}
+                        @test value(dx) == actualval
+                        @test value(dy) == actualval
+                        @test partials(dx, 1) ≈ actualdx nans=true
+                        @test partials(dy, 1) ≈ actualdy nans=true
+                    else
+                        @test dx isa Complex{<:Dual{TestTag()}}
+                        @test dy isa Complex{<:Dual{TestTag()}}
+                        # @test real(value(dx)) == real(actualval)
+                        # @test real(value(dy)) == real(actualval)
+                        # @test imag(value(dx)) == imag(actualval)
+                        # @test imag(value(dy)) == imag(actualval)
+                        @test value(real(dx)) == real(actualval)
+                        @test value(real(dy)) == real(actualval)
+                        @test value(imag(dx)) == imag(actualval)
+                        @test value(imag(dy)) == imag(actualval)
+                        @test partials(real(dx), 1) ≈ real(actualdx) nans=true
+                        @test partials(real(dy), 1) ≈ real(actualdy) nans=true
+                        @test partials(imag(dx), 1) ≈ imag(actualdx) nans=true
+                        @test partials(imag(dy), 1) ≈ imag(actualdy) nans=true
+                    end
                 end
             end
         end
@@ -474,9 +552,7 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
 
     @test all(map(dual_isapprox, ForwardDiff.sincos(FDNUM), (sin(FDNUM), cos(FDNUM))))
 
-    if VERSION >= v"1.6.0-DEV.292"
-        @test all(map(dual_isapprox, sincospi(FDNUM), (sinpi(FDNUM), cospi(FDNUM))))
-    end
+    @test all(map(dual_isapprox, sincospi(FDNUM), (sinpi(FDNUM), cospi(FDNUM))))
 
     if V === Float32
         @test typeof(sqrt(FDNUM)) === typeof(FDNUM)
@@ -492,8 +568,36 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
         @test dual_isapprox(f(FDNUM, PRIMAL2, PRIMAL3), Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL2*PARTIALS))
         @test dual_isapprox(f(PRIMAL, PRIMAL2, FDNUM3), Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PARTIALS3))
     end
+
+    # Functions in Specialfunctions that return tuples and
+    # therefore are not supported by DiffRules
+    @test dual_isapprox(logabsgamma(FDNUM)[1], loggamma(abs(FDNUM)))
+    @test dual_isapprox(logabsgamma(FDNUM)[2], sign(gamma(FDNUM)))
+
+    a = rand(float(V))
+    fdnum = Dual{TestTag()}(1 + PRIMAL, PARTIALS) # 1 + PRIMAL avoids issues with finite differencing close to 0
+    for ind in ((), (0,), (1,), (2,))
+        # Only test if primal method exists
+        # (e.g., older versions of SpecialFunctions don't define `gamma_inc(a, x)` but only `gamma_inc(a, x, ind)`
+        hasmethod(gamma_inc, typeof((a, 1 + PRIMAL, ind...))) || continue
+
+        pq = gamma_inc(a, fdnum, ind...)
+        @test pq isa Tuple{Dual{TestTag()},Dual{TestTag()}}
+        # We have to adjust tolerances if lower accuracy is requested
+        # Therefore we don't use `dual_isapprox`
+        tol = V === Float32 ? 5f-4 : 1e-5
+        tol = tol^(one(tol) / 2^(isempty(ind) ? 0 : first(ind)))
+        for i in 1:2
+            @test value(pq[i]) ≈ gamma_inc(a, 1 + PRIMAL, ind...)[i] rtol=tol
+            @test partials(pq[i]) ≈ PARTIALS * Calculus.derivative(x -> gamma_inc(a, x, ind...)[i], 1 + PRIMAL) rtol=tol
+        end
+    end
 end
 
+#############
+# bug fixes #
+#############
+
 @testset "Exponentiation of zero" begin
     x0 = 0.0
     x1 = Dual{:t1}(x0, 1.0)
@@ -503,6 +607,9 @@ end
     @test pow(x3, 2) === x3^2 === x3 * x3
     @test pow(x2, 1) === x2^1 === x2
     @test pow(x1, 0) === x1^0 === Dual{:t1}(1.0, 0.0)
+    y = Dual{typeof(TestTag())}(1.0, 0.0, 1.0);
+    x = Dual{typeof(OuterTestTag())}(0*y, 0*y);
+    @test iszero(ForwardDiff.partials(ForwardDiff.partials(x^y)[1]))
 end
 
 @testset "Type min/max" begin
@@ -537,13 +644,11 @@ end
     @test length(UnitRange(Dual(1.5,1), Dual(3.5,3))) == 3
 end
 
-if VERSION >= v"1.6.0-rc1"
-    @testset "@printf" begin
-        for T in (Float16, Float32, Float64, BigFloat)
-            d1 = Dual(one(T))
-            @test_nowarn @printf("Testing @printf: %.2e\n", d1)
-            @test @sprintf("Testing @sprintf: %.2e\n", d1) == "Testing @sprintf: 1.00e+00\n"
-        end
+@testset "@printf" begin
+    for T in (Float16, Float32, Float64, BigFloat)
+        d1 = Dual(one(T))
+        @test_nowarn @printf("Testing @printf: %.2e\n", d1)
+        @test @sprintf("Testing @sprintf: %.2e\n", d1) == "Testing @sprintf: 1.00e+00\n"
     end
 end
 
diff --git a/test/GradientTest.jl b/test/GradientTest.jl
index 707fcd68..a74f4304 100644
--- a/test/GradientTest.jl
+++ b/test/GradientTest.jl
@@ -3,6 +3,7 @@ module GradientTest
 import Calculus
 
 using Test
+using LinearAlgebra
 using ForwardDiff
 using ForwardDiff: Dual, Tag
 using StaticArrays
@@ -19,7 +20,7 @@ x = [0.1, 0.2, 0.3]
 v = f(x)
 g = [-9.4, 15.6, 52.0]
 
-for c in (1, 2, 3), tag in (nothing, Tag(f, eltype(x)))
+@testset "Rosenbrock, chunk size = $c and tag = $(repr(tag))" for c in (1, 2, 3), tag in (nothing, Tag(f, eltype(x)))
     println("  ...running hardcoded test with chunk size = $c and tag = $(repr(tag))")
     cfg = ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{c}(), tag)
 
@@ -55,7 +56,7 @@ cfgx = ForwardDiff.GradientConfig(sin, x)
 # test vs. Calculus.jl #
 ########################
 
-for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
+@testset "$f" for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
     v = f(X)
     g = ForwardDiff.gradient(f, X)
     @test isapprox(g, Calculus.gradient(f, X), atol=FINITEDIFF_ERROR)
@@ -83,9 +84,9 @@ end
 
 println("  ...testing specialized StaticArray codepaths")
 
-x = rand(3, 3)
+@testset "$T" for T in (StaticArrays.SArray, StaticArrays.MArray)
+    x = rand(3, 3)
 
-for T in (StaticArrays.SArray, StaticArrays.MArray)
     sx = T{Tuple{3,3}}(x)
 
     cfg = ForwardDiff.GradientConfig(nothing, x)
@@ -148,6 +149,10 @@ end
     @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0,  1.5]), [0.0, 0.0])
 end
 
+#############
+# bug fixes #
+#############
+
 # Issue 399
 @testset "chunk size zero" begin
     f_const(x) = 1.0
@@ -162,6 +167,7 @@ end
     @test_throws DimensionMismatch ForwardDiff.gradient(identity, fill(2pi, 10^6)) # chunk_mode_gradient
 end
 
+# Issue 548
 @testset "ArithmeticStyle" begin
     function f(p)
         sum(collect(0.0:p[1]:p[2]))
@@ -169,4 +175,45 @@ end
     @test ForwardDiff.gradient(f, [0.2,25.0]) == [7875.0, 0.0]
 end
 
+@testset "det with branches" begin
+    # Issue 197
+    det2(A) = return (
+        A[1,1]*(A[2,2]*A[3,3]-A[2,3]*A[3,2]) -
+        A[1,2]*(A[2,1]*A[3,3]-A[2,3]*A[3,1]) +
+        A[1,3]*(A[2,1]*A[3,2]-A[2,2]*A[3,1])
+    )
+
+    A = [1 0 0; 0 2 0; 0 pi 3]
+    @test det2(A) == det(A) == 6
+    @test istril(A)
+
+    ∇A = [6 0 0; 0 3 -pi; 0 0 2]
+    @test ForwardDiff.gradient(det2, A) ≈ ∇A
+    @test ForwardDiff.gradient(det, A) ≈ ∇A
+
+    # And issue 407
+    @test ForwardDiff.hessian(det, A) ≈ ForwardDiff.hessian(det2, A)
+    
+    # https://discourse.julialang.org/t/forwarddiff-and-zygote-return-wrong-jacobian-for-log-det-l/77961
+    S = [1.0 0.8; 0.8 1.0]
+    L = cholesky(S).L
+    @test ForwardDiff.gradient(L -> log(det(L)), Matrix(L)) ≈ [1.0 -1.3333333333333337; 0.0 1.666666666666667]
+    @test ForwardDiff.gradient(L -> logdet(L), Matrix(L)) ≈ [1.0 -1.3333333333333337; 0.0 1.666666666666667]
+end
+
+@testset "branches in mul!" begin
+    a, b = rand(3,3), rand(3,3)
+
+    # Issue 536, version with 3-arg *, Julia 1.7:
+    @test ForwardDiff.derivative(x -> sum(x*a*b), 0.0) ≈ sum(a * b)
+
+    # version with just mul!
+    dx = ForwardDiff.derivative(0.0) do x
+        c = similar(a, typeof(x))
+        mul!(c, a, b, x, false)
+        sum(c)
+    end
+    @test dx ≈ sum(a * b)
+end
+
 end # module
diff --git a/test/HessianTest.jl b/test/HessianTest.jl
index 1df4232b..3f8b08cb 100644
--- a/test/HessianTest.jl
+++ b/test/HessianTest.jl
@@ -3,6 +3,7 @@ module HessianTest
 import Calculus
 
 using Test
+using LinearAlgebra
 using ForwardDiff
 using ForwardDiff: Dual, Tag
 using StaticArrays
@@ -157,4 +158,11 @@ for T in (StaticArrays.SArray, StaticArrays.MArray)
     @test DiffResults.hessian(sresult3) == DiffResults.hessian(result)
 end
 
+@testset "branches in dot" begin
+    # https://github.com/JuliaDiff/ForwardDiff.jl/issues/551
+    H = [1 2 3; 4 5 6; 7 8 9];
+    @test ForwardDiff.hessian(x->dot(x,H,x), fill(0.00001, 3)) ≈ [2 6 10; 6 10 14; 10 14 18]
+    @test ForwardDiff.hessian(x->dot(x,H,x), zeros(3)) ≈ [2 6 10; 6 10 14; 10 14 18]
+end
+
 end # module
diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index 9ea3b0a7..5f926bc7 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -226,6 +226,10 @@ for T in (StaticArrays.SArray, StaticArrays.MArray)
     @test DiffResults.jacobian(sresult3) == DiffResults.jacobian(result)
 end
 
+#########
+# misc. #
+#########
+
 @testset "dimension errors for jacobian" begin
     @test_throws DimensionMismatch ForwardDiff.jacobian(identity, 2pi) # input
     @test_throws DimensionMismatch ForwardDiff.jacobian(sum, fill(2pi, 2)) # vector_mode_jacobian
@@ -235,6 +239,14 @@ end
 @testset "eigen" begin
     @test ForwardDiff.jacobian(x -> eigvals(SymTridiagonal(x, x[1:end-1])), [1.,2.]) ≈ [(1 - 3/sqrt(5))/2 (1 - 1/sqrt(5))/2 ; (1 + 3/sqrt(5))/2 (1 + 1/sqrt(5))/2]
     @test ForwardDiff.jacobian(x -> eigvals(Symmetric(x*x')), [1.,2.]) ≈ [0 0; 2 4]
+    
+    x0 = [1.0, 2.0];
+    ev1(x) = eigen(Symmetric(x*x')).vectors[:,1]
+    @test ForwardDiff.jacobian(ev1, x0) ≈ Calculus.finite_difference_jacobian(ev1, x0)
+    ev2(x) = eigen(SymTridiagonal(x, x[1:end-1])).vectors[:,1]
+    @test ForwardDiff.jacobian(ev2, x0) ≈ Calculus.finite_difference_jacobian(ev2, x0)
+    x0_static = SVector{2}(x0)
+    @test ForwardDiff.jacobian(ev1, x0_static) ≈ Calculus.finite_difference_jacobian(ev1, x0)
 end
 
 @testset "type stability" begin
diff --git a/test/MiscTest.jl b/test/MiscTest.jl
index 2a22cb2f..0ed8039a 100644
--- a/test/MiscTest.jl
+++ b/test/MiscTest.jl
@@ -158,4 +158,7 @@ end
 @test ForwardDiff.derivative(x -> rem2pi(x, RoundUp), rand()) == 1
 @test ForwardDiff.derivative(x -> rem2pi(x, RoundDown), rand()) == 1
 
+# example from https://github.com/JuliaDiff/DiffRules.jl/pull/98#issuecomment-1574420052
+@test only(ForwardDiff.hessian(t -> abs(t[1])^2, [0.0])) == 2
+
 end # module
diff --git a/test/PartialsTest.jl b/test/PartialsTest.jl
index 39fb05d7..591653d4 100644
--- a/test/PartialsTest.jl
+++ b/test/PartialsTest.jl
@@ -7,7 +7,7 @@ using ForwardDiff: Partials
 
 samerng() = MersenneTwister(1)
 
-for N in (0, 3), T in (Int, Float32, Float64)
+@testset "Partials{$N,$T}" for N in (0, 3), T in (Int, Float32, Float64)
     println("  ...testing Partials{$N,$T}")
 
     VALUES = (rand(T,N)...,)
@@ -62,7 +62,6 @@ for N in (0, 3), T in (Int, Float32, Float64)
     @test isequal(PARTIALS, copy(PARTIALS))
     @test isequal(PARTIALS, PARTIALS2) == (N == 0)
 
-    @test hash(PARTIALS) == hash(VALUES, ForwardDiff.PARTIALS_HASH)
     @test hash(PARTIALS) == hash(copy(PARTIALS))
     @test hash(PARTIALS, hash(1)) == hash(copy(PARTIALS), hash(1))
     @test hash(PARTIALS, hash(1)) == hash(copy(PARTIALS), hash(1))
@@ -114,7 +113,7 @@ for N in (0, 3), T in (Int, Float32, Float64)
 
     if N > 0
         @test ForwardDiff._div_partials(PARTIALS, PARTIALS2, X, Y) == ForwardDiff._mul_partials(PARTIALS, PARTIALS2, inv(Y), -X/(Y^2))
-        @test ForwardDiff._mul_partials(PARTIALS, PARTIALS2, X, Y).values == map((a, b) -> (X * a) + (Y * b), VALUES, VALUES2)
+        @test all(isapprox.(ForwardDiff._mul_partials(PARTIALS, PARTIALS2, X, Y).values, map((a, b) -> (X * a) + (Y * b), VALUES, VALUES2)))
         @test ForwardDiff._mul_partials(ZERO_PARTIALS, PARTIALS, X, Y) == Y * PARTIALS
         @test ForwardDiff._mul_partials(PARTIALS, ZERO_PARTIALS, X, Y) == X * PARTIALS
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 0b9b1d8b..e440af65 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,39 +1,55 @@
-using ForwardDiff
-
-println("Testing Partials...")
-t = @elapsed include("PartialsTest.jl")
-println("done (took $t seconds).")
-
-println("Testing Dual...")
-t = @elapsed include("DualTest.jl")
-println("done (took $t seconds).")
-
-println("Testing derivative functionality...")
-t = @elapsed include("DerivativeTest.jl")
-println("done (took $t seconds).")
-
-println("Testing gradient functionality...")
-t = @elapsed include("GradientTest.jl")
-println("done (took $t seconds).")
-
-println("Testing jacobian functionality...")
-t = @elapsed include("JacobianTest.jl")
-println("done (took $t seconds).")
-
-println("Testing hessian functionality...")
-t = @elapsed include("HessianTest.jl")
-println("done (took $t seconds).")
-
-println("Testing perturbation confusion functionality...")
-t = @elapsed include("ConfusionTest.jl")
-println("done (took $t seconds).")
-
-println("Testing miscellaneous functionality...")
-t = @elapsed include("MiscTest.jl")
-println("done (took $t seconds).")
-
-if VERSION >= v"1.5-"
-    println("Testing allocations...")
-    t = @elapsed include("AllocationsTest.jl")
-    println("done (took $t seconds).")
+using ForwardDiff, Test, Random
+
+SEED = trunc(Int, time())
+println("##### Random.seed!($SEED), on VERSION == $VERSION")
+Random.seed!(SEED)
+
+@testset "ForwardDiff.jl" begin
+    t0 = time()
+    @testset "Partials" begin
+        println("##### Testing Partials...")
+        t = @elapsed include("PartialsTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Dual" begin
+        println("##### Testing Dual...")
+        t = @elapsed include("DualTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Derivatives" begin
+        println("##### Testing derivative functionality...")
+        t = @elapsed include("DerivativeTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Gradients" begin
+        println("##### Testing gradient functionality...")
+        t = @elapsed include("GradientTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Jacobians" begin
+        println("##### Testing jacobian functionality...")
+        t = @elapsed include("JacobianTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Hessians" begin
+        println("##### Testing hessian functionality...")
+        t = @elapsed include("HessianTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Perturbation Confusion" begin
+        println("##### Testing perturbation confusion functionality...")
+        t = @elapsed include("ConfusionTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Miscellaneous" begin
+        println("##### Testing miscellaneous functionality...")
+        t = @elapsed include("MiscTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    @testset "Allocations" begin
+        println("##### Testing allocations...")
+        t = @elapsed include("AllocationsTest.jl")
+        println("##### done (took $t seconds).")
+    end
+    println("##### Running all ForwardDiff tests took $(time() - t0) seconds.")
 end