Imports

import Mathlib

35.5 The Subset-Sum Problem

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