Imports

import Mathlib

35.5. The Subset-Sum Problem

The costed execution records the work of the executable approximation scheme.

This section formalizes the subset-sum problem and the fully polynomial-time approximation scheme APPROX-SUBSET-SUM of CLRS §35.5. Given a finite set S = {x₁, ..., xₙ} of positive integers and a target t, the subset-sum problem asks for the subset whose sum is as large as possible without exceeding t. The exact algorithm EXACT-SUBSET-SUM builds the sorted list of all achievable subset sums at most t (its last element is the optimum), but its lists can grow to size 2^n. The approximation scheme instead trims each list so that consecutive kept sums differ by a factor of 1 + δ, retaining a (1 + δ)-representative for every removed value; trimming at each level with δ = ε/(2n) accumulates only a (1 + ε)-factor of error over all n levels.

Main results:

Definition subsetSums: the set of all subset sums of a list of integers, built by the same merge step as EXACT-SUBSET-SUM (Lᵢ = Lᵢ₋₁ ∪ (Lᵢ₋₁ + xᵢ)).
Definition exactLists: the output list of EXACT-SUBSET-SUM — the subset sums of S not exceeding t.
Definition optimalSum: the largest achievable subset sum at most t (the optimum y*).
Definition merge: the merge of two sorted lists of sums.
Definition trim/trimAux: the greedy TRIM of a sorted list with parameter δ, keeping the first element and then every element exceeding (1 + δ) times the last kept element.
Lemma trim_rep (Lemma 35.5): every element of the trimmed list's input is represented within a factor (1 + δ) by a kept element — for every y in a sorted L there is z in trim δ L with z ≤ y ≤ (1 + δ) · z.
Definition approxLists: the trimmed lists Lᵢ of APPROX-SUBSET-SUM.
Lemma approxLists_subset_subsetSums: every value on the trimmed lists is still an achievable subset sum at most t (trimming never introduces or loses validity).
Lemma approxLists_prefix_rep: after i levels of trimming with parameter δ, every achievable sum y ≤ t of the first i elements has a representative z in Lᵢ with y ≤ (1 + δ)^i · z and z ≤ y (the compounded per-level factor, Exercise 35.5-2).
Theorem approxSubsetSum_approx (Theorem 35.7): APPROX-SUBSET-SUM returns a valid subset sum z* ≤ t with the explicit bound y* ≤ (1 + ε/(2n))^n · z*.
Theorem approxSubsetSum_approx_lt (Theorem 35.7): with 0 < ε ≤ 1, the explicit factor is absorbed into (1 + ε) via (1 + ε/(2n))^n ≤ e^{ε/2} ≤ 1 + ε, giving y* ≤ (1 + ε) · z*.
Definition Separated: a list is (1 + δ)-separated when every a before b satisfies (1 + δ) · a < b.
Lemma trim_separated: TRIM returns a (1 + δ)-separated list (the premise of Exercise 35.5-1).
Lemma sep_length_bound (Exercise 35.5-1): a (1 + δ)-separated list of naturals in [1, t] has at most log t / log(1 + δ) + 1 elements.
Lemma approxLists_length_bound: every trimmed list of APPROX-SUBSET-SUM has at most log t / log(1 + δ) + 2 elements (including the least element 0).
Theorem approxSubsetSum_fptas (Theorem 35.8): with δ = ε/(2n), every trimmed list has at most 4n · lg t / ε + 2 elements, so APPROX-SUBSET-SUM runs in time polynomial in the input size and in 1/ε — a fully polynomial-time approximation scheme.

Current gaps: none. The running-time analysis of Theorem 35.8 is formalized at the mathematical level (list-size bounds); a machine-level cost model assigning running times to MERGE-LISTS and TRIM steps is not modeled.

Notation conventions used in this section:

xs/ps : a list of positive integers (the set S, ordered for the algorithm); subset sums are order-independent
s, y : a candidate subset sum
t : the target sum
δ : the trim parameter (0 ≤ δ)
ε : the approximation parameter (0 < ε ≤ 1), with δ = ε/(2n)
L, M : sorted lists of sums
y* : the optimal subset sum (optimalSum)
z* : the value returned by APPROX-SUBSET-SUM (approxSum)

              noncomputable sectionopen scoped BigOperatorsopen scoped Listopen Finsetnamespace CLRSnamespace ApproxSubsetSum

The subset-sum problem and EXACT-SUBSET-SUM

CLRS models the set S = {x₁, ..., xₙ} of positive integers and builds the lists Lᵢ of all sums obtainable from subsets of {x₁, ..., xᵢ}, discarding every sum exceeding the target t. The recurrence Lᵢ = Lᵢ₋₁ ∪ (Lᵢ₋₁ + xᵢ) is captured by the recursive definition of subsetSums; exactLists adds the pruning of sums above t and optimalSum is the largest element of the pruned list.

The set of all subset sums of a list xs of integers. The recursion subsetSums (x :: xs) = subsetSums xs ∪ (subsetSums xs) + x is exactly the merge step of EXACT-SUBSET-SUM (CLRS §35.5, line 2).

def subsetSums : List ℕ → Finset ℕ
  | [] => {0}
  | x :: xs => subsetSums xs ∪ (subsetSums xs).image (fun s => s + x)

The empty subset has sum 0, so 0 is always an achievable sum.

lemma zero_mem_subsetSums (xs : List ℕ) : 0 ∈ subsetSums xs := by
  induction xs with
  | nil => simp [subsetSums]
  | cons x xs ih => simp [subsetSums, ih]

A sum of x :: xs is either a sum of xs alone or a sum of xs plus x. This is the CLRS recurrence Lᵢ = Lᵢ₋₁ ∪ (Lᵢ₋₁ + xᵢ).

lemma mem_subsetSums_cons {x : ℕ} {xs : List ℕ} {s : ℕ} :
    s ∈ subsetSums (x :: xs) ↔ s ∈ subsetSums xs ∨ ∃ t ∈ subsetSums xs, t + x = s := by
  simp [subsetSums]

Every sum of xs is also a sum of x :: xs (drop the new element).

lemma subsetSums_cons_subset {x : ℕ} {xs : List ℕ} : subsetSums xs ⊆ subsetSums (x :: xs) := by
  intro s hs
  simp [subsetSums, hs]

If z is a sum of xs, then z + x is a sum of x :: xs.


lemma mem_subsetSums_add {x : ℕ} {xs : List ℕ} {z : ℕ} :
    z ∈ subsetSums xs → z + x ∈ subsetSums (x :: xs) := by
  intro hz
  change z + x ∈ subsetSums xs ∪ (subsetSums xs).image (fun s => s + x)
  rw [Finset.mem_union]
  right
  rw [Finset.mem_image]
  exact ⟨z, hz, rfl⟩

The output of EXACT-SUBSET-SUM(S, t): the set of all subset sums of S not exceeding the target t (CLRS §35.5, EXACT-SUBSET-SUM line 5).

def exactLists (xs : List ℕ) (t : ℕ) : Finset ℕ :=
  (subsetSums xs).filter (fun s => s ≤ t)

The exact list is nonempty: 0 is always achievable and 0 ≤ t.


lemma exactLists_nonempty (xs : List ℕ) (t : ℕ) : (exactLists xs t).Nonempty := by
  exact ⟨0, by simp [exactLists, zero_mem_subsetSums xs, This simp argument is unused:
  Nat.zero_le

Hint: Omit it from the simp argument list.
  simp [exactLists, zero_mem_subsetSums xs,̵ ̵N̵a̵t̵.̵z̵e̵r̵o̵_̵l̵e̵]

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`Nat.zero_le]⟩

The optimal subset sum y*: the largest achievable subset sum of xs not exceeding t. It is well-defined because 0 is always achievable and 0 ≤ t.

def optimalSum (xs : List ℕ) (t : ℕ) : ℕ :=
  (exactLists xs t).max' (exactLists_nonempty xs t)

The optimum is an achievable subset sum.

lemma optimalSum_mem_subsetSums (xs : List ℕ) (t : ℕ) : optimalSum xs t ∈ subsetSums xs := by
  exact (Finset.mem_filter.mp (Finset.max'_mem (exactLists xs t) (exactLists_nonempty xs t))).1

The optimum does not exceed the target t.

lemma optimalSum_le_t (xs : List ℕ) (t : ℕ) : optimalSum xs t ≤ t := by
  exact (Finset.mem_filter.mp (Finset.max'_mem (exactLists xs t) (exactLists_nonempty xs t))).2

The optimum of the empty list is 0.


lemma optimalSum_nil (t : ℕ) : optimalSum [] t = 0 := by
  rw [optimalSum]
  apply le_antisymm
  · apply Finset.max'_le
    intro y hy
    simp [exactLists, subsetSums] at hy
    exact le_of_eq hy.1
  · apply Finset.le_max'
    simp [exactLists, subsetSums, This simp argument is unused:
  Nat.zero_le

Hint: Omit it from the simp argument list.
  simp [exactLists, subsetSums,̵ ̵N̵a̵t̵.̵z̵e̵r̵o̵_̵l̵e̵]

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`Nat.zero_le]

The merge of two sorted lists

APPROX-SUBSET-SUM keeps the lists sorted so that TRIM can scan them greedily. merge is the merge of two sorted lists of sums (CLRS MERGE-LISTS).

The merge of two sorted lists of natural numbers.

def merge : List ℕ → List ℕ → List ℕ
  | [], ys => ys
  | xs, [] => xs
  | x :: xs, y :: ys => if x ≤ y then x :: merge xs (y :: ys) else y :: merge (x :: xs) ys

An element of the merged list is an element of one of the inputs (and vice versa): merge returns the union of its inputs.


lemma mem_merge {L M : List ℕ} {z : ℕ} : z ∈ merge L M ↔ z ∈ L ∨ z ∈ M := by
  induction L generalizing M with
  | nil => simp [merge]
  | cons x xs ihL =>
      induction M with
      | nil => simp [merge]
      | cons y ys ihM =>
          by_cases hxy : x ≤ y
          · simp [merge, hxy]
            rw [ihL]
            simp
            tauto
          · simp [merge, hxy]
            rw [ihM]
            simp
            tauto

The merge of two sorted lists is sorted.


lemma merge_sorted (L M : List ℕ) (hL : L.Pairwise (· ≤ ·)) (hM : M.Pairwise (· ≤ ·)) :
    (merge L M).Pairwise (· ≤ ·) := by
  induction L generalizing M with
  | nil => simpa [merge] using hM
  | cons x xs ihL =>
      induction M with
      | nil => simpa [merge] using hL
      | cons y ys ihM =>
          by_cases hxy : x ≤ y
          · have hsub : ∀ z ∈ merge xs (y :: ys), x ≤ z := by
              intro z hz
              rw [mem_merge] at hz
              rcases hz with hz | hz
              · exact (List.pairwise_cons.mp hL).1 z hz
              · simp at hz
                rcases hz with hz | hz
                · subst z
                  exact hxy
                · exact le_trans hxy ((List.pairwise_cons.mp hM).1 z hz)
            have hsort : (merge xs (y :: ys)).Pairwise (· ≤ ·) := by
              exact ihL (y :: ys) (List.pairwise_cons.mp hL).2 hM
            simpa [merge, hxy] using (List.pairwise_cons.mpr ⟨hsub, hsort⟩)
          · have hle : y ≤ x := le_of_not_ge hxy
            have hsub : ∀ z ∈ merge (x :: xs) ys, y ≤ z := by
              intro z hz
              rw [mem_merge] at hz
              rcases hz with hz | hz
              · simp at hz
                rcases hz with hz | hz
                · subst z
                  exact hle
                · exact le_trans hle ((List.pairwise_cons.mp hL).1 z hz)
              · exact (List.pairwise_cons.mp hM).1 z hz
            have hsort : (merge (x :: xs) ys).Pairwise (· ≤ ·) := by
              exact ihM (List.pairwise_cons.mp hM).2
            simpa [merge, hxy] using (List.pairwise_cons.mpr ⟨hsub, hsort⟩)

Adding a constant to every element of a sorted list keeps it sorted.

lemma map_add_pairwise (x : ℕ) {L : List ℕ} (hL : L.Pairwise (· ≤ ·)) :
    (L.map (fun s => s + x)).Pairwise (· ≤ ·) := by
  exact List.Pairwise.map (fun s => s + x) (by intro a b h; exact Nat.add_le_add_right h x) hL

TRIM and Lemma 35.5

TRIM scans a sorted list in increasing order and keeps the first element, then keeps an element y only when (1 + δ) · last < y, where last is the last kept element; otherwise y is within a factor (1 + δ) of last and is dropped. Lemma 35.5 states the invariant: every element of the input list has a kept representative within a factor (1 + δ) below it.

The greedy tail scan of TRIM: with the previously kept element last, keeps y iff (1 + δ) · last < y, otherwise recurses keeping last.

def trimAux (δ : ℝ) : ℕ → List ℕ → List ℕ
  | Variable name `last` is not explicitly referenced.

The binding can be removed (if unused) or named `_` (if used implicitly).

Note: This linter can be disabled with `set_option linter.unusedVariables false`last, [] => []
  | last, y :: ys =>
      if (1 + δ) * (last : ℝ) < (y : ℝ) then
        y :: trimAux δ y ys
      else
        trimAux δ last ys

The TRIM of a sorted list L with parameter δ (CLRS §35.5, TRIM(L, δ)): keeps the first element and then every element exceeding (1 + δ) times the last kept element.

def trim (δ : ℝ) : List ℕ → List ℕ
  | [] => []
  | y :: ys => y :: trimAux δ y ys

Every element of a trimmed tail was an element of that tail.

lemma trimAux_mem_subset (δ : ℝ) : ∀ (last : ℕ) {ys : List ℕ} {z : ℕ},
    z ∈ trimAux δ last ys → z ∈ ys := by
  intro last ys z hz
  induction ys generalizing last with
  | nil => simp [trimAux] at hz
  | cons y ys ih =>
      by_cases hkeep : (1 + δ) * (last : ℝ) < (y : ℝ)
      · simp [trimAux, hkeep] at hz
        rcases hz with hz | hz
        · subst z
          simp
        · exact List.mem_cons.mpr (Or.inr (ih y hz))
      · simp [trimAux, hkeep] at hz
        exact List.mem_cons.mpr (Or.inr (ih last hz))

Every kept element of TRIM was an element of the input list.

lemma mem_trim_subset (δ : ℝ) {L : List ℕ} {z : ℕ} (hz : z ∈ trim δ L) : z ∈ L := by
  cases L with
  | nil => simp [trim] at hz
  | cons y ys =>
      simp [trim] at hz
      rcases hz with hz | hz
      · subst z
        simp
      · exact List.mem_cons.mpr (Or.inr (trimAux_mem_subset δ y hz))

The tail scan of TRIM returns a sublist of the scanned tail.

lemma trimAux_sublist (δ : ℝ) : ∀ (last : ℕ) (ys : List ℕ), trimAux δ last ys <+ ys := by
  intro last ys
  induction ys generalizing last with
  | nil => simp [trimAux]
  | cons y ys ih =>
      by_cases hkeep : (1 + δ) * (last : ℝ) < (y : ℝ)
      · have hsub : trimAux δ y ys <+ ys := ih y
        simpa [trimAux, hkeep] using hsub.cons_cons y
      · simpa [trimAux, hkeep] using (ih last).cons y

TRIM returns a sublist of its input list.

lemma trim_sublist (δ : ℝ) : ∀ L : List ℕ, trim δ L <+ L := by
  intro L
  cases L with
  | nil => simp [trim]
  | cons y ys =>
      have hsub : trimAux δ y ys <+ ys := trimAux_sublist δ y ys
      simpa [trim] using hsub.cons_cons y

TRIM keeps a subsequence of a sorted list, so its output is sorted.

lemma trim_sorted {δ : ℝ} {L : List ℕ} (hL : L.Pairwise (· ≤ ·)) : (trim δ L).Pairwise (· ≤ ·) := by
  exact List.Pairwise.sublist (trim_sublist δ L) hL

The tail scan represents every element of the scanned tail: for every y in the tail, some z among the tail scan's output plus the reference last satisfies z ≤ y ≤ (1 + δ) · z. This is the induction core of Lemma 35.5.


lemma trimAux_rep {δ : ℝ} (hδ : 0 ≤ δ) :
    ∀ (last : ℕ) {ys : List ℕ}, ys.Pairwise (· ≤ ·) →
      (∀ y ∈ ys, (last : ℝ) ≤ (y : ℝ)) →
      ∀ y ∈ ys, ∃ z ∈ last :: trimAux δ last ys,
        (z : ℝ) ≤ (y : ℝ) ∧ (y : ℝ) ≤ (1 + δ) * (z : ℝ) := by
  intro last ys
  induction ys generalizing last with
  | nil => intro hys hle y hy; simp at hy
  | cons x xs ih =>
      intro hys hle y hy
      simp at hy
      rcases hy with hy | hy
      · subst y
        by_cases hkeep : (1 + δ) * (last : ℝ) < (x : ℝ)
        · refine ⟨x, by simp [trimAux, hkeep], ?_⟩
          have hx0 : (0 : ℝ) ≤ (x : ℝ) := by exact_mod_cast Nat.zero_le x
          constructor
          · rfl
          · nlinarith
        · refine ⟨last, by simp [trimAux, hkeep], ?_⟩
          have hlastx : (last : ℝ) ≤ (x : ℝ) := by exact_mod_cast (hle x (by simp))
          exact ⟨hlastx, le_of_not_gt hkeep⟩
      · by_cases hkeep : (1 + δ) * (last : ℝ) < (x : ℝ)
        · rcases (ih x (List.pairwise_cons.mp hys).2
            (by intro z hz; exact_mod_cast (List.pairwise_cons.mp hys).1 z hz) y hy) with ⟨z, hz, hlo, hhi⟩
          exact ⟨z, by simp [trimAux, hkeep, hz], hlo, hhi⟩
        · rcases (ih last (List.pairwise_cons.mp hys).2
            (by intro z hz; exact_mod_cast (hle z (by simp [hz]))) y hy) with ⟨z, hz, hlo, hhi⟩
          exact ⟨z, by simp [trimAux, hkeep, hz], hlo, hhi⟩

Lemma 35.5 (TRIM). For a sorted list L and δ ≥ 0, every element y of L is represented in TRIM(L, δ) by some z with z ≤ y ≤ (1 + δ) · z: a kept element represents itself and a dropped element is within a factor (1 + δ) of the last kept element below it (CLRS §35.5, Lemma 35.5).


lemma trim_rep {δ : ℝ} (hδ : 0 ≤ δ) {L : List ℕ} (hL : L.Pairwise (· ≤ ·)) :
    ∀ y ∈ L, ∃ z ∈ trim δ L, (z : ℝ) ≤ (y : ℝ) ∧ (y : ℝ) ≤ (1 + δ) * (z : ℝ) := by
  induction L with
  | nil => intro y hy; simp at hy
  | cons x xs ih =>
      intro y hy
      simp at hy
      rcases hy with hy | hy
      · subst y
        refine ⟨x, by simp [trim], ?_⟩
        have hx0 : (0 : ℝ) ≤ (x : ℝ) := by exact_mod_cast Nat.zero_le x
        constructor
        · rfl
        · nlinarith
      · rcases (trimAux_rep hδ x (List.pairwise_cons.mp hL).2
            (by intro z hz; exact_mod_cast (List.pairwise_cons.mp hL).1 z hz) y hy)
          with ⟨z, hz, hlo, hhi⟩
        exact ⟨z, by simp [trim, hz], hlo, hhi⟩

APPROX-SUBSET-SUM and its approximation guarantee

APPROX-SUBSET-SUM runs the exact algorithm but inserts Lᵢ ← TRIM(Lᵢ, ε/(2n)) after each merge and then drops every element exceeding t; it returns the largest element z* of the final list. The trimmed lists never lose a value beyond a compounded (1 + ε/(2n))^i factor, which gives Theorem 35.7.

The trimmed lists Lᵢ of APPROX-SUBSET-SUM: after processing the whole list, every element exceeding t is dropped. approxLists δ t xs is the final list Lₙ for S = xs with trim parameter δ.

def approxLists (δ : ℝ) (t : ℕ) : List ℕ → List ℕ
  | [] => [0]
  | x :: xs =>
      let L := approxLists δ t xs
      (trim δ (merge L (L.map (fun s => s + x)))).filter (fun s => s ≤ t)

The trimmed lists stay sorted.

lemma approxLists_sorted (δ : ℝ) (t : ℕ) : ∀ xs : List ℕ, (approxLists δ t xs).Pairwise (· ≤ ·) := by
  intro xs
  induction xs with
  | nil => simp [approxLists]
  | cons x xs ih =>
      have hmerge : (merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x))).Pairwise (· ≤ ·) :=
        merge_sorted _ _ ih (map_add_pairwise x ih)
      have htrim : (trim δ (merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x)))).Pairwise (· ≤ ·) :=
        trim_sorted hmerge
      have hfilter : ((trim δ (merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x)))).filter (fun s => s ≤ t)).Pairwise (· ≤ ·) :=
        List.Pairwise.filter (fun s => s ≤ t) htrim
      simpa [approxLists] using hfilter

If 0 belongs to a sorted list, TRIM keeps it (it is the least element, hence the head).


lemma zero_mem_trim {δ : ℝ} {L : List ℕ} (hL : L.Pairwise (· ≤ ·)) (h0 : 0 ∈ L) :
    0 ∈ trim δ L := by
  cases L with
  | nil => simp at h0
  | cons y ys =>
      have hy0 : y = 0 := by
        simp at h0
        rcases h0 with h0 | h0
        · exact h0.symm
        · have hyz : y ≤ 0 := (List.pairwise_cons.mp hL).1 0 h0
          exact le_antisymm hyz (Nat.zero_le y)
      subst y
      simp [trim]

0 is always present on the trimmed lists (the empty subset is never lost and never exceeds the target).


lemma zero_mem_approxLists (δ : ℝ) (t : ℕ) (xs : List ℕ) : 0 ∈ approxLists δ t xs := by
  induction xs with
  | nil => simp [approxLists]
  | cons x xs ih =>
      have h0L : 0 ∈ approxLists δ t xs := ih
      have h0merged : 0 ∈ merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x)) := by
        exact (mem_merge).2 (Or.inl h0L)
      have h0trim : 0 ∈ trim δ (merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x))) := by
        exact zero_mem_trim
          (merge_sorted _ _ (approxLists_sorted δ t xs) (map_add_pairwise x (approxLists_sorted δ t xs))) h0merged
      exact List.mem_filter.mpr ⟨h0trim, by simp⟩

Every value on the trimmed lists is an achievable subset sum of the input list (TRIM only drops values and MERGE only combines existing values), so APPROX-SUBSET-SUM never returns a sum that is not a subset sum.


lemma approxLists_subset_subsetSums (δ : ℝ) (t : ℕ) :
    ∀ xs : List ℕ, ∀ z ∈ approxLists δ t xs, z ∈ subsetSums xs := by
  intro xs
  induction xs with
  | nil => intro z hz; simp [approxLists] at hz; subst z; simp [subsetSums]
  | cons x xs ih =>
      intro z hz
      simp [approxLists] at hz
      rcases hz with ⟨hztrim, hzle⟩
      have hzmerge : z ∈ merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x)) :=
        mem_trim_subset δ hztrim
      rw [mem_merge] at hzmerge
      rcases hzmerge with hzL | hzLx
      · exact subsetSums_cons_subset (ih z hzL)
      · rcases (List.mem_map.mp hzLx) with ⟨s, hsL, hs⟩
        subst z
        exact mem_subsetSums_add (ih s hsL)

Every value on the trimmed lists is at most the target t (they pass the ≤ t filter of APPROX-SUBSET-SUM).

lemma mem_approxLists_le_t (δ : ℝ) (t : ℕ) :
    ∀ xs : List ℕ, ∀ z ∈ approxLists δ t xs, z ≤ t := by
  intro xs
  induction xs with
  | nil => intro z hz; simp [approxLists] at hz; subst z; exact Nat.zero_le t
  | cons x xs ih =>
      intro z hz
      simp [approxLists] at hz
      exact hz.2

After i levels of trimming with parameter δ, every achievable sum y ≤ t of the first i elements has a representative z in Lᵢ with y ≤ (1 + δ)^i · z and z ≤ y. Each level of TRIM contributes one factor (1 + δ); the factor compounds over the i merges (CLRS §35.5, Exercise 35.5-2).


lemma approxLists_prefix_rep {δ : ℝ} (hδ : 0 ≤ δ) (t : ℕ) :
    ∀ ps : List ℕ, ∀ y ∈ subsetSums ps, y ≤ t →
      ∃ z ∈ approxLists δ t ps,
        ((y : ℝ) ≤ (1 + δ) ^ ps.length * (z : ℝ)) ∧ ((z : ℝ) ≤ (y : ℝ)) := by
  intro ps
  induction ps with
  | nil =>
      intro y hy hle
      have hy0 : y = 0 := by
        simpa [subsetSums] using hy
      subst y
      refine ⟨0, by simp [approxLists], ?_⟩
      simp
  | cons x xs ih =>
      intro y hy hle
      rcases (mem_subsetSums_cons.mp hy) with hyxs | ⟨s, hsxs, hs⟩
      · rcases (ih y hyxs hle) with ⟨z0, hz0, hlo0, hhi0⟩
        have hz0merged : z0 ∈ merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x)) := by
          exact (mem_merge).2 (Or.inl hz0)
        rcases (trim_rep hδ
            (merge_sorted _ _ (approxLists_sorted δ t xs) (map_add_pairwise x (approxLists_sorted δ t xs)))
            z0 hz0merged) with ⟨z, hztrim, hzle_z0, hz0_le⟩
        have hzle_t : z ≤ t := by
          exact le_trans (by exact_mod_cast (hzle_z0.trans hhi0)) hle
        refine ⟨z, by simp [approxLists, hztrim, hzle_t], ?_⟩
        constructor
        · rw [List.length_cons]
          have hpow0 : 0 ≤ (1 + δ) ^ xs.length := pow_nonneg (by linarith) _
          have hscale : (1 + δ) ^ xs.length * (z0 : ℝ) ≤ (1 + δ) ^ xs.length * ((1 + δ) * (z : ℝ)) := by
            exact mul_le_mul_of_nonneg_left hz0_le hpow0
          calc
            (y : ℝ) ≤ (1 + δ) ^ xs.length * (z0 : ℝ) := hlo0
            _ ≤ (1 + δ) ^ xs.length * ((1 + δ) * (z : ℝ)) := hscale
            _ = (1 + δ) ^ (xs.length + 1) * (z : ℝ) := by
              rw [pow_succ]
              ring
        · exact hzle_z0.trans hhi0
      · subst y
        have hsle : s ≤ t := by omega
        rcases (ih s hsxs hsle) with ⟨z0, hz0, hlo0, hhi0⟩
        let w : ℕ := z0 + x
        have hwmerged : w ∈ merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x)) := by
          apply (mem_merge).2
          right
          exact List.mem_map.mpr ⟨z0, hz0, rfl⟩
        rcases (trim_rep hδ
            (merge_sorted _ _ (approxLists_sorted δ t xs) (map_add_pairwise x (approxLists_sorted δ t xs)))
            w hwmerged) with ⟨z, hztrim, hzle_w, hw_le⟩
        have hzle_t : z ≤ t := by
          have hwy : w ≤ s + x := by
            dsimp [w]
            exact Nat.add_le_add_right (by exact_mod_cast hhi0) x
          exact le_trans (by exact_mod_cast (hzle_w.trans (by exact_mod_cast hwy : (w : ℝ) ≤ (s + x : ℝ)))) hle
        refine ⟨z, by simp [approxLists, hztrim, hzle_t], ?_⟩
        constructor
        · rw [List.length_cons]
          have hpow0 : 0 ≤ (1 + δ) ^ xs.length := pow_nonneg (by linarith) _
          have hpowge1 : (1 : ℝ) ≤ (1 + δ) ^ xs.length := by
            simpa using (pow_le_pow_left₀ (a := (1 : ℝ)) (b := (1 + δ)) (by norm_num) (by linarith) xs.length)
          have hzx_le : (z0 : ℝ) + (x : ℝ) ≤ (1 + δ) * (z : ℝ) := by
            have hw_eq : (z0 : ℝ) + (x : ℝ) = (w : ℝ) := by
              exact_mod_cast (show z0 + x = w by rfl)
            rw [hw_eq]
            exact hw_le
          have hxscale : (x : ℝ) ≤ (1 + δ) ^ xs.length * (x : ℝ) := by
            exact le_mul_of_one_le_left (by exact_mod_cast Nat.zero_le x) hpowge1
          have hz0x_le : (1 + δ) ^ xs.length * (z0 : ℝ) + (x : ℝ) ≤ (1 + δ) ^ xs.length * ((z0 : ℝ) + (x : ℝ)) := by
            nlinarith [hxscale]
          calc
            ↑(s + x) ≤ (1 + δ) ^ xs.length * (z0 : ℝ) + (x : ℝ) := by
              have hdist : ↑(s + x) = (s : ℝ) + (x : ℝ) := by norm_num
              rw [hdist]
              simpa [add_comm] using add_le_add_right hlo0 (x : ℝ)
            _ ≤ (1 + δ) ^ xs.length * ((z0 : ℝ) + (x : ℝ)) := hz0x_le
            _ ≤ (1 + δ) ^ xs.length * ((1 + δ) * (z : ℝ)) := by
              exact mul_le_mul_of_nonneg_left hzx_le hpow0
            _ = (1 + δ) ^ (xs.length + 1) * (z : ℝ) := by
              rw [pow_succ]
              ring
        · have hwy : w ≤ s + x := by
            dsimp [w]
            exact Nat.add_le_add_right (by exact_mod_cast hhi0) x
          exact (by exact_mod_cast hzle_w : (z : ℝ) ≤ (w : ℝ)).trans (by exact_mod_cast hwy)

The value z* returned by APPROX-SUBSET-SUM(S, t, ε): the largest element of the final trimmed list, computed with trim parameter ε/(2n) (CLRS §35.5, APPROX-SUBSET-SUM line 8).

try 'simp' instead of 'simpa'

Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`
def approxSum (xs : List ℕ) (t : ℕ) (ε : ℝ) : ℕ :=
  ((approxLists (ε / (2 * (xs.length : ℝ))) t xs).toFinset).max' (by
    exact ⟨0, by try 'simp' instead of 'simpa'

Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa [zero_mem_approxLists]⟩)

The value z* belongs to the final trimmed list of APPROX-SUBSET-SUM.


lemma approxSum_mem (xs : List ℕ) (t : ℕ) (ε : ℝ) :
    approxSum xs t ε ∈ approxLists (ε / (2 * (xs.length : ℝ))) t xs := by
  have hm : approxSum xs t ε ∈ (approxLists (ε / (2 * (xs.length : ℝ))) t xs).toFinset := by
    dsimp [approxSum]
    exact Finset.max'_mem (approxLists (ε / (2 * (xs.length : ℝ))) t xs).toFinset
      (by exact ⟨0, by simp [zero_mem_approxLists]⟩)
  simpa using hm

The value returned for the empty list is 0.

lemma approxSum_nil (t : ℕ) (ε : ℝ) : approxSum [] t ε = 0 := by
  simp [approxSum, approxLists]

APPROX-SUBSET-SUM returns a value z* that is an achievable subset sum: z* ∈ subsetSums xs (Theorem 35.7, correctness).

lemma approxSum_mem_subsetSums (xs : List ℕ) (t : ℕ) (ε : ℝ) :
    approxSum xs t ε ∈ subsetSums xs := by
  exact approxLists_subset_subsetSums (ε / (2 * (xs.length : ℝ))) t xs (approxSum xs t ε)
    (approxSum_mem xs t ε)

APPROX-SUBSET-SUM returns a value z* not exceeding the target t (Theorem 35.7, correctness).

lemma approxSum_le_t (xs : List ℕ) (t : ℕ) (ε : ℝ) : approxSum xs t ε ≤ t := by
  exact mem_approxLists_le_t (ε / (2 * (xs.length : ℝ))) t xs (approxSum xs t ε)
    (approxSum_mem xs t ε)

Theorem 35.7 (explicit factor). APPROX-SUBSET-SUM(S, t, ε) returns a value z* such that the optimum y* satisfies y* ≤ (1 + ε/(2n))^n · z*, where n = |S|. This is the compounded (1 + ε/(2n))^n representation error over all n trimming levels (CLRS §35.5, Theorem 35.7).


theorem approxSubsetSum_approx (xs : List ℕ) (t : ℕ) (ε : ℝ) (hε : 0 ≤ ε) :
    (optimalSum xs t : ℝ) ≤
      (1 + ε / (2 * (xs.length : ℝ))) ^ xs.length * (approxSum xs t ε : ℝ) := by
  let δ : ℝ := ε / (2 * (xs.length : ℝ))
  have hδ : 0 ≤ δ := by
    dsimp [δ]
    exact div_nonneg hε (by positivity)
  have hrep := approxLists_prefix_rep (δ := δ) (hδ := hδ) t xs
  rcases hrep (optimalSum xs t) (optimalSum_mem_subsetSums xs t) (optimalSum_le_t xs t)
      with ⟨z, hz, hlo, hhi⟩
  have hzle : (z : ℝ) ≤ (approxSum xs t ε : ℝ) := by
    have hzs : z ∈ (approxLists (ε / (2 * (xs.length : ℝ))) t xs).toFinset := by
      simpa [δ] using hz
    have hzleN : z ≤ approxSum xs t ε := by
      exact Finset.le_max' (approxLists (ε / (2 * (xs.length : ℝ))) t xs).toFinset z hzs
    exact_mod_cast hzleN
  have hpow0 : 0 ≤ (1 + δ) ^ xs.length := pow_nonneg (by linarith) _
  calc
    (optimalSum xs t : ℝ) ≤ (1 + δ) ^ xs.length * (z : ℝ) := hlo
    _ ≤ (1 + δ) ^ xs.length * (approxSum xs t ε : ℝ) := by
      exact mul_le_mul_of_nonneg_left hzle hpow0

The `(1 + ε)` bound (Theorem 35.7, final form)

To reach the clean (1 + ε) factor of CLRS Theorem 35.7, the compounded error (1 + ε/(2n))^n must be absorbed: since (1 + ε/(2n)) ≤ e^{ε/(2n)}, raising to the n-th power gives (1 + ε/(2n))^n ≤ e^{ε/2}, and e^{ε/2} ≤ 1 + ε for 0 ≤ ε ≤ 1.

For 0 ≤ ε ≤ 1, (1 + ε/(2n))^n ≤ 1 + ε: the compounded per-level error (1 + ε/(2n))^n of the n trims is bounded by e^{ε/2} ≤ 1 + ε (CLRS §35.5, inequalities (35.26)-(35.29)).


lemma pow_one_add_half_le_one_add {ε : ℝ} {n : ℕ} (hε0 : 0 ≤ ε) (hε1 : ε ≤ 1) :
    (1 + ε / (2 * (n : ℝ))) ^ n ≤ 1 + ε := by
  by_cases hn : n = 0
  · subst n
    simp [hε0]
  · have hstep1 : (1 + ε / (2 * (n : ℝ))) ^ n ≤ Real.exp (ε / 2) := by
      have hδnz : 0 ≤ 1 + ε / (2 * (n : ℝ)) := by positivity
      have hbase : 1 + ε / (2 * (n : ℝ)) ≤ Real.exp (ε / (2 * (n : ℝ))) := by
        simpa [add_comm] using (Real.add_one_le_exp (ε / (2 * (n : ℝ))))
      have hpow : (1 + ε / (2 * (n : ℝ))) ^ n ≤ (Real.exp (ε / (2 * (n : ℝ)))) ^ n := by
        exact pow_le_pow_left₀ hδnz hbase n
      have htwo : (2 * (n : ℝ)) ≠ 0 := by
        exact_mod_cast (show (2 * n : ℕ) ≠ 0 by omega)
      calc
        (1 + ε / (2 * (n : ℝ))) ^ n ≤ (Real.exp (ε / (2 * (n : ℝ)))) ^ n := hpow
        _ = Real.exp (n * (ε / (2 * (n : ℝ)))) := by rw [Real.exp_nat_mul]
        _ = Real.exp (ε / 2) := by
          congr 1
          field_simp [htwo]
    have hstep2 : Real.exp (ε / 2) ≤ 1 + ε := by
      have hx0 : 0 ≤ ε / 2 := by linarith
      have hx1 : ε / 2 ≤ 1 := by linarith
      have hx2 : |ε / 2| ≤ 1 := by
        rw [abs_of_nonneg hx0]
        exact hx1
      have hb := Real.abs_exp_sub_one_sub_id_le (x := ε / 2) hx2
      have hb' : Real.exp (ε / 2) - 1 - ε / 2 ≤ (ε / 2) ^ 2 := (abs_le.mp hb).2
      have hexp : Real.exp (ε / 2) ≤ 1 + ε / 2 + (ε / 2) ^ 2 := by linarith
      have hsq : (ε / 2) ^ 2 ≤ ε / 2 := by
        have hmul := mul_le_mul_of_nonneg_left hx1 hx0
        nlinarith
      have htail : 1 + ε / 2 + (ε / 2) ^ 2 ≤ 1 + ε := by
        nlinarith [hsq]
      exact le_trans hexp htail
    exact le_trans hstep1 hstep2

Theorem 35.7. For 0 < ε ≤ 1, APPROX-SUBSET-SUM(S, t, ε) is a (1 + ε)-approximation of the subset-sum problem: the optimum y* satisfies y* ≤ (1 + ε) · z* for the returned value z* (CLRS §35.5, Theorem 35.7).


theorem approxSubsetSum_approx_lt (xs : List ℕ) (t : ℕ) (ε : ℝ)
    (hε0 : 0 < ε) (hε1 : ε ≤ 1) :
    (optimalSum xs t : ℝ) ≤ (1 + ε) * (approxSum xs t ε : ℝ) := by
  by_cases hn : xs.length = 0
  · have hlen : xs = [] := by simpa using hn
    subst xs
    rw [optimalSum_nil, approxSum_nil]
    simp
  · have hεnz : 0 ≤ ε := le_of_lt hε0
    have hpow := pow_one_add_half_le_one_add (ε := ε) (n := xs.length) hεnz hε1
    have hmain := approxSubsetSum_approx xs t ε hεnz
    have hz0 : 0 ≤ (approxSum xs t ε : ℝ) := by exact_mod_cast Nat.zero_le (approxSum xs t ε)
    calc
      (optimalSum xs t : ℝ) ≤ (1 + ε / (2 * (xs.length : ℝ))) ^ xs.length * (approxSum xs t ε : ℝ) := hmain
      _ ≤ (1 + ε) * (approxSum xs t ε : ℝ) := by
        exact mul_le_mul_of_nonneg_right hpow hz0

The running-time analysis (Theorem 35.8, FPTAS)

The approximation guarantee of Theorem 35.7 shows APPROX-SUBSET-SUM is a (1 + ε)-approximation; to see that it is a fully polynomial-time scheme the size of the trimmed lists must be bounded. TRIM keeps consecutive kept values more than a multiplicative factor 1 + δ apart, so a (1 + δ)-separated list with every value in [1, t] has at most log t / log(1 + δ) + 1 elements (Exercise 35.5-1). With δ = ε/(2n) and the chord bound log(1 + δ) ≥ δ/2, each list has O(n · lg t / ε) elements; since MERGE-LISTS and TRIM each cost linear time in the list length, the whole algorithm runs in time O(n² · lg t / ε), polynomial in the input size and in 1/ε. This is Theorem 35.8: APPROX-SUBSET-SUM is a fully polynomial-time approximation scheme.

A list is (1 + δ)-separated when every a occurring before b satisfies (1 + δ) · a < b: the kept values of TRIM stay more than a multiplicative factor 1 + δ apart.

def Separated (δ : ℝ) (L : List ℕ) : Prop :=
  L.Pairwise (fun a b : ℕ => (1 + δ) * (a : ℝ) < (b : ℝ))

TRIM's tail scan keeps the previously kept element last followed by the kept tail; every kept value is more than a factor (1 + δ) above every earlier kept value.


lemma trimAux_separated {δ : ℝ} (hδ : 0 ≤ δ) :
    ∀ last ys, Separated δ (last :: trimAux δ last ys) := by
  intro last ys
  induction ys generalizing last with
  | nil => simp [Separated, trimAux]
  | cons y ys ih =>
      by_cases hkeep : (1 + δ) * (last : ℝ) < (y : ℝ)
      · have htail : Separated δ (y :: trimAux δ y ys) := ih y
        rw [trimAux, if_pos hkeep]
        apply List.pairwise_cons.mpr
        constructor
        · intro z hz
          rcases List.mem_cons.mp hz with hz | hz
          · subst z
            exact hkeep
          · have hsep : (1 + δ) * (y : ℝ) < (z : ℝ) := (List.pairwise_cons.mp htail).1 z hz
            have hy0 : (0 : ℝ) ≤ (y : ℝ) := by exact_mod_cast Nat.zero_le y
            have hδy : (0 : ℝ) ≤ δ * (y : ℝ) := mul_nonneg hδ hy0
            have hle : (y : ℝ) ≤ (1 + δ) * (y : ℝ) := by nlinarith
            exact lt_trans hkeep (lt_of_le_of_lt hle hsep)
        · exact htail
      · simpa [trimAux, hkeep] using ih last

TRIM produces a (1 + δ)-separated list: every kept value exceeds every earlier kept value by more than a factor 1 + δ.

lemma trim_separated {δ : ℝ} (hδ : 0 ≤ δ) : ∀ L : List ℕ, Separated δ (trim δ L) := by
  intro L
  cases L with
  | nil => simp [Separated, trim]
  | cons y ys => simpa [trim] using (trimAux_separated hδ y ys)

Filtering a separated list keeps it separated (a sublist of a separated list is separated).

lemma separated_filter {δ : ℝ} {L : List ℕ} (p : ℕ → Prop) [DecidablePred p]
    (h : Separated δ L) : Separated δ (L.filter p) := by
  unfold Separated at *
  exact List.Pairwise.sublist List.filter_sublist h

The trimmed lists of APPROX-SUBSET-SUM are (1 + δ)-separated.

lemma approxLists_separated {δ : ℝ} (hδ : 0 ≤ δ) (t : ℕ) :
    ∀ xs, Separated δ (approxLists δ t xs) := by
  intro xs
  induction xs with
  | nil => simp [approxLists, Separated]
  | cons x xs ih =>
      dsimp [approxLists]
      exact separated_filter (fun s => s ≤ t)
        (trim_separated hδ (merge (approxLists δ t xs) ((approxLists δ t xs).map (fun s => s + x))))

The counting core of Exercise 35.5-1: a (1 + δ)-separated list of naturals whose elements lie in the real interval [c, t] has at most log(t / c) / log(1 + δ) + 1 elements.


lemma sep_length_aux {δ t : ℝ} (hδ : 0 < δ) (c : ℝ) (hc : 1 ≤ c) (hct : c ≤ t) :
    ∀ L : List ℕ, Separated δ L →
      (∀ a ∈ L, c ≤ (a : ℝ)) →
      (∀ a ∈ L, (a : ℝ) ≤ t) →
      (L.length : ℝ) ≤ Real.log (t / c) / Real.log (1 + δ) + 1 := by
  intro L
  induction L generalizing c hc hct with
  | nil =>
      intro hsep hlo hle
      have hc0 : 0 < c := lt_of_lt_of_le zero_lt_one hc
      have h1 : (1 : ℝ) ≤ t / c := by
        rw [le_div_iff₀ hc0]
        simpa using hct
      have hlg : 0 < Real.log (1 + δ) := Real.log_pos (by linarith)
      have hlg0 : (0 : ℝ) ≤ Real.log (t / c) / Real.log (1 + δ) :=
        div_nonneg (Real.log_nonneg h1) (le_of_lt hlg)
      have hz : (0 : ℝ) ≤ Real.log (t / c) / Real.log (1 + δ) + 1 := by linarith [hlg0]
      simpa using hz
  | cons a rest ih =>
      intro hsep hlo hle
      have hsepR : Separated δ rest := (List.pairwise_cons.mp hsep).2
      have hleR : ∀ b ∈ rest, (b : ℝ) ≤ t := by
        intro b hb
        exact hle b (by simp [hb])
      have ha_c : c ≤ (a : ℝ) := hlo a (by simp)
      have hla : (1 : ℝ) ≤ (a : ℝ) := le_trans hc ha_c
      have ha0 : (0 : ℝ) < (a : ℝ) := lt_of_lt_of_le zero_lt_one hla
      have hc0 : 0 < c := lt_of_lt_of_le zero_lt_one hc
      have ht0 : 0 < t := lt_of_lt_of_le hc0 hct
      have hlg : 0 < Real.log (1 + δ) := Real.log_pos (by linarith)
      cases rest with
      | nil =>
          have h1 : (1 : ℝ) ≤ t / c := by
            rw [le_div_iff₀ hc0]
            simpa using hct
          have hlg0 : (0 : ℝ) ≤ Real.log (t / c) / Real.log (1 + δ) :=
            div_nonneg (Real.log_nonneg h1) (le_of_lt hlg)
          have : (1 : ℝ) ≤ Real.log (t / c) / Real.log (1 + δ) + 1 := by linarith [hlg0]
          simpa using this
      | cons b rest' =>
          let c' : ℝ := (1 + δ) * (a : ℝ)
          have hδ0 : 0 ≤ δ := le_of_lt hδ
          have hc'1 : (1 : ℝ) ≤ c' := by
            dsimp [c']
            have hδa : (0 : ℝ) ≤ δ * (a : ℝ) := mul_nonneg hδ0 (by exact_mod_cast Nat.zero_le a)
            nlinarith
          have hc't : c' ≤ t := by
            have hsep_b : (1 + δ) * (a : ℝ) < (b : ℝ) :=
              List.rel_of_pairwise_cons hsep (by simp)
            have hb_t : (b : ℝ) ≤ t := hle b (by simp)
            dsimp [c']
            nlinarith
          have hbc' : ∀ y ∈ b :: rest', c' ≤ (y : ℝ) := by
            intro y hy
            have hsep_y : (1 + δ) * (a : ℝ) < (y : ℝ) :=
              List.rel_of_pairwise_cons hsep (by simpa using hy)
            dsimp [c']
            exact le_of_lt hsep_y
          have hc'0 : 0 < c' := lt_of_lt_of_le zero_lt_one hc'1
          have hloga : Real.log c ≤ Real.log (a : ℝ) := by
            apply Real.log_le_log
            · exact hc0
            · exact ha_c
          have hlogc' : Real.log c' = Real.log (1 + δ) + Real.log (a : ℝ) := by
            dsimp [c']
            rw [Real.log_mul (ne_of_gt (by linarith : (0 : ℝ) < 1 + δ)) (ne_of_gt ha0)]
          have hlt : Real.log (1 + δ) ≤ Real.log c' - Real.log c := by
            nlinarith [hloga, hlogc']
          have hlogtc' : Real.log (t / c') = Real.log t - Real.log c' := by
            rw [Real.log_div (ne_of_gt ht0) (ne_of_gt hc'0)]
          have hlogtc : Real.log (t / c) = Real.log t - Real.log c := by
            rw [Real.log_div (ne_of_gt ht0) (ne_of_gt hc0)]
          have hmain : Real.log (1 + δ) + Real.log (t / c') ≤ Real.log (t / c) := by
            rw [hlogtc', hlogtc]
            nlinarith [hlt]
          have hgoal : (1 : ℝ) + Real.log (t / c') / Real.log (1 + δ) ≤
              Real.log (t / c) / Real.log (1 + δ) := by
            have hmul : (1 + Real.log (t / c') / Real.log (1 + δ)) * Real.log (1 + δ) ≤
                Real.log (t / c) := by
              field_simp [ne_of_gt hlg]
              exact hmain
            exact (le_div_iff₀ hlg).mpr hmul
          have hm : (((b :: rest').length : ℝ) ≤ Real.log (t / c') / Real.log (1 + δ) + 1) :=
            ih c' hc'1 hc't hsepR hbc' hleR
          have hres : ((a :: b :: rest').length : ℝ) ≤
              Real.log (t / c) / Real.log (1 + δ) + 1 := by
            have htail' : ((a :: b :: rest').length : ℝ) =
                ((b :: rest').length : ℝ) + 1 := by
              simp
            rw [htail']
            nlinarith [hm, hgoal]
          simpa using hres

Exercise 35.5-1 (list size). A (1 + δ)-separated list of naturals with every element in [1, t] has at most ⌊log_{1+δ} t⌋ + 1 elements; here the real-valued bound log t / log(1 + δ) + 1 (CLRS §35.5, Exercise 35.5-1).


lemma sep_length_bound {δ : ℝ} (hδ : 0 < δ) (t : ℕ) (ht : 1 ≤ t) {L : List ℕ}
    (hsep : Separated δ L)
    (hpos : ∀ a ∈ L, (1 : ℝ) ≤ (a : ℝ))
    (hle : ∀ a ∈ L, (a : ℝ) ≤ (t : ℝ)) :
    (L.length : ℝ) ≤ Real.log (t : ℝ) / Real.log (1 + δ) + 1 := by
  have h := sep_length_aux hδ 1 (by norm_num) (by exact_mod_cast ht) L hsep hpos hle
  simpa using h

In a sorted list containing 0, the head is 0.


lemma sorted_zero_head {L : List ℕ} (hsorted : L.Pairwise (· ≤ ·)) (h0 : 0 ∈ L) :
    L = 0 :: L.tail := by
  cases L with
  | nil => simp at h0
  | cons a rest =>
      have ha : a = 0 := by
        by_contra hna
        have h0rest : 0 ∈ rest := by
          simp at h0
          rcases h0 with h0 | h0
          · exact False.elim (hna h0.symm)
          · exact h0
        have ha0 : a ≤ 0 := (List.pairwise_cons.mp hsorted).1 0 h0rest
        omega
      subst a
      rfl

In a (1 + δ)-separated list starting with 0, the tail is exactly the filter of the non-zero elements.


lemma cons_zero_filter {δ : ℝ} {rest : List ℕ} (hsep : Separated δ (0 :: rest)) :
    (0 :: rest).filter (fun a => a ≠ 0) = rest := by
  have hne : ∀ z ∈ rest, z ≠ 0 := by
    intro z hz
    have hzsep : (1 + δ) * (0 : ℝ) < (z : ℝ) := by
      simpa [Separated] using List.rel_of_pairwise_cons hsep (by simpa using hz)
    intro hz0
    rw [hz0] at hzsep
    norm_num at hzsep
  rw [List.filter_cons_of_neg]
  · rw [List.filter_eq_self]
    intro z hz
    exact by simp [hne z hz]
  · simp

The positive elements of each trimmed list of APPROX-SUBSET-SUM number at most log t / log(1 + δ) + 1 (Exercise 35.5-1).

lemma approxLists_pos_length_bound {δ : ℝ} (hδ : 0 < δ) (t : ℕ) (ht : 1 ≤ t) :
    ∀ xs, (((approxLists δ t xs).filter (fun a => a ≠ 0)).length : ℝ) ≤
      Real.log (t : ℝ) / Real.log (1 + δ) + 1 := by
  intro xs
  refine sep_length_bound hδ t ht ?_ ?_ ?_
  · exact separated_filter (fun a => a ≠ 0) (approxLists_separated (le_of_lt hδ) t xs)
  · intro a ha
    exact_mod_cast (Nat.succ_le_of_lt (Nat.pos_of_ne_zero (of_decide_eq_true (List.mem_filter.mp ha).2)))
  · intro a ha
    exact_mod_cast (mem_approxLists_le_t δ t xs a (List.mem_filter.mp ha).1)

Each trimmed list of APPROX-SUBSET-SUM has at most log t / log(1 + δ) + 2 elements: the (1 + δ)-separated positive values number at most log t / log(1 + δ) + 1, and 0 contributes one more.


lemma approxLists_length_bound {δ : ℝ} (hδ : 0 < δ) (t : ℕ) (ht : 1 ≤ t) :
    ∀ xs, ((approxLists δ t xs).length : ℝ) ≤ Real.log (t : ℝ) / Real.log (1 + δ) + 2 := by
  intro xs
  have hsep : Separated δ (approxLists δ t xs) := approxLists_separated (le_of_lt hδ) t xs
  have hhead : approxLists δ t xs = 0 :: (approxLists δ t xs).tail :=
    sorted_zero_head (approxLists_sorted δ t xs) (zero_mem_approxLists δ t xs)
  have htail : (approxLists δ t xs).filter (fun a => a ≠ 0) = (approxLists δ t xs).tail := by
    rw [hhead]
    have hsep0 : Separated δ (0 :: (approxLists δ t xs).tail) := by
      rw [← hhead]
      exact hsep
    exact cons_zero_filter hsep0
  have hpos : (((approxLists δ t xs).filter (fun a => a ≠ 0)).length : ℝ) ≤
      Real.log (t : ℝ) / Real.log (1 + δ) + 1 :=
    approxLists_pos_length_bound hδ t ht xs
  rw [htail] at hpos
  rw [hhead]
  have hlen : ((0 :: (approxLists δ t xs).tail).length : ℝ) =
      ((approxLists δ t xs).tail.length : ℝ) + 1 := by
    simp
  rw [hlen]
  linarith

For 0 ≤ δ ≤ 1, δ / 2 ≤ log(1 + δ): the chord bound that turns the log t / log(1 + δ) list size into O(lg t / δ).


lemma log_one_add_ge_half {δ : ℝ} (hδ0 : 0 ≤ δ) (hδ1 : δ ≤ 1) : δ / 2 ≤ Real.log (1 + δ) := by
  have hx0 : 0 < 1 + δ := by linarith
  have hx2 : 0 < 1 - δ / 2 := by linarith
  have hle : Real.log (1 / (1 + δ)) ≤ Real.log (1 - δ / 2) := by
    apply Real.log_le_log
    · positivity
    · rw [div_le_iff₀ hx0]
      nlinarith
  have hsub : Real.log (1 - δ / 2) ≤ -(δ / 2) := by
    have h := Real.log_le_sub_one_of_pos hx2
    nlinarith
  have hloginv : Real.log (1 + δ) = -Real.log (1 / (1 + δ)) := by
    have h : Real.log (1 + δ)⁻¹ = -Real.log (1 + δ) := Real.log_inv (1 + δ)
    have hdiv : 1 / (1 + δ) = (1 + δ)⁻¹ := by ring
    rw [← hdiv] at h
    linarith
  nlinarith [hle, hsub, hloginv]

Theorem 35.8 (FPTAS running time). APPROX-SUBSET-SUM(S, t, ε) is a fully polynomial-time approximation scheme: with δ = ε/(2n) (n = |S|), every trimmed list has at most 4n · lg t / ε + 2 elements — O(n · lg t / ε) — and since MERGE-LISTS and TRIM each cost linear time in the list length, the whole algorithm runs in time O(n² · lg t / ε), polynomial in the input size and in 1/ε (CLRS §35.5, Theorem 35.8; the list-size bound is Exercise 35.5-1).


theorem approxSubsetSum_fptas {xs : List ℕ} {t : ℕ} {ε : ℝ}
    (hε0 : 0 < ε) (hε1 : ε ≤ 1) (ht : 1 ≤ t) :
    ((approxLists (ε / (2 * (xs.length : ℝ))) t xs).length : ℝ) ≤
      (4 * (xs.length : ℝ) * Real.log (t : ℝ)) / ε + 2 := by
  let δ : ℝ := ε / (2 * (xs.length : ℝ))
  by_cases hn0 : xs.length = 0
  · have hxs : xs = [] := List.eq_nil_of_length_eq_zero hn0
    subst xs
    simp [approxLists, This simp argument is unused:
  δ

Hint: Omit it from the simp argument list.
  simp [approxLists,̵ ̵δ̵]

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`δ]
  · have hnpos : 0 < (xs.length : ℝ) := by exact_mod_cast (Nat.pos_of_ne_zero hn0)
    have hn1 : (1 : ℝ) ≤ (xs.length : ℝ) := by exact_mod_cast (Nat.succ_le_of_lt (Nat.pos_of_ne_zero hn0))
    have hδ0 : 0 < δ := by
      dsimp [δ]
      exact div_pos hε0 (by positivity)
    have hδle : δ ≤ 1 := by
      dsimp [δ]
      have h2n : (1 : ℝ) ≤ 2 * (xs.length : ℝ) := by nlinarith [hn1]
      have hεlen : ε ≤ 2 * (xs.length : ℝ) := le_trans hε1 h2n
      exact (div_le_iff₀ (by positivity : (0 : ℝ) < 2 * (xs.length : ℝ))).mpr (by simpa using hεlen)
    have hlogpos : 0 < Real.log (1 + δ) := Real.log_pos (by linarith)
    have hlogt : 0 ≤ Real.log (t : ℝ) := Real.log_nonneg (by exact_mod_cast ht)
    have hlogd : δ / 2 ≤ Real.log (1 + δ) := log_one_add_ge_half (le_of_lt hδ0) hδle
    have hscaled : Real.log (t : ℝ) / Real.log (1 + δ) ≤
        4 * (xs.length : ℝ) * Real.log (t : ℝ) / ε := by
      have hδ2 : 0 < δ / 2 := by positivity
      have hinv : (Real.log (1 + δ))⁻¹ ≤ (δ / 2)⁻¹ := (inv_le_inv₀ hlogpos hδ2).2 hlogd
      have hδ2inv : (δ / 2)⁻¹ = 4 * (xs.length : ℝ) / ε := by
        have h2d : (2 : ℝ) / δ = 4 * (xs.length : ℝ) / ε := by
          dsimp [δ]
          field_simp [hε0.ne', ne_of_gt hnpos]
          ring
        calc
          (δ / 2)⁻¹ = 2 / δ := by field_simp [ne_of_gt hδ0]
          _ = 4 * (xs.length : ℝ) / ε := h2d
      calc
        Real.log (t : ℝ) / Real.log (1 + δ) = Real.log (t : ℝ) * (Real.log (1 + δ))⁻¹ := by ring
        _ ≤ Real.log (t : ℝ) * (δ / 2)⁻¹ := mul_le_mul_of_nonneg_left hinv hlogt
        _ = Real.log (t : ℝ) * (4 * (xs.length : ℝ) / ε) := by rw [hδ2inv]
        _ = 4 * (xs.length : ℝ) * Real.log (t : ℝ) / ε := by ring
    have hlen := approxLists_length_bound (δ := δ) hδ0 t ht xs
    have hfinal : Real.log (t : ℝ) / Real.log (1 + δ) + 2 ≤
        4 * (xs.length : ℝ) * Real.log (t : ℝ) / ε + 2 := by
      nlinarith [hscaled]
    have hgoal' : ((approxLists δ t xs).length : ℝ) ≤
        4 * (xs.length : ℝ) * Real.log (t : ℝ) / ε + 2 := by
      exact le_trans hlen hfinal
    simpa [δ] using hgoal'

end ApproxSubsetSumend CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.Bounds

CLRS Section 35.5 - Work bounds for costed APPROX-SUBSET-SUM

The bounds in this file concern the counters produced by the concrete execution in Costed.Execution. The model charges one unit for each list addition, comparison, or outer-loop composition step.

              noncomputable sectionnamespace CLRSnamespace ApproxSubsetSum

Every intermediate list is bounded using the one trimming parameter selected from the original input length n. In particular, the parameter is not recomputed from the length of the remaining suffix.


theorem approxLists_uniform_length_bound {n t : Nat} {ε : Real}
    (hn : 0 < n) (hε0 : 0 < ε) (hε1 : ε ≤ 1) (ht : 1 ≤ t) (ys : List Nat) :
    ((approxLists (ε / (2 * (n : Real))) t ys).length : Real) ≤
      4 * (n : Real) * Real.log (t : Real) / ε + 2 := by
  let δ : Real := ε / (2 * (n : Real))
  have hnpos : 0 < (n : Real) := by exact_mod_cast hn
  have hn1 : (1 : Real) ≤ (n : Real) := by exact_mod_cast hn
  have hδ0 : 0 < δ := by
    dsimp [δ]
    exact div_pos hε0 (by positivity)
  have hδle : δ ≤ 1 := by
    dsimp [δ]
    have h2n : (1 : Real) ≤ 2 * (n : Real) := by nlinarith [hn1]
    have hεn : ε ≤ 2 * (n : Real) := le_trans hε1 h2n
    exact (div_le_iff₀ (by positivity : (0 : Real) < 2 * (n : Real))).mpr
      (by simpa using hεn)
  have hlogpos : 0 < Real.log (1 + δ) := Real.log_pos (by linarith)
  have hlogt : 0 ≤ Real.log (t : Real) :=
    Real.log_nonneg (by exact_mod_cast ht)
  have hlogd : δ / 2 ≤ Real.log (1 + δ) :=
    log_one_add_ge_half (le_of_lt hδ0) hδle
  have hscaled : Real.log (t : Real) / Real.log (1 + δ) ≤
      4 * (n : Real) * Real.log (t : Real) / ε := by
    have hδ2 : 0 < δ / 2 := by positivity
    have hinv : (Real.log (1 + δ))⁻¹ ≤ (δ / 2)⁻¹ :=
      (inv_le_inv₀ hlogpos hδ2).2 hlogd
    have hδ2inv : (δ / 2)⁻¹ = 4 * (n : Real) / ε := by
      have h2d : (2 : Real) / δ = 4 * (n : Real) / ε := by
        dsimp [δ]
        field_simp [hε0.ne', ne_of_gt hnpos]
        ring
      calc
        (δ / 2)⁻¹ = 2 / δ := by field_simp [ne_of_gt hδ0]
        _ = 4 * (n : Real) / ε := h2d
    calc
      Real.log (t : Real) / Real.log (1 + δ) =
          Real.log (t : Real) * (Real.log (1 + δ))⁻¹ := by ring
      _ ≤ Real.log (t : Real) * (δ / 2)⁻¹ :=
        mul_le_mul_of_nonneg_left hinv hlogt
      _ = Real.log (t : Real) * (4 * (n : Real) / ε) := by rw [hδ2inv]
      _ = 4 * (n : Real) * Real.log (t : Real) / ε := by ring
  have hlen := approxLists_length_bound (δ := δ) hδ0 t ht ys
  have hfinal : Real.log (t : Real) / Real.log (1 + δ) + 2 ≤
      4 * (n : Real) * Real.log (t : Real) / ε + 2 := by
    linarith
  have hgoal : ((approxLists δ t ys).length : Real) ≤
      4 * (n : Real) * Real.log (t : Real) / ε + 2 :=
    le_trans hlen hfinal
  simpa [δ] using hgoal

The five local scans performed in one outer iteration use at most seven times the prior list length.


private theorem localPipeline_work_le (δ : Real) (t x : Nat) (L : List Nat) :
    let shifted := mapAddWithCost x L
    let merged := mergeWithCost L shifted.value
    let trimmed := trimWithCost δ merged.value
    let kept := filterAtMostWithCost t trimmed.value
    shifted.work + merged.work + trimmed.work + kept.work ≤ 7 * L.length := by
  dsimp only
  have hshiftWork := mapAddWithCost_work x L
  have hshiftLen : (mapAddWithCost x L).value.length = L.length := by
    rw [mapAddWithCost_value]
    simp
  have hmergeWork := mergeWithCost_work_le L (mapAddWithCost x L).value
  have hmergeLen := mergeWithCost_length L (mapAddWithCost x L).value
  have htrimWork :=
    trimWithCost_work_le δ (mergeWithCost L (mapAddWithCost x L).value).value
  have htrimLen :=
    trimWithCost_length_le δ (mergeWithCost L (mapAddWithCost x L).value).value
  have hfilterWork := filterAtMostWithCost_work t
    (trimWithCost δ (mergeWithCost L (mapAddWithCost x L).value).value).value
  omega

If all semantic intermediate lists are bounded by B, the actual outer execution counter is at most length * (7B + 1).


theorem approxListsWithCost_work_le_of_length_bound {δ B : Real} {t : Nat}
    (_hB : 0 ≤ B)
    (hbound : ∀ ys, ((approxLists δ t ys).length : Real) ≤ B) (xs : List Nat) :
    ((approxListsWithCost δ t xs).work : Real) ≤
      (xs.length : Real) * (7 * B + 1) := by
  induction xs with
  | nil => simp [approxListsWithCost]
  | cons x xs ih =>
      let prior := approxListsWithCost δ t xs
      let shifted := mapAddWithCost x prior.value
      let merged := mergeWithCost prior.value shifted.value
      let trimmed := trimWithCost δ merged.value
      let kept := filterAtMostWithCost t trimmed.value
      have hprior : (prior.work : Real) ≤ (xs.length : Real) * (7 * B + 1) := by
        simpa [prior] using ih
      have hlen : (prior.value.length : Real) ≤ B := by
        rw [approxListsWithCost_value]
        exact hbound xs
      have hpipeline :
          shifted.work + merged.work + trimmed.work + kept.work ≤
            7 * prior.value.length := by
        simpa [shifted, merged, trimmed, kept] using
          localPipeline_work_le δ t x prior.value
      have hpipelineReal :
          ((shifted.work + merged.work + trimmed.work + kept.work : Nat) : Real) ≤
            7 * (prior.value.length : Real) := by
        exact_mod_cast hpipeline
      push_cast at hpipelineReal
      have hstep :
          ((shifted.work + merged.work + trimmed.work + kept.work + 1 : Nat) : Real) ≤
            7 * B + 1 := by
        push_cast
        nlinarith
      rw [approxListsWithCost]
      dsimp only
      change ((prior.work + shifted.work + merged.work + trimmed.work + kept.work + 1 : Nat) : Real) ≤
        (((xs.length + 1 : Nat) : Real) * (7 * B + 1))
      push_cast
      push_cast at hstep
      ring_nf at hprior hstep ⊢
      nlinarith

The counter returned by the complete execution is polynomial in the input length, log t, and 1 / ε. The bound includes construction of every intermediate list and the final maximum scan.


theorem approxSubsetSumWithCost_work_le {xs : List Nat} {t : Nat} {ε : Real}
    (hε0 : 0 < ε) (hε1 : ε ≤ 1) (ht : 1 ≤ t) :
    ((approxSubsetSumWithCost xs t ε).work : Real) ≤
      48 * ((xs.length : Real) + 1) ^ 2 *
        (Real.log (t : Real) + 1) / ε := by
  have hlog : 0 ≤ Real.log (t : Real) :=
    Real.log_nonneg (by exact_mod_cast ht)
  by_cases hnil : xs = []
  · subst xs
    have hratio : 1 ≤ (Real.log (t : Real) + 1) / ε := by
      apply (le_div_iff₀ hε0).2
      linarith
    have h48 : (0 : Real) ≤ 48 := by norm_num
    simp only [List.length_nil, Nat.cast_zero, zero_add, one_pow]
    rw [approxSubsetSumWithCost]
    simp [approxListsWithCost, maximumWithCost, maximumAuxWithCost]
    calc
      (1 : Real) ≤ 48 := by norm_num
      _ ≤ 48 * ((Real.log (t : Real) + 1) / ε) :=
        by simpa using mul_le_mul_of_nonneg_left hratio h48
      _ = 48 * (Real.log (t : Real) + 1) / ε := by ring
  · have hn : 0 < xs.length := by
      exact Nat.pos_of_ne_zero (fun hlen => hnil (List.eq_nil_of_length_eq_zero hlen))
    let n : Real := xs.length
    let B : Real := 4 * n * Real.log (t : Real) / ε + 2
    have hn0 : 0 ≤ n := by positivity
    have hn1 : 1 ≤ n := by
      dsimp [n]
      exact_mod_cast hn
    have hB0 : 0 ≤ B := by
      dsimp [B]
      have : 0 ≤ 4 * n * Real.log (t : Real) / ε := by positivity
      linarith
    have hbound : ∀ ys,
        ((approxLists (ε / (2 * (xs.length : Real))) t ys).length : Real) ≤ B := by
      intro ys
      simpa [B, n] using
        (approxLists_uniform_length_bound (n := xs.length) hn hε0 hε1 ht ys)
    have hlists := approxListsWithCost_work_le_of_length_bound
      (δ := ε / (2 * (xs.length : Real))) (B := B) hB0 hbound xs
    have hvalueLen :
        ((approxListsWithCost (ε / (2 * (xs.length : Real))) t xs).value.length : Real) ≤ B := by
      rw [approxListsWithCost_value]
      exact hbound xs
    have htotal :
        ((approxSubsetSumWithCost xs t ε).work : Real) ≤
          n * (7 * B + 1) + B := by
      rw [approxSubsetSumWithCost]
      simp only
      rw [maximumWithCost_work]
      push_cast
      simpa [n] using add_le_add hlists hvalueLen
    have hBstep : B ≤ 7 * B + 1 := by nlinarith
    have hcoarse :
        ((approxSubsetSumWithCost xs t ε).work : Real) ≤
          (n + 1) * (7 * B + 1) := by
      calc
        ((approxSubsetSumWithCost xs t ε).work : Real)
            ≤ n * (7 * B + 1) + B := htotal
        _ ≤ n * (7 * B + 1) + (7 * B + 1) := add_le_add_right hBstep _
        _ = (n + 1) * (7 * B + 1) := by ring
    have hnlog : n * Real.log (t : Real) ≤
        (n + 1) * (Real.log (t : Real) + 1) := by
      nlinarith [mul_nonneg hn0 hlog]
    have hεprod : ε ≤ (n + 1) * (Real.log (t : Real) + 1) := by
      have hone : 1 ≤ (n + 1) * (Real.log (t : Real) + 1) := by
        nlinarith [mul_nonneg hn0 hlog]
      exact le_trans hε1 hone
    have hinside :
        28 * n * Real.log (t : Real) + 15 * ε ≤
          48 * (n + 1) * (Real.log (t : Real) + 1) := by
      nlinarith
    have hpoly : (n + 1) * (7 * B + 1) ≤
        48 * (n + 1) ^ 2 * (Real.log (t : Real) + 1) / ε := by
      have hleft : (n + 1) * (7 * B + 1) =
          ((n + 1) * (28 * n * Real.log (t : Real) + 15 * ε)) / ε := by
        dsimp [B]
        field_simp [hε0.ne']
        ring
      rw [hleft]
      apply (div_le_div_iff_of_pos_right hε0).2
      have hmul := mul_le_mul_of_nonneg_left hinside (by linarith : 0 ≤ n + 1)
      nlinarith
    simpa [n] using le_trans hcoarse hpoly

Kernel-checked FPTAS bundle for the costed execution: feasibility, approximation quality, and polynomial work all refer to the same run.


theorem approxSubsetSumWithCost_fptas {xs : List Nat} {t : Nat} {ε : Real}
    (hε0 : 0 < ε) (hε1 : ε ≤ 1) (ht : 1 ≤ t) :
    let run := approxSubsetSumWithCost xs t ε
    run.value ∈ subsetSums xs ∧
      run.value ≤ t ∧
      (optimalSum xs t : Real) ≤ (1 + ε) * (run.value : Real) ∧
      (run.work : Real) ≤
        48 * ((xs.length : Real) + 1) ^ 2 *
          (Real.log (t : Real) + 1) / ε := by
  dsimp only
  refine ⟨?_, ?_, ?_, approxSubsetSumWithCost_work_le hε0 hε1 ht⟩
  · rw [approxSubsetSumWithCost_value]
    exact approxSum_mem_subsetSums xs t ε
  · rw [approxSubsetSumWithCost_value]
    exact approxSum_le_t xs t ε
  · rw [approxSubsetSumWithCost_value]
    exact approxSubsetSum_approx_lt xs t ε hε0 hε1

end ApproxSubsetSumend CLRS

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.Definitions

CLRS Section 35.5 - Costed local scans

Execution records for the concrete list passes used by APPROX-SUBSET-SUM. Each counter is produced by the same recursion as its returned value.

              noncomputable sectionnamespace CLRSnamespace ApproxSubsetSum

Result and unit-operation count of a list-producing scan.

structure ListExecution where
  value : List Nat
  work : Nat
deriving Repr

Result and unit-operation count of a natural-number-producing scan.

structure NatExecution where
  value : Nat
  work : Nat
deriving Repr

Add x to every element, charging one addition per element.

def mapAddWithCost (x : Nat) : List Nat → ListExecution
  | [] => ⟨[], 0⟩
  | y :: ys =>
      let rest := mapAddWithCost x ys
      ⟨(y + x) :: rest.value, rest.work + 1⟩

CLRS MERGE-LISTS, charging one comparison whenever both inputs are nonempty.

def mergeWithCost : (L M : List Nat) → ListExecution
  | [], ys => ⟨ys, 0⟩
  | xs, [] => ⟨xs, 0⟩
  | x :: xs, y :: ys =>
      if x ≤ y then
        let rest := mergeWithCost xs (y :: ys)
        ⟨x :: rest.value, rest.work + 1⟩
      else
        let rest := mergeWithCost (x :: xs) ys
        ⟨y :: rest.value, rest.work + 1⟩
termination_by L M => L.length + M.length

Tail scan of TRIM, charging one threshold comparison per scanned value.

def trimAuxWithCost (δ : Real) : Nat → List Nat → ListExecution
  | _last, [] => ⟨[], 0⟩
  | last, y :: ys =>
      if (1 + δ) * (last : Real) < (y : Real) then
        let rest := trimAuxWithCost δ y ys
        ⟨y :: rest.value, rest.work + 1⟩
      else
        let rest := trimAuxWithCost δ last ys
        ⟨rest.value, rest.work + 1⟩

CLRS TRIM with a counter for its tail comparisons.

def trimWithCost (δ : Real) : List Nat → ListExecution
  | [] => ⟨[], 0⟩
  | y :: ys =>
      let rest := trimAuxWithCost δ y ys
      ⟨y :: rest.value, rest.work⟩

Keep values at most t, charging one target comparison per value.

def filterAtMostWithCost (t : Nat) : List Nat → ListExecution
  | [] => ⟨[], 0⟩
  | y :: ys =>
      let rest := filterAtMostWithCost t ys
      if y ≤ t then
        ⟨y :: rest.value, rest.work + 1⟩
      else
        ⟨rest.value, rest.work + 1⟩

Maximum scan with an explicit accumulator and one comparison per value.

def maximumAuxWithCost (best : Nat) : List Nat → NatExecution
  | [] => ⟨best, 0⟩
  | y :: ys =>
      let rest := maximumAuxWithCost (max best y) ys
      ⟨rest.value, rest.work + 1⟩

Maximum of a natural-number list, using 0 for the empty case.

def maximumWithCost (xs : List Nat) : NatExecution :=
  maximumAuxWithCost 0 xs

end ApproxSubsetSumend CLRS

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.Execution

CLRS Section 35.5 - Costed APPROX-SUBSET-SUM execution

The outer recursion composes the costed map, merge, trim, and filter scans. The final wrapper scans the resulting list for its maximum.

              noncomputable sectionnamespace CLRSnamespace ApproxSubsetSum

Costed construction of the trimmed lists. The stored work is exactly the sum of the recursively executed local scan counters plus one outer-loop unit.

def approxListsWithCost (δ : Real) (t : Nat) : List Nat → ListExecution
  | [] => ⟨[0], 0⟩
  | x :: xs =>
      let prior := approxListsWithCost δ t xs
      let shifted := mapAddWithCost x prior.value
      let merged := mergeWithCost prior.value shifted.value
      let trimmed := trimWithCost δ merged.value
      let kept := filterAtMostWithCost t trimmed.value
      ⟨kept.value,
        prior.work + shifted.work + merged.work + trimmed.work + kept.work + 1⟩

Erasing the outer counter gives the existing semantic trimmed-list construction.

theorem approxListsWithCost_value (δ : Real) (t : Nat) (xs : List Nat) :
    (approxListsWithCost δ t xs).value = approxLists δ t xs := by
  induction xs with
  | nil => simp [approxListsWithCost, approxLists]
  | cons x xs ih =>
      simp [approxListsWithCost, approxLists, ih, mapAddWithCost_value,
        mergeWithCost_value, trimWithCost_value, filterAtMostWithCost_value]

A fold by max never drops its initial accumulator.

theorem le_foldl_max (best : Nat) (L : List Nat) :
    best ≤ L.foldl max best := by
  induction L generalizing best with
  | nil => simp
  | cons y ys ih =>
      exact (Nat.le_max_left best y).trans (ih (max best y))

Every list member is bounded by a fold by max.

theorem mem_le_foldl_max (best : Nat) {L : List Nat} {x : Nat}
    (hx : x ∈ L) : x ≤ L.foldl max best := by
  induction L generalizing best with
  | nil => simp at hx
  | cons y ys ih =>
      simp only [List.mem_cons] at hx
      simp only [List.foldl_cons]
      rcases hx with hxy | hx
      · subst x
        exact (Nat.le_max_right best y).trans (le_foldl_max (max best y) ys)
      · exact ih (max best y) hx

A fold by max is below every common upper bound for its accumulator and list elements.

theorem foldl_max_le (L : List Nat) (best bound : Nat)
    (hbest : best ≤ bound) (hmem : ∀ x ∈ L, x ≤ bound) :
    L.foldl max best ≤ bound := by
  induction L generalizing best with
  | nil => simpa using hbest
  | cons y ys ih =>
      simp only [List.foldl_cons]
      apply ih (max best y)
      · exact max_le hbest (hmem y (by simp))
      · intro x hx
        exact hmem x (by simp [hx])

On a list containing 0, the maximum scan agrees with the nonempty Finset maximum used by approxSum.


theorem foldl_max_eq_toFinset_max' {L : List Nat} (h0 : 0 ∈ L) :
    L.foldl max 0 = L.toFinset.max' ⟨0, by simpa using h0⟩ := by
  apply le_antisymm
  · apply foldl_max_le
    · exact L.toFinset.le_max' 0 (by simpa using h0)
    · intro x hx
      exact L.toFinset.le_max' x (by simpa using hx)
  · have hmax : L.toFinset.max' ⟨0, by simpa using h0⟩ ∈ L := by
      simpa using L.toFinset.max'_mem ⟨0, by simpa using h0⟩
    exact mem_le_foldl_max 0 hmax

Costed APPROX-SUBSET-SUM, including the final maximum scan.

def approxSubsetSumWithCost (xs : List Nat) (t : Nat) (ε : Real) : NatExecution :=
  let lists := approxListsWithCost (ε / (2 * (xs.length : Real))) t xs
  let answer := maximumWithCost lists.value
  ⟨answer.value, lists.work + answer.work⟩

Erasing the complete execution counter gives the existing approxSum result.


theorem approxSubsetSumWithCost_value (xs : List Nat) (t : Nat) (ε : Real) :
    (approxSubsetSumWithCost xs t ε).value = approxSum xs t ε := by
  rw [approxSubsetSumWithCost]
  simp only
  rw [maximumWithCost_value, approxListsWithCost_value]
  unfold approxSum
  apply foldl_max_eq_toFinset_max'
  exact zero_mem_approxLists _ _ _

end ApproxSubsetSumend CLRS

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.LocalCorrectness

CLRS Section 35.5 - Local scan correctness and work

Erasure and linear counter bounds for the concrete list scans used by the costed APPROX-SUBSET-SUM execution.

              noncomputable sectionnamespace CLRSnamespace ApproxSubsetSumtheorem mapAddWithCost_value (x : Nat) (L : List Nat) :
    (mapAddWithCost x L).value = L.map (fun y => y + x) := by
  induction L with
  | nil => simp [mapAddWithCost]
  | cons y ys ih => simp [mapAddWithCost, ih]theorem mapAddWithCost_work (x : Nat) (L : List Nat) :
    (mapAddWithCost x L).work = L.length := by
  induction L with
  | nil => simp [mapAddWithCost]
  | cons y ys ih => simp [mapAddWithCost, ih]theorem mergeWithCost_value (L M : List Nat) :
    (mergeWithCost L M).value = merge L M := by
  induction L generalizing M with
  | nil => simp [mergeWithCost, merge]
  | cons x xs ihL =>
      induction M with
      | nil => simp [mergeWithCost, merge]
      | cons y ys ihM =>
          by_cases hxy : x ≤ y
          · simp [mergeWithCost, merge, hxy, ihL]
          · simp [mergeWithCost, merge, hxy, ihM]

theorem mergeWithCost_length (L M : List Nat) :
    (mergeWithCost L M).value.length = L.length + M.length := by
  induction L generalizing M with
  | nil => simp [mergeWithCost]
  | cons x xs ihL =>
      induction M with
      | nil => simp [mergeWithCost]
      | cons y ys ihM =>
          by_cases hxy : x ≤ y
          · simp only [mergeWithCost, hxy, ↓reduceIte, List.length_cons]
            rw [ihL]
            simp only [List.length_cons]
            omega
          · simp only [mergeWithCost, hxy, ↓reduceIte, List.length_cons]
            rw [ihM]
            simp only [List.length_cons]
            omegatheorem mergeWithCost_work_le (L M : List Nat) :
    (mergeWithCost L M).work ≤ L.length + M.length := by
  induction L generalizing M with
  | nil => simp [mergeWithCost]
  | cons x xs ihL =>
      induction M with
      | nil => simp [mergeWithCost]
      | cons y ys ihM =>
          by_cases hxy : x ≤ y
          · simp only [mergeWithCost, hxy, ↓reduceIte, List.length_cons]
            have h := ihL (y :: ys)
            simp only [List.length_cons] at h
            omega
          · simp only [mergeWithCost, hxy, ↓reduceIte, List.length_cons]
            have h := ihM
            simp only [List.length_cons] at h
            omegatheorem trimAuxWithCost_value (δ : Real) (last : Nat) (ys : List Nat) :
    (trimAuxWithCost δ last ys).value = trimAux δ last ys := by
  induction ys generalizing last with
  | nil => simp [trimAuxWithCost, trimAux]
  | cons y ys ih =>
      by_cases hkeep : (1 + δ) * (last : Real) < (y : Real)
      · simp [trimAuxWithCost, trimAux, hkeep, ih]
      · simp [trimAuxWithCost, trimAux, hkeep, ih]theorem trimAuxWithCost_work (δ : Real) (last : Nat) (ys : List Nat) :
    (trimAuxWithCost δ last ys).work = ys.length := by
  induction ys generalizing last with
  | nil => simp [trimAuxWithCost]
  | cons y ys ih =>
      by_cases hkeep : (1 + δ) * (last : Real) < (y : Real)
      · simp [trimAuxWithCost, hkeep, ih]
      · simp [trimAuxWithCost, hkeep, ih]theorem trimWithCost_value (δ : Real) (L : List Nat) :
    (trimWithCost δ L).value = trim δ L := by
  cases L with
  | nil => simp [trimWithCost, trim]
  | cons y ys => simp [trimWithCost, trim, trimAuxWithCost_value]theorem trimWithCost_work_le (δ : Real) (L : List Nat) :
    (trimWithCost δ L).work ≤ L.length := by
  cases L with
  | nil => simp [trimWithCost]
  | cons y ys => simp [trimWithCost, trimAuxWithCost_work]

theorem trimWithCost_length_le (δ : Real) (L : List Nat) :
    (trimWithCost δ L).value.length ≤ L.length := by
  rw [trimWithCost_value]
  exact (trim_sublist δ L).length_letheorem filterAtMostWithCost_value (t : Nat) (L : List Nat) :
    (filterAtMostWithCost t L).value = L.filter (fun y => y ≤ t) := by
  induction L with
  | nil => simp [filterAtMostWithCost]
  | cons y ys ih =>
      by_cases hy : y ≤ t
      · simp [filterAtMostWithCost, hy, ih]
      · simp [filterAtMostWithCost, hy, ih]theorem filterAtMostWithCost_work (t : Nat) (L : List Nat) :
    (filterAtMostWithCost t L).work = L.length := by
  induction L with
  | nil => simp [filterAtMostWithCost]
  | cons y ys ih =>
      by_cases hy : y ≤ t
      · simp [filterAtMostWithCost, hy, ih]
      · simp [filterAtMostWithCost, hy, ih]

theorem filterAtMostWithCost_length_le (t : Nat) (L : List Nat) :
    (filterAtMostWithCost t L).value.length ≤ L.length := by
  rw [filterAtMostWithCost_value]
  exact List.length_filter_le (fun y => y ≤ t) Ltheorem maximumAuxWithCost_value (best : Nat) (L : List Nat) :
    (maximumAuxWithCost best L).value = L.foldl max best := by
  induction L generalizing best with
  | nil => simp [maximumAuxWithCost]
  | cons y ys ih => simp [maximumAuxWithCost, ih]theorem maximumAuxWithCost_work (best : Nat) (L : List Nat) :
    (maximumAuxWithCost best L).work = L.length := by
  induction L generalizing best with
  | nil => simp [maximumAuxWithCost]
  | cons y ys ih => simp [maximumAuxWithCost, ih]theorem maximumWithCost_value (L : List Nat) :
    (maximumWithCost L).value = L.foldl max 0 := by
  exact maximumAuxWithCost_value 0 Ltheorem maximumWithCost_work (L : List Nat) :
    (maximumWithCost L).work = L.length := by
  exact maximumAuxWithCost_work 0 Lend ApproxSubsetSumend CLRS

35.5. The Subset-Sum Problem

The subset-sum problem and EXACT-SUBSET-SUM

The merge of two sorted lists

TRIM and Lemma 35.5

APPROX-SUBSET-SUM and its approximation guarantee

The (1 + ε) bound (Theorem 35.7, final form)

The running-time analysis (Theorem 35.8, FPTAS)

Definitions and proofs

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.Bounds

CLRS Section 35.5 - Work bounds for costed APPROX-SUBSET-SUM

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.Definitions

CLRS Section 35.5 - Costed local scans

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.Execution

CLRS Section 35.5 - Costed APPROX-SUBSET-SUM execution

CLRSLean.FourthEdition.Chapter_35.Section_35_5_The_Subset_Sum_Problem.Costed.LocalCorrectness

CLRS Section 35.5 - Local scan correctness and work

The `(1 + ε)` bound (Theorem 35.7, final form)