@@ -80,10 +80,15 @@ is
 
   `series_out[k-1] = f^(k)(x₀) / k!`.
 
-To recover the unscaled `k`-th derivative, multiply by `k!`. (This differs
-from `jax.experimental.jet`, which returns derivative coefficients without
-the factorial scaling. The relationship is documented and tested in
-`tests/test_against_jax_jet.py`.)
+To recover the unscaled `k`-th derivative, multiply by `k!`. This is the
+convention used in *Evaluating Derivatives* (Griewank & Walther, §13).
+
+`jax.experimental.jet` exposes both conventions through its
+`factorial_scaled` keyword argument (defaults to `True`, returning
+derivative coefficients `f^(k)(x₀)`). Calling
+`jax.experimental.jet.jet(..., factorial_scaled=False)` returns the same
+Taylor coefficients as `jet_array`. The equivalence is tested
+coefficient-by-coefficient in `tests/test_against_jax_jet.py`.
 
 ## Examples
 
@@ -227,13 +232,27 @@ combine with `jax.vmap`.
 
 ### Caveats
 
-- Only primitives in `_rules_with_effective_order` honor the hint. Other
-  primitives compute the full series. So for a function like `sin(x)`
-  alone, `effective_order` has no effect — `sin` always computes all `K`
-  coefficients.
-- Output entries with index `≥ effective_order` are not guaranteed to
-  contain anything specific (zero, garbage, or a partial answer). Slice
-  to `series_out[:effective_order]` before using them.
+- The primitives a JAX program uses fall into three categories:
+  - **Linear / zero-derivative** (`add`, `sub`, multiplication by a
+    constant, `neg`, broadcasting, slicing, comparisons, …): the operation
+    applies elementwise along the order axis — `(a + b)[k] = a[k] + b[k]`
+    — so it is already a single `O(K)` XLA op over the series array.
+    There is no convolution loop to short-circuit, and `effective_order`
+    has no effect. (Output entries above `effective_order` still hold
+    correct values for these primitives, because the linear op computes
+    every entry uniformly.)
+  - **Nonlinear without an `effective_order`-aware rule** (`sin`, `cos`,
+    `sinh`, `cosh`, `square`, `sqrt`, `erf`, …): the rule computes the
+    full `O(K²)` convolution regardless. `effective_order` is silently
+    ignored on this path.
+  - **Nonlinear with an `effective_order`-aware rule** (`exp`, `expm1`,
+    `log`, `log1p`, `pow`, `logistic`, `tanh`, `erf_inv`, `div`): the
+    convolution scan short-circuits iterations beyond `effective_order`.
+    Output entries above `effective_order` are unspecified — slice to
+    `series_out[:effective_order]` before using them.
+- Because output entries above `effective_order` are unspecified for the
+  third category, treat the whole `series_out[effective_order:]` slice as
+  garbage even if upstream linear ops would otherwise leave it intact.
 - The win is largest when the dominant cost is in the convolution-style
   loops that the hint short-circuits. For shallow programs the overhead of
   the conditional may exceed the savings.
 
@@ -1,19 +1,12 @@
 """Correctness tests against ``jax.experimental.jet``.
 
-``jax.experimental.jet`` returns derivative coefficients
-``[f(x), f'(x), f''(x), f'''(x), ...]``; ``jet_array`` returns Taylor
-coefficients ``[f(x), f'(x), f''(x)/2!, f'''(x)/3!, ...]``. So the
-relationship is::
-
-    jet_array_series[k] == standard_jet_series[k] / (k+1)!  for k = 0, 1, ...
-
-(Note the index shift: ``standard_jet`` returns the primal as
-``series[0]`` only when the user requests it; in our framing
-``primal_out`` is returned separately and ``series_out[k]`` is the
-``(k+1)``-th coefficient.)
-
-Each test scales the standard-jet output by 1/k! and compares against
-``jet_array`` to within tight tolerance.
+Both libraries return Taylor coefficients ``[f'(x)/1!, f''(x)/2!, ...]``
+when invoked with matching conventions. ``jet_array.jet`` always uses
+this convention; ``jax.experimental.jet.jet`` returns it when called
+with ``factorial_scaled=False`` (with default ``True`` it returns the
+unscaled derivative coefficients ``[f'(x), f''(x), ...]``). These tests
+exercise the matching path and assert coefficient-by-coefficient
+equivalence.
 """
 
 import math
@@ -35,26 +28,22 @@
 ATOL = 1e-5
 
 
-def _scale_to_taylor(series_in_deriv):
-    """Convert standard-jet derivative coefficients to Taylor coefficients."""
-    return [coeff / math.factorial(k + 1) for k, coeff in enumerate(series_in_deriv)]
-
-
 def _compare(name, x0, fn, order=5):
     """Compare jet_array.jet against jax.experimental.jet at one point."""
     x0_arr = jnp.asarray(x0, dtype=jnp.float64)
     series_in_arr = jnp.zeros(order, dtype=jnp.float64).at[0].set(1.0)
     p_arr, s_arr = jet(fn, (x0_arr,), (series_in_arr,))
 
     series_in_std = [1.0] + [0.0] * (order - 1)
-    p_std, s_std = standard_jet.jet(fn, (x0_arr,), (series_in_std,))
-    s_std_taylor = _scale_to_taylor(s_std)
+    p_std, s_std = standard_jet.jet(
+        fn, (x0_arr,), (series_in_std,), factorial_scaled=False
+    )
 
     np.testing.assert_allclose(
         p_arr, p_std, rtol=RTOL, atol=ATOL,
         err_msg=f"{name}@{x0}: primal mismatch",
     )
-    for k, (a, b) in enumerate(zip(s_arr, s_std_taylor)):
+    for k, (a, b) in enumerate(zip(s_arr, s_std)):
         np.testing.assert_allclose(
             a, b, rtol=RTOL, atol=ATOL,
             err_msg=f"{name}@{x0}: coefficient {k+1} mismatch",
@@ -92,11 +81,12 @@ def test_multivariate_matches_standard_jet(name, fn, order):
 
     sx_std = [1.0] + [0.0] * (order - 1)
     sy_std = [1.0] + [0.0] * (order - 1)
-    p_std, s_std = standard_jet.jet(fn, (x0, y0), (sx_std, sy_std))
-    s_std_taylor = _scale_to_taylor(s_std)
+    p_std, s_std = standard_jet.jet(
+        fn, (x0, y0), (sx_std, sy_std), factorial_scaled=False
+    )
 
     np.testing.assert_allclose(p_arr, p_std, rtol=RTOL, atol=ATOL)
-    for k, (a, b) in enumerate(zip(s_arr, s_std_taylor)):
+    for k, (a, b) in enumerate(zip(s_arr, s_std)):
         np.testing.assert_allclose(
             a, b, rtol=RTOL, atol=ATOL,
             err_msg=f"{name}: coefficient {k+1} mismatch",