diff --git a/CMakeLists.txt b/CMakeLists.txt
index 33f7861..7e35fd0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -43,7 +43,7 @@ if(EXISTS "${KIMODO_GGML_SOURCE_DIR}/CMakeLists.txt")
   endif()
   target_link_libraries(kmd-inspect PRIVATE ggml)
   target_compile_definitions(kmd-inspect PRIVATE KIMODO_HAVE_GGML=1)
-  target_sources(kimodo PRIVATE src/ggml_weights.cpp src/denoiser.cpp)
+  target_sources(kimodo PRIVATE src/ggml_weights.cpp src/denoiser.cpp src/sequence.cpp)
   target_sources(kimodo PRIVATE src/llm_tokenizer.cpp src/llm_text_encoder.cpp)
   add_executable(kmd-sample-fixture src/sample_fixture.cpp)
   target_sources(kmd-sample-fixture PRIVATE
@@ -72,7 +72,7 @@ if(EXISTS "${KIMODO_GGML_SOURCE_DIR}/CMakeLists.txt")
   add_executable(kimodo-llm-embedding-parity tests/llm_embedding_parity.cpp)
   target_link_libraries(kimodo-llm-embedding-parity PRIVATE ggml ggml-vulkan)
   add_executable(kimodo-multiprompt-fixture-parity tests/multiprompt_fixture_parity.cpp)
-  target_sources(kimodo-multiprompt-fixture-parity PRIVATE src/gguf.cpp src/ggml_weights.cpp src/denoiser.cpp src/motion_rep.cpp src/diffusion.cpp)
+  target_sources(kimodo-multiprompt-fixture-parity PRIVATE src/gguf.cpp src/ggml_weights.cpp src/denoiser.cpp src/sequence.cpp src/motion_rep.cpp src/diffusion.cpp)
   target_include_directories(kimodo-multiprompt-fixture-parity PRIVATE src)
   target_link_libraries(kimodo-multiprompt-fixture-parity PRIVATE ggml ggml-vulkan)
   target_compile_definitions(kimodo-multiprompt-fixture-parity PRIVATE KIMODO_HAVE_GGML=1 KIMODO_HAVE_GGML_VULKAN=1)
diff --git a/reference/dump_kimodo_multiprompt_reference.py b/reference/dump_kimodo_multiprompt_reference.py
index f2c5351..f479cc6 100644
--- a/reference/dump_kimodo_multiprompt_reference.py
+++ b/reference/dump_kimodo_multiprompt_reference.py
@@ -106,6 +106,9 @@ def main() -> None:
                                  text_encoder=encoder, return_resolved_name=True)
     calls: list[dict[str, list[torch.Tensor] | torch.Tensor]] = []
     inverse_inputs: list[torch.Tensor] = []
+    # The first inverse call decodes the tail of the preceding segment.  Keep
+    # its exact FK products as a transition-encoder oracle for native ports.
+    inverse_outputs: list[dict[str, torch.Tensor]] = []
     root_calls: list[tuple[tuple[torch.Tensor, ...], torch.Tensor]] = []
     body_calls: list[tuple[tuple[torch.Tensor, ...], torch.Tensor]] = []
     original_step = model.denoising_step
@@ -143,7 +146,13 @@ def main() -> None:
 
     def capture_inverse(motion: torch.Tensor, *values: Any, **kwargs: Any):
         inverse_inputs.append(motion.detach().clone())
-        return original_inverse(motion, *values, **kwargs)
+        result = original_inverse(motion, *values, **kwargs)
+        inverse_outputs.append({
+            key: value.detach().clone()
+            for key, value in result.items()
+            if isinstance(value, torch.Tensor)
+        })
+        return result
 
     model.motion_rep.inverse = capture_inverse
     model.denoising_step = capture_step
@@ -164,6 +173,14 @@ def main() -> None:
         raise RuntimeError(f"expected {len(args.prompt)} segment trajectories, got {len(calls)}")
 
     arrays: dict[str, np.ndarray] = {"stitched_motion_rep": as_f32(inverse_inputs[-1])}
+    # `_multiprompt` calls inverse on the preceding tail before creating the
+    # continuation constraints.  This is deliberately separate from the final
+    # output decode below, whose first five frames have already been blended.
+    if len(inverse_outputs) < 2:
+        raise RuntimeError("expected transition-tail and final inverse calls")
+    arrays["transition_source_motion"] = as_f32(inverse_inputs[0])
+    for key, value in inverse_outputs[0].items():
+        arrays["transition_source_" + key] = as_f32(value)
     for stage, stage_calls in (("root", root_calls), ("body", body_calls)):
         values, stage_output = stage_calls[0]
         arrays[f"{stage}_output"] = as_f32(stage_output)
diff --git a/scripts/convert_motion_to_gguf.py b/scripts/convert_motion_to_gguf.py
index 685f7a7..6ee6e07 100755
--- a/scripts/convert_motion_to_gguf.py
+++ b/scripts/convert_motion_to_gguf.py
@@ -19,7 +19,7 @@ from pathlib import Path
 
 ALIGNMENT = 32
 GGUF_MAGIC, GGUF_VERSION, GGML_TYPE_F32 = 0x46554747, 3, 0
-TYPE_UINT64, TYPE_STRING = 10, 8
+TYPE_UINT64, TYPE_STRING, TYPE_FLOAT32 = 10, 8, 6
 TYPE_UINT32 = 4
 
 @dataclass(frozen=True)
@@ -155,6 +155,9 @@ def metadata_uint(key: str, value: int) -> bytes:
 def metadata_uint32(key: str, value: int) -> bytes:
     return string(key) + struct.pack("<I", TYPE_UINT32) + struct.pack("<I", value)
 
+def metadata_float32(key: str, value: float) -> bytes:
+    return string(key) + struct.pack("<I", TYPE_FLOAT32) + struct.pack("<f", value)
+
 def tensor_info(tensor: Tensor, offset: int) -> bytes:
     # GGML stores dim 0 as the contiguous dimension. PyTorch F32 storage is
     # row-major, so reverse dimensions without changing the underlying bytes.
@@ -215,6 +218,8 @@ def main() -> None:
         metadata_uint("kimodo.num_text_tokens", 50),
         metadata_uint("kimodo.base_diffusion_steps", 1000),
         metadata_uint("kimodo.fps", 30),
+        # kimodo.motion_rep.stats.Stats uses sqrt(std**2 + eps).
+        metadata_float32("kimodo.normalization_epsilon", 1.0e-5),
     ]
     offsets, cursor = [], 0
     for tensor in tensors:
diff --git a/src/denoiser.hpp b/src/denoiser.hpp
index b3d3f1b..504b16f 100644
--- a/src/denoiser.hpp
+++ b/src/denoiser.hpp
@@ -9,6 +9,22 @@
 namespace kimodo::detail {
 class ggml_motion_weights;
 
+// One sequence segment with an already-encoded text condition and caller-
+// supplied initial noise. Continuation noise includes its transition prefix.
+struct sampled_sequence_segment {
+    std::span<const float> embedding;
+    std::span<const float> initial_noise;
+    std::size_t frames;
+};
+
+struct sequence_transition {
+    std::vector<float> observed;
+    std::vector<float> observed_mask;
+    float first_heading = 0.F;
+    float origin_x = 0.F;
+    float origin_z = 0.F;
+};
+
 // F32 TransformerEncoderBlock. Inputs/outputs are row-major [B,T,D], while
 // the implementation creates GGML [D,T,B] views over the same byte order.
 std::expected<std::vector<float>, std::string> run_motion_transformer(
@@ -46,4 +62,17 @@ std::expected<std::vector<float>, std::string> sample_motion_from_noise_conditio
     std::span<const float> embedding, std::span<const float> observed,
     std::span<const float> observed_mask, float first_heading, std::size_t frames,
     unsigned steps, float text_weight, float constraint_weight);
+
+// End-to-end upstream `_multiprompt` orchestration.  DDIM operates in
+// normalized motion space; the returned joined representation is raw so its
+// translated roots and blended tail preserve upstream semantics.
+std::expected<std::vector<float>, std::string> sample_motion_sequence_from_noise(
+    const ggml_motion_weights &weights, std::span<const sampled_sequence_segment> segments,
+    unsigned transition_frames, unsigned steps, float text_weight, float constraint_weight);
+
+// Build the exact condition consumed by the next `_multiprompt` DDIM run.
+// Exposed for the raw fixture test as well as the runtime orchestrator.
+std::expected<sequence_transition, std::string> prepare_sequence_transition(
+    const ggml_motion_weights &weights, std::span<const float> previous,
+    std::size_t continuation_frames, unsigned transition_frames);
 }
diff --git a/src/model.cpp b/src/model.cpp
index 45af2a3..70de857 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -9,7 +9,6 @@
 
 #include <cmath>
 #include <algorithm>
-#include <array>
 #include <random>
 
 namespace kimodo {
@@ -110,80 +109,38 @@ std::expected<motion_data, std::string> model::generate_text_sequence(
         if (!loaded) return std::unexpected(loaded.error());
         impl_->weights = std::move(*loaded);
     }
-    auto gm=impl_->weights->f32_values("stats.global_root.mean"), gs=impl_->weights->f32_values("stats.global_root.std");
     auto bm=impl_->weights->f32_values("stats.body.mean"), bs=impl_->weights->f32_values("stats.body.std");
+    auto gm=impl_->weights->f32_values("stats.global_root.mean"), gs=impl_->weights->f32_values("stats.global_root.std");
     if (!gm || !gs || !bm || !bs) return std::unexpected("motion GGUF lacks normalization statistics");
     std::mt19937_64 rng(seed); std::normal_distribution<float> normal(0.f, 1.f);
-    std::vector<float> joined, previous;
+    std::vector<std::array<float, embedding_width>> embeddings;
+    std::vector<std::vector<float>> noise;
+    std::vector<detail::sampled_sequence_segment> sampled;
+    embeddings.reserve(segments.size()); noise.reserve(segments.size()); sampled.reserve(segments.size());
     for (size_t index=0; index<segments.size(); ++index) {
         const auto &segment=segments[index];
         if (segment.prompt.empty() || segment.frames < 2 || segment.frames > 300)
             return std::unexpected("each sequence segment must contain a prompt and have 2..300 frames");
         if (index && transition_frames >= segment.frames) return std::unexpected("transition must be shorter than every following segment");
-        auto embedding=impl_->text->encode(segment.prompt); if (!embedding) return std::unexpected(embedding.error());
-        // NVIDIA's _multiprompt samples an additional conditioned prefix for
-        // every continuation.  That prefix replaces the old tail, leaving the
-        // caller-requested number of new frames after it is discarded.
+        auto embedding=impl_->text->encode(segment.prompt);
+        if (!embedding) return std::unexpected(embedding.error());
         const auto sampled_frames = static_cast<size_t>(segment.frames) +
             (index == 0 ? 0 : transition_frames);
-        std::vector<float> noise(sampled_frames*273); for (float &value : noise) value=normal(rng);
-        std::vector<float> current;
-        if (index == 0) {
-            auto sampled=detail::sample_motion_from_noise(*impl_->weights,noise,*embedding,segment.frames,steps,text_cfg,constraint_cfg);
-            if (!sampled) return std::unexpected(sampled.error());
-            current=std::move(*sampled);
-        } else {
-            // Derived from NVIDIA's Apache-2.0 `_multiprompt` sampler:
-            // https://github.com/nv-tlabs/kimodo/blob/main/kimodo/model/kimodo_model.py
-            // Preserve the prior tail as observed motion for the next DDIM
-            // run, then use its transition frames to replace the old tail.
-            std::vector<float> observed(noise.size()), observed_mask(noise.size());
-            const auto overlap=static_cast<size_t>(transition_frames);
-            const auto previous_start=previous.size()-overlap*273;
-            // FullBodyConstraintSet's captured mask is deliberately sparse:
-            // global root/posed joints [0,71), smooth root [113,125), and
-            // global rotations [191,203).  In particular, the gaps contain
-            // generated velocities and must not be treated as observed.
-            constexpr std::array<std::pair<size_t, size_t>, 3> constrained = {{
-                {0, 71}, {113, 125}, {191, 203},
-            }};
-            for (size_t frame=0; frame<overlap; ++frame) {
-                std::copy_n(previous.data()+previous_start+frame*273,203,observed.data()+frame*273);
-                for (const auto &[first, last] : constrained)
-                    std::fill(observed_mask.begin()+static_cast<std::ptrdiff_t>(frame*273+first),
-                              observed_mask.begin()+static_cast<std::ptrdiff_t>(frame*273+last), 1.f);
-            }
-            const float origin_x=observed[0]*(*gs)[0]+(*gm)[0];
-            const float origin_z=observed[2]*(*gs)[2]+(*gm)[2];
-            for (size_t frame=0; frame<overlap; ++frame) {
-                auto *row=observed.data()+frame*273;
-                row[0]=((row[0]*(*gs)[0]+(*gm)[0])-origin_x)/(*gs)[0];
-                row[2]=((row[2]*(*gs)[2]+(*gm)[2])-origin_z)/(*gs)[2];
-            }
-            const auto p=(previous.size()/273-overlap)*273;
-            const float heading=std::atan2(previous[p+4]*(*gs)[4]+(*gm)[4],previous[p+3]*(*gs)[3]+(*gm)[3]);
-            auto sampled=detail::sample_motion_from_noise_conditioned(*impl_->weights,noise,*embedding,observed,observed_mask,heading,sampled_frames,steps,text_cfg,constraint_cfg);
-            if (!sampled) return std::unexpected(sampled.error());
-            current=std::move(*sampled);
-            // `_multiprompt` samples in the translated local coordinates, then
-            // restores the prior segment's planar smooth-root origin.
-            for (size_t frame=0; frame<sampled_frames; ++frame) {
-                auto *row=current.data()+frame*273;
-                row[0]=((row[0]*(*gs)[0]+(*gm)[0])+origin_x)/(*gs)[0];
-                row[2]=((row[2]*(*gs)[2]+(*gm)[2])+origin_z)/(*gs)[2];
-            }
-            const auto start=joined.size()-overlap*273;
-            for (size_t frame=0; frame<overlap; ++frame) {
-                const float alpha=overlap==1?.5f:1.f-float(frame)/float(overlap-1);
-                for (size_t d=0; d<273; ++d) joined[start+frame*273+d]=alpha*joined[start+frame*273+d]+(1.f-alpha)*current[frame*273+d];
-            }
-            joined.insert(joined.end(),current.begin()+static_cast<std::ptrdiff_t>(overlap*273),current.end());
-        }
-        if (index == 0) joined=current;
-        previous=std::move(current);
+        embeddings.push_back(*embedding);
+        noise.emplace_back(sampled_frames*273);
+        for (float &value : noise.back()) value=normal(rng);
+        sampled.push_back({embeddings.back(), noise.back(), segment.frames});
     }
-    const auto frames=static_cast<unsigned>(joined.size()/273);
-    auto decoded=detail::decode_smplx22(joined,frames,*gm,*gs,*bm,*bs);
+    auto joined=detail::sample_motion_sequence_from_noise(*impl_->weights,sampled,transition_frames,steps,text_cfg,constraint_cfg);
+    if (!joined) return std::unexpected(joined.error());
+    const auto frames=static_cast<unsigned>(joined->size()/273);
+    auto normalized=*joined;
+    for (size_t row=0; row<frames; ++row) {
+        auto *value=normalized.data()+row*273;
+        for (size_t d=0; d<5; ++d) value[d]=(value[d]-(*gm)[d])/std::sqrt((*gs)[d]*(*gs)[d]+1.e-5F);
+        for (size_t d=0; d<268; ++d) value[5+d]=(value[5+d]-(*bm)[d])/std::sqrt((*bs)[d]*(*bs)[d]+1.e-5F);
+    }
+    auto decoded=detail::decode_smplx22(normalized,frames,*gm,*gs,*bm,*bs);
     if (!decoded) return std::unexpected(decoded.error());
     motion_data result; result.frames=frames; result.joints=22;
     result.local_rotations_xyzw=std::move(decoded->local_xyzw); result.root_positions=std::move(decoded->root_positions);
diff --git a/src/motion_decode.cpp b/src/motion_decode.cpp
index 53c0fdf..31edd04 100644
--- a/src/motion_decode.cpp
+++ b/src/motion_decode.cpp
@@ -9,5 +9,5 @@ M tr(const M&a){M r{};for(int i=0;i<3;++i)for(int j=0;j<3;++j)r.v[i*3+j]=a.v[j*3
 M six(const float*x){float n=std::sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);float a[3]={x[0]/n,x[1]/n,x[2]/n};float z[3]={a[1]*x[5]-a[2]*x[4],a[2]*x[3]-a[0]*x[5],a[0]*x[4]-a[1]*x[3]};n=std::sqrt(z[0]*z[0]+z[1]*z[1]+z[2]*z[2]);for(float&v:z)v/=n;float b[3]={z[1]*a[2]-z[2]*a[1],z[2]*a[0]-z[0]*a[2],z[0]*a[1]-z[1]*a[0]};return M{{a[0],b[0],z[0],a[1],b[1],z[1],a[2],b[2],z[2]}};}
 void quat(const M&m,float*q){float w,x,y,z,t=m.v[0]+m.v[4]+m.v[8];if(t>0){float s=2*std::sqrt(t+1);w=.25f*s;x=(m.v[7]-m.v[5])/s;y=(m.v[2]-m.v[6])/s;z=(m.v[3]-m.v[1])/s;}else if(m.v[0]>m.v[4]&&m.v[0]>m.v[8]){float s=2*std::sqrt(1+m.v[0]-m.v[4]-m.v[8]);w=(m.v[7]-m.v[5])/s;x=.25f*s;y=(m.v[1]+m.v[3])/s;z=(m.v[2]+m.v[6])/s;}else if(m.v[4]>m.v[8]){float s=2*std::sqrt(1+m.v[4]-m.v[0]-m.v[8]);w=(m.v[2]-m.v[6])/s;x=(m.v[1]+m.v[3])/s;y=.25f*s;z=(m.v[5]+m.v[7])/s;}else{float s=2*std::sqrt(1+m.v[8]-m.v[0]-m.v[4]);w=(m.v[3]-m.v[1])/s;x=(m.v[2]+m.v[6])/s;y=(m.v[5]+m.v[7])/s;z=.25f*s;}q[0]=x;q[1]=y;q[2]=z;q[3]=w;}
 }
-std::expected<decoded_motion,std::string> decode_smplx22(std::span<const float>x,size_t T,std::span<const float>gm,std::span<const float>gs,std::span<const float>bm,std::span<const float>bs){if(x.size()!=T*273||gm.size()!=5||gs.size()!=5||bm.size()!=268||bs.size()!=268)return std::unexpected("invalid SMPL-X decode inputs");decoded_motion o;o.local_xyzw.resize(T*22*4);o.root_positions.resize(T*3);for(size_t t=0;t<T;++t){const float*in=x.data()+t*273;std::array<float,273> f{};for(int i=0;i<5;++i)f[i]=in[i]*gs[i]+gm[i];for(int i=0;i<268;++i)f[5+i]=in[5+i]*bs[i]+bm[i];o.root_positions[t*3]=f[0]+f[5];o.root_positions[t*3+1]=f[6];o.root_positions[t*3+2]=f[2]+f[7];M g[22],l[22];for(int j=0;j<22;++j)g[j]=six(f.data()+71+j*6);for(int j=0;j<22;++j)l[j]=parent[j]<0?g[j]:mul(tr(g[parent[j]]),g[j]);for(int j=0;j<22;++j)quat(l[j],o.local_xyzw.data()+(t*22+j)*4);}return o;}
+std::expected<decoded_motion,std::string> decode_smplx22(std::span<const float>x,size_t T,std::span<const float>gm,std::span<const float>gs,std::span<const float>bm,std::span<const float>bs){if(x.size()!=T*273||gm.size()!=5||gs.size()!=5||bm.size()!=268||bs.size()!=268)return std::unexpected("invalid SMPL-X decode inputs");decoded_motion o;o.local_xyzw.resize(T*22*4);o.root_positions.resize(T*3);auto scale=[](float s){return std::sqrt(s*s+1.e-5f);};for(size_t t=0;t<T;++t){const float*in=x.data()+t*273;std::array<float,273> f{};for(int i=0;i<5;++i)f[i]=in[i]*scale(gs[i])+gm[i];for(int i=0;i<268;++i)f[5+i]=in[5+i]*scale(bs[i])+bm[i];o.root_positions[t*3]=f[0]+f[5];o.root_positions[t*3+1]=f[6];o.root_positions[t*3+2]=f[2]+f[7];M g[22],l[22];for(int j=0;j<22;++j)g[j]=six(f.data()+71+j*6);for(int j=0;j<22;++j)l[j]=parent[j]<0?g[j]:mul(tr(g[parent[j]]),g[j]);for(int j=0;j<22;++j)quat(l[j],o.local_xyzw.data()+(t*22+j)*4);}return o;}
 }
diff --git a/src/motion_rep.cpp b/src/motion_rep.cpp
index d27b770..6e87a7c 100644
--- a/src/motion_rep.cpp
+++ b/src/motion_rep.cpp
@@ -13,6 +13,8 @@ std::expected<std::vector<float>, std::string> global_root_to_local_root(
         !std::isfinite(fps) || fps <= 0.f) return std::unexpected("invalid global-root conversion input");
     for (float x : global_std) if (!std::isfinite(x) || x == 0.f) return std::unexpected("invalid global root standard deviation");
     for (float x : local_std) if (!std::isfinite(x) || x == 0.f) return std::unexpected("invalid local root standard deviation");
+    // Match kimodo.motion_rep.stats.Stats: sqrt(std**2 + eps), eps=1e-5.
+    auto scale=[](float stddev) { return std::sqrt(stddev*stddev + 1.e-5f); };
     std::vector<float> result(batch*frames*4);
     for (std::size_t b=0;b<batch;++b) {
         std::size_t length=0; for(std::size_t t=0;t<frames;++t) length += mask[b*frames+t] > .5f;
@@ -20,15 +22,15 @@ std::expected<std::vector<float>, std::string> global_root_to_local_root(
         std::vector<float> angle(frames), x(frames), y(frames), z(frames);
         for (std::size_t t=0;t<frames;++t) {
             const auto p=(b*frames+t)*5;
-            x[t]=root[p]*global_std[0]+global_mean[0]; y[t]=root[p+1]*global_std[1]+global_mean[1]; z[t]=root[p+2]*global_std[2]+global_mean[2];
-            angle[t]=std::atan2(root[p+4]*global_std[4]+global_mean[4], root[p+3]*global_std[3]+global_mean[3]);
+            x[t]=root[p]*scale(global_std[0])+global_mean[0]; y[t]=root[p+1]*scale(global_std[1])+global_mean[1]; z[t]=root[p+2]*scale(global_std[2])+global_mean[2];
+            angle[t]=std::atan2(root[p+4]*scale(global_std[4])+global_mean[4], root[p+3]*scale(global_std[3])+global_mean[3]);
         }
         for(std::size_t t=0;t<frames;++t) {
             const std::size_t next=t+1<length?t+1:length-1, previous=t+1<length?t:length-2;
             const float cos_diff=std::cos(angle[next])*std::cos(angle[previous])+std::sin(angle[next])*std::sin(angle[previous]);
             const float sin_diff=std::sin(angle[next])*std::cos(angle[previous])-std::cos(angle[next])*std::sin(angle[previous]);
             const float raw[]{fps*std::atan2(sin_diff,cos_diff), fps*(x[next]-x[previous]), fps*(z[next]-z[previous]), y[t]};
-            for(std::size_t d=0;d<4;++d) result[(b*frames+t)*4+d]=(raw[d]-local_mean[d])/local_std[d];
+            for(std::size_t d=0;d<4;++d) result[(b*frames+t)*4+d]=(raw[d]-local_mean[d])/scale(local_std[d]);
         }
     }
     return result;
diff --git a/src/sequence.cpp b/src/sequence.cpp
new file mode 100644
index 0000000..031c253
--- /dev/null
+++ b/src/sequence.cpp
@@ -0,0 +1,132 @@
+#include "denoiser.hpp"
+#include "ggml_weights.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+
+namespace kimodo::detail {
+namespace {
+constexpr int parent[22]={-1,0,0,0,1,2,3,4,5,6,7,8,9,9,9,12,13,14,16,17,18,19};
+constexpr float offset[22][3]={{0,0,0},{.052299179F,-.093935639F,-.027606763F},{-.057192899F,-.106548190F,-.022217851F},{-.001495834F,.112929940F,-.024981268F},{.058866613F,-.416441321F,-.006556974F},{-.048074268F,-.397559673F,-.014061437F},{.006900469F,.145636231F,-.006858510F},{-.041737989F,-.437583506F,-.029511765F},{.014489345F,-.446852267F,-.018029511F},{-.010334037F,.056081813F,.021115851F},{.049293540F,-.065279245F,.126259089F},{-.040575184F,-.065286517F,.127075911F},{-.011025756F,.171365142F,-.028827066F},{.047724526F,.087643057F,-.008375450F},{-.046636276F,.086612143F,-.014864366F},{.024654359F,.175390735F,.024463326F},{.126284808F,.057680372F,-.013885141F},{-.109341696F,.053674292F,-.009117880F},{.272907287F,-.069853373F,-.039094493F},{-.292028785F,-.035440356F,-.024564851F},{.276173830F,.021254137F,-.002478220F},{-.271878421F,-.004834589F,-.016445294F}};
+struct mat { double v[9]; };
+mat mul(const mat&a,const mat&b){mat r{};for(int i=0;i<3;++i)for(int j=0;j<3;++j)for(int k=0;k<3;++k)r.v[i*3+j]+=a.v[i*3+k]*b.v[k*3+j];return r;}
+mat trans(const mat&a){mat r{};for(int i=0;i<3;++i)for(int j=0;j<3;++j)r.v[i*3+j]=a.v[j*3+i];return r;}
+mat cont6(const float *x) { double a[3]={x[0],x[1],x[2]}, n=std::sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]); for(double &q:a)q/=n; double z[3]={a[1]*x[5]-a[2]*x[4],a[2]*x[3]-a[0]*x[5],a[0]*x[4]-a[1]*x[3]}; n=std::sqrt(z[0]*z[0]+z[1]*z[1]+z[2]*z[2]);for(double&q:z)q/=n; double b[3]={z[1]*a[2]-z[2]*a[1],z[2]*a[0]-z[0]*a[2],z[0]*a[1]-z[1]*a[0]};return {{a[0],b[0],z[0],a[1],b[1],z[1],a[2],b[2],z[2]}}; }
+void rotate(const mat&m,const float *x,float *o){for(int i=0;i<3;++i)o[i]=static_cast<float>(m.v[i*3]*x[0]+m.v[i*3+1]*x[1]+m.v[i*3+2]*x[2]);}
+}
+std::expected<sequence_transition, std::string> prepare_sequence_transition(
+    const ggml_motion_weights &, std::span<const float> previous,
+    std::size_t continuation_frames, unsigned transition_frames) {
+    constexpr size_t features = 273;
+    const size_t overlap=transition_frames;
+    if (!overlap || overlap>=continuation_frames || previous.size()<=(overlap*features))
+        return std::unexpected("invalid sequence transition");
+    sequence_transition result;
+    result.observed.resize((continuation_frames+overlap)*features);
+    result.observed_mask.resize(result.observed.size());
+    const size_t previous_start=previous.size()-overlap*features;
+    constexpr std::array<std::pair<size_t, size_t>, 3> constrained = {{{0, 71}, {113, 125}, {191, 203}}};
+    for (size_t frame=0; frame<overlap; ++frame) {
+        const size_t base=frame*features;
+        std::array<float, features> value{};
+        const float *raw=previous.data()+previous_start+base;
+        std::copy_n(raw,features,value.data());
+        mat decoded[22],local[22],global[22]; for(int j=0;j<22;++j)decoded[j]=cont6(value.data()+71+j*6);
+        for(int j=0;j<22;++j) local[j]=parent[j]<0?decoded[j]:mul(trans(decoded[parent[j]]),decoded[j]);
+        float root[3]={value[0]+value[5],value[6],value[2]+value[7]}, posed[22][3]{};
+        for(int j=0;j<22;++j){if(parent[j]<0){global[j]=local[j];posed[j][0]=root[0];posed[j][1]=root[1];posed[j][2]=root[2];}else{global[j]=mul(global[parent[j]],local[j]);float d[3];rotate(global[parent[j]],offset[j],d);for(int k=0;k<3;++k)posed[j][k]=posed[parent[j]][k]+d[k];}}
+        // FullBodyConstraintSet: smooth root, root Y, heading, all joint
+        // positions; EndEffectorConstraintSet adds its four rotation blocks.
+        value[0]=value[0]; value[1]=root[1]; value[2]=value[2];
+        // `compute_heading_angle`: right hip minus left hip.
+        const float dx=posed[2][0]-posed[1][0], dz=posed[2][2]-posed[1][2], angle=std::atan2(dz,-dx);
+        value[3]=std::cos(angle); value[4]=std::sin(angle);
+        for(int j=0;j<22;++j){value[5+j*3]=posed[j][0]-value[0];value[6+j*3]=posed[j][1];value[7+j*3]=posed[j][2]-value[2];}
+        for(int j: {7,8,20,21}) for(int d=0;d<6;++d)
+            value[71+j*6+d]=static_cast<float>(global[j].v[(d%3)*3+d/3]);
+        std::copy_n(value.data(),203,result.observed.data()+base);
+        for (const auto &[first,last] : constrained)
+            std::fill(result.observed_mask.begin()+static_cast<std::ptrdiff_t>(base+first),
+                      result.observed_mask.begin()+static_cast<std::ptrdiff_t>(base+last),1.F);
+    }
+    result.origin_x=result.observed[0];
+    result.origin_z=result.observed[2];
+    for (size_t frame=0; frame<overlap; ++frame) {
+        auto *row=result.observed.data()+frame*features;
+        row[0]-=result.origin_x;
+        row[2]-=result.origin_z;
+    }
+    // First heading comes from the first retained full-body constraint.
+    const size_t first=(previous.size()-overlap*features);
+    const float *raw=previous.data()+first;
+    std::array<float, features> value{}; std::copy_n(raw,features,value.data());
+    mat decoded[22],local[22],global[22]; for(int j=0;j<22;++j)decoded[j]=cont6(value.data()+71+j*6);
+    for(int j=0;j<22;++j)local[j]=parent[j]<0?decoded[j]:mul(trans(decoded[parent[j]]),decoded[j]);
+    float root[3]={value[0]+value[5],value[6],value[2]+value[7]}, posed[22][3]{};
+    for(int j=0;j<22;++j){if(parent[j]<0){global[j]=local[j];for(int k=0;k<3;++k)posed[j][k]=root[k];}else{global[j]=mul(global[parent[j]],local[j]);float d[3];rotate(global[parent[j]],offset[j],d);for(int k=0;k<3;++k)posed[j][k]=posed[parent[j]][k]+d[k];}}
+    result.first_heading=std::atan2(posed[2][2]-posed[1][2],-(posed[2][0]-posed[1][0]));
+    return result;
+}
+
+std::expected<std::vector<float>, std::string> sample_motion_sequence_from_noise(
+    const ggml_motion_weights &weights, std::span<const sampled_sequence_segment> segments,
+    unsigned transition_frames, unsigned steps, float text_weight, float constraint_weight) {
+    constexpr size_t features = 273;
+    if (segments.empty() || !transition_frames)
+        return std::unexpected("sequence requires segments and a transition");
+    auto gm=weights.f32_values("stats.global_root.mean"), gs=weights.f32_values("stats.global_root.std");
+    auto bm=weights.f32_values("stats.body.mean"), bs=weights.f32_values("stats.body.std");
+    if (!gm || !gs || !bm || !bs) return std::unexpected("motion GGUF lacks motion statistics");
+    // Upstream Stats normalizes with sqrt(std^2 + 1e-5), rather than raw std.
+    auto scale = [](float stddev) { return std::sqrt(stddev * stddev + 1.e-5F); };
+    auto unnormalize = [&](std::vector<float> &motion) { for(size_t row=0;row<motion.size()/features;++row) { auto *v=motion.data()+row*features; for(size_t d=0;d<5;++d)v[d]=v[d]*scale((*gs)[d])+(*gm)[d]; for(size_t d=0;d<268;++d)v[5+d]=v[5+d]*scale((*bs)[d])+(*bm)[d]; } };
+    auto normalize = [&](std::vector<float> &motion) { for(size_t row=0;row<motion.size()/features;++row) { auto *v=motion.data()+row*features; for(size_t d=0;d<5;++d)v[d]=(v[d]-(*gm)[d])/scale((*gs)[d]); for(size_t d=0;d<268;++d)v[5+d]=(v[5+d]-(*bm)[d])/scale((*bs)[d]); } };
+
+    std::vector<float> joined, previous;
+    for (size_t index=0; index<segments.size(); ++index) {
+        const auto &segment=segments[index];
+        const size_t sampled_frames=segment.frames+(index ? transition_frames : 0);
+        if (segment.frames < 2 || segment.embedding.size()!=4096 ||
+            segment.initial_noise.size()!=sampled_frames*features)
+            return std::unexpected("invalid sampled sequence segment");
+        std::vector<float> current;
+        if (!index) {
+            auto sampled=sample_motion_from_noise(weights,segment.initial_noise,segment.embedding,
+                                                  sampled_frames,steps,text_weight,constraint_weight);
+            if (!sampled) return std::unexpected(sampled.error());
+            current=std::move(*sampled);
+            unnormalize(current);
+        } else {
+            const size_t overlap=transition_frames;
+            if (overlap >= segment.frames || previous.size()<overlap*features)
+                return std::unexpected("transition must be shorter than every following segment");
+            auto transition=prepare_sequence_transition(weights,previous,segment.frames,transition_frames);
+            if (!transition) return std::unexpected(transition.error());
+            const float origin_x=transition->origin_x;
+            const float origin_z=transition->origin_z;
+            normalize(transition->observed);
+            auto sampled=sample_motion_from_noise_conditioned(weights,segment.initial_noise,segment.embedding,
+                transition->observed,transition->observed_mask,transition->first_heading,sampled_frames,steps,text_weight,constraint_weight);
+            if (!sampled) return std::unexpected(sampled.error());
+            current=std::move(*sampled);
+            unnormalize(current);
+            for (size_t frame=0; frame<sampled_frames; ++frame) {
+                auto *row=current.data()+frame*features;
+                row[0]+=origin_x;
+                row[2]+=origin_z;
+            }
+            const size_t start=joined.size()-overlap*features;
+            for (size_t frame=0; frame<overlap; ++frame) {
+                const float alpha=overlap==1?.5F:1.F-float(frame)/float(overlap-1);
+                for (size_t d=0; d<features; ++d)
+                    joined[start+frame*features+d]=alpha*joined[start+frame*features+d]+(1.F-alpha)*current[frame*features+d];
+            }
+            joined.insert(joined.end(),current.begin()+static_cast<std::ptrdiff_t>(overlap*features),current.end());
+        }
+        if (!index) joined=current;
+        previous=std::move(current);
+    }
+    return joined;
+}
+} // namespace kimodo::detail
diff --git a/tests/multiprompt_fixture_parity.cpp b/tests/multiprompt_fixture_parity.cpp
index 95a3836..dae9b53 100644
--- a/tests/multiprompt_fixture_parity.cpp
+++ b/tests/multiprompt_fixture_parity.cpp
@@ -5,6 +5,7 @@
 #include "ggml_weights.hpp"
 
 #include <algorithm>
+#include <array>
 #include <cmath>
 #include <cstdio>
 #include <cstring>
@@ -42,15 +43,23 @@ error compare(const std::vector<float> &actual, const std::vector<float> &expect
     return result;
 }
 
-void unnormalize(std::vector<float> &motion, const std::vector<float> &global_mean,
-                 const std::vector<float> &global_std, const std::vector<float> &body_mean,
-                 const std::vector<float> &body_std) {
-    for (std::size_t row = 0; row < motion.size() / features; ++row) {
-        auto *value = motion.data() + row * features;
-        for (std::size_t d = 0; d < 5; ++d) value[d] = value[d] * global_std[d] + global_mean[d];
-        for (std::size_t d = 0; d < 268; ++d) value[5 + d] = value[5 + d] * body_std[d] + body_mean[d];
+error compare_masked(const std::vector<float> &actual, const std::vector<float> &expected,
+                     const std::vector<float> &mask) {
+    if (actual.size() != expected.size() || actual.size() != mask.size())
+        throw std::runtime_error("masked fixture shape mismatch");
+    double squared_error = 0, squared_reference = 0;
+    error result;
+    for (std::size_t i = 0; i < actual.size(); ++i) {
+        if (mask[i] == 0.F) continue;
+        const float difference = actual[i] - expected[i];
+        result.max_abs = std::max(result.max_abs, std::abs(difference));
+        squared_error += static_cast<double>(difference) * difference;
+        squared_reference += static_cast<double>(expected[i]) * expected[i];
     }
+    result.relative_l2 = std::sqrt(squared_error / squared_reference);
+    return result;
 }
+
 }
 
 int main(int argc, char **argv) try {
@@ -76,37 +85,50 @@ int main(int argc, char **argv) try {
     if (!second) throw std::runtime_error(second.error());
     const auto second_error = compare(*second, read(directory + "segment_01_sampling_output_001.f32"));
 
-    auto global_mean = (*weights)->f32_values("stats.global_root.mean");
-    auto global_std = (*weights)->f32_values("stats.global_root.std");
-    auto body_mean = (*weights)->f32_values("stats.body.mean");
-    auto body_std = (*weights)->f32_values("stats.body.std");
-    if (!global_mean || !global_std || !body_mean || !body_std) throw std::runtime_error("missing motion statistics");
-    auto stitched_first = *first;
-    auto stitched_second = *second;
-    unnormalize(stitched_first, *global_mean, *global_std, *body_mean, *body_std);
-    unnormalize(stitched_second, *global_mean, *global_std, *body_mean, *body_std);
-    constexpr std::size_t overlap = 5;
-    // The captured observed tensor is already translated to local origin;
-    // recover the world origin from the first segment's retained tail.
-    const float origin_x = stitched_first[(30 - overlap) * features];
-    const float origin_z = stitched_first[(30 - overlap) * features + 2];
-    for (std::size_t frame = 0; frame < 35; ++frame) {
-        stitched_second[frame * features] += origin_x;
-        stitched_second[frame * features + 2] += origin_z;
+    auto prior=read(directory + "segment_00_sampling_output_001.f32");
+    auto gm=(*weights)->f32_values("stats.global_root.mean"), gs=(*weights)->f32_values("stats.global_root.std");
+    auto bm=(*weights)->f32_values("stats.body.mean"), bs=(*weights)->f32_values("stats.body.std");
+    if (!gm || !gs || !bm || !bs) throw std::runtime_error("missing motion statistics");
+    auto scale=[](float stddev) { return std::sqrt(stddev*stddev+1.e-5F); };
+    for (std::size_t row=0;row<30;++row) { auto *v=prior.data()+row*features; for(std::size_t d=0;d<5;++d)v[d]=v[d]*scale((*gs)[d])+(*gm)[d]; for(std::size_t d=0;d<268;++d)v[5+d]=v[5+d]*scale((*bs)[d])+(*bm)[d]; }
+    const auto transition=kimodo::detail::prepare_sequence_transition(**weights,prior,30,5);
+    if (!transition) throw std::runtime_error(transition.error());
+    auto actual_observed=transition->observed;
+    for (std::size_t row=0;row<35;++row) { auto *v=actual_observed.data()+row*features; for(std::size_t d=0;d<5;++d)v[d]=(v[d]-(*gm)[d])/scale((*gs)[d]); for(std::size_t d=0;d<268;++d)v[5+d]=(v[5+d]-(*bm)[d])/scale((*bs)[d]); }
+    const auto expected_observed=read(directory + "segment_01_observed_motion.f32");
+    const auto expected_mask=read(directory + "segment_01_motion_mask.f32");
+    const auto observed_error=compare_masked(actual_observed,expected_observed,expected_mask);
+    std::size_t worst=0; float worst_value=0.F;
+    for (std::size_t i=0;i<expected_mask.size();++i) if (expected_mask[i] != 0.F) {
+        const float difference=std::abs(actual_observed[i]-expected_observed[i]);
+        if (difference>worst_value) { worst_value=difference; worst=i; }
     }
-    for (std::size_t frame = 0; frame < overlap; ++frame) {
-        const float alpha = 1.F - static_cast<float>(frame) / static_cast<float>(overlap - 1);
-        for (std::size_t d = 0; d < features; ++d)
-            stitched_first[(30 - overlap + frame) * features + d] =
-                alpha * stitched_first[(30 - overlap + frame) * features + d] +
-                (1.F - alpha) * stitched_second[frame * features + d];
-    }
-    stitched_first.insert(stitched_first.end(), stitched_second.begin() + static_cast<std::ptrdiff_t>(overlap * features), stitched_second.end());
-    const auto stitched_error = compare(stitched_first, read(directory + "stitched_motion_rep.f32"));
-    std::printf("segment0 max_abs=%g rel_l2=%g\nsegment1 max_abs=%g rel_l2=%g\nstitched max_abs=%g rel_l2=%g\n",
+    const auto mask_error=compare(transition->observed_mask,expected_mask);
+    const float heading_error=std::abs(transition->first_heading-heading.at(0));
+    const auto constructed_second=kimodo::detail::sample_motion_from_noise_conditioned(
+        **weights, read(directory + "segment_01_sampling_input_000.f32"),
+        read(directory + "segment_01_text_features.f32"), actual_observed, transition->observed_mask,
+        transition->first_heading, 35, 2, 2.F, 2.F);
+    if (!constructed_second) throw std::runtime_error(constructed_second.error());
+    const auto constructed_error=compare(*constructed_second,read(directory + "segment_01_sampling_output_001.f32"));
+    const auto first_noise=read(directory + "segment_00_sampling_input_000.f32");
+    const auto first_text=read(directory + "segment_00_text_features.f32");
+    const auto second_noise=read(directory + "segment_01_sampling_input_000.f32");
+    const auto second_text=read(directory + "segment_01_text_features.f32");
+    const std::array<kimodo::detail::sampled_sequence_segment, 2> segments{{
+        {first_text, first_noise, 30}, {second_text, second_noise, 30},
+    }};
+    const auto joined=kimodo::detail::sample_motion_sequence_from_noise(
+        **weights, segments, 5, 2, 2.F, 2.F);
+    if (!joined) throw std::runtime_error(joined.error());
+    const auto joined_error=compare(*joined,read(directory + "stitched_motion_rep.f32"));
+    std::printf("segment0 max_abs=%g rel_l2=%g\nsegment1 max_abs=%g rel_l2=%g\ntransition observed max_abs=%g worst=%zu mask max_abs=%g heading_abs=%g constructed max_abs=%g stitched max_abs=%g\n",
                 first_error.max_abs, first_error.relative_l2, second_error.max_abs, second_error.relative_l2,
-                stitched_error.max_abs, stitched_error.relative_l2);
-    return (first_error.max_abs <= 3.e-3F && second_error.max_abs <= 3.e-3F && stitched_error.max_abs <= 3.e-3F) ? 0 : 1;
+                observed_error.max_abs, worst, mask_error.max_abs, heading_error, constructed_error.max_abs,
+                joined_error.max_abs);
+    return (first_error.max_abs <= 3.e-3F && second_error.max_abs <= 3.e-3F &&
+            observed_error.max_abs <= 3.e-5F && mask_error.max_abs == 0.F && heading_error <= 2.e-3F &&
+            joined_error.max_abs <= 3.e-3F) ? 0 : 1;
 } catch (const std::exception &error) {
     std::fprintf(stderr, "multi-prompt fixture parity error: %s\n", error.what());
     return 1;