diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml
index 81caf5536..b479b3970 100644
--- a/.github/workflows/Test.yml
+++ b/.github/workflows/Test.yml
@@ -98,6 +98,7 @@ jobs:
           - GTPSA
           - HyperHessians
           - Mooncake
+          - Piste
           - PolyesterForwardDiff
           - ReverseDiff
           - SparsityDetector
diff --git a/DifferentiationInterface/CHANGELOG.md b/DifferentiationInterface/CHANGELOG.md
index 575d9cff7..8f064eff0 100644
--- a/DifferentiationInterface/CHANGELOG.md
+++ b/DifferentiationInterface/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Support for [Piste.jl](https://github.com/EvoArt/Piste.jl) with `AutoPiste`, in forward and reverse mode
 - New operators `value_gradient_and_hvp` and `value_gradient_and_hvp!`, which also return the primal value ([#1086](https://github.com/JuliaDiff/DifferentiationInterface.jl/pull/1086))
 
 ### Deprecated
diff --git a/DifferentiationInterface/Project.toml b/DifferentiationInterface/Project.toml
index 594a3f743..ae107ac3c 100644
--- a/DifferentiationInterface/Project.toml
+++ b/DifferentiationInterface/Project.toml
@@ -22,6 +22,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
 Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
+Piste = "f69a1a94-130c-4fa4-b8d2-e36d656cbe66"
 PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -44,6 +45,7 @@ DifferentiationInterfaceGPUArraysCoreExt = ["GPUArraysCore", "Adapt"]
 DifferentiationInterfaceHyperHessiansExt = "HyperHessians"
 DifferentiationInterfaceGTPSAExt = "GTPSA"
 DifferentiationInterfaceMooncakeExt = "Mooncake"
+DifferentiationInterfacePisteExt = "Piste"
 DifferentiationInterfacePolyesterForwardDiffExt = [
     "PolyesterForwardDiff",
     "ForwardDiff",
@@ -60,7 +62,7 @@ DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"]
 
 [compat]
 Adapt = "4.5.0"
-ADTypes = "1.23.0"
+ADTypes = "1.25.0"
 ChainRulesCore = "1.23.0"
 DiffResults = "1.1.0"
 Diffractor = "=0.2.6"
@@ -75,6 +77,7 @@ GTPSA = "1.4.0"
 HyperHessians = "0.3"
 LinearAlgebra = "1"
 Mooncake = "0.5.25"
+Piste = "0.4"
 PolyesterForwardDiff = "0.1.2"
 ReverseDiff = "1.15.1"
 SparseArrays = "1"
diff --git a/DifferentiationInterface/README.md b/DifferentiationInterface/README.md
index b6c2c6a59..d1b9f8a59 100644
--- a/DifferentiationInterface/README.md
+++ b/DifferentiationInterface/README.md
@@ -38,6 +38,7 @@ We support the following backends defined by [ADTypes.jl](https://github.com/Sci
 - [HyperHessians.jl](https://github.com/KristofferC/HyperHessians.jl)
 - [GTPSA.jl](https://github.com/bmad-sim/GTPSA.jl)
 - [Mooncake.jl](https://github.com/chalk-lab/Mooncake.jl)
+- [Piste.jl](https://github.com/EvoArt/Piste.jl)
 - [PolyesterForwardDiff.jl](https://github.com/JuliaDiff/PolyesterForwardDiff.jl)
 - [ReverseDiff.jl](https://github.com/JuliaDiff/ReverseDiff.jl)
 - [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl)
diff --git a/DifferentiationInterface/docs/src/explanation/backends.md b/DifferentiationInterface/docs/src/explanation/backends.md
index 36193920d..c21da3e17 100644
--- a/DifferentiationInterface/docs/src/explanation/backends.md
+++ b/DifferentiationInterface/docs/src/explanation/backends.md
@@ -14,6 +14,7 @@ We support the following dense backend choices from [ADTypes.jl](https://github.
 - [`AutoGTPSA`](@extref ADTypes.AutoGTPSA)
 - [`AutoHyperHessians`](@extref ADTypes.AutoHyperHessians)
 - [`AutoMooncake`](@extref ADTypes.AutoMooncake) and [`AutoMooncakeForward`](@extref ADTypes.AutoMooncake) (the latter is experimental)
+- [`AutoPiste`](@extref ADTypes.AutoPiste)
 - [`AutoPolyesterForwardDiff`](@extref ADTypes.AutoPolyesterForwardDiff)
 - [`AutoReverseDiff`](@extref ADTypes.AutoReverseDiff)
 - [`AutoSymbolics`](@extref ADTypes.AutoSymbolics)
@@ -52,6 +53,8 @@ In practice, many AD backends have custom implementations for high-level operato
     | `AutoHyperHessians`        | ❌    | ❌    | ❌     | ❌      | ❌     | ✅      | ✅     | ✅      |
     | `AutoMooncake`             | ❌    | ✅    | ❌     | ❌      | ❌     | ❌      | ❌     | ❌      |
     | `AutoMooncakeForward`      | ✅    | ❌    | ❌     | ❌      | ❌     | ❌      | ❌     | ❌      |
+    | `AutoPiste` (forward)      | ✅    | ❌    | ❌     | ✅      | ❌     | ❌      | ❌     | ❌      |
+    | `AutoPiste` (reverse)      | ❌    | ✅    | ❌     | ✅      | ❌     | ❌      | ❌     | ❌      |
     | `AutoPolyesterForwardDiff` | 🔀    | ❌    | 🔀     | ✅      | ✅     | 🔀      | 🔀     | 🔀      |
     | `AutoReverseDiff`          | ❌    | 🔀    | ❌     | ✅      | ✅     | ✅      | ❌     | ❌      |
     | `AutoSymbolics`            | ✅    | ❌    | ✅     | ✅      | ✅     | ✅      | ✅     | ✅      |
@@ -74,6 +77,8 @@ Moreover, each context type is supported by a specific subset of backends:
 | `AutoHyperHessians`        | ✅                  | ✅               |
 | `AutoMooncake`             | ✅                  | ✅               |
 | `AutoMooncakeForward`      | ✅                  | ✅               |
+| `AutoPiste` (forward)      | ✅                  | ✅               |
+| `AutoPiste` (reverse)      | ✅                  | ✅               |
 | `AutoPolyesterForwardDiff` | ✅                  | ✅               |
 | `AutoReverseDiff`          | ✅                  | ❌               |
 | `AutoSymbolics`            | ✅                  | ✅               |
@@ -197,6 +202,13 @@ For all operators, preparation generates an [executable function](https://docs.s
 
 For `pullback`, preparation [builds the reverse rule](https://chalk-lab.github.io/Mooncake.jl/stable/understanding_mooncake/rule_system/) of the function.
 
+### Piste
+
+In forward mode, `pushforward` propagates `Piste.Dual` numbers carrying a whole batch of tangents, and `gradient` (for a `Vector` input) calls Piste's own chunked driver, with its workspace allocated during preparation.
+In reverse mode, `pullback` and `gradient` record the function on a `Piste.Tape` and sweep it backwards, reusing a workspace allocated during preparation.
+
+Piste is designed to survive `juliac --trim`, and prepared `gradient` calls through DifferentiationInterface do too, provided that the `chunksize` of a forward-mode backend is set explicitly.
+
 ### Tracker
 
 We implement `pullback` based on `Tracker.back`.
diff --git a/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/DifferentiationInterfacePisteExt.jl b/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/DifferentiationInterfacePisteExt.jl
new file mode 100644
index 000000000..851ffa4bc
--- /dev/null
+++ b/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/DifferentiationInterfacePisteExt.jl
@@ -0,0 +1,24 @@
+module DifferentiationInterfacePisteExt
+
+using ADTypes: AutoPiste, ForwardMode, ReverseMode
+import DifferentiationInterface as DI
+using Piste:
+    Piste,
+    Dual,
+    GradientWorkspace,
+    ReverseWorkspace,
+    TapedReal,
+    gradient!,
+    pickchunk,
+    rev_value_and_pullback,
+    rev_value_and_pullback!,
+    value,
+    value_and_gradient!
+
+DI.check_available(::AutoPiste) = true
+DI.inplace_support(::AutoPiste) = DI.InPlaceSupported()
+
+include("forward.jl")
+include("reverse.jl")
+
+end # module
diff --git a/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/forward.jl b/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/forward.jl
new file mode 100644
index 000000000..0a6f15339
--- /dev/null
+++ b/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/forward.jl
@@ -0,0 +1,339 @@
+## Forward mode: `Piste.Dual{B,T}` carries `B` directional derivatives at once
+
+const AutoPisteForward{chunksize} = AutoPiste{ForwardMode, chunksize}
+
+function DI.pick_batchsize(::AutoPisteForward{nothing}, N::Integer)
+    return DI.BatchSizeSettings(DI.reasonable_batchsize(N, 8), N)
+end
+
+function DI.pick_batchsize(::AutoPisteForward{chunksize}, N::Integer) where {chunksize}
+    return DI.BatchSizeSettings(chunksize, N)
+end
+
+function DI.threshold_batchsize(
+        ::AutoPisteForward{chunksize1}, chunksize2::Integer
+    ) where {chunksize1}
+    chunksize = isnothing(chunksize1) ? nothing : min(chunksize1, chunksize2)
+    return AutoPiste(; mode = ForwardMode(), chunksize)
+end
+
+### Utilities
+
+dual_eltype(x) = float(eltype(x))
+
+make_dual(::Type{T}, x::Number, tx::NTuple{B}) where {T, B} = Dual{B, T}(T(x), map(T, tx))
+
+function make_dual(::Type{T}, x::AbstractArray, tx::NTuple{B}) where {T, B}
+    return map((xi, txi...) -> Dual{B, T}(T(xi), map(T, txi)), x, tx...)
+end
+
+partial(y::Dual, b::Integer) = Piste.partials(y, b)
+partial(y::Number, ::Integer) = zero(y)  # an output that does not depend on `x`
+
+mypartials(::Val{B}, ydual::Number) where {B} = ntuple(b -> partial(ydual, b), Val(B))
+
+function mypartials(::Val{B}, ydual::AbstractArray) where {B}
+    return ntuple(b -> map(yi -> partial(yi, b), ydual), Val(B))
+end
+
+function mypartials!(ty::NTuple{B}, ydual::AbstractArray) where {B}
+    for b in 1:B
+        map!(yi -> partial(yi, b), ty[b], ydual)
+    end
+    return ty
+end
+
+myvalue(ydual::Number) = value(ydual)
+myvalue(ydual::AbstractArray) = map(value, ydual)
+myvalue!(y::AbstractArray, ydual::AbstractArray) = map!(value, y, ydual)
+
+translate(::Type{D}, c::DI.Cache) where {D} = DI.recursive_similar(DI.unwrap(c), D)
+translate(::Type{D}, c::DI.Context) where {D} = DI.unwrap(c)
+
+function translate(::Type{D}, contexts::NTuple{C, DI.Context}) where {D, C}
+    return map(c -> translate(D, c), contexts)
+end
+
+### Pushforward, one argument
+
+struct PisteOneArgPushforwardPrep{SIG} <: DI.PushforwardPrep{SIG}
+    _sig::Val{SIG}
+end
+
+function DI.prepare_pushforward_nokwarg(
+        strict::Val,
+        f::F,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    _sig = DI.signature(f, backend, x, tx, contexts...; strict)
+    return PisteOneArgPushforwardPrep(_sig)
+end
+
+function compute_ydual(
+        f::F, x, tx::NTuple{B}, contexts::Vararg{DI.Context, C}
+    ) where {F, B, C}
+    T = dual_eltype(x)
+    xdual = make_dual(T, x, tx)
+    contexts_dual = translate(Dual{B, T}, contexts)
+    return f(xdual, contexts_dual...)
+end
+
+function DI.value_and_pushforward(
+        f::F,
+        prep::PisteOneArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple{B},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, B, C}
+    DI.check_prep(f, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual(f, x, tx, contexts...)
+    return myvalue(ydual), mypartials(Val(B), ydual)
+end
+
+function DI.value_and_pushforward!(
+        f::F,
+        ty::NTuple,
+        prep::PisteOneArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual(f, x, tx, contexts...)
+    mypartials!(ty, ydual)
+    return myvalue(ydual), ty
+end
+
+function DI.pushforward(
+        f::F,
+        prep::PisteOneArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple{B},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, B, C}
+    DI.check_prep(f, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual(f, x, tx, contexts...)
+    return mypartials(Val(B), ydual)
+end
+
+function DI.pushforward!(
+        f::F,
+        ty::NTuple,
+        prep::PisteOneArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual(f, x, tx, contexts...)
+    return mypartials!(ty, ydual)
+end
+
+### Pushforward, two arguments
+
+struct PisteTwoArgPushforwardPrep{SIG, YD} <: DI.PushforwardPrep{SIG}
+    _sig::Val{SIG}
+    ydual::YD
+end
+
+function DI.prepare_pushforward_nokwarg(
+        strict::Val,
+        f!::F,
+        y,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple{B},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, B, C}
+    _sig = DI.signature(f!, y, backend, x, tx, contexts...; strict)
+    ydual = similar(y, Dual{B, dual_eltype(x)})
+    return PisteTwoArgPushforwardPrep(_sig, ydual)
+end
+
+function compute_ydual!(
+        f!::F,
+        prep::PisteTwoArgPushforwardPrep,
+        x,
+        tx::NTuple{B},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, B, C}
+    (; ydual) = prep
+    T = dual_eltype(x)
+    xdual = make_dual(T, x, tx)
+    contexts_dual = translate(Dual{B, T}, contexts)
+    f!(ydual, xdual, contexts_dual...)
+    return ydual
+end
+
+function DI.value_and_pushforward(
+        f!::F,
+        y,
+        prep::PisteTwoArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple{B},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, B, C}
+    DI.check_prep(f!, y, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual!(f!, prep, x, tx, contexts...)
+    myvalue!(y, ydual)
+    return y, mypartials(Val(B), ydual)
+end
+
+function DI.value_and_pushforward!(
+        f!::F,
+        y,
+        ty::NTuple,
+        prep::PisteTwoArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f!, y, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual!(f!, prep, x, tx, contexts...)
+    myvalue!(y, ydual)
+    mypartials!(ty, ydual)
+    return y, ty
+end
+
+function DI.pushforward(
+        f!::F,
+        y,
+        prep::PisteTwoArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple{B},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, B, C}
+    DI.check_prep(f!, y, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual!(f!, prep, x, tx, contexts...)
+    myvalue!(y, ydual)
+    return mypartials(Val(B), ydual)
+end
+
+function DI.pushforward!(
+        f!::F,
+        y,
+        ty::NTuple,
+        prep::PisteTwoArgPushforwardPrep,
+        backend::AutoPisteForward,
+        x,
+        tx::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f!, y, prep, backend, x, tx, contexts...)
+    ydual = compute_ydual!(f!, prep, x, tx, contexts...)
+    myvalue!(y, ydual)
+    return mypartials!(ty, ydual)
+end
+
+### Gradient: Piste's own chunked driver, with a reusable workspace
+
+# `Piste.gradient!` needs `x` and the gradient to be `Vector`s of one float type.
+# Any other input falls back to DI's default (batched pushforwards, above).
+
+struct PisteForwardGradientPrep{SIG, N, T} <: DI.GradientPrep{SIG}
+    _sig::Val{SIG}
+    ws::GradientWorkspace{N, T}
+end
+
+# A fixed `chunksize` stays a compile-time constant, which is what lets this
+# path compile under `juliac --trim`. With `nothing` the width is picked from
+# `length(x)`, which is fine in the JIT but not statically known.
+chunk(::AutoPisteForward{nothing}, K::Int) = pickchunk(K)
+chunk(::AutoPisteForward{chunksize}, ::Int) where {chunksize} = Val(chunksize)
+
+function DI.prepare_gradient_nokwarg(
+        strict::Val,
+        f::F,
+        backend::AutoPisteForward,
+        x::Vector{T},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, T <: AbstractFloat, C}
+    _sig = DI.signature(f, backend, x, contexts...; strict)
+    ws = GradientWorkspace(x, chunk(backend, length(x)))
+    return PisteForwardGradientPrep(_sig, ws)
+end
+
+function fix_contexts(
+        f::F, ::GradientWorkspace{N, T}, contexts::Vararg{DI.Context, C}
+    ) where {F, N, T, C}
+    return DI.fix_tail(f, translate(Dual{N, T}, contexts)...)
+end
+
+function DI.value_and_gradient!(
+        f::F,
+        grad::Vector{T},
+        prep::PisteForwardGradientPrep{SIG, N, T},
+        backend::AutoPisteForward,
+        x::Vector{T},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, SIG, N, T, C}
+    DI.check_prep(f, prep, backend, x, contexts...)
+    fc = fix_contexts(f, prep.ws, contexts...)
+    return value_and_gradient!(grad, fc, x, Val(N), prep.ws)
+end
+
+function DI.gradient!(
+        f::F,
+        grad::Vector{T},
+        prep::PisteForwardGradientPrep{SIG, N, T},
+        backend::AutoPisteForward,
+        x::Vector{T},
+        contexts::Vararg{DI.Context, C},
+    ) where {F, SIG, N, T, C}
+    DI.check_prep(f, prep, backend, x, contexts...)
+    fc = fix_contexts(f, prep.ws, contexts...)
+    return gradient!(grad, fc, x, Val(N), prep.ws)
+end
+
+function DI.value_and_gradient!(
+        f::F,
+        grad,
+        prep::PisteForwardGradientPrep,
+        backend::AutoPisteForward,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    y, new_grad = DI.value_and_gradient(f, prep, backend, x, contexts...)
+    return y, copyto!(grad, new_grad)
+end
+
+function DI.gradient!(
+        f::F,
+        grad,
+        prep::PisteForwardGradientPrep,
+        backend::AutoPisteForward,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return copyto!(grad, DI.gradient(f, prep, backend, x, contexts...))
+end
+
+function DI.value_and_gradient(
+        f::F,
+        prep::PisteForwardGradientPrep,
+        backend::AutoPisteForward,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return DI.value_and_gradient!(f, similar(x), prep, backend, x, contexts...)
+end
+
+function DI.gradient(
+        f::F,
+        prep::PisteForwardGradientPrep,
+        backend::AutoPisteForward,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return DI.gradient!(f, similar(x), prep, backend, x, contexts...)
+end
diff --git a/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/reverse.jl b/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/reverse.jl
new file mode 100644
index 000000000..6f97f431d
--- /dev/null
+++ b/DifferentiationInterface/ext/DifferentiationInterfacePisteExt/reverse.jl
@@ -0,0 +1,258 @@
+## Reverse mode: one recording of `f` on a `Piste.Tape`, one sweep per cotangent
+
+const AutoPisteReverse{chunksize} = AutoPiste{ReverseMode, chunksize}
+
+### Utilities
+
+translate_rev(c::DI.Cache) = DI.recursive_similar(DI.unwrap(c), TapedReal)
+translate_rev(c::DI.Context) = DI.unwrap(c)
+
+function fix_contexts_rev(f::F, contexts::Vararg{DI.Context, C}) where {F, C}
+    return DI.fix_tail(f, map(translate_rev, contexts)...)
+end
+
+# An in-place `f!(y, x)` is recorded through an out-of-place closure that
+# writes into a tracked copy of `y`. Starting from detached zeros means an
+# entry `f!` never writes is a constant, not uninitialized memory.
+struct InPlaceAsOutOfPlace{F, Y}
+    f!::F
+    y::Y
+end
+
+function (g::InPlaceAsOutOfPlace)(x, contexts::Vararg{Any, C}) where {C}
+    ytape = similar(g.y, TapedReal)
+    fill!(ytape, zero(TapedReal))
+    g.f!(ytape, x, contexts...)
+    return ytape
+end
+
+workspace(x::Number) = ReverseWorkspace(1)
+workspace(x::AbstractArray) = ReverseWorkspace(length(x))
+
+vjp(fc::F, ws, x::Number, dy) where {F} = rev_value_and_pullback(fc, x, dy, ws)
+
+function vjp(fc::F, ws, x::AbstractArray, dy) where {F}
+    return rev_value_and_pullback!(similar(x, float(eltype(x))), fc, x, dy, ws)
+end
+
+vjp!(dx, fc::F, ws, x::AbstractArray, dy) where {F} = rev_value_and_pullback!(dx, fc, x, dy, ws)
+
+# One recording per cotangent: the first also supplies the primal.
+function value_and_vjps(fc::F, ws, x, ty::NTuple{B}) where {F, B}
+    y, dx1 = vjp(fc, ws, x, first(ty))
+    rest = map(dy -> last(vjp(fc, ws, x, dy)), Base.tail(ty))
+    return y, (dx1, rest...)
+end
+
+function value_and_vjps!(tx::NTuple{B}, fc::F, ws, x, ty::NTuple{B}) where {F, B}
+    y, _ = vjp!(first(tx), fc, ws, x, first(ty))
+    for b in 2:B
+        vjp!(tx[b], fc, ws, x, ty[b])
+    end
+    return y, tx
+end
+
+### Pullback, one argument
+
+struct PistePullbackPrep{SIG} <: DI.PullbackPrep{SIG}
+    _sig::Val{SIG}
+    ws::ReverseWorkspace
+end
+
+function DI.prepare_pullback_nokwarg(
+        strict::Val,
+        f::F,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    _sig = DI.signature(f, backend, x, ty, contexts...; strict)
+    return PistePullbackPrep(_sig, workspace(x))
+end
+
+function DI.value_and_pullback(
+        f::F,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f, prep, backend, x, ty, contexts...)
+    fc = fix_contexts_rev(f, contexts...)
+    return value_and_vjps(fc, prep.ws, x, ty)
+end
+
+function DI.value_and_pullback!(
+        f::F,
+        tx::NTuple,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f, prep, backend, x, ty, contexts...)
+    fc = fix_contexts_rev(f, contexts...)
+    return value_and_vjps!(tx, fc, prep.ws, x, ty)
+end
+
+function DI.pullback(
+        f::F,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return last(DI.value_and_pullback(f, prep, backend, x, ty, contexts...))
+end
+
+function DI.pullback!(
+        f::F,
+        tx::NTuple,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return last(DI.value_and_pullback!(f, tx, prep, backend, x, ty, contexts...))
+end
+
+### Pullback, two arguments
+
+function DI.prepare_pullback_nokwarg(
+        strict::Val,
+        f!::F,
+        y,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    _sig = DI.signature(f!, y, backend, x, ty, contexts...; strict)
+    return PistePullbackPrep(_sig, workspace(x))
+end
+
+function DI.value_and_pullback(
+        f!::F,
+        y,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f!, y, prep, backend, x, ty, contexts...)
+    fc = fix_contexts_rev(InPlaceAsOutOfPlace(f!, y), contexts...)
+    yval, tx = value_and_vjps(fc, prep.ws, x, ty)
+    copyto!(y, yval)
+    return y, tx
+end
+
+function DI.value_and_pullback!(
+        f!::F,
+        y,
+        tx::NTuple,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f!, y, prep, backend, x, ty, contexts...)
+    fc = fix_contexts_rev(InPlaceAsOutOfPlace(f!, y), contexts...)
+    yval, _ = value_and_vjps!(tx, fc, prep.ws, x, ty)
+    copyto!(y, yval)
+    return y, tx
+end
+
+function DI.pullback(
+        f!::F,
+        y,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return last(DI.value_and_pullback(f!, y, prep, backend, x, ty, contexts...))
+end
+
+function DI.pullback!(
+        f!::F,
+        y,
+        tx::NTuple,
+        prep::PistePullbackPrep,
+        backend::AutoPisteReverse,
+        x,
+        ty::NTuple,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return last(DI.value_and_pullback!(f!, y, tx, prep, backend, x, ty, contexts...))
+end
+
+### Gradient: a single pullback with cotangent 1, reusing the workspace
+
+struct PisteReverseGradientPrep{SIG} <: DI.GradientPrep{SIG}
+    _sig::Val{SIG}
+    ws::ReverseWorkspace
+end
+
+function DI.prepare_gradient_nokwarg(
+        strict::Val,
+        f::F,
+        backend::AutoPisteReverse,
+        x::AbstractArray,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    _sig = DI.signature(f, backend, x, contexts...; strict)
+    return PisteReverseGradientPrep(_sig, workspace(x))
+end
+
+function DI.value_and_gradient!(
+        f::F,
+        grad,
+        prep::PisteReverseGradientPrep,
+        backend::AutoPisteReverse,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    DI.check_prep(f, prep, backend, x, contexts...)
+    fc = fix_contexts_rev(f, contexts...)
+    return vjp!(grad, fc, prep.ws, x, 1.0)
+end
+
+function DI.gradient!(
+        f::F,
+        grad,
+        prep::PisteReverseGradientPrep,
+        backend::AutoPisteReverse,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return last(DI.value_and_gradient!(f, grad, prep, backend, x, contexts...))
+end
+
+function DI.value_and_gradient(
+        f::F,
+        prep::PisteReverseGradientPrep,
+        backend::AutoPisteReverse,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    grad = similar(x, float(eltype(x)))
+    return DI.value_and_gradient!(f, grad, prep, backend, x, contexts...)
+end
+
+function DI.gradient(
+        f::F,
+        prep::PisteReverseGradientPrep,
+        backend::AutoPisteReverse,
+        x,
+        contexts::Vararg{DI.Context, C},
+    ) where {F, C}
+    return last(DI.value_and_gradient(f, prep, backend, x, contexts...))
+end
diff --git a/DifferentiationInterface/test/Back/Piste/Project.toml b/DifferentiationInterface/test/Back/Piste/Project.toml
new file mode 100644
index 000000000..ac93c6936
--- /dev/null
+++ b/DifferentiationInterface/test/Back/Piste/Project.toml
@@ -0,0 +1,13 @@
+[deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+DifferentiationInterfaceTest = "a82114a7-5aa3-49a8-9643-716bb13727a3"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+Piste = "f69a1a94-130c-4fa4-b8d2-e36d656cbe66"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[sources]
+DifferentiationInterface = { path = "../../.." }
diff --git a/DifferentiationInterface/test/Back/Piste/test.jl b/DifferentiationInterface/test/Back/Piste/test.jl
new file mode 100644
index 000000000..227f515a4
--- /dev/null
+++ b/DifferentiationInterface/test/Back/Piste/test.jl
@@ -0,0 +1,63 @@
+include("../../testutils.jl")
+
+using ADTypes: ADTypes
+using DifferentiationInterface, DifferentiationInterfaceTest
+using Piste: Piste
+using Test
+import JET
+
+using ExplicitImports
+check_no_implicit_imports(DifferentiationInterface)
+
+backends = [
+    AutoPiste(),
+    AutoPiste(; chunksize = 2),
+    AutoPiste(; mode = ADTypes.ReverseMode()),
+]
+
+for backend in backends
+    @test check_available(backend)
+    @test check_inplace(backend)
+    test_counterparts(backend)
+end
+
+@testset "Dense" begin
+    test_differentiation(
+        backends,
+        default_scenarios(; include_constantified = true);
+        excluded = SECOND_ORDER,
+        logging = LOGGING,
+    )
+
+    test_differentiation(
+        backends,
+        default_scenarios(;
+            include_normal = false,
+            include_batchified = false,
+            include_cachified = true,
+            include_constantorcachified = true,
+            use_tuples = true,
+            include_smaller = true,
+        );
+        excluded = SECOND_ORDER,
+        logging = LOGGING,
+    )
+
+    test_differentiation(
+        backends;
+        correctness = false,
+        type_stability = safetypestab(:prepared),
+        excluded = SECOND_ORDER,
+        logging = LOGGING,
+    )
+end
+
+@testset "Sparse" begin
+    test_differentiation(
+        MyAutoSparse.(backends),
+        sparse_scenarios(; include_constantified = true);
+        sparsity = true,
+        excluded = SECOND_ORDER,
+        logging = LOGGING,
+    )
+end