-- Module: Mathlib.NumberTheory.BernoulliPolynomials
Goal:
A : Type u_1
inst : CommRing A
inst_1 : Algebra ℚ A
t : A
n : ℕ
h :
  Eq (α := A)
    ((PowerSeries.coeff 0 : (a : A⟦X⟧) → A)
        ((PowerSeries.mk fun (n : ℕ) ↦ (aeval t : (a : ℚ[X]) → A) (HDiv.hDiv (α := ℚ) 1 (↑n ! : ℚ) • bernoulli n)) *
          (exp A - 1)) :
      A)
    ((PowerSeries.coeff 0 : (a : A⟦X⟧) → A)
        (HMul.hMul (α := A⟦X⟧) (β := A⟦X⟧) PowerSeries.X ((rescale t : (a : A⟦X⟧) → A⟦X⟧) (exp A) : A⟦X⟧)) :
      A)
h₁ :
  ∀ (n : ℕ),
    Eq (α := A)
      ((PowerSeries.coeff (n + 1) : (a : A⟦X⟧) → A)
          ((PowerSeries.mk fun (n : ℕ) ↦ (aeval t : (a : ℚ[X]) → A) (HDiv.hDiv (α := ℚ) 1 (↑n ! : ℚ) • bernoulli n)) *
            (exp A - 1)) :
        A)
      ((PowerSeries.coeff (n + 1) : (a : A⟦X⟧) → A)
          (HMul.hMul (α := A⟦X⟧) (β := A⟦X⟧) PowerSeries.X ((rescale t : (a : A⟦X⟧) → A⟦X⟧) (exp A) : A⟦X⟧)) :
        A)
⊢ Eq (α := A)
    ((PowerSeries.coeff n : (a : A⟦X⟧) → A)
        ((PowerSeries.mk fun (n : ℕ) ↦ (aeval t : (a : ℚ[X]) → A) (HDiv.hDiv (α := ℚ) 1 (↑n ! : ℚ) • bernoulli n)) *
          (exp A - 1)) :
      A)
    ((PowerSeries.coeff n : (a : A⟦X⟧) → A)
        (HMul.hMul (α := A⟦X⟧) (β := A⟦X⟧) PowerSeries.X ((rescale t : (a : A⟦X⟧) → A⟦X⟧) (exp A) : A⟦X⟧)) :
      A)

Tactic:
cases n with | zero => simp | succ n =>
  -- n ≥ 1, the coefficients is a sum to n+2, so use `sum_range_succ` to write as
  -- last term plus sum to n+1
  rw [coeff_succ_X_mul, coeff_rescale, coeff_exp, PowerSeries.coeff_mul,
    Nat.sum_antidiagonal_eq_sum_range_succ_mk, sum_range_succ]
  -- last term is zero so kill with `add_zero`
  simp only [map_sub, tsub_self, constantCoeff_one, constantCoeff_exp,
    coeff_zero_eq_constantCoeff, mul_zero, sub_self, add_zero]
  -- Let's multiply both sides by (n+1)! (OK because it's a unit)
  have hnp1 : IsUnit ((n + 1)! : ℚ) := IsUnit.mk0 _ (mod_cast factorial_ne_zero (n + 1))
  rw [← (hnp1.map (algebraMap ℚ A)).mul_right_inj]
  -- do trivial rearrangements to make RHS (n+1)*t^n
  rw [mul_left_comm, ← map_mul]
  change _ = t ^ n * algebraMap ℚ A (((n + 1) * n ! : ℕ) * (1 / n !))
  rw [cast_mul, mul_assoc,
    mul_one_div_cancel (show (n ! : ℚ) ≠ 0 from cast_ne_zero.2 (factorial_ne_zero n)), mul_one,
    mul_comm (t ^ n), ← aeval_monomial, cast_add, cast_one]
  -- But this is the RHS of `Polynomial.sum_bernoulli`
  rw [← sum_bernoulli, Finset.mul_sum, map_sum]
  -- and now we have to prove a sum is a sum, but all the terms are equal.
  apply Finset.sum_congr rfl
  -- The rest is just trivialities, hampered by the fact that we're coercing
  -- factorials and binomial coefficients between ℕ and ℚ and A.
  intro i hi
  -- deal with coefficients of e^X-1
  simp only [Nat.cast_choose ℚ (mem_range_le hi), coeff_mk, if_neg (mem_range_sub_ne_zero hi),
    PowerSeries.coeff_one, coeff_exp, sub_zero, Algebra.smul_def,
    mul_right_comm _ ((aeval t) _), ← mul_assoc, ← map_mul, ← Polynomial.C_eq_algebraMap,
    Polynomial.aeval_mul, Polynomial.aeval_C]
  -- finally cancel the Bernoulli polynomial and the algebra_map
  field_simp

Premises used: Iff.mpr, Finset.mul_sum, PowerSeries.coeff_rescale, Polynomial.sum_bernoulli, HMul.hMul, map_sub, Algebra.algebraMap, sub_self, sub_zero, Finset.Nat.sum_antidiagonal_eq_sum_range_succ_mk, Rat, mul_assoc, tsub_self, PowerSeries.constantCoeff_one, HDiv.hDiv, PowerSeries.coeff_mk, Polynomial.C_eq_algebraMap, mul_left_comm, PowerSeries.coeff_exp, map_sum, Polynomial.aeval_mul, Nat.cast, mul_one_div_cancel, MulZeroClass.mul_zero, Nat.cast_add, PowerSeries.coeff_mul, map_mul, Nat.factorial_ne_zero, Nat.cast_mul, HPow.hPow, HAdd.hAdd, Nat.factorial, PowerSeries.coeff_one, mul_right_comm, Polynomial.aeval_C, PowerSeries.coeff_zero_eq_constantCoeff, mul_comm, PowerSeries.constantCoeff_exp, Nat.zero, Polynomial.aeval_monomial, Algebra.smul_def, PowerSeries.coeff_succ_X_mul, Polynomial.aeval, mul_one, add_zero, OfNat.ofNat, Nat.succ, Nat.cast_one, DFunLike.coe, Finset.sum_range_succ, Finset.sum_congr, rfl, Nat.cast_ne_zero