Verify and test higher-order grad

Nested calls such as grad(grad(f, x), x) already yield higher-order
derivatives: the plan rewrite walks expressions bottom-up (transform_up), so
the inner grad is differentiated to a plain expression first and the outer
grad differentiates that result. No code change was needed; this adds tests
and documents the behavior.

- Rust: a unit test that differentiation composes (d2/dx2 sin = -sin).
- Python: second derivatives of sin (-sin) and x^3 (6x) and the third
  derivative of sin (-cos), against numpy.
- Doc: note higher-order support in the module overview.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017mDoFJgsm9kS7SicGoCVF6

Author: claude

src/autograd.rs

@@ -36,6 +36,10 @@
 //! full gradient or Jacobian is expressed as several scalar columns (e.g.
 //! `grad(f, x) AS dfdx, grad(f, y) AS dfdy`) rather than a nested array, which
 //! would break the one-value-per-coordinate model.
+//!
+//! Calls nest, giving higher-order derivatives for free: the rewrite walks
+//! bottom-up, so the inner call in `grad(grad(f, x), x)` is differentiated
+//! first and the outer call differentiates that result.
 
 #![allow(dead_code)]
 
@@ -596,6 +600,14 @@ mod tests {
         assert!(differentiate(&e, "x").is_err());
     }
 
+    #[test]
+    fn higher_order_derivative() {
+        // Differentiation composes: d2/dx2 sin(x) = -sin(x).
+        let d1 = differentiate(&expr_fn::sin(col("x")), "x").unwrap();
+        let d2 = differentiate(&d1, "x").unwrap();
+        assert_eq!(d2, neg(expr_fn::sin(col("x"))));
+    }
+
     #[test]
     fn jvp_seeds_a_tangent_on_one_input() {
         // jvp(x*y, {x: dx}) = product rule with tangent(x)=dx, tangent(y)=0

tests/test_autograd.py

@@ -80,6 +80,31 @@ def test_grad_quotient_and_power(ctx):
     np.testing.assert_allclose(res["dcube"], 3.0 * val**2)
 
 
+def test_higher_order_grad(ctx):
+    # Nested grad() differentiates repeatedly: the inner call is rewritten
+    # first, then the outer differentiates its result.
+    val = np.linspace(0.1, 3.0, 16)
+    res = _ordered(
+        ctx.sql(
+            "SELECT i, "
+            "grad(grad(sin(val), val), val) AS d2_sin, "
+            "grad(grad(power(val, 3), val), val) AS d2_cube FROM t"
+        )
+    )
+    np.testing.assert_allclose(res["d2_sin"], -np.sin(val))  # -sin
+    np.testing.assert_allclose(res["d2_cube"], 6.0 * val)  # d2/dx2 x^3 = 6x
+
+
+def test_third_order_grad(ctx):
+    val = np.linspace(0.1, 3.0, 16)
+    res = _ordered(
+        ctx.sql(
+            "SELECT i, grad(grad(grad(sin(val), val), val), val) AS d3 FROM t"
+        )
+    )
+    np.testing.assert_allclose(res["d3"], -np.cos(val))  # d3/dx3 sin = -cos
+
+
 def test_non_grad_query_is_unaffected(ctx):
     # Queries without grad() bypass the rewrite and behave normally.
     res = _ordered(ctx.sql("SELECT i, val FROM t"))