diff --git a/src/api/python/expression.rs b/src/api/python/expression.rs
index 936b8a27..a7fc3769 100644
--- a/src/api/python/expression.rs
+++ b/src/api/python/expression.rs
@@ -8816,9 +8816,9 @@ impl PythonExpression {
     /// functions: Sequence[FunctionDefinition] = []
     ///     Function definitions available to the evaluator.
     /// iterations: int, optional
-    ///     The number of optimization iterations to perform.
+    ///     The number of Horner schemes to try.
     /// cpe_iterations: Optional[int], optional
-    ///     The number of common subexpression elimination iterations to perform.
+    ///     The maximum number of common-pair elimination iterations. If `None`, continue until no more pairs are eliminated.
     /// n_cores: int, optional
     ///     The number of cores to use for the optimization.
     /// verbose: bool, optional
@@ -8839,6 +8839,10 @@ impl PythonExpression {
     ///     The maximum number of entries in the common pair cache.
     /// max_common_pair_distance: int, optional
     ///     The maximum distance between common pairs. Used when clearing cache entries.
+    /// horner_use_cpe: bool, optional
+    ///     If `False` (the default), rank Horner candidates using common-subexpression elimination only.
+    ///     If `True`, also run common-pair elimination for each candidate; this can improve the selected scheme but is more expensive.
+    ///     The selected scheme receives final common-pair elimination in either mode, controlled by `cpe_iterations`.
     #[pyo3(signature =
         (
         params,
@@ -8854,7 +8858,8 @@ impl PythonExpression {
         jit_options = HashMap::default(),
         max_horner_scheme_variables = 500,
         max_common_pair_cache_entries = 1_000_000,
-        max_common_pair_distance = 100),
+        max_common_pair_distance = 100,
+        horner_use_cpe = false),
         )]
     pub fn evaluator(
         &self,
@@ -8872,6 +8877,7 @@ impl PythonExpression {
         max_horner_scheme_variables: usize,
         max_common_pair_cache_entries: usize,
         max_common_pair_distance: usize,
+        horner_use_cpe: bool,
         py: Python,
     ) -> PyResult<PythonExpressionEvaluator> {
         #[cfg(not(feature = "native_code_generation"))]
@@ -8905,6 +8911,7 @@ impl PythonExpression {
 
         let settings = OptimizationSettings {
             horner_iterations: iterations,
+            horner_use_cpe,
             cpe_iterations,
             n_cores,
             verbose: verbose.into(),
@@ -8979,9 +8986,9 @@ impl PythonExpression {
     /// functions: Sequence[FunctionDefinition]
     ///     Function definitions available to the evaluator.
     /// iterations: int
-    ///     The number of optimization passes to run.
+    ///     The number of Horner schemes to try.
     /// cpe_iterations: int | None
-    ///     The number of common subexpression elimination iterations to perform.
+    ///     The maximum number of common-pair elimination iterations. If `None`, continue until no more pairs are eliminated.
     /// n_cores: int
     ///     The number of CPU cores used for parallel optimization.
     /// verbose: bool
@@ -9002,6 +9009,10 @@ impl PythonExpression {
     ///     The maximum number of common-subexpression pairs to cache.
     /// max_common_pair_distance: int
     ///     The maximum distance between factors when searching for common pairs.
+    /// horner_use_cpe: bool
+    ///     If `False` (the default), rank Horner candidates using common-subexpression elimination only.
+    ///     If `True`, also run common-pair elimination for each candidate; this can improve the selected scheme but is more expensive.
+    ///     The selected scheme receives final common-pair elimination in either mode, controlled by `cpe_iterations`.
     #[classmethod]
     #[pyo3(signature =
         (exprs,
@@ -9018,7 +9029,8 @@ impl PythonExpression {
         jit_options = HashMap::default(),
         max_horner_scheme_variables = 500,
         max_common_pair_cache_entries = 1_000_000,
-        max_common_pair_distance = 100)
+        max_common_pair_distance = 100,
+        horner_use_cpe = false)
     )]
     pub fn evaluator_multiple(
         _cls: &Bound<'_, PyType>,
@@ -9037,6 +9049,7 @@ impl PythonExpression {
         max_horner_scheme_variables: usize,
         max_common_pair_cache_entries: usize,
         max_common_pair_distance: usize,
+        horner_use_cpe: bool,
         py: Python,
     ) -> PyResult<PythonExpressionEvaluator> {
         #[cfg(not(feature = "native_code_generation"))]
@@ -9070,6 +9083,7 @@ impl PythonExpression {
 
         let settings = OptimizationSettings {
             horner_iterations: iterations,
+            horner_use_cpe,
             cpe_iterations,
             n_cores,
             verbose: verbose.into(),
diff --git a/src/atom/core.rs b/src/atom/core.rs
index 6675cb7b..8bac5ccf 100644
--- a/src/atom/core.rs
+++ b/src/atom/core.rs
@@ -331,7 +331,7 @@ pub trait AtomCore: private::Sealed + Sized {
         let vs = variables.map(|v| v.iter().map(|x| x.clone().into()).collect::<Vec<_>>());
 
         self.as_atom_view()
-            .horner_scheme(vs.as_deref(), false, false)
+            .horner_scheme_local(vs.as_deref(), false)
             .wrap(self)
     }
 
diff --git a/src/collect.rs b/src/collect.rs
index 69e048c1..2e85967d 100644
--- a/src/collect.rs
+++ b/src/collect.rs
@@ -274,7 +274,10 @@ impl<'a> AtomView<'a> {
             let f = AlgebraicExtension::complex(Q);
             let f2 = FloatField::from_rep(Complex::new(Rational::zero(), Rational::one()));
             let poly = self.try_to_factorized_rational_polynomial::<_, _, u32>(&f, &f, None);
-            if matches!(&poly, Err(PolynomialConversionError::ZeroDenominator { .. })) {
+            if matches!(
+                &poly,
+                Err(PolynomialConversionError::ZeroDenominator { .. })
+            ) {
                 out.set_from_view(self);
                 return;
             }
@@ -323,7 +326,10 @@ impl<'a> AtomView<'a> {
             }
         } else {
             let poly = self.try_to_factorized_rational_polynomial::<_, _, u32>(&Q, &Z, None);
-            if matches!(&poly, Err(PolynomialConversionError::ZeroDenominator { .. })) {
+            if matches!(
+                &poly,
+                Err(PolynomialConversionError::ZeroDenominator { .. })
+            ) {
                 out.set_from_view(self);
                 return;
             }
@@ -1435,16 +1441,17 @@ impl<'a> AtomView<'a> {
         }
     }
 
-    /// Construct a Horner scheme for the given variables. If no variables are provided,
-    /// a heuristically determined near-optimal ordering is used.
-    pub(crate) fn horner_scheme<'b>(
+    /// Construct a Horner scheme within this algebraic region. Function atoms
+    /// and powers are terminal values; their arguments are handled separately.
+    /// Positive integer powers can still be factored as powers of an indeterminate.
+    /// If no variables are provided, their occurrence counts determine the order.
+    pub(crate) fn horner_scheme_local<T: AtomCore>(
         &self,
-        xs: Option<&[Indeterminate]>,
-        enter_functions: bool,
+        xs: Option<&[T]>,
         collect_by_coefficient: bool,
     ) -> Atom {
         let r = if let Some(xs) = xs {
-            self.horner_scheme_impl(xs, enter_functions)
+            self.horner_scheme_impl(xs)
         } else {
             // sort variables by their occurrence count
             let mut v = HashMap::default();
@@ -1457,29 +1464,175 @@ impl<'a> AtomView<'a> {
                 .map(|(x, _)| Indeterminate::try_from(x.to_owned()).unwrap())
                 .collect::<Vec<_>>();
 
-            self.horner_scheme_impl(&res, enter_functions)
+            self.horner_scheme_impl(&res)
         };
 
         if collect_by_coefficient {
-            r.as_view().collect_by_coefficient()
+            Workspace::get_local().with(|ws| {
+                let grouped = r.as_view().cleanup_horner_region(ws, false);
+                grouped
+                    .as_view()
+                    .cleanup_horner_region(ws, true)
+                    .into_owned()
+            })
         } else {
             r
         }
     }
 
-    pub(crate) fn horner_scheme_impl<'b>(
-        &self,
-        xs: &[Indeterminate],
-        enter_functions: bool,
-    ) -> Atom {
+    /// Clean coefficients bottom-up within an Add/Mul region, retaining unchanged
+    /// children by reference. Function and power children must remain untouched.
+    /// Keep grouping and content extraction as separate passes so parent grouping
+    /// sees the same coefficients as the general coefficient collector.
+    fn cleanup_horner_region(&self, ws: &Workspace, canonicalize_content: bool) -> AtomOrView<'a> {
+        let (is_sum, children) = match self {
+            AtomView::Add(sum) => (true, sum.iter()),
+            AtomView::Mul(product) => (false, product.iter()),
+            _ => return self.into(),
+        };
+        let mut rebuilt: Option<Atom> = None;
+        for (i, child) in children.enumerate() {
+            let cleaned = child.cleanup_horner_region(ws, canonicalize_content);
+            if rebuilt.is_none() && matches!(cleaned, AtomOrView::Atom(_)) {
+                let mut result = Atom::new();
+                if is_sum {
+                    let sum = result.to_add();
+                    for previous in children.take(i) {
+                        sum.extend(previous);
+                    }
+                } else {
+                    let product = result.to_mul();
+                    for previous in children.take(i) {
+                        product.extend(previous);
+                    }
+                }
+                rebuilt = Some(result);
+            }
+            if let Some(result) = &mut rebuilt {
+                match result {
+                    Atom::Add(sum) => sum.extend(cleaned.as_view()),
+                    Atom::Mul(product) => product.extend(cleaned.as_view()),
+                    _ => unreachable!(),
+                }
+            }
+        }
+        let current: AtomOrView<'a> = if let Some(rebuilt) = rebuilt {
+            let mut normalized = Atom::new();
+            rebuilt.as_view().normalize(ws, &mut normalized);
+            normalized.into()
+        } else {
+            self.into()
+        };
+        let mut result = Atom::new();
+        let mut out = Settable::from(&mut result);
+        if canonicalize_content {
+            current
+                .as_view()
+                .canonicalize_horner_content_into(ws, &mut out);
+        } else {
+            current
+                .as_view()
+                .collect_num_exact_no_norm_impl(ws, &mut out);
+        }
+        if out.is_set() {
+            if !canonicalize_content {
+                let mut normalized = Atom::new();
+                result.as_view().normalize(ws, &mut normalized);
+                return normalized.into();
+            }
+            result.into()
+        } else {
+            current
+        }
+    }
+
+    /// Expose proportional sums to CSE, including sums with implicit unit
+    /// coefficients. Restrict this to exact real rationals: extracting content
+    /// from floating point, finite-field or symbolic coefficients needs a
+    /// different cost and normalization policy.
+    fn canonicalize_horner_content(&self, ws: &Workspace) -> Atom {
+        let mut result = Atom::new();
+        let mut out = Settable::from(&mut result);
+        self.canonicalize_horner_content_into(ws, &mut out);
+        if !out.is_set() {
+            out.set_from_view(self);
+        }
+        result
+    }
+
+    /// Leave `out` unset when an already-normalized expression has no content
+    /// to extract. Bottom-up maps can then retain its parents without copying
+    /// and normalizing their entire subtrees again.
+    fn canonicalize_horner_content_into(&self, ws: &Workspace, out: &mut Settable<'_, Atom>) {
+        let normalized = if self.needs_normalization() {
+            self.normalize(ws, out);
+            out.as_view()
+        } else {
+            *self
+        };
+        let AtomView::Add(sum) = normalized else {
+            return;
+        };
+        let mut content: Option<Rational> = None;
+        let mut negative = false;
+        for term in sum {
+            let coefficient = match term {
+                AtomView::Num(_) => Some(term),
+                AtomView::Mul(product) => product.get_coefficient(),
+                _ => None,
+            };
+            let value = if let Some(AtomView::Num(number)) = coefficient {
+                let Coefficient::Complex(value) = number.get_coeff_view().to_owned() else {
+                    return;
+                };
+                if !value.is_real() {
+                    return;
+                }
+                value.re
+            } else {
+                Rational::one()
+            };
+            if let Some(current) = &mut content {
+                *current = current.gcd(&value);
+            } else {
+                negative = value.is_negative();
+                content = Some(value.gcd(&value));
+            }
+        }
+        let mut content = content.unwrap();
+        if negative {
+            content = -content;
+        }
+        if content.is_one() {
+            return;
+        }
+
+        let coefficient = ws.new_num(content);
+        let mut primitive = ws.new_atom();
+        let addition = primitive.to_add();
+        let mut quotient = ws.new_atom();
+        for term in sum {
+            term.div_with_ws_into(ws, coefficient.as_view(), &mut quotient);
+            addition.extend(quotient.as_view());
+        }
+        let mut result = Atom::new();
+        primitive
+            .as_view()
+            .mul_no_norm(ws, coefficient.as_view())
+            .as_view()
+            .normalize(ws, &mut result);
+        **out = result;
+    }
+
+    pub(crate) fn horner_scheme_impl<T: AtomCore>(&self, xs: &[T]) -> Atom {
         // map every indeterminate to its position in the ordering, the first position wins
         let mut index = HashMap::default();
         for (i, x) in xs.iter().enumerate() {
-            index.entry(x.as_view()).or_insert(i);
+            index.entry(x.as_atom_view()).or_insert(i);
         }
 
         Workspace::get_local().with(|ws| {
-            let r = self.horner_scheme_impl_no_norm(ws, xs, &index, 0, enter_functions);
+            let r = self.horner_scheme_impl_no_norm(ws, xs, &index, 0);
 
             let mut out = Atom::new();
             r.as_view().normalize(ws, &mut out);
@@ -1513,14 +1666,19 @@ impl<'a> AtomView<'a> {
         index: &HashMap<AtomView<'b>, usize>,
         start: usize,
     ) -> Option<(usize, i64)> {
+        // Evaluator candidates can be whole non-polynomial expressions, such
+        // as exp(x), stored as a non-integer power. Factor only exact matches:
+        // decomposing their exponents would require additional algebraic laws.
+        if let Some(&i) = index.get(&atom)
+            && i >= start
+        {
+            return Some((i, 1));
+        }
+
         match atom {
-            AtomView::Var(_) | AtomView::Fun(_) => {
-                index.get(&atom).filter(|&&i| i >= start).map(|&i| (i, 1))
-            }
             AtomView::Pow(p) => {
                 let (b, e) = p.get_base_exp();
-                if matches!(b, AtomView::Var(_) | AtomView::Fun(_))
-                    && let AtomView::Num(n) = e
+                if let AtomView::Num(n) = e
                     && let CoefficientView::Natural(n, 1, 0, _) = n.get_coeff_view()
                     && n > 0
                     && let Some(&i) = index.get(&b)
@@ -1559,49 +1717,20 @@ impl<'a> AtomView<'a> {
 
     /// Apply the Horner scheme for the indeterminates `xs[start..]`, where `index` maps every
     /// indeterminate to its position in `xs`. The result is not normalized.
-    fn horner_scheme_impl_no_norm<'b>(
+    fn horner_scheme_impl_no_norm<'b, T: AtomCore>(
         &self,
         ws: &Workspace,
-        xs: &'b [Indeterminate],
+        xs: &'b [T],
         index: &HashMap<AtomView<'b>, usize>,
         start: usize,
-        enter_functions: bool,
     ) -> AtomOrView<'a> {
         if start >= xs.len() {
             return self.into();
         }
 
         match self {
-            AtomView::Num(_) | AtomView::Var(_) => self.into(),
-            AtomView::Fun(f) => {
-                if enter_functions {
-                    let mut tmp = ws.new_atom();
-                    let fun = tmp.to_fun(f.get_symbol());
-
-                    for arg in f {
-                        let r =
-                            arg.horner_scheme_impl_no_norm(ws, xs, index, start, enter_functions);
-                        fun.add_arg(r.as_view());
-                    }
-
-                    tmp.into_inner().into()
-                } else {
-                    self.into()
-                }
-            }
-            AtomView::Pow(p) => {
-                let (b, e) = p.get_base_exp();
-
-                let bb = b.horner_scheme_impl_no_norm(ws, xs, index, start, enter_functions);
-                let ee = e.horner_scheme_impl_no_norm(ws, xs, index, start, enter_functions);
-
-                if matches!(bb, AtomOrView::Atom(_)) || matches!(ee, AtomOrView::Atom(_)) {
-                    let mut pow = Atom::new();
-                    pow.to_pow(bb.as_view(), ee.as_view());
-                    pow.into()
-                } else {
-                    self.into()
-                }
+            AtomView::Num(_) | AtomView::Var(_) | AtomView::Fun(_) | AtomView::Pow(_) => {
+                self.into()
             }
             AtomView::Mul(m) => {
                 let mut tmp = ws.new_atom();
@@ -1609,7 +1738,7 @@ impl<'a> AtomView<'a> {
 
                 let mut changed = false;
                 for arg in m {
-                    let r = arg.horner_scheme_impl_no_norm(ws, xs, index, start, enter_functions);
+                    let r = arg.horner_scheme_impl_no_norm(ws, xs, index, start);
                     changed |= matches!(r, AtomOrView::Atom(_));
                     mul.extend(r.as_view());
                 }
@@ -1653,7 +1782,7 @@ impl<'a> AtomView<'a> {
                 let mut pow = ws.new_atom();
                 for group in grouped.chunk_by(|a, b| a.0 == b.0) {
                     let i = group[0].0;
-                    let x = xs[i].as_view();
+                    let x = xs[i].as_atom_view();
                     let min_power = group.iter().map(|(_, p, _)| *p).min().unwrap();
 
                     let mut coeff = ws.new_atom();
@@ -1689,7 +1818,7 @@ impl<'a> AtomView<'a> {
                 drop(grouped);
 
                 for (i, min_power, coeff) in coeffs {
-                    let x = xs[i].as_view();
+                    let x = xs[i].as_atom_view();
 
                     // higher powers of x may remain in the coefficient, so start from x again
                     let new_key = (if let AtomView::Add(c) = coeff.as_view()
@@ -1699,13 +1828,7 @@ impl<'a> AtomView<'a> {
                     } else {
                         coeff.as_view()
                     })
-                    .horner_scheme_impl_no_norm(
-                        ws,
-                        xs,
-                        index,
-                        i,
-                        enter_functions,
-                    );
+                    .horner_scheme_impl_no_norm(ws, xs, index, i);
 
                     let v = if min_power == 1 {
                         new_key.as_view().mul_no_norm(ws, x)
@@ -1722,7 +1845,9 @@ impl<'a> AtomView<'a> {
                     return new_arg.into_inner().into();
                 }
 
-                res.into_inner().into()
+                // Preserve numerical content at this remainder boundary, before
+                // a parent addition can flatten it and hide proportional sums.
+                res.as_view().canonicalize_horner_content(ws).into()
             }
         }
     }
@@ -1732,7 +1857,10 @@ impl<'a> AtomView<'a> {
 mod test {
     use crate::{
         atom::{Atom, AtomCore, representation::InlineVar},
-        function, parse, symbol,
+        function, parse,
+        state::Workspace,
+        symbol,
+        utils::Settable,
     };
 
     #[test]
@@ -1751,6 +1879,101 @@ mod test {
         assert_eq!(collected.expand(), expr);
     }
 
+    #[test]
+    fn horner_content_exposes_proportional_remainders() {
+        let input = parse!("2*a*x+4*a*y+3*b*x+6*b*y-5*c*x-10*c*y");
+        let collected = input.collect_horner(Some(&[
+            symbol!("a"),
+            symbol!("b"),
+            symbol!("c"),
+            symbol!("x"),
+            symbol!("y"),
+        ]));
+        assert_eq!(collected.expand(), input);
+        assert_eq!(collected, parse!("2*a*(x+2*y)+3*b*(x+2*y)-5*c*(x+2*y)"));
+    }
+
+    #[test]
+    fn horner_content_handles_unit_fractional_and_negative_coefficients() {
+        Workspace::get_local().with(|ws| {
+            for (source, expected) in [
+                ("2*x+y", "2*x+y"),
+                ("-2*x-4*y", "-2*(x+2*y)"),
+                ("2/3*x+4/3*y", "2/3*(x+2*y)"),
+                ("-x-y", "-(x+y)"),
+            ] {
+                let source = parse!(source);
+                let result = source.as_view().canonicalize_horner_content(ws);
+                assert_eq!(result.expand(), source);
+                assert_eq!(result, parse!(expected));
+            }
+        });
+    }
+
+    #[test]
+    fn horner_content_leaves_normalized_primitive_sums_unset() {
+        Workspace::get_local().with(|ws| {
+            for source in ["2*x+y", "1+2*x+4*y", "x*(a+b)+y*(c+d)"] {
+                let source = parse!(source);
+                let mut result = Atom::new();
+                let mut out = Settable::from(&mut result);
+                source
+                    .as_view()
+                    .canonicalize_horner_content_into(ws, &mut out);
+                assert!(!out.is_set(), "Unchanged content in {source}");
+            }
+
+            let source = parse!("x*(a+b)+y*(c+d)");
+            let mut changed_children = 0;
+            let result = source.as_view().replace_map_bottom_up(
+                |term, context, out| {
+                    changed_children += usize::from(context.child_changed);
+                    if matches!(term, crate::atom::AtomView::Add(_)) {
+                        term.canonicalize_horner_content_into(ws, out);
+                    }
+                },
+                true,
+            );
+            assert_eq!(result, source);
+            assert_eq!(changed_children, 0);
+        });
+    }
+
+    #[test]
+    fn horner_content_normalizes_before_extracting_and_is_idempotent() {
+        Workspace::get_local().with(|ws| {
+            let mut source = Atom::new();
+            let sum = source.to_add();
+            for term in ["2*x", "4*y", "2*x", "4"] {
+                sum.extend(parse!(term).as_view());
+            }
+            assert!(source.as_view().needs_normalization());
+            let result = source.as_view().canonicalize_horner_content(ws);
+            assert_eq!(result, parse!("4*(1+x+y)"));
+            assert_eq!(result.expand(), parse!("4+4*x+4*y"));
+
+            for source in [result, parse!("-2*x-4*y"), parse!("2/3*x+4/3*y")] {
+                let canonical = source.as_view().canonicalize_horner_content(ws);
+                assert_eq!(canonical.expand(), source.expand());
+                let mut second = Atom::new();
+                let mut out = Settable::from(&mut second);
+                canonical
+                    .as_view()
+                    .canonicalize_horner_content_into(ws, &mut out);
+                assert!(!out.is_set(), "Canonical content in {canonical}");
+            }
+
+            let mut source = Atom::new();
+            let sum = source.to_add();
+            sum.extend(parse!("x").as_view());
+            sum.extend(parse!("-x").as_view());
+            assert_eq!(
+                source.as_view().canonicalize_horner_content(ws),
+                Atom::new()
+            );
+        });
+    }
+
     #[test]
     fn collect_horner() {
         let expr = parse!("1 + v1*v2 + 2 v1*v2*v3 + v1^2 + v1^3*y + v1^4*z + v1^10");
@@ -1758,6 +1981,135 @@ mod test {
         assert_eq!(collected.expand(), expr);
     }
 
+    #[test]
+    fn local_horner_preserves_function_and_power_children() {
+        let variables = [
+            parse!("a"),
+            parse!("b"),
+            parse!("x"),
+            parse!("y"),
+            parse!("z"),
+        ];
+        for terminal in [
+            "f(2*x+2*x*y)",
+            "(2*x+2*x*y)^(z+z*x)",
+            "exp(2*x+2*x*y)",
+            "(2*x+2*x*y)^3",
+        ] {
+            let terminal = parse!(terminal);
+            let source = parse!("a+a*b") * &terminal;
+            for cleanup in [false, true] {
+                let result = source
+                    .as_view()
+                    .horner_scheme_local(Some(&variables), cleanup);
+                assert_eq!(result, parse!("a*(1+b)") * &terminal);
+                assert_eq!(result.expand(), source.expand());
+                assert_eq!(
+                    terminal
+                        .as_view()
+                        .horner_scheme_local(Some(&variables), cleanup),
+                    terminal
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn local_horner_retains_positive_integer_power_factoring() {
+        let x = parse!("x");
+        let source = parse!("x+x^2+x^3");
+        let result = source.as_view().horner_scheme_local(Some(&[x]), true);
+        assert_eq!(result, parse!("x*(1+x*(1+x))"));
+        assert_eq!(result.expand(), source);
+
+        let base = parse!("a+a*b");
+        let source = parse!("(a+a*b)^2+(a+a*b)^3");
+        let result = source.as_view().horner_scheme_local(Some(&[base]), true);
+        assert_eq!(result, parse!("(a+a*b)^2*(1+a+a*b)"));
+        assert_eq!(result.expand(), source.expand());
+    }
+
+    #[test]
+    fn local_horner_coefficient_cleanup_remains_bottom_up() {
+        for (source, expected) in [
+            ("z*(2*a+2*b)*(3*c+3*d)", "6*z*(a+b)*(c+d)"),
+            ("u*(2*a+2*b)+v*(2*a+2*b)", "2*(u*(a+b)+v*(a+b))"),
+        ] {
+            let source = parse!(source);
+            let result = source
+                .as_view()
+                .horner_scheme_local(Some(&[] as &[Atom]), true);
+            assert_eq!(result, parse!(expected));
+            assert_eq!(result.expand(), source.expand());
+            assert_eq!(
+                result
+                    .as_view()
+                    .horner_scheme_local(Some(&[] as &[Atom]), true),
+                result
+            );
+        }
+
+        Workspace::get_local().with(|ws| {
+            let source = parse!("x*(a+b)+y*f(2*a+2*b)+(2*a+2*b)^(c+c*d)");
+            for canonicalize_content in [false, true] {
+                assert!(matches!(
+                    source
+                        .as_view()
+                        .cleanup_horner_region(ws, canonicalize_content),
+                    crate::atom::AtomOrView::View(_)
+                ));
+            }
+        });
+    }
+
+    #[test]
+    fn evaluator_horner_candidates_can_include_whole_exponentials() {
+        let exponential = parse!("exp(x+y)");
+        assert!(matches!(
+            exponential.as_view(),
+            crate::atom::AtomView::Pow(_)
+        ));
+        // The evaluator's internal candidate set is broader than the public
+        // polynomial indeterminate API, whose restrictions remain unchanged.
+        assert!(crate::atom::Indeterminate::try_from(exponential.clone()).is_err());
+
+        let source = parse!("a*exp(x+y)+b*exp(x+y)+c");
+        let collected = source
+            .as_view()
+            .horner_scheme_local(Some(&[exponential]), true);
+        assert_eq!(collected, parse!("(a+b)*exp(x+y)+c"));
+        assert_eq!(collected.expand(), source);
+    }
+
+    #[test]
+    fn evaluator_horner_candidates_do_not_infer_power_relations() {
+        for (source, candidate, expected) in [
+            (
+                "a*x^(1/3)+b*x^(1/3)+c*x^(2/3)",
+                "x^(1/3)",
+                "(a+b)*x^(1/3)+c*x^(2/3)",
+            ),
+            ("a*x^y+b*x^y+c*x^(2*y)", "x^y", "(a+b)*x^y+c*x^(2*y)"),
+            (
+                "a*exp(x)+b*exp(x)+c*exp(2*x)",
+                "exp(x)",
+                "(a+b)*exp(x)+c*exp(2*x)",
+            ),
+            (
+                "a*(1+x)^(1/2)+b*(1+x)^(1/2)+c*(1+x)^(-1/2)",
+                "(1+x)^(1/2)",
+                "(a+b)*(1+x)^(1/2)+c*(1+x)^(-1/2)",
+            ),
+        ] {
+            let source = parse!(source);
+            let collected = source
+                .as_view()
+                .horner_scheme_local(Some(&[parse!(candidate)]), true);
+            assert_eq!(collected, parse!(expected));
+            assert_eq!(collected.expand(), source);
+        }
+    }
+
     #[test]
     fn collect_factors() {
         let input = parse!("v1*(v1+v2*v1+v1^2+v2*(v1+v1^2))");
diff --git a/src/evaluate.rs b/src/evaluate.rs
index 20a72c0f..9091aca9 100644
--- a/src/evaluate.rs
+++ b/src/evaluate.rs
@@ -13,7 +13,7 @@ use std::{
     panic,
     sync::{
         Arc, Mutex,
-        atomic::{AtomicBool, AtomicUsize, Ordering},
+        atomic::{AtomicBool, Ordering},
     },
 };
 #[cfg(feature = "native_code_generation")]
diff --git a/src/evaluate/evaluator.rs b/src/evaluate/evaluator.rs
index 2a6f155a..aedf4ce1 100644
--- a/src/evaluate/evaluator.rs
+++ b/src/evaluate/evaluator.rs
@@ -1706,6 +1706,42 @@ mod tests {
     }
 }
 
+#[cfg(test)]
+mod merge_operand_remapping_tests {
+    use super::*;
+
+    #[test]
+    fn merge_normalizes_remapped_operands_before_cpe() {
+        let exact = |x| Complex::new(Rational::from(x), Rational::from(0));
+        let mut evaluator = ExpressionEvaluator {
+            stack: [1i64, 2, 0].map(exact).to_vec(),
+            param_count: 0,
+            reserved_indices: 2,
+            instructions: vec![(Instr::Pow(2, 1, 2), ComplexPhase::Any)],
+            result_indices: vec![2],
+            external_fns: vec![],
+            // An existing abort skips the first CSE pass while leaving CPE enabled.
+            settings: OptimizationSettings::new().abort_level(1),
+        };
+        let other = ExpressionEvaluator {
+            stack: [2i64, 1, 2, 0].map(exact).to_vec(),
+            param_count: 0,
+            reserved_indices: 3,
+            instructions: vec![(Instr::Add(3, vec![0, 0, 1, 2, 2, 2]), ComplexPhase::Any)],
+            result_indices: vec![3],
+            external_fns: vec![],
+            settings: OptimizationSettings::new(),
+        };
+        // Remapping the incoming constants produces [1, 1, 0, 1, 1, 1].
+        // CPE must see sorted operands to compact repeated arguments correctly.
+        evaluator.merge(other, None).unwrap();
+        let mut evaluator = evaluator.map_coeff(&|x| x.re.to_f64());
+        let mut output = [0.; 2];
+        evaluator.evaluate(&[], &mut output);
+        assert_eq!(output, [4., 11.]);
+    }
+}
+
 #[cfg(test)]
 mod cse_tests {
     use super::*;
diff --git a/src/evaluate/function_map.rs b/src/evaluate/function_map.rs
index 25ba9685..a62ba4ea 100644
--- a/src/evaluate/function_map.rs
+++ b/src/evaluate/function_map.rs
@@ -399,6 +399,18 @@ impl<'a> EvaluatorBuilder<'a> {
         self
     }
 
+    /// Run common-pair elimination when scoring each Horner scheme.
+    ///
+    /// Disabled by default: schemes are ranked after common-subexpression
+    /// elimination, and common-pair elimination runs only on the selected
+    /// scheme. Enabling this can improve the chosen scheme but substantially
+    /// increases optimization time. The `cpe_iterations` setting controls
+    /// common-pair elimination in either mode.
+    pub fn horner_use_cpe(mut self, horner_use_cpe: bool) -> Self {
+        self.optimization_settings.horner_use_cpe = horner_use_cpe;
+        self
+    }
+
     /// Set the number of CPU cores to use during optimization.
     pub fn cores(mut self, n_cores: usize) -> Self {
         self.optimization_settings.n_cores = n_cores;
@@ -467,7 +479,7 @@ impl<'a> EvaluatorBuilder<'a> {
         {
             AtomView::to_evaluator(
                 &self.exprs,
-                &self.fn_map,
+                self.fn_map,
                 &self.params,
                 self.optimization_settings,
             )
@@ -500,6 +512,8 @@ pub(super) struct Expr {
 #[derive(Clone)]
 pub struct OptimizationSettings {
     pub(crate) horner_iterations: usize,
+    #[cfg_attr(feature = "serde", serde(default))]
+    pub(crate) horner_use_cpe: bool,
     pub(crate) n_cores: usize,
     pub(crate) cpe_iterations: Option<usize>,
     pub(crate) hot_start: Option<Vec<Expression<Complex<Rational>>>>,
@@ -530,7 +544,10 @@ impl bincode::Encode for OptimizationSettings {
         bincode::Encode::encode(&self.max_common_pair_cache_entries, encoder)?;
         bincode::Encode::encode(&self.max_common_pair_distance, encoder)?;
         bincode::Encode::encode(&self.verbose, encoder)?;
-        bincode::Encode::encode(&self.direct_translation, encoder)?;
+        // Preserve the legacy trailing bool byte when Horner CPE is disabled.
+        // Values 2 and 3 were invalid booleans and now encode the extra flag.
+        let flags = u8::from(self.direct_translation) | (u8::from(self.horner_use_cpe) << 1);
+        bincode::Encode::encode(&flags, encoder)?;
         Ok(())
     }
 }
@@ -540,8 +557,9 @@ impl<Context> bincode::Decode<Context> for OptimizationSettings {
     fn decode<D: bincode::de::Decoder<Context = Context>>(
         decoder: &mut D,
     ) -> core::result::Result<Self, bincode::error::DecodeError> {
-        Ok(Self {
+        let mut settings = Self {
             horner_iterations: bincode::Decode::decode(decoder)?,
+            horner_use_cpe: false,
             n_cores: bincode::Decode::decode(decoder)?,
             cpe_iterations: bincode::Decode::decode(decoder)?,
             hot_start: bincode::Decode::decode(decoder)?,
@@ -551,8 +569,17 @@ impl<Context> bincode::Decode<Context> for OptimizationSettings {
             max_common_pair_cache_entries: bincode::Decode::decode(decoder)?,
             max_common_pair_distance: bincode::Decode::decode(decoder)?,
             verbose: bincode::Decode::decode(decoder)?,
-            direct_translation: bincode::Decode::decode(decoder)?,
-        })
+            direct_translation: false,
+        };
+        let flags: u8 = bincode::Decode::decode(decoder)?;
+        if flags & !3 != 0 {
+            return Err(bincode::error::DecodeError::Other(
+                "Invalid evaluator optimization flags",
+            ));
+        }
+        settings.direct_translation = flags & 1 != 0;
+        settings.horner_use_cpe = flags & 2 != 0;
+        Ok(settings)
     }
 }
 
@@ -560,6 +587,7 @@ impl std::fmt::Debug for OptimizationSettings {
     fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
         f.debug_struct("OptimizationSettings")
             .field("horner_iterations", &self.horner_iterations)
+            .field("horner_use_cpe", &self.horner_use_cpe)
             .field("n_cores", &self.n_cores)
             .field("cpe_iterations", &self.cpe_iterations)
             .field("hot_start", &self.hot_start)
@@ -574,6 +602,7 @@ impl Default for OptimizationSettings {
     fn default() -> Self {
         OptimizationSettings {
             horner_iterations: 10,
+            horner_use_cpe: false,
             n_cores: 1,
             cpe_iterations: None,
             hot_start: None,
@@ -600,6 +629,18 @@ impl OptimizationSettings {
         self
     }
 
+    /// Run common-pair elimination when scoring each Horner scheme.
+    ///
+    /// Disabled by default: schemes are ranked after common-subexpression
+    /// elimination, and common-pair elimination runs only on the selected
+    /// scheme. Enabling this can improve the chosen scheme but substantially
+    /// increases optimization time. The `cpe_iterations` setting controls
+    /// common-pair elimination in either mode.
+    pub fn horner_use_cpe(mut self, horner_use_cpe: bool) -> Self {
+        self.horner_use_cpe = horner_use_cpe;
+        self
+    }
+
     /// Set the number of CPU cores to use during optimization.
     pub fn cores(mut self, n_cores: usize) -> Self {
         self.n_cores = n_cores;
@@ -660,3 +701,72 @@ impl OptimizationSettings {
         self
     }
 }
+
+#[cfg(all(test, feature = "bincode"))]
+mod optimization_settings_serialization_tests {
+    use super::*;
+
+    fn check_encoding<C: bincode::config::Config>(config: C) {
+        for direct_translation in [false, true] {
+            let settings = OptimizationSettings::new()
+                .horner_iterations(321)
+                .cores(3)
+                .cpe_iterations(Some(7))
+                .direct_translation(direct_translation);
+            // The old schema ended with two bools and had no Horner CPE flag.
+            let legacy = (
+                settings.horner_iterations,
+                settings.n_cores,
+                settings.cpe_iterations,
+                &settings.hot_start,
+                settings.max_horner_scheme_variables,
+                settings.max_common_pair_cache_entries,
+                settings.max_common_pair_distance,
+                settings.verbose,
+                settings.direct_translation,
+            );
+            let bytes = bincode::encode_to_vec((legacy, 0x12345678_u32), config).unwrap();
+            let ((decoded, following), consumed): ((OptimizationSettings, u32), _) =
+                bincode::decode_from_slice(&bytes, config).unwrap();
+            assert_eq!(consumed, bytes.len());
+            assert_eq!(following, 0x12345678);
+            assert_eq!(decoded.horner_iterations, 321);
+            assert_eq!(decoded.n_cores, 3);
+            assert_eq!(decoded.cpe_iterations, Some(7));
+            assert_eq!(decoded.direct_translation, direct_translation);
+            assert!(!decoded.horner_use_cpe);
+            assert_eq!(
+                bincode::encode_to_vec((&settings, 0x12345678_u32), config).unwrap(),
+                bytes
+            );
+
+            for horner_use_cpe in [false, true] {
+                let settings = settings.clone().horner_use_cpe(horner_use_cpe);
+                let bytes = bincode::encode_to_vec((&settings, 0x12345678_u32), config).unwrap();
+                let ((decoded, following), consumed): ((OptimizationSettings, u32), _) =
+                    bincode::decode_from_slice(&bytes, config).unwrap();
+                assert_eq!(consumed, bytes.len());
+                assert_eq!(following, 0x12345678);
+                assert_eq!(decoded.direct_translation, direct_translation);
+                assert_eq!(decoded.horner_use_cpe, horner_use_cpe);
+                assert_eq!(decoded.horner_iterations, settings.horner_iterations);
+            }
+
+            let mut invalid = bincode::encode_to_vec(&settings, config).unwrap();
+            *invalid.last_mut().unwrap() = 4;
+            assert!(
+                bincode::decode_from_slice::<OptimizationSettings, _>(&invalid, config).is_err()
+            );
+        }
+    }
+
+    #[test]
+    fn horner_cpe_flag_preserves_legacy_settings_and_nested_streams() {
+        check_encoding(bincode::config::standard());
+        check_encoding(
+            bincode::config::standard()
+                .with_fixed_int_encoding()
+                .with_big_endian(),
+        );
+    }
+}
diff --git a/src/evaluate/optimize.rs b/src/evaluate/optimize.rs
index fa17e3fb..6e274995 100644
--- a/src/evaluate/optimize.rs
+++ b/src/evaluate/optimize.rs
@@ -478,6 +478,13 @@ impl<T: Default + Clone + Eq + Hash> ExpressionEvaluator<T> {
 
         self.undo_stack_optimization();
 
+        // Constant and temporary remapping can change operand order. Keep
+        // normalization independent of whether an elimination pass runs.
+        for (instruction, _) in &mut self.instructions {
+            if let Instr::Add(_, args) | Instr::Mul(_, args) = instruction {
+                args.sort_unstable();
+            }
+        }
         loop {
             if self.settings.abort_level > 0 || self.remove_common_instructions() == 0 {
                 self.settings.abort_level = 0;
@@ -490,6 +497,10 @@ impl<T: Default + Clone + Eq + Hash> ExpressionEvaluator<T> {
                 self.settings.abort_level = 0;
                 break;
             }
+
+            // Reuse functions and other instructions whose operands became
+            // identical after extracting common pairs.
+            self.remove_common_instructions();
         }
 
         self.optimize_stack();
diff --git a/src/evaluate/tree.rs b/src/evaluate/tree.rs
index 378d96db..afee419d 100644
--- a/src/evaluate/tree.rs
+++ b/src/evaluate/tree.rs
@@ -1,9 +1,355 @@
 use std::{collections::HashMap as StdHashMap, hash::BuildHasher};
 
-use crate::{atom::InlineVar, coefficient::ConvertToRing};
+use crate::{
+    atom::{AtomOrView, InlineVar},
+    coefficient::ConvertToRing,
+    state::Workspace,
+};
 
 use super::*;
 
+/// Keep multiplication count as the primary objective, then minimize additions.
+fn horner_cost(ops: OperationCount) -> (usize, usize) {
+    (ops.multiplications, ops.additions)
+}
+
+/// Once search cancellation is requested, also stop optimization of the final
+/// candidate even when the caller's callback only returns true once.
+fn latched_horner_settings(
+    settings: &OptimizationSettings,
+) -> (OptimizationSettings, Arc<AtomicBool>) {
+    let cancelled = Arc::new(AtomicBool::new(false));
+    let mut search_settings = settings.clone();
+    let requested_abort = settings.abort_check.clone();
+    let abort_latch = cancelled.clone();
+    search_settings.abort_check = Some(Box::new(move || {
+        if abort_latch.load(Ordering::Relaxed) {
+            return true;
+        }
+        if requested_abort.as_ref().is_some_and(|check| check()) {
+            abort_latch.store(true, Ordering::Relaxed);
+            true
+        } else {
+            false
+        }
+    }));
+    (search_settings, cancelled)
+}
+
+/// Mix changes of different sizes: moving a variable through a block can expose
+/// a factorization that no improving single transposition reaches.
+fn mutate_horner_order<T>(order: &mut [T], rng: &mut impl Rng, mutation: usize) {
+    let size = order.len();
+    if mutation % 4 == 0 {
+        let first = rng.random_range(0..size - 1);
+        order.swap(first, first + 1);
+        return;
+    }
+    let first = rng.random_range(0..size);
+    let mut second = rng.random_range(0..size - 1);
+    if second >= first {
+        second += 1;
+    }
+    match mutation % 4 {
+        1 => order.swap(first, second),
+        2 if first < second => order[first..=second].rotate_left(1),
+        2 => order[second..=first].rotate_right(1),
+        _ => order[first.min(second)..=first.max(second)].reverse(),
+    }
+}
+
+/// Score the instructions that will be executed, including binary powering.
+/// Keep this local to the search: the public operation-count convention is
+/// unchanged.
+fn evaluator_horner_cost<T>(evaluator: &ExpressionEvaluator<T>) -> OperationCount {
+    let mut count = OperationCount::default();
+    for (instruction, _) in &evaluator.instructions {
+        match instruction {
+            Instr::Add(_, args) => count.additions += args.len() - 1,
+            Instr::Mul(_, args) => count.multiplications += args.len() - 1,
+            Instr::Pow(_, _, exponent) => {
+                count.inversions += usize::from(*exponent < 0);
+                let power = exponent.unsigned_abs();
+                if power > 1 {
+                    count.multiplications += (power.ilog2() + power.count_ones() - 1) as usize;
+                }
+            }
+            Instr::Powf(..) | Instr::BuiltinFun(..) | Instr::ExternalFun(..) => {
+                count.add_function_call();
+            }
+            _ => {}
+        }
+    }
+    count
+}
+
+fn search_horner_order<T: Clone + Ord + Send + Sync>(
+    vars: &[T],
+    settings: &OptimizationSettings,
+    score: impl Fn(&[T]) -> OperationCount + Sync,
+) -> Vec<T> {
+    search_horner_order_with_seed(vars, None, settings, score)
+}
+
+fn search_horner_order_with_seed<T: Clone + Ord + Send + Sync>(
+    vars: &[T],
+    alternative: Option<&[T]>,
+    settings: &OptimizationSettings,
+    score: impl Fn(&[T]) -> OperationCount + Sync,
+) -> Vec<T> {
+    let alternatives = alternative.into_iter().collect::<Vec<_>>();
+    search_horner_order_with_seeds(vars, &alternatives, settings, score)
+}
+
+fn search_horner_order_with_seeds<T: Clone + Ord + Send + Sync>(
+    vars: &[T],
+    alternatives: &[&[T]],
+    settings: &OptimizationSettings,
+    score: impl Fn(&[T]) -> OperationCount + Sync,
+) -> Vec<T> {
+    if settings.horner_iterations == 0 {
+        return vars.to_vec();
+    }
+    if vars.len() < 2 || settings.horner_iterations == 1 {
+        return vars.to_vec();
+    }
+
+    // Exhaustive searches may use fewer candidates than the requested
+    // iteration budget.
+    let mut scored_orders = vec![vars.to_vec()];
+    for &alternative in alternatives {
+        if scored_orders.len() >= settings.horner_iterations {
+            break;
+        }
+        if scored_orders.iter().any(|order| order == alternative) {
+            continue;
+        }
+        scored_orders.push(alternative.to_vec());
+    }
+    let permutations =
+        if vars.len() < 10 && Integer::factorial(vars.len() as u32) <= settings.horner_iterations {
+            let indices: Vec<_> = (0..vars.len()).collect();
+            // Seed orders were already scored. Visit every other permutation once.
+            Some(
+                unique_permutations(&indices)
+                    .1
+                    .into_iter()
+                    .filter(|p| {
+                        !scored_orders
+                            .iter()
+                            .any(|order| p.iter().zip(order).all(|(&i, value)| &vars[i] == value))
+                    })
+                    .collect::<Vec<_>>(),
+            )
+        } else {
+            None
+        };
+    let candidates = permutations
+        .as_ref()
+        .map_or(settings.horner_iterations - scored_orders.len(), Vec::len);
+    let mut initial = score(vars);
+
+    let mut initial_order = vars.to_vec();
+    for alternative in scored_orders.iter().skip(1) {
+        let alternative_cost = score(alternative);
+
+        if horner_cost(alternative_cost) < horner_cost(initial)
+            || horner_cost(alternative_cost) == horner_cost(initial) && alternative < &initial_order
+        {
+            initial = alternative_cost;
+            initial_order = alternative.clone();
+        }
+    }
+    let best = Mutex::new((initial, initial_order.clone()));
+    if candidates == 0 {
+        return initial_order;
+    }
+    let workers = LicenseManager::max_threads(settings.n_cores)
+        .max(1)
+        .min(candidates);
+    let abort = AtomicBool::new(false);
+    let inherited_unlock = crate::license::InheritedLibraryUnlock::capture();
+
+    std::thread::scope(|scope| {
+        for worker in 0..workers {
+            let best = &best;
+            let abort = &abort;
+            let permutations = &permutations;
+            let score = &score;
+            let initial_order = &initial_order;
+            let run = move || {
+                let _unlock = inherited_unlock.activate();
+                let mut rng = MonteCarloRng::new(0, worker);
+                let mut order = initial_order.clone();
+                let mut last = initial;
+                let worker_candidates = (candidates - 1 - worker) / workers + 1;
+                // Keep short searches in one basin. Longer searches also try
+                // reversed, shuffled, and perturbed-best starts, all charged to
+                // the same candidate budget.
+                let chain_length = if worker_candidates >= 128 {
+                    worker_candidates.div_ceil(4)
+                } else {
+                    worker_candidates
+                };
+                let mut mutation = 0;
+                // Striding includes the remainder when the budget is not divisible
+                // by the number of workers.
+                for (step, candidate) in (worker..candidates).step_by(workers).enumerate() {
+                    if abort.load(Ordering::Relaxed) {
+                        break;
+                    }
+                    if worker + 1 == workers
+                        && let Some(check) = &settings.abort_check
+                        && check()
+                    {
+                        abort.store(true, Ordering::Relaxed);
+                        break;
+                    }
+                    let previous = order.clone();
+                    let restart = permutations.is_none() && step > 0 && step % chain_length == 0;
+                    if let Some(permutations) = &permutations {
+                        order = permutations[candidate]
+                            .iter()
+                            .map(|&i| vars[i].clone())
+                            .collect();
+                    } else if restart {
+                        match step / chain_length {
+                            1 => order = vars.iter().rev().cloned().collect(),
+                            2 => {
+                                order = vars.to_vec();
+                                for i in (1..order.len()).rev() {
+                                    order.swap(i, rng.random_range(0..=i));
+                                }
+                            }
+                            _ => {
+                                order = best.lock().unwrap().1.clone();
+                                for _ in 0..3 {
+                                    mutate_horner_order(&mut order, &mut rng, mutation);
+                                    mutation += 1;
+                                }
+                            }
+                        }
+                    } else {
+                        mutate_horner_order(&mut order, &mut rng, mutation);
+                        mutation += 1;
+                    }
+                    let current = score(&order);
+
+                    if permutations.is_some()
+                        || restart
+                        || horner_cost(current) <= horner_cost(last)
+                    {
+                        last = current;
+                        // Keep the score and order under one lock: independent
+                        // atomics can expose a score assembled from different winners.
+                        let mut winner = best.lock().unwrap();
+                        if horner_cost(current) < horner_cost(winner.0)
+                            || horner_cost(current) == horner_cost(winner.0) && order < winner.1
+                        {
+                            *winner = (current, order.clone());
+                        }
+                    } else {
+                        order = previous;
+                    }
+                }
+            };
+            if worker + 1 == workers {
+                run();
+            } else {
+                scope.spawn(run);
+            }
+        }
+    });
+    let (ops, order) = best.into_inner().unwrap();
+    if settings.verbose {
+        info!("Final Horner scheme ops: {}", ops);
+    }
+
+    order
+}
+
+#[cfg(test)]
+mod horner_search_tests {
+    use super::*;
+
+    #[test]
+    fn complete_permutations_with_worker_remainder_and_deterministic_ties() {
+        let visited = Mutex::new(Vec::new());
+        let settings = OptimizationSettings {
+            horner_iterations: 24,
+            n_cores: 3,
+            ..Default::default()
+        };
+        let result = search_horner_order(&[3, 2, 1, 0], &settings, |order| {
+            visited.lock().unwrap().push(order.to_vec());
+            OperationCount::default()
+        });
+        let mut visited = visited.into_inner().unwrap();
+        assert_eq!(visited.len(), 24);
+        visited.sort();
+        visited.dedup();
+        assert_eq!(visited.len(), 24);
+        assert_eq!(result, [0, 1, 2, 3]);
+    }
+
+    #[test]
+    fn multiplication_priority_then_addition_and_initial_score() {
+        let settings = OptimizationSettings {
+            horner_iterations: 6,
+            n_cores: 1,
+            ..Default::default()
+        };
+        let result = search_horner_order(&[0, 1, 2], &settings, |order| {
+            if order == [1, 0, 2] {
+                OperationCount::new(50, 1, 0, 0)
+            } else if order == [2, 1, 0] {
+                OperationCount::new(20, 1, 0, 0)
+            } else {
+                OperationCount::new(0, 2, 0, 0)
+            }
+        });
+        assert_eq!(result, [2, 1, 0]);
+
+        let result = search_horner_order(
+            &[0, 1, 2],
+            &OptimizationSettings {
+                horner_iterations: 2,
+                ..settings
+            },
+            |order| OperationCount::new(0, usize::from(order != [0, 1, 2]), 0, 0),
+        );
+        assert_eq!(result, [0, 1, 2]);
+    }
+
+    #[test]
+    fn random_search_spends_entire_budget_without_no_op_swaps() {
+        let visited = Mutex::new(Vec::new());
+        search_horner_order(
+            &[0, 1, 2, 3],
+            &OptimizationSettings {
+                horner_iterations: 11,
+                n_cores: 1,
+                ..Default::default()
+            },
+            |order| {
+                visited.lock().unwrap().push(order.to_vec());
+                OperationCount::default()
+            },
+        );
+        let visited = visited.into_inner().unwrap();
+        assert_eq!(visited.len(), 11);
+        assert!(visited.windows(2).all(|orders| orders[0] != orders[1]));
+    }
+
+    #[test]
+    fn one_variable_needs_no_search() {
+        assert_eq!(
+            search_horner_order(&[0], &OptimizationSettings::default(), |_| panic!()),
+            [0]
+        );
+    }
+}
+
 /// Errors that can occur while building or performing numerical expression evaluation.
 #[derive(Clone, Debug, PartialEq, Eq, Hash)]
 pub enum EvaluationError {
@@ -272,6 +618,292 @@ struct LinearizationContext<'a> {
     expanding_functions: HashSet<usize>,
 }
 
+/// Optimize nested algebra regions once, bottom-up, while keeping function calls
+/// and power contents terminal to each Horner search. The cache belongs to one
+/// binding scope; exact evaluator inputs are never traversed or renamed.
+struct HornerRegions<'a, 's> {
+    tags: &'s HashMap<Symbol, usize>,
+    settings: &'s OptimizationSettings,
+    cancelled: Arc<AtomicBool>,
+    parameters: HashSet<AtomView<'a>>,
+    repeated: HashSet<AtomView<'a>>,
+    optimized: HashMap<AtomView<'a>, AtomOrView<'a>>,
+}
+
+impl<'a, 's> HornerRegions<'a, 's> {
+    fn new(
+        tags: &'s HashMap<Symbol, usize>,
+        settings: &'s OptimizationSettings,
+        cancelled: Arc<AtomicBool>,
+        parameters: impl IntoIterator<Item = AtomView<'a>>,
+        roots: &[AtomView<'a>],
+    ) -> Self {
+        let mut regions = Self {
+            tags,
+            settings,
+            cancelled,
+            parameters: parameters.into_iter().collect(),
+            repeated: HashSet::default(),
+            optimized: HashMap::default(),
+        };
+        let mut counts = HashMap::default();
+        for &root in roots {
+            regions.count_children(root, &mut counts);
+        }
+        regions.repeated = counts
+            .into_iter()
+            .filter_map(|(region, count)| (count > 1).then_some(region))
+            .collect();
+        regions
+    }
+
+    fn tag_count(&self, symbol: Symbol) -> usize {
+        self.tags.get(&symbol).copied().unwrap_or_else(|| {
+            symbol
+                .get_evaluation_info()
+                .map_or(0, |info| info.get_tag_count())
+        })
+    }
+
+    /// Count calls and power boundaries reached while preparing an algebra
+    /// region. Ordinary Add/Mul subtrees are not separate cached regions.
+    fn count_children(&self, expression: AtomView<'a>, counts: &mut HashMap<AtomView<'a>, usize>) {
+        if self.parameters.contains(&expression) {
+            return;
+        }
+        match expression {
+            AtomView::Fun(_) | AtomView::Pow(_) => self.count_region(expression, counts),
+            AtomView::Add(sum) => {
+                for term in sum {
+                    self.count_children(term, counts);
+                }
+            }
+            AtomView::Mul(product) => {
+                for factor in product {
+                    self.count_children(factor, counts);
+                }
+            }
+            _ => {}
+        }
+    }
+
+    /// A repeated whole region will reuse its first result, so its descendants
+    /// are only prepared once. Borrow all keys from the original scope roots.
+    fn count_region(&self, expression: AtomView<'a>, counts: &mut HashMap<AtomView<'a>, usize>) {
+        if matches!(expression, AtomView::Num(_) | AtomView::Var(_))
+            || self.parameters.contains(&expression)
+        {
+            return;
+        }
+        let count = counts.entry(expression).or_insert(0);
+        *count += 1;
+        if *count > 1 {
+            return;
+        }
+        match expression {
+            AtomView::Fun(function) => {
+                for argument in function.iter().skip(self.tag_count(function.get_symbol())) {
+                    self.count_region(argument, counts);
+                }
+            }
+            AtomView::Pow(power) => {
+                let (base, exponent) = power.get_base_exp();
+                self.count_region(base, counts);
+                self.count_region(exponent, counts);
+            }
+            _ => self.count_children(expression, counts),
+        }
+    }
+
+    fn cache_result(&mut self, original: AtomView<'a>, result: &AtomOrView<'a>) {
+        if self.repeated.contains(&original) {
+            self.optimized.insert(original, result.clone());
+        }
+    }
+
+    /// Borrow unchanged regions. Rebuild an Add/Mul only after its first
+    /// changed child, so large unaffected siblings are copied just once.
+    fn children(&mut self, expression: AtomView<'a>) -> AtomOrView<'a> {
+        if self.cancelled.load(Ordering::Relaxed) || self.parameters.contains(&expression) {
+            return expression.into();
+        }
+        let (is_sum, children) = match expression {
+            AtomView::Num(_) | AtomView::Var(_) => return expression.into(),
+            AtomView::Fun(function) => {
+                if let Some(cached) = self.optimized.get(&expression) {
+                    return cached.clone();
+                }
+                let tag_count = self.tag_count(function.get_symbol());
+                let arguments = function
+                    .iter()
+                    .enumerate()
+                    .map(|(i, argument)| {
+                        if i < tag_count {
+                            argument.into()
+                        } else {
+                            self.region(argument)
+                        }
+                    })
+                    .collect::<Vec<_>>();
+                let result = if function
+                    .iter()
+                    .zip(&arguments)
+                    .all(|(before, after)| before == after.as_view())
+                {
+                    expression.into()
+                } else {
+                    let views = arguments
+                        .iter()
+                        .map(AtomOrView::as_view)
+                        .collect::<Vec<_>>();
+                    function.get_symbol().call(&views).into()
+                };
+                self.cache_result(expression, &result);
+                return result;
+            }
+            AtomView::Pow(power) => {
+                if let Some(cached) = self.optimized.get(&expression) {
+                    return cached.clone();
+                }
+                let (base, exponent) = power.get_base_exp();
+                let new_base = self.region(base);
+                let new_exponent = self.region(exponent);
+                let result = if base == new_base.as_view() && exponent == new_exponent.as_view() {
+                    expression.into()
+                } else {
+                    new_base.as_view().pow(new_exponent.as_view()).into()
+                };
+                self.cache_result(expression, &result);
+                return result;
+            }
+            AtomView::Add(sum) => (true, sum.iter()),
+            AtomView::Mul(product) => (false, product.iter()),
+        };
+        let mut rebuilt: Option<Atom> = None;
+        for (i, child) in children.enumerate() {
+            let prepared = self.children(child);
+            if rebuilt.is_none() && prepared.as_view() != child {
+                let mut result = Atom::new();
+                if is_sum {
+                    let sum = result.to_add();
+                    for previous in children.take(i) {
+                        sum.extend(previous);
+                    }
+                } else {
+                    let product = result.to_mul();
+                    for previous in children.take(i) {
+                        product.extend(previous);
+                    }
+                }
+                rebuilt = Some(result);
+            }
+            if let Some(result) = &mut rebuilt {
+                match result {
+                    Atom::Add(sum) => sum.extend(prepared.as_view()),
+                    Atom::Mul(product) => product.extend(prepared.as_view()),
+                    _ => unreachable!(),
+                }
+            }
+        }
+        if let Some(rebuilt) = rebuilt {
+            Workspace::get_local().with(|ws| {
+                let mut normalized = Atom::new();
+                rebuilt.as_view().normalize(ws, &mut normalized);
+                normalized.into()
+            })
+        } else {
+            expression.into()
+        }
+    }
+
+    fn region(&mut self, expression: AtomView<'a>) -> AtomOrView<'a> {
+        if matches!(expression, AtomView::Num(_) | AtomView::Var(_))
+            || self.parameters.contains(&expression)
+            || self.cancelled.load(Ordering::Relaxed)
+        {
+            return expression.into();
+        }
+        if let Some(result) = self.optimized.get::<[u8]>(expression.get_data()) {
+            return result.clone();
+        }
+        let prepared = self.children(expression);
+        if self.cancelled.load(Ordering::Relaxed) {
+            return prepared;
+        }
+        let result = if matches!(prepared.as_view(), AtomView::Add(_) | AtomView::Mul(_)) {
+            let orders =
+                AtomView::horner_orders(&[prepared.as_view()], &self.parameters, self.settings);
+            let mut score_params = HashSet::default();
+            if self.settings.horner_use_cpe {
+                // Functions are already optimized, fixed values for this region.
+                // Use them directly as numerical inputs when scoring its CSE/CPE.
+                prepared.as_view().visitor(&mut |part| match part {
+                    AtomView::Fun(_) => {
+                        score_params.insert(part);
+                        false
+                    }
+                    AtomView::Var(variable) => {
+                        if self.parameters.contains(&part)
+                            || !matches!(
+                                variable.get_symbol().get_id(),
+                                Symbol::E_ID | Symbol::PI_ID
+                            )
+                        {
+                            score_params.insert(part);
+                        }
+                        false
+                    }
+                    _ => true,
+                });
+            }
+            let mut score_params: Vec<_> = score_params.into_iter().collect();
+            score_params.sort();
+            let order = search_horner_order_with_seeds(
+                &orders[0],
+                &[&orders[1], &orders[2]],
+                self.settings,
+                |order| {
+                    let candidate = prepared.as_view().horner_scheme_local(Some(order), true);
+                    if self.settings.horner_use_cpe {
+                        if let Ok(evaluator) = AtomView::linearize_multiple(
+                            &[candidate.as_view()],
+                            &FunctionMap::new(),
+                            &score_params,
+                            self.settings.clone(),
+                        ) {
+                            return evaluator_horner_cost(
+                                &AtomView::optimize_evaluator_instructions(evaluator),
+                            );
+                        }
+                    }
+                    // Construction errors belong to the final evaluator build,
+                    // including errors in branches which may be eliminated.
+                    candidate
+                        .as_view()
+                        .count_operations_with_subexpressions(&mut HashSet::default())
+                },
+            );
+            if order.is_empty() && !orders[0].is_empty() || self.cancelled.load(Ordering::Relaxed) {
+                prepared
+            } else {
+                prepared
+                    .as_view()
+                    .horner_scheme_local(Some(&order), true)
+                    .into()
+            }
+        } else {
+            prepared
+        };
+        let result = if result.as_view() == expression {
+            expression.into()
+        } else {
+            result
+        };
+        self.cache_result(expression, &result);
+        result
+    }
+}
+
 impl LinearizationContext<'_> {
     fn constant_slot(&mut self, constant: Complex<Rational>) -> Slot {
         let index = *self
@@ -305,9 +937,81 @@ impl LinearizationContext<'_> {
 }
 
 impl<'a> AtomView<'a> {
+    /// Candidate orders for one algebra region, without looking inside calls
+    /// or the children of powers. Each nested region gets its own search.
+    fn horner_orders<'e>(
+        expressions: &[AtomView<'e>],
+        parameters: &HashSet<AtomView<'_>>,
+        settings: &OptimizationSettings,
+    ) -> [Vec<AtomView<'e>>; 3] {
+        let mut counts = HashMap::default();
+        for expression in expressions {
+            expression.count_horner_candidates(parameters, &mut counts);
+        }
+        let mut counts: Vec<_> = counts.into_iter().filter(|(_, count)| *count > 1).collect();
+        counts.sort_by(|a, b| b.1.cmp(&a.1).then_with(|| a.0.cmp(&b.0)));
+        counts.truncate(settings.max_horner_scheme_variables);
+        let occurrence: Vec<_> = counts.into_iter().map(|(atom, _)| atom).collect();
+        let mut computed_first = occurrence.clone();
+        computed_first.sort_by_key(|atom| matches!(atom, AtomView::Var(_)));
+        let mut lexical = occurrence.clone();
+        lexical.sort_by_cached_key(|atom| atom.to_canonical_string());
+        [computed_first, occurrence, lexical]
+    }
+
+    /// Evaluation treats a non-integer power as an atomic computed value, just
+    /// like a function call. In particular exp(x) is stored as E^x, so looking
+    /// only for variables and function atoms misses a useful common factor.
+    fn count_horner_candidates(
+        self,
+        parameters: &HashSet<AtomView<'_>>,
+        counts: &mut HashMap<AtomView<'a>, usize>,
+    ) {
+        if parameters.contains(&self) {
+            *counts.entry(self).or_insert(0) += 1;
+            return;
+        }
+        match self {
+            AtomView::Num(_) => {}
+            AtomView::Var(variable) => {
+                if !matches!(variable.get_symbol().get_id(), Symbol::E_ID | Symbol::PI_ID) {
+                    *counts.entry(self).or_insert(0) += 1;
+                }
+            }
+            AtomView::Fun(_) => {
+                *counts.entry(self).or_insert(0) += 1;
+            }
+            AtomView::Pow(power) => {
+                let (base, exponent) = power.get_base_exp();
+                // The base is one indeterminate at this level. Retain positive
+                // integer degree without descending into the base's algebra.
+                if i64::try_from(exponent).is_ok_and(|power| power > 0) {
+                    if parameters.contains(&base)
+                        || base != InlineVar::new(Symbol::E).as_view()
+                            && base != InlineVar::new(Symbol::PI).as_view()
+                    {
+                        *counts.entry(base).or_insert(0) += 1;
+                    }
+                } else {
+                    *counts.entry(self).or_insert(0) += 1;
+                }
+            }
+            AtomView::Mul(product) => {
+                for factor in product {
+                    factor.count_horner_candidates(parameters, counts);
+                }
+            }
+            AtomView::Add(sum) => {
+                for term in sum {
+                    term.count_horner_candidates(parameters, counts);
+                }
+            }
+        }
+    }
+
     pub(crate) fn to_evaluator(
         expressions: &[Self],
-        fn_map: &FunctionMap,
+        mut fn_map: FunctionMap,
         params: &[Atom],
         settings: OptimizationSettings,
     ) -> Result<ExpressionEvaluator<Complex<Rational>>, EvaluationError> {
@@ -321,40 +1025,220 @@ impl<'a> AtomView<'a> {
         }
 
         if settings.horner_iterations == 0 {
-            return Self::linearize_multiple(expressions, fn_map, params, settings);
+            return Self::linearize_multiple(expressions, &fn_map, params, settings);
         }
 
-        let v = match &settings.hot_start {
-            Some(_) => {
-                return Err(EvaluationError::UnsupportedHotStart);
+        if settings.hot_start.is_some() {
+            return Err(EvaluationError::UnsupportedHotStart);
+        }
+
+        // All argument, function-body, and root searches share a one-shot abort
+        // latch, so cancellation of an inner region also stops later searches.
+        let (search_settings, cancelled) = latched_horner_settings(&settings);
+
+        // Function bodies have independent bindings and independent orders.
+        // Alias bodies inherit their callers, so retain every possible formal
+        // input identity when preparing their shared definition.
+        // Split each owned definition into its immutable bindings and mutable
+        // body. Alias bindings can then borrow arguments from all definitions
+        // while the bodies are updated in place.
+        let mut definitions: Vec<_> = fn_map
+            .tagged_fn_map
+            .values_mut()
+            .map(|function| {
+                (
+                    function.args.as_slice(),
+                    &mut function.body,
+                    function.is_alias,
+                )
+            })
+            .collect();
+        let alias_parameters: Vec<_> = params
+            .iter()
+            .map(Atom::as_view)
+            .chain(
+                definitions
+                    .iter()
+                    .flat_map(|(args, _, _)| (*args).iter().map(Indeterminate::as_view)),
+            )
+            .collect();
+        for (args, body, is_alias) in &mut definitions {
+            let parameters = if *is_alias {
+                alias_parameters.clone()
+            } else {
+                params
+                    .iter()
+                    .map(Atom::as_view)
+                    .chain(args.iter().map(Indeterminate::as_view))
+                    .collect()
+            };
+            let mut regions = HornerRegions::new(
+                &fn_map.tag,
+                &search_settings,
+                cancelled.clone(),
+                parameters,
+                &[body.as_view()],
+            );
+            let prepared = regions.region(body.as_view());
+            drop(regions);
+            if let AtomOrView::Atom(prepared) = prepared {
+                **body = prepared;
             }
-            None => {
-                // start with an occurence order Horner scheme
-                let mut v = HashMap::default();
+        }
+        drop(definitions);
+        let fn_map = &fn_map;
+        let mut regions = HornerRegions::new(
+            &fn_map.tag,
+            &search_settings,
+            cancelled.clone(),
+            params.iter().map(Atom::as_view),
+            expressions,
+        );
+        let prepared = expressions
+            .iter()
+            .map(|expression| regions.children(*expression))
+            .collect::<Vec<_>>();
+        drop(regions);
+        let expression_views = prepared.iter().map(AtomOrView::as_view).collect::<Vec<_>>();
+        let expressions = expression_views.as_slice();
+        if cancelled.load(Ordering::Relaxed) {
+            let mut evaluator =
+                Self::linearize_multiple(expressions, fn_map, params, search_settings)?;
+            evaluator.fix_labels();
+            evaluator.optimize_stack();
+            evaluator.settings = settings;
+            evaluator.settings.abort_level = 0;
+            return Ok(evaluator);
+        }
+        let parameters = params.iter().map(Atom::as_view).collect();
+        let [computed_first, v, lexical] = Self::horner_orders(expressions, &parameters, &settings);
+        if settings.horner_iterations == 1 || computed_first.len() < 2 {
+            let result = Self::build_horner_evaluator(
+                expressions,
+                fn_map,
+                params,
+                &computed_first,
+                search_settings,
+            );
 
-                for t in expressions {
-                    t.count_indeterminates(true, &mut v);
-                }
+            let mut evaluator = result?;
+            evaluator.optimize_stack();
+            evaluator.settings = settings;
+            evaluator.settings.abort_level = 0;
+            return Ok(evaluator);
+        }
 
-                let mut v: Vec<_> = v.into_iter().collect();
-                v.retain(|(_, vv)| *vv > 1);
-                v.sort_by(|a, b| b.1.cmp(&a.1).then_with(|| a.0.cmp(&b.0)));
-                v.truncate(settings.max_horner_scheme_variables);
-                v.into_iter()
-                    .map(|(k, _)| Indeterminate::try_from(k.to_owned()).unwrap())
-                    .collect::<Vec<_>>()
-            }
-        };
+        if !settings.horner_use_cpe {
+            let scheme = search_horner_order_with_seeds(
+                &computed_first,
+                &[&v, &lexical],
+                &search_settings,
+                |order| Self::horner_scheme_cost(expressions, order),
+            );
+            // The inexpensive search still receives the complete CSE/CPE
+            // pipeline once, after selecting its winning Horner scheme.
+            let mut evaluator = if scheme.is_empty() {
+                Self::build_unhornered_evaluator(expressions, fn_map, params, search_settings)?
+            } else {
+                Self::build_horner_evaluator(expressions, fn_map, params, &scheme, search_settings)?
+            };
+            evaluator.settings = settings;
+            evaluator.settings.abort_level = 0;
+            evaluator.optimize_stack();
 
-        let scheme = if settings.horner_iterations > 1 {
-            Self::optimize_horner_scheme_multiple(expressions, &v, &settings)
-        } else {
-            v
+            return Ok(evaluator);
+        }
+
+        // Build and score the final CSE/CPE program for every candidate. A cheap
+        // tree score can reject a scheme whose new common pairs make the final
+        // evaluator smaller. Retain only the winning program, avoiding a final
+        // rebuild and keeping memory bounded by one candidate per worker.
+        let best_evaluator = Mutex::new(
+            None::<(
+                OperationCount,
+                Vec<AtomView<'_>>,
+                ExpressionEvaluator<Complex<Rational>>,
+            )>,
+        );
+        let first_error = Mutex::new(None);
+        let mut candidate_settings = search_settings.clone();
+        candidate_settings.verbose = false;
+        search_horner_order_with_seeds(
+            &computed_first,
+            &[&v, &lexical],
+            &search_settings,
+            |order| {
+                let unavailable = OperationCount::new(0, usize::MAX, 0, 0);
+                // The first candidate is built even if cancellation was requested,
+                // so an interrupted optimization still returns a usable evaluator.
+                if first_error.lock().unwrap().is_some()
+                    || cancelled.load(Ordering::Relaxed) && best_evaluator.lock().unwrap().is_some()
+                {
+                    return unavailable;
+                }
+                let evaluator = match Self::build_horner_evaluator(
+                    expressions,
+                    fn_map,
+                    params,
+                    order,
+                    candidate_settings.clone(),
+                ) {
+                    Ok(evaluator) => evaluator,
+                    Err(error) => {
+                        first_error.lock().unwrap().get_or_insert(error);
+                        cancelled.store(true, Ordering::Relaxed);
+                        return unavailable;
+                    }
+                };
+                let cost = evaluator_horner_cost(&evaluator);
+                let mut winner = best_evaluator.lock().unwrap();
+                if winner.as_ref().is_none_or(|(best_cost, best_order, _)| {
+                    horner_cost(cost) < horner_cost(*best_cost)
+                        || horner_cost(cost) == horner_cost(*best_cost)
+                            && order < best_order.as_slice()
+                }) {
+                    *winner = Some((cost, order.to_vec(), evaluator));
+                }
+                cost
+            },
+        );
+        if let Some(error) = first_error.into_inner().unwrap() {
+            return Err(error);
+        }
+        let mut evaluator = match best_evaluator.into_inner().unwrap() {
+            Some((_, _, evaluator)) => evaluator,
+            None => Self::build_unhornered_evaluator(expressions, fn_map, params, search_settings)?,
         };
+        evaluator.settings = settings;
+        evaluator.settings.abort_level = 0;
+        evaluator.optimize_stack();
+
+        Ok(evaluator)
+    }
+
+    /// Skipping the Horner stage still permits later optimization stages.
+    fn build_unhornered_evaluator(
+        expressions: &[AtomView<'_>],
+        fn_map: &FunctionMap,
+        params: &[Atom],
+        settings: OptimizationSettings,
+    ) -> Result<ExpressionEvaluator<Complex<Rational>>, EvaluationError> {
+        let evaluator = Self::linearize_multiple(expressions, fn_map, params, settings)?;
+        Ok(Self::optimize_evaluator_instructions(evaluator))
+    }
 
+    /// Construct one candidate without changing its variable order. Stack
+    /// allocation is delayed until the winning candidate is selected.
+    fn build_horner_evaluator<T: AtomCore>(
+        expressions: &[AtomView<'_>],
+        fn_map: &FunctionMap,
+        params: &[Atom],
+        scheme: &[T],
+        settings: OptimizationSettings,
+    ) -> Result<ExpressionEvaluator<Complex<Rational>>, EvaluationError> {
         let hornered_expressions = expressions
             .iter()
-            .map(|x| x.horner_scheme(Some(&scheme), true, true))
+            .map(|x| x.horner_scheme_local(Some(scheme), true))
             .collect::<Vec<_>>();
 
         if settings.horner_iterations == 1 && settings.verbose {
@@ -366,16 +1250,16 @@ impl<'a> AtomView<'a> {
             info!("Horner scheme ops: {}", count);
         }
 
-        let mut f = fn_map.clone();
-        for Expr { body, .. } in f.tagged_fn_map.values_mut() {
-            *body = body.as_view().horner_scheme(Some(&scheme), true, true);
-        }
-
-        let mut e = Self::linearize_multiple(&hornered_expressions, &f, params, settings)?;
+        let e = Self::linearize_multiple(&hornered_expressions, fn_map, params, settings)?;
 
-        drop(f);
         drop(hornered_expressions);
 
+        Ok(Self::optimize_evaluator_instructions(e))
+    }
+
+    fn optimize_evaluator_instructions(
+        mut e: ExpressionEvaluator<Complex<Rational>>,
+    ) -> ExpressionEvaluator<Complex<Rational>> {
         loop {
             let r = e.remove_common_instructions();
 
@@ -397,14 +1281,17 @@ impl<'a> AtomView<'a> {
                 break;
             }
 
+            // Pair extraction can make previously different instruction DAGs
+            // identical, including functions and powers that CPE cannot merge.
+            e.remove_common_instructions();
+
             if e.settings.verbose {
                 let count = e.count_operations();
                 info!("Removed {} common pairs: {}", r, count);
             }
         }
 
-        e.optimize_stack();
-        Ok(e)
+        e
     }
 
     /// Search for a Horner variable order that reduces the combined operation
@@ -415,198 +1302,40 @@ impl<'a> AtomView<'a> {
         vars: &[Indeterminate],
         settings: &OptimizationSettings,
     ) -> Vec<Indeterminate> {
-        if vars.is_empty() {
-            return vars.to_vec();
-        }
+        search_horner_order(vars, settings, |order| {
+            Self::horner_scheme_cost(expressions, order)
+        })
+    }
 
-        let horner: Vec<_> = expressions
+    fn horner_scheme_cost<T: AtomCore>(
+        expressions: &[AtomView<'_>],
+        order: &[T],
+    ) -> OperationCount {
+        let hornered: Vec<_> = expressions
             .iter()
-            .map(|x| x.horner_scheme(Some(vars), true, true))
+            .map(|expression| expression.horner_scheme_local(Some(order), true))
             .collect();
         let mut subexpr = HashSet::default();
-        let mut best_ops = OperationCount::default();
-        for h in &horner {
-            let ops = h
+        let mut operations = OperationCount::default();
+        for expression in &hornered {
+            operations += expression
                 .as_view()
                 .count_operations_with_subexpressions(&mut subexpr);
-            best_ops += ops;
-        }
-
-        if settings.verbose {
-            info!("Initial Horner scheme ops: {}", best_ops);
-        }
-
-        let best_mul = Arc::new(AtomicUsize::new(best_ops.multiplications));
-        let best_add = Arc::new(AtomicUsize::new(best_ops.additions));
-        let best_scheme = Arc::new(Mutex::new(vars.to_vec()));
-
-        let n_iterations = settings.horner_iterations.max(1) - 1;
-
-        let permutations = if vars.len() < 10
-            && Integer::factorial(vars.len() as u32) <= settings.horner_iterations.max(1)
-        {
-            let v: Vec<_> = (0..vars.len()).collect();
-            Some(unique_permutations(&v).1)
-        } else {
-            None
-        };
-        let p_ref = &permutations;
-
-        let n_cores = LicenseManager::max_threads(settings.n_cores).min(n_iterations);
-        let inherited_unlock = crate::license::InheritedLibraryUnlock::capture();
-
-        std::thread::scope(|s| {
-            let abort = Arc::new(AtomicBool::new(false));
-
-            for i in 0..n_cores {
-                let mut rng = MonteCarloRng::new(0, i);
-
-                let mut cvars = vars.to_vec();
-                let best_scheme = best_scheme.clone();
-                let best_mul = best_mul.clone();
-                let best_add = best_add.clone();
-                let mut last_mul = usize::MAX;
-                let mut last_add = usize::MAX;
-                let abort = abort.clone();
-
-                let mut op = move || {
-                    let _unlock = inherited_unlock.activate();
-                    for j in 0..n_iterations / n_cores {
-                        if abort.load(Ordering::Relaxed) {
-                            return;
-                        }
-
-                        if i == n_cores - 1
-                            && let Some(a) = &settings.abort_check
-                            && a()
-                        {
-                            abort.store(true, Ordering::Relaxed);
-
-                            if settings.verbose {
-                                info!(
-                                    "Aborting Horner optimization at step {}/{}.",
-                                    j,
-                                    settings.horner_iterations / n_cores
-                                );
-                            }
-
-                            return;
-                        }
-
-                        // try a random swap
-                        let mut t1 = 0;
-                        let mut t2 = 0;
-
-                        if let Some(p) = p_ref {
-                            if j >= p.len() / n_cores {
-                                break;
-                            }
-
-                            let perm = &p[i * (p.len() / n_cores) + j];
-                            cvars = perm.iter().map(|x| vars[*x].clone()).collect();
-                        } else {
-                            t1 = rng.random_range(0..cvars.len());
-                            t2 = rng.random_range(0..cvars.len() - 1);
-
-                            cvars.swap(t1, t2);
-                        }
-
-                        let horner: Vec<_> = expressions
-                            .iter()
-                            .map(|x| x.horner_scheme(Some(&cvars), true, true))
-                            .collect();
-                        let mut subexpr = HashSet::default();
-                        let mut cur_ops = OperationCount::default();
-
-                        for h in &horner {
-                            let ops = h
-                                .as_view()
-                                .count_operations_with_subexpressions(&mut subexpr);
-                            cur_ops += ops;
-                        }
-
-                        // prefer fewer multiplications
-                        if cur_ops.multiplications <= last_mul
-                            || cur_ops.multiplications == last_mul && cur_ops.additions <= last_add
-                        {
-                            if settings.verbose {
-                                info!(
-                                    "Accept move at step {}/{}: {}",
-                                    j,
-                                    settings.horner_iterations / n_cores,
-                                    cur_ops
-                                );
-                            }
-
-                            last_add = cur_ops.additions;
-                            last_mul = cur_ops.multiplications;
-
-                            if cur_ops.multiplications <= best_mul.load(Ordering::Relaxed)
-                                || cur_ops.multiplications == best_mul.load(Ordering::Relaxed)
-                                    && cur_ops.additions <= best_add.load(Ordering::Relaxed)
-                            {
-                                let mut best_scheme = best_scheme.lock().unwrap();
-
-                                // check again if it is the best now that we have locked
-                                let best_mul_l = best_mul.load(Ordering::Relaxed);
-                                let best_add_l = best_add.load(Ordering::Relaxed);
-                                if cur_ops.multiplications <= best_mul_l
-                                    || cur_ops.multiplications == best_mul_l
-                                        && cur_ops.additions <= best_add_l
-                                {
-                                    if cur_ops.additions == best_add_l
-                                        && cur_ops.multiplications == best_mul_l
-                                    {
-                                        if *best_scheme < cvars {
-                                            // on a draw, accept the lexicographical minimum
-                                            // to get a deterministic scheme
-                                            *best_scheme = cvars.clone();
-                                        }
-                                    } else {
-                                        best_mul.store(cur_ops.multiplications, Ordering::Relaxed);
-                                        best_add.store(cur_ops.additions, Ordering::Relaxed);
-                                        *best_scheme = cvars.clone();
-                                    }
-                                }
-                            }
-                        } else {
-                            cvars.swap(t1, t2);
-                        }
-                    }
-                };
-
-                if i + 1 < n_cores {
-                    s.spawn(op);
-                } else {
-                    // execute in the main thread and do the abort check on the main thread
-                    // this helps with catching ctrl-c
-                    op()
-                }
-            }
-        });
-
-        if settings.verbose {
-            info!(
-                "Final scheme: {} + and {} ×",
-                best_add.load(Ordering::Relaxed),
-                best_mul.load(Ordering::Relaxed)
-            );
         }
-
-        Arc::try_unwrap(best_scheme).unwrap().into_inner().unwrap()
+        operations
     }
 
-    pub(crate) fn linearize_multiple<T: AtomCore>(
+    pub(crate) fn linearize_multiple<T: AtomCore, P: AtomCore>(
         expressions: &[T],
         fn_map: &FunctionMap,
-        params: &[Atom],
+        params: &[P],
         settings: OptimizationSettings,
     ) -> Result<ExpressionEvaluator<Complex<Rational>>, EvaluationError> {
         let expressions = expressions
             .iter()
             .map(|e| e.as_atom_view())
             .collect::<Vec<_>>();
-        let params = params.iter().map(|p| p.as_view()).collect::<Vec<_>>();
+        let params = params.iter().map(|p| p.as_atom_view()).collect::<Vec<_>>();
         let mut context = LinearizationContext {
             external_functions: Vec::new(),
             functions: HashMap::default(),
@@ -2135,6 +2864,16 @@ impl EvalTree<Complex<Rational>> {
     pub fn linearize(
         &mut self,
         settings: &OptimizationSettings,
+    ) -> ExpressionEvaluator<Complex<Rational>> {
+        let mut evaluator = self.linearize_without_stack_optimization(&settings);
+        evaluator.optimize_stack();
+
+        evaluator
+    }
+
+    fn linearize_without_stack_optimization(
+        &mut self,
+        settings: &OptimizationSettings,
     ) -> ExpressionEvaluator<Complex<Rational>> {
         // strip every constant and move them into the stack after the params
         let mut constants = vec![];
@@ -2173,35 +2912,10 @@ impl EvalTree<Complex<Rational>> {
             settings: settings.clone(),
         };
 
-        loop {
-            let r = e.remove_common_instructions();
-
-            if r == 0 || e.settings.abort_level > 0 {
-                e.settings.abort_level = 0;
-                break;
-            }
-
-            if settings.verbose {
-                let count = e.count_operations();
-                info!("Removed {} common instructions: {}", r, count);
-            }
-        }
-
-        for _ in 0..settings.cpe_iterations.unwrap_or(usize::MAX) {
-            let r = e.remove_common_pairs();
-            if r == 0 || e.settings.abort_level > 0 {
-                e.settings.abort_level = 0;
-                break;
-            }
-
-            if settings.verbose {
-                let count = e.count_operations();
-                info!("Removed {} common pairs: {}", r, count);
-            }
-        }
-
-        e.optimize_stack();
-        e
+        // Resolving branch destinations is required lowering, even when all
+        // optional CSE/CPE stages are interrupted by cancellation.
+        e.fix_labels();
+        AtomView::optimize_evaluator_instructions(e)
     }
 
     fn strip_constants(&mut self, constants: &mut Vec<Complex<Rational>>, param_count: usize) {
@@ -2627,12 +3341,26 @@ impl EvalTree<Complex<Rational>> {
             info!("Initial ops: {}", count);
         }
 
+        if settings.horner_use_cpe && settings.horner_iterations > 1 {
+            let (_, mut evaluator) = self.optimize_horner_scheme_with_cpe(settings);
+            evaluator.settings = settings.clone();
+            evaluator.settings.abort_level = 0;
+            evaluator.optimize_stack();
+
+            return evaluator;
+        }
+
+        let (search_settings, _) = latched_horner_settings(settings);
         if settings.horner_iterations > 0 {
-            let _ = self.optimize_horner_scheme(settings);
+            let _ = self.optimize_horner_scheme(&search_settings);
         }
 
         self.common_subexpression_elimination();
-        self.linearize(settings)
+        let mut evaluator = self.linearize(&search_settings);
+        evaluator.settings = settings.clone();
+        evaluator.settings.abort_level = 0;
+
+        evaluator
     }
 
     /// Write the expressions in a Horner scheme where the variables
@@ -2664,6 +3392,10 @@ impl EvalTree<Complex<Rational>> {
         &mut self,
         settings: &OptimizationSettings,
     ) -> Vec<Expression<Complex<Rational>>> {
+        if settings.horner_use_cpe && settings.horner_iterations > 1 {
+            return self.optimize_horner_scheme_with_cpe(settings).0;
+        }
+
         let v = match &settings.hot_start {
             Some(a) => a.clone(),
             None => {
@@ -2731,7 +3463,121 @@ impl EvalTree<Complex<Rational>> {
             }
         }
 
-        scheme
+        scheme
+    }
+
+    /// Preserve the legacy path's independent root/function searches, but score
+    /// each candidate in the complete program where its function calls occur.
+    fn optimize_horner_scheme_with_cpe(
+        &mut self,
+        settings: &OptimizationSettings,
+    ) -> (
+        Vec<Expression<Complex<Rational>>>,
+        ExpressionEvaluator<Complex<Rational>>,
+    ) {
+        let (search_settings, cancelled) = latched_horner_settings(settings);
+        let mut candidate_settings = search_settings.clone();
+        candidate_settings.verbose = false;
+        let mut root_scheme = Vec::new();
+        let mut final_evaluator = None;
+
+        for scope in 0..=self.functions.len() {
+            if cancelled.load(Ordering::Relaxed) && final_evaluator.is_some() {
+                break;
+            }
+            let expressions = if scope == 0 {
+                &self.expressions
+            } else {
+                &self.functions[scope - 1].2
+            };
+            let vars = if scope == 0 && settings.hot_start.is_some() {
+                settings.hot_start.clone().unwrap()
+            } else {
+                let mut counts = HashMap::default();
+                for expression in expressions.tree.iter().chain(&expressions.subexpressions) {
+                    expression.find_all_variables(&mut counts);
+                }
+                let mut counts: Vec<_> = counts.into_iter().collect();
+                counts.retain(|(_, count)| *count > 1);
+                counts.sort_by(|a, b| b.1.cmp(&a.1).then_with(|| a.0.cmp(&b.0)));
+                counts.truncate(settings.max_horner_scheme_variables);
+                counts
+                    .into_iter()
+                    .map(|(variable, _)| variable)
+                    .collect::<Vec<_>>()
+            };
+            // An empty scheme leaves this scope unchanged, so a previously
+            // compiled winner is still valid.
+            if vars.is_empty() && final_evaluator.is_some() {
+                continue;
+            }
+            let best = Mutex::new(
+                None::<(
+                    OperationCount,
+                    Vec<Expression<Complex<Rational>>>,
+                    ExpressionEvaluator<Complex<Rational>>,
+                )>,
+            );
+            let score = |order: &[Expression<Complex<Rational>>]| {
+                if cancelled.load(Ordering::Relaxed) && best.lock().unwrap().is_some() {
+                    return OperationCount::new(0, usize::MAX, 0, 0);
+                }
+                let mut candidate = self.clone();
+                candidate.apply_horner_scheme_to_scope(scope, order);
+                candidate.common_subexpression_elimination();
+                let evaluator = candidate.linearize_without_stack_optimization(&candidate_settings);
+                let cost = evaluator_horner_cost(&evaluator);
+                let mut winner = best.lock().unwrap();
+                if winner.as_ref().is_none_or(|(best_cost, best_order, _)| {
+                    horner_cost(cost) < horner_cost(*best_cost)
+                        || horner_cost(cost) == horner_cost(*best_cost)
+                            && order < best_order.as_slice()
+                }) {
+                    *winner = Some((cost, order.to_vec(), evaluator));
+                }
+                cost
+            };
+            let scheme = if vars.len() < 2 {
+                score(&vars);
+                vars
+            } else {
+                search_horner_order(&vars, &search_settings, score)
+            };
+            if let Some((_, _, evaluator)) = best.into_inner().unwrap() {
+                self.apply_horner_scheme_to_scope(scope, &scheme);
+                if scope == 0 {
+                    root_scheme = scheme;
+                }
+                final_evaluator = Some(evaluator);
+            }
+        }
+        (
+            root_scheme,
+            final_evaluator.unwrap_or_else(|| {
+                let mut candidate = self.clone();
+                candidate.common_subexpression_elimination();
+                candidate.linearize_without_stack_optimization(&search_settings)
+            }),
+        )
+    }
+
+    fn apply_horner_scheme_to_scope(
+        &mut self,
+        scope: usize,
+        scheme: &[Expression<Complex<Rational>>],
+    ) {
+        let expressions = if scope == 0 {
+            &mut self.expressions
+        } else {
+            &mut self.functions[scope - 1].2
+        };
+        for expression in expressions
+            .tree
+            .iter_mut()
+            .chain(&mut expressions.subexpressions)
+        {
+            expression.apply_horner_scheme(scheme);
+        }
     }
 }
 
@@ -3128,190 +3974,23 @@ impl Expression<Complex<Rational>> {
         vars: &[Self],
         settings: &OptimizationSettings,
     ) -> Vec<Self> {
-        if vars.is_empty() {
-            return vars.to_vec();
-        }
-
-        let horner: Vec<_> = expressions
-            .iter()
-            .map(|x| {
-                let mut h = x.clone();
-                h.apply_horner_scheme(vars);
-                h.rehashed(true)
-            })
-            .collect();
-        let mut subexpr = HashMap::default();
-        let mut best_ops = OperationCount::default();
-        for h in &horner {
-            let ops = h.count_operations_with_subexpression(&mut subexpr);
-            best_ops += ops;
-        }
-
-        if settings.verbose {
-            info!("Initial Horner scheme ops: {}", best_ops);
-        }
-
-        let best_mul = Arc::new(AtomicUsize::new(best_ops.multiplications));
-        let best_add = Arc::new(AtomicUsize::new(best_ops.additions));
-        let best_scheme = Arc::new(Mutex::new(vars.to_vec()));
-
-        let n_iterations = settings.horner_iterations.max(1) - 1;
-
-        let permutations = if vars.len() < 10
-            && Integer::factorial(vars.len() as u32) <= settings.horner_iterations.max(1)
-        {
-            let v: Vec<_> = (0..vars.len()).collect();
-            Some(unique_permutations(&v).1)
-        } else {
-            None
-        };
-        let p_ref = &permutations;
-
-        let n_cores = LicenseManager::max_threads(settings.n_cores).min(n_iterations);
-        let inherited_unlock = crate::license::InheritedLibraryUnlock::capture();
-
-        std::thread::scope(|s| {
-            let abort = Arc::new(AtomicBool::new(false));
-
-            for i in 0..n_cores {
-                let mut rng = MonteCarloRng::new(0, i);
-
-                let mut cvars = vars.to_vec();
-                let best_scheme = best_scheme.clone();
-                let best_mul = best_mul.clone();
-                let best_add = best_add.clone();
-                let mut last_mul = usize::MAX;
-                let mut last_add = usize::MAX;
-                let abort = abort.clone();
-
-                let mut op = move || {
-                    let _unlock = inherited_unlock.activate();
-                    for j in 0..n_iterations / n_cores {
-                        if abort.load(Ordering::Relaxed) {
-                            return;
-                        }
-
-                        if i == n_cores - 1
-                            && let Some(a) = &settings.abort_check
-                            && a()
-                        {
-                            abort.store(true, Ordering::Relaxed);
-
-                            if settings.verbose {
-                                info!(
-                                    "Aborting Horner optimization at step {}/{}.",
-                                    j,
-                                    settings.horner_iterations / n_cores
-                                );
-                            }
-
-                            return;
-                        }
-
-                        // try a random swap
-                        let mut t1 = 0;
-                        let mut t2 = 0;
-
-                        if let Some(p) = p_ref {
-                            if j >= p.len() / n_cores {
-                                break;
-                            }
-
-                            let perm = &p[i * (p.len() / n_cores) + j];
-                            cvars = perm.iter().map(|x| vars[*x].clone()).collect();
-                        } else {
-                            t1 = rng.random_range(0..cvars.len());
-                            t2 = rng.random_range(0..cvars.len() - 1);
-
-                            cvars.swap(t1, t2);
-                        }
-
-                        let horner: Vec<_> = expressions
-                            .iter()
-                            .map(|x| {
-                                let mut h = x.clone();
-                                h.apply_horner_scheme(&cvars);
-                                h.rehash(true);
-                                h
-                            })
-                            .collect();
-                        let mut subexpr = HashMap::default();
-                        let mut cur_ops = OperationCount::default();
-
-                        for h in &horner {
-                            let ops = h.count_operations_with_subexpression(&mut subexpr);
-                            cur_ops += ops;
-                        }
-
-                        // prefer fewer multiplications
-                        if cur_ops.multiplications <= last_mul
-                            || cur_ops.multiplications == last_mul && cur_ops.additions <= last_add
-                        {
-                            if settings.verbose {
-                                info!(
-                                    "Accept move at step {}/{}: {}",
-                                    j,
-                                    settings.horner_iterations / n_cores,
-                                    cur_ops
-                                );
-                            }
-
-                            last_add = cur_ops.additions;
-                            last_mul = cur_ops.multiplications;
-
-                            if cur_ops.multiplications <= best_mul.load(Ordering::Relaxed)
-                                || cur_ops.multiplications == best_mul.load(Ordering::Relaxed)
-                                    && cur_ops.additions <= best_add.load(Ordering::Relaxed)
-                            {
-                                let mut best_scheme = best_scheme.lock().unwrap();
-
-                                // check again if it is the best now that we have locked
-                                let best_mul_l = best_mul.load(Ordering::Relaxed);
-                                let best_add_l = best_add.load(Ordering::Relaxed);
-                                if cur_ops.multiplications <= best_mul_l
-                                    || cur_ops.multiplications == best_mul_l
-                                        && cur_ops.additions <= best_add_l
-                                {
-                                    if cur_ops.additions == best_add_l
-                                        && cur_ops.multiplications == best_mul_l
-                                    {
-                                        if *best_scheme < cvars {
-                                            // on a draw, accept the lexicographical minimum
-                                            // to get a deterministic scheme
-                                            *best_scheme = cvars.clone();
-                                        }
-                                    } else {
-                                        best_mul.store(cur_ops.multiplications, Ordering::Relaxed);
-                                        best_add.store(cur_ops.additions, Ordering::Relaxed);
-                                        *best_scheme = cvars.clone();
-                                    }
-                                }
-                            }
-                        } else {
-                            cvars.swap(t1, t2);
-                        }
-                    }
-                };
-
-                if i + 1 < n_cores {
-                    s.spawn(op);
-                } else {
-                    // execute in the main thread and do the abort check on the main thread
-                    // this helps with catching ctrl-c
-                    op()
-                }
+        search_horner_order(vars, settings, |order| {
+            let hornered: Vec<_> = expressions
+                .iter()
+                .map(|expression| {
+                    let mut hornered = expression.clone();
+                    hornered.apply_horner_scheme(order);
+                    hornered.rehash(true);
+                    hornered
+                })
+                .collect();
+            let mut subexpr = HashMap::default();
+            let mut operations = OperationCount::default();
+            for expression in &hornered {
+                operations += expression.count_operations_with_subexpression(&mut subexpr);
             }
-        });
-
-        if settings.verbose {
-            info!(
-                "Final scheme: {} + and {} ×",
-                best_add.load(Ordering::Relaxed),
-                best_mul.load(Ordering::Relaxed)
-            );
-        }
-
-        Arc::try_unwrap(best_scheme).unwrap().into_inner().unwrap()
+            operations
+        })
     }
 
     fn find_all_variables(&self, vars: &mut HashMap<Expression<Complex<Rational>>, usize>) {
@@ -4580,6 +5259,110 @@ impl<'a> AtomView<'a> {
     }
 }
 
+#[cfg(test)]
+mod post_cpe_cse_tests {
+    use super::*;
+
+    type Expr = Expression<Complex<Rational>>;
+
+    fn parameter(index: usize) -> Expr {
+        Expression::Parameter(0, index)
+    }
+
+    fn sum(terms: Vec<Expr>) -> Expr {
+        Expression::Add(0, terms)
+    }
+
+    fn sine(argument: Expr) -> Expr {
+        Expression::BuiltinFun(0, Symbol::SIN, Box::new(argument))
+    }
+
+    fn equivalent_sums() -> [Expr; 2] {
+        [
+            sum(vec![parameter(0), parameter(1), parameter(2)]),
+            sum(vec![sum(vec![parameter(0), parameter(1)]), parameter(2)]),
+        ]
+    }
+
+    fn build(
+        expressions: Vec<Expr>,
+        param_count: usize,
+        cpe_iterations: Option<usize>,
+    ) -> ExpressionEvaluator<Complex<Rational>> {
+        EvalTree {
+            functions: vec![],
+            external_functions: vec![],
+            expressions: SplitExpression {
+                tree: expressions,
+                subexpressions: vec![],
+            },
+            param_count,
+        }
+        .linearize(&OptimizationSettings {
+            cpe_iterations,
+            ..Default::default()
+        })
+    }
+
+    #[test]
+    fn pair_extraction_exposes_duplicate_functions_in_tree() {
+        // Different addition trees initially hide the identical sine calls.
+        // Extracting (a+b) makes both argument DAGs identical, even with one
+        // CPE round. CSE must then propagate that identity through the sine.
+        for rounds in [Some(1), None] {
+            let mut evaluator = build(equivalent_sums().map(sine).to_vec(), 3, rounds)
+                .map_coeff(&|coefficient| coefficient.re.to_f64());
+            assert_eq!(evaluator.count_operations().function_calls, 1);
+            assert_eq!(evaluator.count_operations().additions, 2);
+            for point in [[0.2, 0.3, 0.4], [-0.7, 0.1, -0.2], [0., 0., 0.]] {
+                let mut outputs = [0.; 2];
+                evaluator.evaluate(&point, &mut outputs);
+                let expected = (point[0] + point[1] + point[2]).sin();
+                assert_eq!(outputs[0], outputs[1]);
+                assert!((outputs[0] - expected).abs() < 1e-14);
+            }
+        }
+    }
+
+    #[test]
+    fn pair_extraction_exposes_duplicate_functions_after_merge() {
+        let [flat, nested] = equivalent_sums();
+        let mut evaluator = build(vec![sine(flat)], 3, Some(0));
+        let other = build(vec![sine(nested)], 3, Some(0));
+        evaluator.merge(other, Some(1)).unwrap();
+        assert_eq!(evaluator.count_operations().function_calls, 1);
+        let mut evaluator = evaluator.map_coeff(&|coefficient| coefficient.re.to_f64());
+        let mut outputs = [0.; 2];
+        evaluator.evaluate(&[0.2, 0.3, 0.4], &mut outputs);
+        assert_eq!(outputs[0], outputs[1]);
+        assert!((outputs[0] - 0.9_f64.sin()).abs() < 1e-14);
+    }
+
+    #[test]
+    fn pair_extraction_merges_functions_only_in_dominating_branches() {
+        let branch = || sum(equivalent_sums().map(sine).to_vec());
+        let conditional = Expression::IfElse(
+            0,
+            Box::new((parameter(3), branch(), sum(vec![branch(), parameter(4)]))),
+        );
+        let mut evaluator =
+            build(vec![conditional], 5, None).map_coeff(&|coefficient| coefficient.re.to_f64());
+        // Each branch can share its two sine calls, but a sine evaluated only
+        // in the true branch cannot supply the false branch (or vice versa).
+        assert_eq!(evaluator.count_operations().function_calls, 2);
+        for point in [
+            [0.2, 0.3, 0.4, 0., 1.25],
+            [-0.7, 0.1, -0.2, 1., 1.25],
+            [0.6, -0.2, 0.1, 0., -0.5],
+            [0., 0., 0., 1., 7.],
+        ] {
+            let extra = if point[3] == 0. { point[4] } else { 0. };
+            let expected = 2. * (point[0] + point[1] + point[2]).sin() + extra;
+            assert!((evaluator.evaluate_single(&point) - expected).abs() < 1e-14);
+        }
+    }
+}
+
 #[cfg(test)]
 mod external_cache_tests {
     use crate::{
@@ -4647,3 +5430,662 @@ mod external_cache_tests {
         assert_eq!(calls.load(Ordering::SeqCst), 1);
     }
 }
+
+#[cfg(test)]
+mod evaluator_horner_candidate_tests {
+    use super::*;
+
+    #[test]
+    fn horner_orders_borrow_the_input_storage() {
+        let source = crate::parse!("a*f(x^4+x^3+x^2+x+1)+b*f(x^4+x^3+x^2+x+1)+c*sin(x)+d*sin(x)");
+        let view = source.as_view();
+        let storage = view.get_data();
+        let start = storage.as_ptr() as usize;
+        let end = start + storage.len();
+        let orders = AtomView::horner_orders(
+            &[view],
+            &HashSet::default(),
+            &OptimizationSettings::default(),
+        );
+        assert_eq!(orders[0].len(), 2);
+        for order in orders {
+            for atom in order {
+                let data = atom.get_data();
+                let address = data.as_ptr() as usize;
+                assert!(address >= start && address + data.len() <= end);
+            }
+        }
+    }
+
+    #[test]
+    fn unchanged_regions_and_repeated_arguments_remain_borrowed() {
+        let source = crate::parse!("sin(x+1)+cos(x+1)+x^2*y+z^5");
+        let tags = HashMap::default();
+        let settings = OptimizationSettings::default();
+        let mut regions = HornerRegions::new(
+            &tags,
+            &settings,
+            Arc::new(AtomicBool::new(false)),
+            [],
+            &[source.as_view()],
+        );
+        assert!(matches!(
+            regions.children(source.as_view()),
+            AtomOrView::View(view) if view.get_data().as_ptr() == source.as_view().get_data().as_ptr()
+        ));
+        let argument = crate::parse!("x+1");
+        assert!(matches!(
+            regions.optimized.get(&argument.as_view()),
+            Some(AtomOrView::View(_))
+        ));
+    }
+
+    #[test]
+    fn two_candidate_budget_compares_seeds_with_multiplication_priority() {
+        let visited = Mutex::new(Vec::new());
+        let original = [0, 1, 2];
+        let alternative = [2, 1, 0];
+        let result = search_horner_order_with_seed(
+            &original,
+            Some(&alternative),
+            &OptimizationSettings {
+                horner_iterations: 2,
+                ..Default::default()
+            },
+            |order| {
+                visited.lock().unwrap().push(order.to_vec());
+                if order == alternative {
+                    OperationCount::new(100, 1, 0, 0)
+                } else {
+                    OperationCount::new(0, 2, 0, 0)
+                }
+            },
+        );
+        assert_eq!(result, alternative);
+        assert_eq!(
+            visited.into_inner().unwrap(),
+            [original.to_vec(), alternative.to_vec()]
+        );
+    }
+
+    #[test]
+    fn duplicate_seed_does_not_consume_an_extra_candidate() {
+        let visited = Mutex::new(Vec::new());
+        let original = [0, 1, 2];
+        search_horner_order_with_seed(
+            &original,
+            Some(&original),
+            &OptimizationSettings {
+                horner_iterations: 2,
+                n_cores: 1,
+                ..Default::default()
+            },
+            |order| {
+                visited.lock().unwrap().push(order.to_vec());
+                OperationCount::default()
+            },
+        );
+        let visited = visited.into_inner().unwrap();
+        assert_eq!(visited.len(), 2);
+        assert_eq!(visited[0], original);
+        assert_ne!(visited[1], original);
+    }
+
+    #[test]
+    fn exhaustive_search_with_two_seeds_visits_each_order_once() {
+        let visited = Mutex::new(Vec::new());
+        let original = [3, 2, 1, 0];
+        let alternative = [0, 2, 3, 1];
+        let result = search_horner_order_with_seed(
+            &original,
+            Some(&alternative),
+            &OptimizationSettings {
+                horner_iterations: 24,
+                n_cores: 3,
+                ..Default::default()
+            },
+            |order| {
+                visited.lock().unwrap().push(order.to_vec());
+                OperationCount::default()
+            },
+        );
+        let mut visited = visited.into_inner().unwrap();
+        assert_eq!(visited.len(), 24);
+        visited.sort();
+        visited.dedup();
+        assert_eq!(visited.len(), 24);
+        assert!(visited.contains(&original.to_vec()));
+        assert!(visited.contains(&alternative.to_vec()));
+        assert_eq!(result, [0, 1, 2, 3]);
+    }
+
+    #[test]
+    fn equal_cost_seeds_choose_same_lexicographic_winner() {
+        let first = [2, 1, 0];
+        let second = [0, 2, 1];
+        for (original, alternative) in [(&first, &second), (&second, &first)] {
+            let result = search_horner_order_with_seed(
+                original,
+                Some(alternative),
+                &OptimizationSettings {
+                    horner_iterations: 2,
+                    ..Default::default()
+                },
+                |_| OperationCount::new(3, 5, 0, 0),
+            );
+            assert_eq!(result, second);
+        }
+    }
+
+    #[test]
+    fn discovery_counts_computed_values_but_not_builtin_constants() {
+        let source = crate::parse!("a*exp(x+y)+b*exp(x+y)+c*sin(x)+d*sin(x)+x^4+pi");
+        let exponential = crate::parse!("exp(x+y)");
+        let sine = crate::parse!("sin(x)");
+        let integer_power = crate::parse!("x^4");
+        let x = crate::parse!("x");
+        let y = crate::parse!("y");
+        let e = Atom::var(Symbol::E);
+        let pi = Atom::var(Symbol::PI);
+        let mut candidates = HashMap::default();
+        source
+            .as_view()
+            .count_horner_candidates(&HashSet::default(), &mut candidates);
+        assert_eq!(candidates.get(&exponential.as_view()), Some(&2));
+        assert_eq!(candidates.get(&sine.as_view()), Some(&2));
+        assert_eq!(candidates.get(&x.as_view()), Some(&1));
+        assert!(!candidates.contains_key(&y.as_view()));
+        assert!(!candidates.contains_key(&integer_power.as_view()));
+        assert!(!candidates.contains_key(&e.as_view()));
+        assert!(!candidates.contains_key(&pi.as_view()));
+    }
+
+    #[test]
+    fn discovery_keeps_explicit_parameters_opaque_even_for_builtin_constants() {
+        let parameter = crate::parse!("f(x+x*y)");
+        let source = crate::parse!("a*f(x+x*y)+b*f(x+x*y)");
+        let parameters = [parameter.as_view()].into_iter().collect();
+        let mut candidates = HashMap::default();
+        source
+            .as_view()
+            .count_horner_candidates(&parameters, &mut candidates);
+        assert_eq!(candidates.get(&parameter.as_view()), Some(&2));
+        assert!(!candidates.contains_key(&crate::parse!("x").as_view()));
+        assert!(!candidates.contains_key(&crate::parse!("y").as_view()));
+
+        let e = Atom::var(Symbol::E);
+        let pi = Atom::var(Symbol::PI);
+        let source = &e + &pi;
+        let parameters = [e.as_view(), pi.as_view()].into_iter().collect();
+        let mut candidates = HashMap::default();
+        source
+            .as_view()
+            .count_horner_candidates(&parameters, &mut candidates);
+        assert_eq!(candidates.get(&e.as_view()), Some(&1));
+        assert_eq!(candidates.get(&pi.as_view()), Some(&1));
+    }
+
+    #[test]
+    fn direct_builder_preserves_single_occurrence_function_parameter() {
+        let parameter = crate::parse!("f(x+x*y)");
+        let source = crate::parse!("(x+x*y)*f(x+x*y)");
+        let parameters = [parameter, crate::parse!("x"), crate::parse!("y")];
+        for iterations in [1, 10] {
+            let mut evaluator = source
+                .evaluator(&parameters)
+                .optimization_settings(OptimizationSettings {
+                    direct_translation: true,
+                    horner_iterations: iterations,
+                    ..Default::default()
+                })
+                .build()
+                .unwrap()
+                .map_coeff(&|coefficient| coefficient.re.to_f64());
+            for point in [[2., 0.5, -0.25], [3., -2., 0.25], [0., 7., 2.]] {
+                let expected = point[0] * (point[1] + point[1] * point[2]);
+                assert!((evaluator.evaluate_single(&point) - expected).abs() < 1e-14);
+            }
+        }
+    }
+
+    #[test]
+    fn direct_builder_factors_exponential_and_preserves_numerical_values() {
+        let source = crate::parse!("a*exp(x)+b*exp(x)+c");
+        let parameters = [
+            crate::parse!("x"),
+            crate::parse!("a"),
+            crate::parse!("b"),
+            crate::parse!("c"),
+        ];
+        let original = source
+            .evaluator(&parameters)
+            .optimization_settings(OptimizationSettings {
+                direct_translation: true,
+                horner_iterations: 0,
+                ..Default::default()
+            })
+            .build()
+            .unwrap();
+        assert_eq!(original.count_operations().multiplications, 2);
+        for iterations in [1, 2, 10] {
+            let mut evaluator = source
+                .evaluator(&parameters)
+                .optimization_settings(OptimizationSettings {
+                    direct_translation: true,
+                    horner_iterations: iterations,
+                    ..Default::default()
+                })
+                .build()
+                .unwrap()
+                .map_coeff(&|coefficient| coefficient.re.to_f64());
+            let operations = evaluator.count_operations();
+            assert_eq!(operations.multiplications, 1);
+            assert_eq!(operations.function_calls, 1);
+            for point in [
+                [-2., 2., -3., 0.7],
+                [0., -1., 4., 2.],
+                [0.5, 3., 5., -2.],
+                [2., 0.25, -0.5, -1.],
+            ] {
+                let expected = (point[1] + point[2]) * point[0].exp() + point[3];
+                assert!((evaluator.evaluate_single(&point) - expected).abs() < 1e-13);
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod local_function_parameter_scope_tests {
+    use super::*;
+
+    fn check_local_function_parameter(inherited_alias: bool) {
+        let source = crate::parse!("g(u,x,y)");
+        let root_parameters = [crate::parse!("u"), crate::parse!("x"), crate::parse!("y")];
+        let declared_parameters: Vec<Indeterminate> = [
+            crate::parse!("h(2*x+2*y)"),
+            crate::parse!("x"),
+            crate::parse!("y"),
+        ]
+        .into_iter()
+        .map(|argument| argument.try_into().unwrap())
+        .collect();
+        let body = crate::parse!("(x+x*y)*h(2*x+2*y)");
+        let alias = crate::parse!("local_parameter_scope_alias");
+        for policy in [InliningPolicy::Always, InliningPolicy::Never] {
+            for iterations in [0, 1, 10] {
+                let function_body = if inherited_alias { &alias } else { &body };
+                let mut builder = source
+                    .evaluator(&root_parameters)
+                    .add_function_with_options(
+                        crate::symbol!("g"),
+                        declared_parameters.clone(),
+                        function_body.clone(),
+                        FunctionRegistrationOptions::new().inlining(policy),
+                    )
+                    .unwrap()
+                    .optimization_settings(OptimizationSettings {
+                        direct_translation: true,
+                        horner_iterations: iterations,
+                        ..Default::default()
+                    });
+                if inherited_alias {
+                    builder = builder
+                        .add_aliases([(alias.clone(), body.clone())])
+                        .unwrap();
+                }
+                let mut evaluator = builder
+                    .build()
+                    .unwrap()
+                    .map_coeff(&|coefficient| coefficient.re.to_f64());
+                for point in [[2., 3., 4.], [-3., 0.5, -0.25], [0., 7., 2.]] {
+                    let expected = point[0] * (point[1] + point[1] * point[2]);
+                    assert!((evaluator.evaluate_single(&point) - expected).abs() < 1e-14);
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn inline_and_deferred_function_arguments_keep_their_atom_identity() {
+        check_local_function_parameter(false);
+    }
+
+    #[test]
+    fn inherited_aliases_preserve_inline_and_deferred_function_arguments() {
+        check_local_function_parameter(true);
+    }
+}
+
+#[cfg(test)]
+mod optimized_horner_search_tests {
+    use super::*;
+    use std::sync::atomic::AtomicUsize;
+
+    #[test]
+    fn mixed_mutations_preserve_permutations_and_include_larger_moves() {
+        let original: Vec<_> = (0..12).collect();
+        let mut rng = MonteCarloRng::new(0, 0);
+        let mut large_insert = false;
+        let mut large_reversal = false;
+        for mutation in 0..128 {
+            let mut order = original.clone();
+            mutate_horner_order(&mut order, &mut rng, mutation);
+            assert_ne!(order, original);
+            let changed = order.iter().zip(&original).filter(|(a, b)| a != b).count();
+            large_insert |= mutation % 4 == 2 && changed > 2;
+            large_reversal |= mutation % 4 == 3 && changed > 2;
+            order.sort_unstable();
+            assert_eq!(order, original);
+        }
+        assert!(large_insert && large_reversal);
+    }
+
+    #[test]
+    fn multiple_seeds_are_unique_and_charged_against_the_budget() {
+        for budget in [2, 3] {
+            let visited = Mutex::new(Vec::new());
+            let initial = [0, 1, 2, 3];
+            let second = [1, 0, 2, 3];
+            let third = [3, 2, 1, 0];
+            let selected = search_horner_order_with_seeds(
+                &initial,
+                &[&initial, &second, &second, &third],
+                &OptimizationSettings {
+                    horner_iterations: budget,
+                    ..Default::default()
+                },
+                |order| {
+                    visited.lock().unwrap().push(order.to_vec());
+                    OperationCount::new(0, if order == third { 1 } else { 2 }, 0, 0)
+                },
+            );
+            let visited = visited.into_inner().unwrap();
+            assert_eq!(visited.len(), budget);
+            assert_eq!(visited[0], initial);
+            assert_eq!(visited[1], second);
+            if budget == 3 {
+                assert_eq!(visited[2], third);
+                assert_eq!(selected, third);
+            }
+        }
+    }
+
+    #[test]
+    fn long_search_restarts_without_exceeding_budget_or_losing_the_winner() {
+        let initial: Vec<_> = (0..10).collect();
+        let visited = Mutex::new(Vec::new());
+        let selected = search_horner_order(
+            &initial,
+            &OptimizationSettings {
+                horner_iterations: 130,
+                n_cores: 1,
+                ..Default::default()
+            },
+            |order| {
+                visited.lock().unwrap().push(order.to_vec());
+                OperationCount::new(0, usize::from(order != initial), 0, 0)
+            },
+        );
+        let visited = visited.into_inner().unwrap();
+        assert_eq!(visited.len(), 130);
+        assert_eq!(
+            visited[34],
+            initial.iter().rev().copied().collect::<Vec<_>>()
+        );
+        assert_eq!(selected, initial);
+    }
+
+    #[test]
+    fn finished_program_score_counts_binary_powering() {
+        // A lowered instruction fixture keeps powers explicit: the Atom
+        // linearizer expands positive integer powers into multiplication lists.
+        let mut evaluator = ExpressionEvaluator {
+            stack: vec![0_f64; 5],
+            param_count: 2,
+            reserved_indices: 2,
+            instructions: vec![
+                (Instr::Pow(2, 0, 7), ComplexPhase::default()),
+                (Instr::Pow(3, 1, -4), ComplexPhase::default()),
+                (Instr::Add(4, vec![2, 3]), ComplexPhase::default()),
+            ],
+            result_indices: vec![4],
+            external_fns: vec![],
+            settings: OptimizationSettings::default(),
+        };
+        let cost = evaluator_horner_cost(&evaluator);
+        assert_eq!(cost.multiplications, 6);
+        assert_eq!(cost.additions, 1);
+        assert_eq!(cost.inversions, 1);
+        assert_eq!(evaluator.evaluate_single(&[2., 2.]), 128.0625);
+    }
+
+    #[test]
+    fn search_selects_fewer_multiplications_after_pair_elimination() {
+        let parameters = [
+            crate::parse!("horner_cost_a"),
+            crate::parse!("horner_cost_b"),
+            crate::parse!("horner_cost_c"),
+        ];
+        let permutations = unique_permutations(&[0_usize, 1, 2]).1;
+        let settings = OptimizationSettings {
+            horner_iterations: 6,
+            horner_use_cpe: true,
+            ..Default::default()
+        };
+        let mut strict_misranking = false;
+        // Coefficient sign normalization depends on the canonical variable
+        // order. Check every relabeling of the same small polynomial instead
+        // of relying on symbols registered by another parallel test.
+        for relabeling in &permutations {
+            let [a, b, c] =
+                std::array::from_fn::<_, 3, _>(|i| parameters[relabeling[i]].to_canonical_string());
+            let source = crate::parse!(&format!(
+                "{b}*{c}^3-3*{b}^2-{a}*{c}^2+2*{a}^2*{c}-2*{a}^2*{b}^2-3*{a}^3*{c}+2*{a}^3*{c}^2+3*{a}^3*{b}^2*{c}-2*{a}^3*{b}^2*{c}^2"
+            ));
+            let expressions = [source.as_view()];
+            let mut outcomes = Vec::new();
+            for permutation in &permutations {
+                let order = permutation
+                    .iter()
+                    .map(|&i| parameters[i].clone())
+                    .collect::<Vec<_>>();
+                let before = horner_cost(AtomView::horner_scheme_cost(&expressions, &order));
+                let evaluator = AtomView::build_horner_evaluator(
+                    &expressions,
+                    &FunctionMap::new(),
+                    &parameters,
+                    &order,
+                    settings.clone(),
+                )
+                .unwrap();
+                outcomes.push((
+                    before,
+                    horner_cost(evaluator_horner_cost(&evaluator)),
+                    order,
+                ));
+            }
+            let cheap_winner = outcomes
+                .iter()
+                .min_by(|left, right| left.0.cmp(&right.0).then_with(|| left.2.cmp(&right.2)))
+                .unwrap();
+            let finished_minimum = outcomes.iter().map(|outcome| outcome.1).min().unwrap();
+            strict_misranking |= cheap_winner.1.0 > finished_minimum.0;
+            for horner_use_cpe in [false, true] {
+                let evaluator = source
+                    .evaluator(&parameters)
+                    .optimization_settings(settings.clone().horner_use_cpe(horner_use_cpe))
+                    .build()
+                    .unwrap();
+                assert_eq!(
+                    horner_cost(evaluator_horner_cost(&evaluator)),
+                    if horner_use_cpe {
+                        finished_minimum
+                    } else {
+                        cheap_winner.1
+                    }
+                );
+                let mut evaluator = evaluator.map_coeff(&|coefficient| coefficient.re.to_f64());
+                for point in [[1.2_f64, -0.3, 0.7], [-2., 1.5, -0.2], [0., 3., 2.]] {
+                    let [a, b, c] = std::array::from_fn::<_, 3, _>(|i| point[relabeling[i]]);
+                    let expected = b * c.powi(3) - 3. * b.powi(2) - a * c.powi(2)
+                        + 2. * a.powi(2) * c
+                        - 2. * a.powi(2) * b.powi(2)
+                        - 3. * a.powi(3) * c
+                        + 2. * a.powi(3) * c.powi(2)
+                        + 3. * a.powi(3) * b.powi(2) * c
+                        - 2. * a.powi(3) * b.powi(2) * c.powi(2);
+                    assert!((evaluator.evaluate_single(&point) - expected).abs() < 1e-12);
+                }
+            }
+        }
+        assert!(
+            strict_misranking,
+            "the fixture must expose a strictly worse pre-CPE winner"
+        );
+    }
+
+    #[test]
+    fn legacy_function_search_scores_all_calls_and_reuses_the_winner() {
+        let parameters = [
+            crate::parse!("legacy_score_x"),
+            crate::parse!("legacy_score_y"),
+            crate::parse!("legacy_score_z"),
+        ];
+        let body = crate::parse!(
+            "legacy_score_b*legacy_score_c^3-3*legacy_score_b^2\
+            -legacy_score_a*legacy_score_c^2+2*legacy_score_a^2*legacy_score_c\
+            -2*legacy_score_a^2*legacy_score_b^2-3*legacy_score_a^3*legacy_score_c\
+            +2*legacy_score_a^3*legacy_score_c^2\
+            +3*legacy_score_a^3*legacy_score_b^2*legacy_score_c\
+            -2*legacy_score_a^3*legacy_score_b^2*legacy_score_c^2"
+        );
+        let source = crate::parse!(
+            "legacy_score_f(legacy_score_x,legacy_score_y,legacy_score_z)\
+            +legacy_score_f(legacy_score_y,legacy_score_x,legacy_score_z)"
+        );
+        let mut functions = FunctionMap::new();
+        functions
+            .add_function(
+                crate::symbol!("legacy_score_f"),
+                vec![
+                    crate::symbol!("legacy_score_a"),
+                    crate::symbol!("legacy_score_b"),
+                    crate::symbol!("legacy_score_c"),
+                ],
+                body,
+            )
+            .unwrap();
+        let tree = source
+            .as_view()
+            .to_evaluation_tree(&functions, &parameters)
+            .unwrap();
+        assert_eq!(tree.functions.len(), 1);
+        let mut variables = HashMap::default();
+        for expression in &tree.functions[0].2.tree {
+            expression.find_all_variables(&mut variables);
+        }
+        let mut variables = variables
+            .into_iter()
+            .filter(|(_, count)| *count > 1)
+            .map(|(variable, _)| variable)
+            .collect::<Vec<_>>();
+        variables.sort();
+        assert_eq!(variables.len(), 3);
+        let settings = OptimizationSettings {
+            direct_translation: false,
+            horner_use_cpe: true,
+            horner_iterations: 6,
+            ..Default::default()
+        };
+        let minimum = unique_permutations(&variables)
+            .1
+            .into_iter()
+            .map(|order| {
+                let mut candidate = tree.clone();
+                candidate.apply_horner_scheme_to_scope(1, &order);
+                candidate.common_subexpression_elimination();
+                horner_cost(evaluator_horner_cost(&candidate.linearize(&settings)))
+            })
+            .min()
+            .unwrap();
+        let mut selected_tree = tree.clone();
+        let (_, cached) = selected_tree.optimize_horner_scheme_with_cpe(&settings);
+        assert_eq!(horner_cost(evaluator_horner_cost(&cached)), minimum);
+        // The tree returned by the public scheme API must agree with the cached
+        // full program even though only the latter has undergone CSE and CPE.
+        selected_tree.common_subexpression_elimination();
+        assert_eq!(
+            horner_cost(evaluator_horner_cost(&selected_tree.linearize(&settings))),
+            minimum
+        );
+        let evaluator = tree.optimize(&settings);
+        assert_eq!(horner_cost(evaluator_horner_cost(&evaluator)), minimum);
+        let mut evaluator = evaluator.map_coeff(&|coefficient| coefficient.re.to_f64());
+        let expected_body = |a: f64, b: f64, c: f64| {
+            b * c.powi(3) - 3. * b.powi(2) - a * c.powi(2) + 2. * a.powi(2) * c
+                - 2. * a.powi(2) * b.powi(2)
+                - 3. * a.powi(3) * c
+                + 2. * a.powi(3) * c.powi(2)
+                + 3. * a.powi(3) * b.powi(2) * c
+                - 2. * a.powi(3) * b.powi(2) * c.powi(2)
+        };
+        for point in [[1.2, -0.3, 0.7], [-2., 1.5, -0.2], [0., 3., 2.]] {
+            let [a, b, c] = point;
+            let expected = expected_body(a, b, c) + expected_body(b, a, c);
+            assert!((evaluator.evaluate_single(&point) - expected).abs() < 1e-11);
+        }
+    }
+
+    #[test]
+    fn candidate_build_errors_are_returned() {
+        let source = crate::parse!("unregistered_horner_function(x+y)+x*y");
+        for direct_translation in [false, true] {
+            for horner_use_cpe in [false, true] {
+                let result = source
+                    .evaluator(&[crate::parse!("x"), crate::parse!("y")])
+                    .optimization_settings(OptimizationSettings {
+                        horner_iterations: 1000,
+                        horner_use_cpe,
+                        direct_translation,
+                        n_cores: 3,
+                        ..Default::default()
+                    })
+                    .build();
+                assert!(matches!(
+                    result,
+                    Err(EvaluationError::UndefinedFunction { .. })
+                ));
+            }
+        }
+    }
+
+    #[test]
+    fn cancellation_is_latched_and_returns_a_working_candidate() {
+        let source = crate::parse!("x*y+x*z+y*z+x*y*z");
+        for direct_translation in [false, true] {
+            for horner_use_cpe in [false, true] {
+                let polls = Arc::new(AtomicUsize::new(0));
+                let callback_polls = polls.clone();
+                let mut evaluator = source
+                    .evaluator(&[crate::parse!("x"), crate::parse!("y"), crate::parse!("z")])
+                    .optimization_settings(OptimizationSettings {
+                        horner_iterations: 1000,
+                        horner_use_cpe,
+                        direct_translation,
+                        n_cores: 3,
+                        // A one-shot callback tests the latch rather than relying on
+                        // a caller that keeps returning true after cancellation.
+                        abort_check: Some(Box::new(move || {
+                            callback_polls.fetch_add(1, Ordering::Relaxed) == 0
+                        })),
+                        ..Default::default()
+                    })
+                    .build()
+                    .unwrap()
+                    .map_coeff(&|coefficient| coefficient.re.to_f64());
+                assert_eq!(polls.load(Ordering::Relaxed), 1);
+                assert_eq!(evaluator.evaluate_single(&[2., 3., 4.]), 50.);
+            }
+        }
+    }
+}
diff --git a/src/id.rs b/src/id.rs
index 373173b3..09e20521 100644
--- a/src/id.rs
+++ b/src/id.rs
@@ -7100,7 +7100,7 @@ mod test {
     fn nested() {
         let res = parse!("f(x+x*y+f(x+x^2),x,f(x+x^2))").replace_map_bottom_up(|a, c, o| {
             if c.parent_type == Some(AtomType::Fun) {
-                let r = a.horner_scheme(None, false, false);
+                let r = a.horner_scheme_local::<crate::atom::Indeterminate>(None, false);
                 if r.as_view() != a {
                     **o = r;
                 }
diff --git a/symbolica.pyi b/symbolica.pyi
index f8c92721..b2b1628a 100644
--- a/symbolica.pyi
+++ b/symbolica.pyi
@@ -5457,6 +5457,7 @@ class Expression:
         max_horner_scheme_variables: int = 500,
         max_common_pair_cache_entries: int = 1000000,
         max_common_pair_distance: int = 100,
+        horner_use_cpe: bool = False,
     ) -> Evaluator:
         """
         Create an evaluator that can evaluate (nested) expressions in an optimized fashion.
@@ -5498,7 +5499,7 @@ class Expression:
         iterations: int, optional
             The number of Horner schemes to try.
         cpe_iterations: int | None, optional
-            The number of CPE iterations to perform. The number if unbounded if `None`.
+            The maximum number of common-pair elimination iterations. If `None`, continue until no more pairs are eliminated.
         n_cores: int, optional
             The number of CPU cores used for the optimization.
         verbose: bool, optional
@@ -5520,6 +5521,10 @@ class Expression:
             The maximum number of entries in the common pair cache.
         max_common_pair_distance: int, optional
             The maximum distance between common pairs. Used when clearing cache entries.
+        horner_use_cpe: bool, optional
+            If `False` (the default), rank Horner candidates using common-subexpression elimination only.
+            If `True`, also run common-pair elimination for each candidate; this can improve the selected scheme but is more expensive.
+            The selected scheme receives final common-pair elimination in either mode, controlled by `cpe_iterations`.
         """
 
     @classmethod
@@ -5540,6 +5545,7 @@ class Expression:
         max_horner_scheme_variables: int = 500,
         max_common_pair_cache_entries: int = 1000000,
         max_common_pair_distance: int = 100,
+        horner_use_cpe: bool = False,
     ) -> Evaluator:
         """
         Create an evaluator that can jointly evaluate (nested) expressions in an optimized fashion.
@@ -5564,9 +5570,9 @@ class Expression:
         functions: Sequence[FunctionDefinition]
             Function definitions available to the evaluator.
         iterations: int, optional
-            The number of optimization passes to run.
+            The number of Horner schemes to try.
         cpe_iterations: int | None, optional
-            The number of common subexpression elimination iterations to perform.
+            The maximum number of common-pair elimination iterations. If `None`, continue until no more pairs are eliminated.
         n_cores: int, optional
             The number of CPU cores used for parallel optimization.
         verbose: bool, optional
@@ -5587,6 +5593,10 @@ class Expression:
             The maximum number of common-subexpression pairs to cache.
         max_common_pair_distance: int, optional
             The maximum distance between factors when searching for common pairs.
+        horner_use_cpe: bool, optional
+            If `False` (the default), rank Horner candidates using common-subexpression elimination only.
+            If `True`, also run common-pair elimination for each candidate; this can improve the selected scheme but is more expensive.
+            The selected scheme receives final common-pair elimination in either mode, controlled by `cpe_iterations`.
         """
 
     def canonize_tensors(
diff --git a/tests/horner_function_parameters.rs b/tests/horner_function_parameters.rs
new file mode 100644
index 00000000..3423aeb8
--- /dev/null
+++ b/tests/horner_function_parameters.rs
@@ -0,0 +1,241 @@
+use symbolica::evaluate::{
+    FunctionRegistrationOptions, InliningPolicy, Instruction, OperationCount, OptimizationSettings,
+};
+use symbolica::prelude::*;
+
+fn configurations() -> impl Iterator<Item = (bool, bool, usize)> {
+    [false, true].into_iter().flat_map(|direct| {
+        [false, true].into_iter().flat_map(move |use_cpe| {
+            [0, 1, 10].map(move |iterations| (direct, use_cpe, iterations))
+        })
+    })
+}
+
+fn settings(direct: bool, use_cpe: bool, iterations: usize) -> OptimizationSettings {
+    OptimizationSettings::new()
+        .direct_translation(direct)
+        .horner_use_cpe(use_cpe)
+        .horner_iterations(iterations)
+        .cores(1)
+}
+
+#[test]
+fn single_function_inputs_survive_factoring_and_coefficient_cleanup() {
+    let parameters = [
+        parse!("horner_inputs::f(x+x*y)"),
+        parse!("horner_inputs::h(2*x+2*y)"),
+        parse!("x"),
+        parse!("y"),
+    ];
+    // Each function occurs only once in its output. Its argument must remain
+    // the input's exact identity even when the surrounding algebra is factored.
+    let outputs = [
+        parse!("(x+x*y)*horner_inputs::f(x+x*y)"),
+        parse!("(6*x+10*y)*horner_inputs::h(2*x+2*y)"),
+    ];
+    for (direct, use_cpe, iterations) in configurations() {
+        let mut evaluator = Atom::evaluator_multiple(&outputs, &parameters)
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        assert_eq!(evaluator.count_operations().function_calls, 0);
+        for [f, h, x, y] in [[2., 7., 3., 4.], [-3., 0.5, -2., 0.25], [0., 2., 7., 0.]] {
+            let mut actual = [0.; 2];
+            evaluator.evaluate(&[f, h, x, y], &mut actual);
+            assert_eq!(actual, [f * (x + x * y), h * (6. * x + 10. * y)]);
+        }
+    }
+}
+
+#[test]
+fn nested_function_inputs_remain_independent_in_either_parameter_order() {
+    let f = parse!("horner_nested_inputs::f(x+x*y)");
+    let g = parse!("horner_nested_inputs::g(horner_nested_inputs::f(x+x*y))");
+    let source = parse!(
+        "(x+x*y)*horner_nested_inputs::f(x+x*y)\
+         +(x+x*y)*horner_nested_inputs::g(horner_nested_inputs::f(x+x*y))"
+    );
+    for reverse in [false, true] {
+        let parameters = if reverse {
+            [g.clone(), f.clone(), parse!("x"), parse!("y")]
+        } else {
+            [f.clone(), g.clone(), parse!("x"), parse!("y")]
+        };
+        for (direct, use_cpe, iterations) in configurations() {
+            let mut evaluator = source
+                .evaluator(&parameters)
+                .optimization_settings(settings(direct, use_cpe, iterations))
+                .build()
+                .unwrap()
+                .map_coeff(&|c| c.re.to_f64());
+            assert_eq!(evaluator.count_operations().function_calls, 0);
+            for [f, g, x, y] in [
+                [2., 7., 3., 4.],
+                [-3., 0.5, -2., 0.25],
+                [0., 2., 0.25, -0.5],
+            ] {
+                let point = if reverse { [g, f, x, y] } else { [f, g, x, y] };
+                assert_eq!(evaluator.evaluate_single(&point), (f + g) * (x + x * y));
+            }
+        }
+    }
+}
+
+#[test]
+fn normalization_hooks_cannot_observe_temporary_parameter_symbols() {
+    let guard = symbol!(
+        "horner_guard_inputs::guard",
+        norm = |view, out| {
+            if let AtomView::Fun(function) = view
+                && matches!(function.iter().next(), Some(AtomView::Var(_)))
+            {
+                out.to_num(0);
+            }
+        }
+    );
+    assert!(
+        guard
+            .call((parse!("horner_guard_inputs::probe"),))
+            .is_zero()
+    );
+    let parameter = parse!("horner_guard_inputs::p(x+x*y)");
+    let guarded = guard.call((&parameter,));
+    let nested_guarded = guard.call((&guarded,));
+    let factor = parse!("x+x*x");
+    let outputs = [&guarded * &factor, &nested_guarded * &factor];
+    let parameters = [parameter, parse!("x")];
+    for (direct, use_cpe, iterations) in configurations() {
+        let mut evaluator = Atom::evaluator_multiple(&outputs, &parameters)
+            .add_function(
+                guard,
+                vec![symbol!("horner_guard_inputs::z")],
+                parse!("horner_guard_inputs::z+1"),
+            )
+            .unwrap()
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        // Both enclosing custom-normalized calls must keep their original
+        // function arguments until lowering; a placeholder Var would erase them.
+        for ([p, x], expected) in [
+            ([3., 2.], [24., 30.]),
+            ([-2., 0.5], [-0.75, 0.]),
+            ([7., -2.], [16., 18.]),
+        ] {
+            let mut actual = [0.; 2];
+            evaluator.evaluate(&[p, x], &mut actual);
+            assert_eq!(actual, expected);
+        }
+    }
+}
+
+#[test]
+fn compound_local_arguments_and_aliases_keep_their_binding_scopes() {
+    let parameter = parse!("horner_local_inputs::h(2*x+2*y)");
+    let parameters = [parse!("u"), parse!("x"), parse!("y"), parameter.clone()];
+    let declared: Vec<Indeterminate> = [parameter, parse!("x"), parse!("y")]
+        .into_iter()
+        .map(|p| p.try_into().unwrap())
+        .collect();
+    let alias = parse!("horner_local_inputs::value");
+    let source = parse!("horner_local_inputs::outer(u,x+1,y+2)+horner_local_inputs::value");
+    for policy in [InliningPolicy::Always, InliningPolicy::Never] {
+        for (direct, use_cpe, iterations) in configurations() {
+            let mut evaluator = source
+                .evaluator(&parameters)
+                .add_aliases([(
+                    alias.clone(),
+                    parse!("(x+x*y)*horner_local_inputs::h(2*x+2*y)"),
+                )])
+                .unwrap()
+                .add_function_with_options(
+                    symbol!("horner_local_inputs::inner"),
+                    Vec::<Symbol>::new(),
+                    alias.clone(),
+                    FunctionRegistrationOptions::new().inlining(InliningPolicy::Never),
+                )
+                .unwrap()
+                .add_function_with_options(
+                    symbol!("horner_local_inputs::outer"),
+                    declared.clone(),
+                    parse!("horner_local_inputs::value+horner_local_inputs::inner()"),
+                    FunctionRegistrationOptions::new().inlining(policy),
+                )
+                .unwrap()
+                .optimization_settings(settings(direct, use_cpe, iterations))
+                .build()
+                .unwrap()
+                .map_coeff(&|c| c.re.to_f64());
+            // The alias in outer inherits its arguments; the alias in inner
+            // and at the root sees the global inputs instead.
+            for [u, x, y, h] in [[2., 3., 4., 11.], [-3., 0.5, -0.25, 7.]] {
+                let expected = u * ((x + 1.) + (x + 1.) * (y + 2.)) + 2. * h * (x + x * y);
+                assert_eq!(evaluator.evaluate_single(&[u, x, y, h]), expected);
+            }
+        }
+    }
+}
+
+#[test]
+fn temporary_parameter_names_do_not_capture_existing_symbols() {
+    let source = parse!(
+        "(x+x*y)*horner_collision::f(x+x*y)\
+         +symbolica::horner_parameter_0*horner_collision::f(x+x*y)\
+         +symbolica::horner_parameter_0"
+    );
+    let parameters = [
+        parse!("horner_collision::f(x+x*y)"),
+        parse!("x"),
+        parse!("y"),
+        parse!("symbolica::horner_parameter_0"),
+    ];
+    for (direct, use_cpe, iterations) in configurations() {
+        let mut evaluator = source
+            .evaluator(&parameters)
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        for [f, x, y, existing] in [[2., 3., 4., 11.], [-3., 0.5, -0.25, 7.]] {
+            assert_eq!(
+                evaluator.evaluate_single(&[f, x, y, existing]),
+                f * (x + x * y + existing) + existing
+            );
+        }
+    }
+}
+
+#[test]
+fn ordinary_symbol_inputs_retain_exported_arithmetic_quality() {
+    let source = parse!("a*x+b*x+c");
+    let parameters = [parse!("a"), parse!("b"), parse!("c"), parse!("x")];
+    for (direct, use_cpe, iterations) in configurations() {
+        let mut evaluator = source
+            .evaluator(&parameters)
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        let expected = OperationCount::new(2, if iterations == 0 { 2 } else { 1 }, 0, 0);
+        assert_eq!(evaluator.count_operations(), expected);
+        let exported = evaluator.export_instructions();
+        assert_eq!((exported.input_count, exported.output_count), (4, 1));
+        assert!(exported.sub_evaluators.is_empty());
+        assert!(exported.constant_functions.is_empty());
+        let mut exported_cost = OperationCount::default();
+        for instruction in &exported.instructions {
+            match instruction {
+                Instruction::Add(_, inputs, _) => exported_cost.additions += inputs.len() - 1,
+                Instruction::Mul(_, inputs, _) => exported_cost.multiplications += inputs.len() - 1,
+                Instruction::Assign(..) => {}
+                other => panic!("unexpected instruction in polynomial evaluator: {other:?}"),
+            }
+        }
+        assert_eq!(exported_cost, expected);
+        for [a, b, c, x] in [[2., 7., 3., 4.], [-3., 0.5, -2., 0.25]] {
+            assert_eq!(evaluator.evaluate_single(&[a, b, c, x]), (a + b) * x + c);
+        }
+    }
+}
diff --git a/tests/horner_regions.rs b/tests/horner_regions.rs
new file mode 100644
index 00000000..b67446a3
--- /dev/null
+++ b/tests/horner_regions.rs
@@ -0,0 +1,298 @@
+use std::sync::{
+    Arc,
+    atomic::{AtomicUsize, Ordering},
+};
+
+use symbolica::evaluate::{
+    FunctionRegistrationOptions, InliningPolicy, Instruction, OperationCount, OptimizationSettings,
+};
+use symbolica::prelude::*;
+
+fn configurations() -> impl Iterator<Item = (bool, bool, usize)> {
+    [false, true].into_iter().flat_map(|direct| {
+        [false, true].into_iter().flat_map(move |use_cpe| {
+            [0, 1, 10].map(move |iterations| (direct, use_cpe, iterations))
+        })
+    })
+}
+
+fn settings(direct: bool, use_cpe: bool, iterations: usize) -> OptimizationSettings {
+    OptimizationSettings::new()
+        .direct_translation(direct)
+        .horner_use_cpe(use_cpe)
+        .horner_iterations(iterations)
+        .cores(1)
+}
+
+fn close(actual: f64, expected: f64) {
+    assert!(
+        (actual - expected).abs() <= 1e-12 * (1. + expected.abs()),
+        "{actual} != {expected}"
+    );
+}
+
+#[test]
+fn repeated_and_nested_function_arguments_are_factored_and_shared() {
+    let parameters = [parse!("a"), parse!("b"), parse!("x")];
+    let outputs = [
+        parse!("sin(a*x+b*x)"),
+        parse!("cos(a*x+b*x)"),
+        parse!("sin(sin(a*x+b*x)+cos(a*x+b*x))"),
+    ];
+    for (direct, use_cpe, iterations) in configurations() {
+        let mut evaluator = Atom::evaluator_multiple(&outputs, &parameters)
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        // The common argument must use one multiplication after Horner and
+        // appear only once even though it is used at different nesting depths.
+        if direct {
+            assert_eq!(
+                evaluator.count_operations(),
+                OperationCount::new(2, if iterations == 0 { 2 } else { 1 }, 0, 3),
+                "direct={direct}, use_cpe={use_cpe}, iterations={iterations}"
+            );
+        }
+        for [a, b, x] in [[2., 3., 0.25], [-3., 0.5, -2.], [7., -7., 3.]] {
+            let argument: f64 = (a + b) * x;
+            let expected = [
+                argument.sin(),
+                argument.cos(),
+                (argument.sin() + argument.cos()).sin(),
+            ];
+            let mut actual = [0.; 3];
+            evaluator.evaluate(&[a, b, x], &mut actual);
+            for (actual, expected) in actual.into_iter().zip(expected) {
+                close(actual, expected);
+            }
+        }
+    }
+}
+
+#[test]
+fn tagged_calls_keep_symbolic_tags_and_optimize_only_runtime_arguments() {
+    let tag = parse!("horner_tag::a+horner_tag::a*horner_tag::b");
+    let function = symbol!("horner_tag::evaluate");
+    let source = function.call((&tag, parse!("a*x+b*x")));
+    let parameters = [parse!("a"), parse!("b"), parse!("x")];
+    let mut functions = FunctionMap::new();
+    functions
+        .add_tagged_function_with_options(
+            function,
+            vec![tag.clone()],
+            vec![symbol!("horner_tag::z")],
+            parse!("horner_tag::z^2+1"),
+            FunctionRegistrationOptions::new().inlining(InliningPolicy::Never),
+        )
+        .unwrap();
+    for (direct, use_cpe, iterations) in configurations() {
+        let mut evaluator = source
+            .evaluator(&parameters)
+            .function_map(functions.clone())
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        // The tag's symbols are intentionally absent from the input list:
+        // they identify the definition and must never become runtime inputs.
+        assert_eq!(
+            evaluator.count_operations(),
+            OperationCount::new(1, if iterations == 0 { 2 } else { 1 }, 0, 1),
+            "direct={direct}, use_cpe={use_cpe}, iterations={iterations}"
+        );
+        let exported = evaluator.export_instructions();
+        let body = exported
+            .sub_evaluators
+            .iter()
+            .find(|body| body.symbol == function)
+            .unwrap();
+        assert_eq!(body.tags, [tag.to_canonical_string()]);
+        for [a, b, x] in [[2., 3., 0.25], [-3., 0.5, -2.], [7., -7., 3.]] {
+            let argument = (a + b) * x;
+            assert_eq!(
+                evaluator.evaluate_single(&[a, b, x]),
+                argument * argument + 1.
+            );
+        }
+    }
+}
+
+#[test]
+fn power_children_are_optimized_without_losing_integer_power_factoring() {
+    let parameters = [
+        parse!("a"),
+        parse!("b"),
+        parse!("x"),
+        parse!("c"),
+        parse!("d"),
+        parse!("y"),
+    ];
+    for (direct, use_cpe, iterations) in configurations() {
+        let mut evaluator = parse!("(a*x+b*x)^(c*y+d*y)")
+            .evaluator(&parameters)
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        if direct {
+            assert_eq!(
+                evaluator.count_operations(),
+                OperationCount::new(2, if iterations == 0 { 4 } else { 2 }, 0, 1),
+                "direct={direct}, use_cpe={use_cpe}, iterations={iterations}"
+            );
+        }
+        for [a, b, x, c, d, y] in [[2., 3., 0.5, 1., 2., 0.25], [1., 0.5, 2., -3., 1., 0.5]] {
+            let base: f64 = (a + b) * x;
+            let exponent = (c + d) * y;
+            close(
+                evaluator.evaluate_single(&[a, b, x, c, d, y]),
+                base.powf(exponent),
+            );
+        }
+
+        let mut squared = parse!("(a*x+b*x)^2")
+            .evaluator(&parameters[..3])
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        if direct {
+            // Atom normalization distributes an integer power over the
+            // factored product: (x*(a+b))^2 becomes x^2*(a+b)^2.
+            assert_eq!(
+                squared.count_operations(),
+                OperationCount::new(1, 3, 0, 0),
+                "direct={direct}, use_cpe={use_cpe}, iterations={iterations}"
+            );
+        }
+        assert_eq!(squared.evaluate_single(&[2., 3., 0.5]), 6.25);
+
+        let mut affine_squared = parse!("(a*x+b*x+c)^2")
+            .evaluator(&parameters[..4])
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        if direct {
+            assert_eq!(
+                affine_squared.count_operations(),
+                OperationCount::new(2, if iterations == 0 { 3 } else { 2 }, 0, 0),
+                "direct={direct}, use_cpe={use_cpe}, iterations={iterations}"
+            );
+        }
+        assert_eq!(affine_squared.evaluate_single(&[2., 3., 0.5, 1.]), 12.25);
+
+        let mut polynomial = parse!("x^2+x")
+            .evaluator(&[parse!("x")])
+            .optimization_settings(settings(direct, use_cpe, iterations))
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        assert_eq!(
+            polynomial.count_operations(),
+            OperationCount::new(1, 1, 0, 0),
+            "direct={direct}, use_cpe={use_cpe}, iterations={iterations}"
+        );
+        if direct && iterations > 0 {
+            let exported = polynomial.export_instructions();
+            assert!(
+                exported.instructions.iter().any(|instruction| {
+                    matches!(instruction, Instruction::Mul(_, inputs, _)
+                        if inputs.len() == 2 && inputs[0] != inputs[1])
+                }),
+                "x^2+x should become x*(x+1): {:?}",
+                exported.instructions
+            );
+        }
+        for x in [-2., 0., 0.5, 3.] {
+            assert_eq!(polynomial.evaluate_single(&[x]), x * (x + 1.));
+        }
+    }
+}
+
+#[test]
+fn rewritten_function_calls_can_match_a_supplied_input() {
+    let source = parse!("sin(x+x*y)");
+    let parameters = [parse!("sin(x*(1+y))"), parse!("x"), parse!("y")];
+    for use_cpe in [false, true] {
+        for iterations in [0, 1, 10] {
+            let mut evaluator = source
+                .evaluator(&parameters)
+                .optimization_settings(settings(true, use_cpe, iterations))
+                .build()
+                .unwrap()
+                .map_coeff(&|c| c.re.to_f64());
+            let actual = evaluator.evaluate_single(&[99., 1., 1.]);
+            if iterations == 0 {
+                close(actual, 2_f64.sin());
+            } else {
+                // Rebuilding sin after optimizing its argument makes it
+                // exactly equal to the supplied function-shaped input.
+                assert_eq!(actual, 99., "use_cpe={use_cpe}, iterations={iterations}");
+            }
+        }
+    }
+}
+
+#[test]
+fn nested_search_cancellation_is_latched_across_later_regions() {
+    let source = parse!("sin(x*y+x*z+y*z+x*y*z)+cos(x*y+x*z+y*z+2*x*y*z)");
+    let parameters = [parse!("x"), parse!("y"), parse!("z")];
+    for use_cpe in [false, true] {
+        let polls = Arc::new(AtomicUsize::new(0));
+        let callback_polls = polls.clone();
+        let settings = settings(true, use_cpe, 10).abort_check(Some(Box::new(move || {
+            callback_polls.fetch_add(1, Ordering::Relaxed) == 0
+        })));
+        let mut evaluator = source
+            .evaluator(&parameters)
+            .optimization_settings(settings)
+            .build()
+            .unwrap()
+            .map_coeff(&|c| c.re.to_f64());
+        assert_eq!(polls.load(Ordering::Relaxed), 1, "use_cpe={use_cpe}");
+        for [x, y, z] in [[2., 3., 4.], [-2., 0.5, 0.25]] {
+            let common: f64 = x * y + x * z + y * z;
+            let product = x * y * z;
+            close(
+                evaluator.evaluate_single(&[x, y, z]),
+                (common + product).sin() + (common + 2. * product).cos(),
+            );
+        }
+    }
+}
+
+#[test]
+fn shared_alias_arguments_keep_root_and_local_bindings_distinct() {
+    let parameters = [parse!("a"), parse!("b"), parse!("x"), parse!("u")];
+    let alias = parse!("horner_region_scope::value");
+    let source = parse!("horner_region_scope::value+horner_region_scope::outer(u)");
+    for policy in [InliningPolicy::Always, InliningPolicy::Never] {
+        for (direct, use_cpe, iterations) in configurations() {
+            let mut evaluator = source
+                .evaluator(&parameters)
+                .add_aliases([(alias.clone(), parse!("sin(a*x+b*x)"))])
+                .unwrap()
+                .add_function_with_options(
+                    symbol!("horner_region_scope::outer"),
+                    vec![symbol!("x")],
+                    alias.clone(),
+                    FunctionRegistrationOptions::new().inlining(policy),
+                )
+                .unwrap()
+                .optimization_settings(settings(direct, use_cpe, iterations))
+                .build()
+                .unwrap()
+                .map_coeff(&|c| c.re.to_f64());
+            for [a, b, x, u] in [[2., 3., 0.25, 0.5], [-3., 0.5, -2., 3.]] {
+                let root: f64 = (a + b) * x;
+                let local: f64 = (a + b) * u;
+                close(
+                    evaluator.evaluate_single(&[a, b, x, u]),
+                    root.sin() + local.sin(),
+                );
+            }
+        }
+    }
+}
diff --git a/tests/python/test_horner_cost.py b/tests/python/test_horner_cost.py
new file mode 100644
index 00000000..0c51d8d2
--- /dev/null
+++ b/tests/python/test_horner_cost.py
@@ -0,0 +1,106 @@
+"""Horner scoring options preserve results and final pair elimination."""
+from itertools import permutations
+
+import pytest
+
+from symbolica import Expression, S
+
+
+def _polynomial(a, b, c):
+    return (b * c**3 - 3 * b**2 - a * c**2 + 2 * a**2 * c
+            - 2 * a**2 * b**2 - 3 * a**3 * c + 2 * a**3 * c**2
+            + 3 * a**3 * b**2 * c - 2 * a**3 * b**2 * c**2)
+
+
+def _build(expression, params, multiple, **settings):
+    if multiple:
+        return Expression.evaluator_multiple(
+            [expression, expression + 1], params, **settings)
+    return expression.evaluator(params, **settings)
+
+
+def _cost(evaluator):
+    additions = multiplications = 0
+    for instruction in evaluator.get_instructions().instructions:
+        if instruction[0] == "add":
+            additions += len(instruction[2]) - 1
+        elif instruction[0] == "mul":
+            multiplications += len(instruction[2]) - 1
+        elif instruction[0] == "pow":
+            exponent = instruction[3]
+            assert exponent > 0
+            multiplications += exponent.bit_length() + exponent.bit_count() - 2
+        else:
+            raise AssertionError(f"Unexpected polynomial instruction: {instruction}")
+    return multiplications, additions
+
+
+@pytest.mark.parametrize("multiple", [False, True])
+@pytest.mark.parametrize("horner_use_cpe", [False, True])
+@pytest.mark.parametrize("direct_translation", [False, True])
+def test_horner_cost_option_preserves_numerical_outputs(multiple, horner_use_cpe, direct_translation):
+    a, b, c = S("horner_option_a", "horner_option_b", "horner_option_c")
+    evaluator = _build(
+        _polynomial(a, b, c), [a, b, c], multiple,
+        iterations=6, n_cores=1, jit_compile=False,
+        direct_translation=direct_translation,
+        horner_use_cpe=horner_use_cpe)
+    points = [[2., 3., 4.], [-1., 2., -3.], [0.25, -0.5, 0.75]]
+    expected = [_polynomial(*point) for point in points]
+    rows = [[value, value + 1] if multiple else [value] for value in expected]
+    assert evaluator.evaluate(points).tolist() == rows
+
+
+@pytest.mark.parametrize("multiple", [False, True])
+def test_horner_cost_default_and_existing_positional_arguments(multiple):
+    a, b, c = S("horner_position_a", "horner_position_b", "horner_position_c")
+    expression = _polynomial(a, b, c)
+    params = [a, b, c]
+    # Every existing optional positional argument retains its original slot.
+    arguments = (params, [], 6, None, 1, False, False, True, False, 3,
+                 {}, 500, 1000000, 100)
+    if multiple:
+        constructor = Expression.evaluator_multiple
+        arguments = ([expression, expression + 1], *arguments)
+    else:
+        constructor = expression.evaluator
+    default = constructor(*arguments)
+    explicit_false = constructor(*arguments, False)
+    explicit_true = constructor(*arguments, True)
+    assert default.get_instructions().instructions == explicit_false.get_instructions().instructions
+    points = [[2., 3., 4.], [0.25, -0.5, 0.75]]
+    assert default.evaluate(points).tolist() == explicit_true.evaluate(points).tolist()
+
+
+@pytest.mark.parametrize("multiple", [False, True])
+def test_full_cpe_candidate_cost_can_select_fewer_multiplications(multiple):
+    params = list(S("horner_rank_a", "horner_rank_b", "horner_rank_c"))
+    improved = False
+    # Relabel the same polynomial to avoid depending on symbol interning order.
+    for a, b, c in permutations(params):
+        expression = _polynomial(a, b, c)
+        settings = dict(iterations=6, n_cores=1, jit_compile=False)
+        cheap = _build(expression, params, multiple, horner_use_cpe=False, **settings)
+        full = _build(expression, params, multiple, horner_use_cpe=True, **settings)
+        cheap_cost, full_cost = _cost(cheap), _cost(full)
+        # Six candidates exhaust every order of these three variables.
+        assert full_cost <= cheap_cost
+        improved |= full_cost[0] < cheap_cost[0]
+        assert full.evaluate([[2., 3., 4.]]).tolist() == cheap.evaluate([[2., 3., 4.]]).tolist()
+    assert improved
+
+
+@pytest.mark.parametrize("multiple", [False, True])
+@pytest.mark.parametrize("horner_use_cpe", [False, True])
+@pytest.mark.parametrize("direct_translation", [False, True])
+def test_final_pair_elimination_is_independent_of_horner_cost(multiple, horner_use_cpe, direct_translation):
+    a, b, c, d = S("horner_pair_a", "horner_pair_b", "horner_pair_c", "horner_pair_d")
+    expression = (a + b + c) * (a + b + d)
+    settings = dict(iterations=1, n_cores=1, jit_compile=False,
+                    direct_translation=direct_translation,
+                    horner_use_cpe=horner_use_cpe)
+    without_pairs = _build(expression, [a, b, c, d], multiple, cpe_iterations=0, **settings)
+    with_pairs = _build(expression, [a, b, c, d], multiple, cpe_iterations=None, **settings)
+    assert _cost(with_pairs)[1] < _cost(without_pairs)[1]
+    assert with_pairs.evaluate([[1., 2., 3., 4.]]).tolist() == (
+        [[42., 43.]] if multiple else [[42.]])