diff --git a/thys/Hyperdual/AnalyticTestFunction.thy b/thys/Hyperdual/AnalyticTestFunction.thy
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/AnalyticTestFunction.thy
@@ -0,0 +1,441 @@
+(*  Title:   AnalyticTestFunction.thy
+    Authors: Filip Smola and Jacques D. Fleuriot, University of Edinburgh, 2019-2021
+*)
+
+theory AnalyticTestFunction
+  imports HyperdualFunctionExtension "HOL-Decision_Procs.Approximation"
+begin
+
+subsection\<open>Analytic Test Function\<close>
+text\<open>
+  We investigate the analytic test function used by Fike and Alonso in their
+  2011 paper~@{cite "fike_alonso-2011"} as a relatively non-trivial example.
+  The function is defined as: @{term "\<lambda>x. exp x / (sqrt (sin x ^ 3 + cos x ^ 3))"}.
+\<close>
+definition fa_test :: "real \<Rightarrow> real"
+  where "fa_test x = exp x / (sqrt (sin x ^ 3 + cos x ^ 3))"
+
+text\<open>
+  We define the same composition of functions but using the relevant hyperdual versions.
+  Note that we implicitly use the facts that hyperdual extensions of plus, times and inverse are the
+  same operations on hyperduals.
+\<close>
+definition hyp_fa_test :: "real hyperdual \<Rightarrow> real hyperdual"
+  where "hyp_fa_test x = ((*h* exp) x) / ((*h* sqrt) (((*h* sin) x) ^ 3 + ((*h* cos) x) ^ 3))"
+
+text\<open>We prove lemmas useful to show when this function is well defined.\<close>
+lemma sin_cube_plus_cos_cube:
+  "sin x ^ 3 + cos x ^ 3 = 1/2 * (sin x + cos x) * (2 - sin (2 * x))"
+  for x :: "'a::{real_normed_field,banach}"
+proof -
+  have "sin x ^ 3 + cos x ^ 3 = (sin x + cos x) * (cos x ^ 2 - cos x * sin x + sin x ^ 2)"
+    by (smt (z3) add.commute combine_common_factor diff_add_cancel distrib_left
+          mult.commute mult.left_commute power2_eq_square power3_eq_cube)
+  also have "\<dots> = (sin x + cos x) * (1 - cos x * sin x)"
+    by simp
+  finally show ?thesis
+    by (smt (z3) mult.commute mult.left_neutral mult_2 nonzero_mult_div_cancel_left
+          right_diff_distrib' sin_add times_divide_eq_left zero_neq_numeral)
+qed
+
+lemma sin_cube_plus_cos_cube_gt_zero_iff:
+  "(sin x ^ 3 + cos x ^ 3 > 0) = (sin x + cos x > 0)"
+  for x ::real
+  by (smt (verit, best) cos_zero power3_eq_cube power_zero_numeral sin_cube_plus_cos_cube
+       sin_le_one sin_zero zero_less_mult_iff)
+
+lemma sin_plus_cos_eq_45:
+  "sin x + cos x = sqrt 2 * sin (x + pi/4)"
+  apply (simp add: sin_add sin_45 cos_45 )
+  by (simp add: field_simps)
+
+lemma sin_cube_plus_cos_cube_gt_zero_iff':
+  "(sin x ^ 3 + cos x ^ 3 > 0) = (sin (x + pi/4) > 0)"
+  by (smt (verit, best) mult_pos_pos real_sqrt_gt_0_iff
+      sin_cube_plus_cos_cube_gt_zero_iff sin_plus_cos_eq_45 zero_less_mult_pos)
+
+lemma sin_less_zero_pi:
+  "\<lbrakk>- pi < x; x < 0\<rbrakk> \<Longrightarrow> sin x < 0"
+  by (metis add.inverse_inverse add.inverse_neutral neg_less_iff_less sin_gt_zero sin_minus)
+
+lemma sin_45_positive_intervals:
+  "(sin (x + pi/4) > 0) = (x \<in> (\<Union>n::int. {-pi/4 + 2*pi*n <..< 3*pi/4 + 2*pi*n}))"
+proof (standard ; (elim UnionE rangeE | -))
+  obtain y :: real and n :: int
+    where "x = y + 2*pi*n" and "- pi \<le> y" and "y \<le> pi"
+    using sincos_principal_value sin_cos_eq_iff less_eq_real_def by metis
+  note yn = this
+
+  assume "0 < sin (x + pi / 4)"
+  then have a: "0 < sin (y + pi / 4)"
+    using yn by (metis sin_add sin_cos_eq_iff)
+  then have "y \<in> {- pi / 4<..<3 * pi / 4}"
+  proof (unfold greaterThanLessThan_iff, safe)
+    show "- pi / 4 < y"
+    proof (rule ccontr)
+      assume "\<not> - pi / 4 < y"
+      then have "y < - pi / 4 \<or> y = - pi / 4"
+        by (simp add: not_less le_less)
+      then show False
+        using a sin_less_zero_pi[where x = "y + pi/4"] yn sin_zero
+        by force
+    qed
+    show "y < 3 * pi / 4"
+    proof (rule ccontr)
+      assume "\<not> y < 3 * pi / 4"
+      then have "pi \<le> y + pi/4"
+        by (simp add: not_less)
+      then show False
+        using a sin_le_zero yn pi_ge_two by fastforce
+    qed
+  qed
+  then have "x \<in> {- pi / 4 + 2*pi*n<..<3 * pi / 4 + 2*pi*n}"
+    using yn greaterThanLessThan_iff by simp
+  then show "x \<in> (\<Union>n::int. {- pi / 4 + 2*pi*n<..<3 * pi / 4 + 2*pi*n})"
+    by blast
+next
+  fix X and n :: int
+  assume "x \<in> X"
+     and "X = {- pi / 4 + 2*pi*n<..<3 * pi / 4 + 2*pi*n}"
+  then have "x \<in> {- pi / 4 + 2*pi*n<..<3 * pi / 4 + 2*pi*n}"
+    by simp
+  then obtain y :: real and n :: int
+    where "x = y + 2*pi*n" and "- pi / 4 < y" and "y < 3 * pi / 4"
+    by (smt (z3) greaterThanLessThan_iff)
+  note yn = this
+
+  have "0 < sin (y + pi / 4)"
+    using sin_gt_zero yn by force
+  then show "0 < sin (x + pi / 4)"
+    using yn sin_cos_eq_iff[of "x + pi / 4" "y + pi / 4"] by simp
+qed
+
+text\<open>When the function is well defined our hyperdual definition is equal to the hyperdual extension.\<close>
+lemma hypext_fa_test:
+  assumes "Base x \<in> (\<Union>n::int. {-pi/4 + 2*pi*n <..< 3*pi/4 + 2*pi*n})"
+  shows "(*h* fa_test) x = hyp_fa_test x"
+proof -
+  have inverse_sqrt_valid: "0 < sin (Base x) ^ 3 + cos (Base x) ^ 3"
+    using assms sin_45_positive_intervals sin_cube_plus_cos_cube_gt_zero_iff' by force
+
+  have "\<And>f. (\<lambda>x. (sin x) ^ 3) twice_field_differentiable_at Base x"
+   and "\<And>f. (\<lambda>x. (cos x) ^ 3) twice_field_differentiable_at Base x"
+    by (simp_all add: twice_field_differentiable_at_compose[OF _ twice_field_differentiable_at_power])
+  then have "(*h* (\<lambda>x. sin x ^ 3 + cos x ^ 3)) x = (*h* sin) x ^ 3 + (*h* cos) x ^ 3"
+        and d_sincos: "(\<lambda>x. sin x ^ 3 + cos x ^ 3) twice_field_differentiable_at Base x"
+    using hypext_fun_add[of "\<lambda>x. sin x ^ 3" x "\<lambda>x. cos x ^ 3"]
+    by (simp_all add: hypext_fun_power twice_field_differentiable_at_add)
+  then have "(*h* (\<lambda>x. sqrt (sin x ^ 3 + cos x ^ 3))) x = (*h* sqrt) ((*h* sin) x ^ 3 + (*h* cos) x ^ 3)"
+    using inverse_sqrt_valid hypext_compose[of "\<lambda>x. sin x ^ 3 + cos x ^ 3" x] by simp
+  moreover have d_sqrt: "(\<lambda>x. sqrt (sin x ^ 3 + cos x ^ 3)) twice_field_differentiable_at Base x"
+    using inverse_sqrt_valid d_sincos twice_field_differentiable_at_compose twice_field_differentiable_at_sqrt
+    by blast
+  ultimately have "(*h* (\<lambda>x. inverse (sqrt (sin x ^ 3 + cos x ^ 3)))) x = inverse ((*h* sqrt) ((*h* sin) x ^ 3 + (*h* cos) x ^ 3))"
+    using inverse_sqrt_valid hypext_fun_inverse[of "\<lambda>x. sqrt (sin x ^ 3 + cos x ^ 3)" x]
+    by simp
+  moreover have "(\<lambda>x. inverse (sqrt (sin x ^ 3 + cos x ^ 3))) twice_field_differentiable_at Base x"
+    using inverse_sqrt_valid d_sqrt real_sqrt_eq_zero_cancel_iff
+          twice_field_differentiable_at_compose twice_field_differentiable_at_inverse
+          less_numeral_extra(3)
+    by force
+  ultimately have
+    "(*h* (\<lambda>x. exp x * inverse (sqrt (sin x ^ 3 + cos x ^ 3)))) x =
+     (*h* exp) x * inverse ((*h* sqrt) ((*h* sin) x ^ 3 + (*h* cos) x ^ 3))"
+    by (simp add: hypext_fun_mult)
+  then have
+    "(*h* (\<lambda>x. exp x / sqrt (sin x ^ 3 + cos x ^ 3))) x =
+     (*h* exp) x / (*h* sqrt) ((*h* sin) x ^ 3 + (*h* cos) x ^ 3)"
+    by (simp add: inverse_eq_divide hyp_divide_inverse)
+  then show ?thesis
+    unfolding fa_test_def hyp_fa_test_def .
+qed
+
+text\<open>
+  We can show that our hyperdual extension gives (approximately) the same values as those found by
+  Fike and Alonso when evaluated at @{term "1.5"}.
+\<close>
+lemma
+  assumes "x = hyp_fa_test (\<beta> 1.5)"
+  shows "\<bar>Base x - 4.4978\<bar> \<le> 0.00005"
+    and "\<bar>Eps1 x - 4.0534\<bar> \<le> 0.00005"
+    and "\<bar>Eps12 x - 9.4631\<bar> \<le> 0.00005"
+proof -
+  show "\<bar>Base x - 4.4978\<bar> \<le> 0.00005"
+    by (simp add: assms hyp_fa_test_def, approximation 20)
+  have d: "0 < sin (3/2 :: real) ^ 3 + cos (3/2 :: real) ^ 3"
+    using sin_cube_plus_cos_cube_gt_zero_iff' sin_gt_zero pos_add_strict pi_gt3 by force
+  show "\<bar>Eps1 x - 4.0534\<bar> \<le> 0.00005"
+    by (simp add: assms hyp_fa_test_def d, approximation 22)
+  show "\<bar>Eps12 x - 9.4631\<bar> \<le> 0.00005"
+    by (simp add: assms hyp_fa_test_def d, approximation 24)
+qed
+
+text\<open>A number of additional lemmas that will be required to prove the derivatives:\<close>
+lemma hypext_sqrt_hyperdual_parts:
+      "a > 0 \<Longrightarrow>  (*h* sqrt) (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) =
+       sqrt a *\<^sub>H ba + (b * inverse (sqrt a) / 2) *\<^sub>H e1 + (c * inverse (sqrt a) / 2) *\<^sub>H e2 +
+       (d * inverse (sqrt a) / 2 - b * c * inverse (sqrt a ^ 3) / 4) *\<^sub>H e12"
+  by (metis Hyperdual_eq hypext_sqrt_Hyperdual)
+
+lemma cos_multiple: "cos (n * x) = 2 * cos x * cos ((n - 1) * x) - cos ((n - 2) * x)"
+  for x :: "'a :: {banach,real_normed_field}"
+proof -
+  have "cos ((n - 1) * x + x) + cos ((n - 1) * x - x) = 2 * cos ((n - 1) * x) * cos x"
+    by (simp add: cos_add cos_diff)
+  then show ?thesis
+    by (simp add: left_diff_distrib' eq_diff_eq)
+qed
+
+lemma sin_multiple: "sin (n * x) = 2 * cos x * sin ((n - 1) * x) - sin ((n - 2) * x)"
+  for x :: "'a :: {banach,real_normed_field}"
+proof -
+  have "sin ((n - 1) * x + x) + sin ((n - 1) * x - x) = 2 * cos x * sin ((n - 1) * x)"
+    by (simp add: sin_add sin_diff)
+  then show ?thesis
+    by (simp add: left_diff_distrib' eq_diff_eq)
+qed
+
+lemma power5:
+  fixes z :: "'a :: monoid_mult"
+  shows "z ^ 5 = z * z * z * z * z"
+  by (simp add: mult.assoc power2_eq_square power_numeral_odd)
+
+lemma power6:
+  fixes z :: "'a :: monoid_mult"
+  shows "z ^ 6 = z * z * z * z * z * z"
+  by (simp add: mult.assoc power3_eq_cube power_numeral_even)
+
+text\<open>
+  We find the derivatives of @{const fa_test} by applying a Wengert list approach, as done by Fike
+  and Alonso, to make the same composition but in hyperduals.
+  We know that this is equal to the hyperdual extension which in turn gives us the derivatives.
+\<close>
+lemma Wengert_autodiff_fa_test:
+  assumes "x \<in> (\<Union>n::int. {-pi/4 + 2*pi*n <..< 3*pi/4 + 2*pi*n})"
+  shows "First (autodiff fa_test x) =
+    (exp x * (3 * cos x + 5 * cos (3 * x) + 9 * sin x + sin (3 * x))) /
+        (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)"
+    and "Second (autodiff fa_test x) =
+    (exp x * (130 - 12 * cos (2 * x) +
+              30 * cos (4 * x) + 12 * cos (6 * x) -
+              111 * sin (2 * x) +
+              48 * sin (4 * x) + 5 * sin (6 * x))) /
+                (64 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 5)"
+proof -
+  have s3_c3_gt_zero: "(sin x) ^ 3 + (cos x) ^ 3 > 0"
+    using assms sin_45_positive_intervals sin_cube_plus_cos_cube_gt_zero_iff' sin_gt_zero
+    by force
+   \<comment> \<open>Work out @{const hyp_fa_test} as Wengert list of basic operations\<close>
+  let ?w0 = "\<beta> x"
+  have w0: "?w0 = x *\<^sub>H ba + 1 *\<^sub>H e1 + 1 *\<^sub>H e2 + 0 *\<^sub>H e12"
+    by (simp add: Hyperdual_eq hyperdualx_def)
+
+  let ?w1 = "(*h* exp) ?w0"
+  have w1: "?w1 = exp x *\<^sub>H ba + exp x *\<^sub>H e1 + exp x *\<^sub>H e2 + exp x *\<^sub>H e12"
+    using hypext_exp_extract by blast
+
+  let ?w2 = "(*h* sin) ?w0"
+  have w2: "?w2 = sin x *\<^sub>H ba + cos x *\<^sub>H e1 + cos x *\<^sub>H e2 + - sin x *\<^sub>H e12"
+    by (simp add: hypext_sin_extract of_comp_minus scaleH_times)
+
+  let ?w3 = "(*h* (\<lambda>x. x ^ 3)) ?w2"
+  have w3: "?w3 = (sin x) ^ 3 *\<^sub>H ba + (3 * cos x * (sin x)\<^sup>2) *\<^sub>H e1 + (3 * cos x * (sin x)\<^sup>2) *\<^sub>H e2 + - (3/4 * (sin x - 3 * sin (3 * x))) *\<^sub>H e12"
+    by (simp add: w2 hypext_power_Hyperdual_parts power2_eq_square cos_times_cos sin_times_sin
+          sin_times_cos distrib_left right_diff_distrib' divide_simps)
+
+  let ?w4 = "(*h* cos) ?w0"
+  have w4: "?w4 = cos x *\<^sub>H ba + - sin x *\<^sub>H e1 + - sin x *\<^sub>H e2 + - cos x *\<^sub>H e12"
+    by (simp add: hypext_cos_extract of_comp_minus scaleH_times)
+
+  let ?w5 = "(*h* (\<lambda>x. x ^ 3)) ?w4"
+  have w5: "?w5 = (cos x) ^ 3 *\<^sub>H ba + - (3 * sin x * (cos x)\<^sup>2) *\<^sub>H e1 + - (3 * sin x * (cos x)\<^sup>2) *\<^sub>H e2 + - (3/4 * (cos x + 3 * cos (3 * x))) *\<^sub>H e12"
+    by (simp add: w4 hypext_power_Hyperdual_parts sin_times_sin right_diff_distrib' cos_times_cos
+              power2_eq_square distrib_left sin_times_cos divide_simps)
+
+  let ?w6 = "?w3 + ?w5"
+  have sqrt_pos: "Base ?w6 > 0"
+      using s3_c3_gt_zero by auto
+  have w6: "?w6 = (sin x ^ 3 + cos x ^ 3) *\<^sub>H ba +
+                  (3 * cos x * sin x * (sin x - cos x)) *\<^sub>H e1 +
+                  (3 * cos x * sin x * (sin x - cos x)) *\<^sub>H e2 +
+                  - (3/4 * (sin x + cos x + 3 * cos (3 * x) - 3 * sin (3 * x))) *\<^sub>H e12"
+    by (auto simp add: w3 w5 add_hyperdual_parts power2_eq_square right_diff_distrib' divide_simps)
+
+  let ?w7 = "inverse ((*h* sqrt) ?w6)"
+  have w7: "?w7 = inverse(sqrt(sin x ^ 3 + cos x ^ 3)) *\<^sub>H ba +
+                  - ((3 * cos x * sin x * (sin x - cos x))/(2 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)) *\<^sub>H e1 +
+                  - ((3 * cos x * sin x * (sin x - cos x))/(2 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)) *\<^sub>H e2 +
+                  ((3 * (30 + 2 * cos (4 * x) - 41 * sin (2 * x) + 3 * sin (6 * x)))/(64 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 5)) *\<^sub>H e12"
+    \<comment> \<open>Apply the functions in two steps, then simplify the result to match\<close>
+  proof -
+    let ?w7a = "(*h* sqrt) ?w6"
+    have w7a: "?w7a =
+      sqrt (sin x ^ 3 + cos x ^ 3) *\<^sub>H ba +
+      ((3 * (cos x * (sin x * (sin x - cos x)))) * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) / 2) *\<^sub>H e1 +
+      ((3 * (cos x * (sin x * (sin x - cos x)))) * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) / 2) *\<^sub>H e2 +
+      (- ((3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) / 8) +
+       - 9 * (cos x * (sin x * ((sin x - cos x) * (cos x * (sin x * (sin x - cos x)))))) * inverse (sqrt (sin x ^ 3 + cos x ^ 3) ^ 3) / 4) *\<^sub>H e12"
+      unfolding w6
+      using sqrt_pos by (simp add: hypext_sqrt_hyperdual_parts mult.assoc)
+
+    let ?w7b = "inverse ?w7a"
+    have "?w7b =
+      (1 / sqrt (sin x ^ 3 + cos x ^ 3)) *\<^sub>H ba +
+      - (3 * (cos x * (sin x * (sin x - cos x))) * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) / (2 * (sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2)) *\<^sub>H e1 +
+      - (3 * (cos x * (sin x * (sin x - cos x))) * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) / (2 * (sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2)) *\<^sub>H e2 +
+      (9 * (cos x * (cos x * (sin x * (sin x * ((sin x - cos x) * ((sin x - cos x) * (inverse (sqrt (sin x ^ 3 + cos x ^ 3)) * inverse (sqrt (sin x ^ 3 + cos x ^ 3))))))))) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3) -
+        (- ((3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) / 8) -
+        9 * (cos x * (sin x * ((sin x - cos x) * (cos x * (sin x * (sin x - cos x)))))) * inverse (sqrt (sin x ^ 3 + cos x ^ 3) ^ 3) / 4) /
+     (sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2) *\<^sub>H e12"
+      \<comment> \<open>Push the inverse in\<close>
+      by (simp add: w7a inverse_hyperdual_parts)
+    then have w7b: "?w7b =
+      (inverse (sqrt (sin x ^ 3 + cos x ^ 3))) *\<^sub>H ba +
+      - (3 * cos x * sin x * (sin x - cos x) / (2 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)) *\<^sub>H e1 +
+      - (3 * cos x * sin x * (sin x - cos x) / (2 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)) *\<^sub>H e2 +
+      (9 * cos x * cos x * sin x * sin x * ((sin x - cos x) * (sin x - cos x) / ((sqrt (sin x ^ 3 + cos x ^ 3)) ^ 2)) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3) -
+        (- ((3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) / (8 * sqrt (sin x ^ 3 + cos x ^ 3))) -
+         9 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x) / (4 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3))
+        / (sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2) *\<^sub>H e12"
+      \<comment> \<open>Simplify powers and parents\<close>
+      by (simp add: power2_eq_square power3_eq_cube field_simps)
+
+    have
+      "9 * cos x * cos x * sin x * sin x * ((sin x - cos x) * (sin x - cos x)) / ((sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2 * (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3)) -
+       (- ((3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) / (8 * sqrt (sin x ^ 3 + cos x ^ 3))) -
+          9 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x) / (4 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3)) /
+       (sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2 =
+       (90 + 6 * cos (4 * x) - 123 * sin (2 * x) + 9 * sin (6 * x)) / (64 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5)"
+      \<comment> \<open>Equate the last component's coefficients\<close>
+    proof -
+      have
+        "9 * cos x * cos x * sin x * sin x * ((sin x - cos x) * (sin x - cos x)) / ((sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2 * (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3)) -
+         (- ((3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) / (8 * sqrt (sin x ^ 3 + cos x ^ 3))) -
+           9 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x) / (4 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3)) /
+         (sqrt (sin x ^ 3 + cos x ^ 3))\<^sup>2 =
+         9 * cos x * cos x * sin x * sin x * (sin x - cos x) * (sin x - cos x) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5) +
+         ( (3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) / (8 * sqrt (sin x ^ 3 + cos x ^ 3)) +
+           9 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x) / (4 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3)) /
+         (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 2"
+          by (simp add: divide_simps)
+      also have "... =
+         9 * cos x * cos x * sin x * sin x * (sin x - cos x) * (sin x - cos x) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5) +
+         ( (3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) / (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3) +
+           9 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x) / (4 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 2))"
+        by (simp add: add_divide_distrib power2_eq_square power3_eq_cube mult.commute mult.left_commute)
+      also have "... =
+         9 * cos x * cos x * sin x * sin x * (sin x - cos x) * (sin x - cos x) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5) +
+         ( (3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) / (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3) +
+           9 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x) / (4 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5))"
+        by (simp add: mult.commute mult.left_commute)
+      also have "... =
+         9 * cos x * cos x * sin x * sin x * (sin x - cos x) * (sin x - cos x) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5) +
+         ( (3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) * (sin x ^ 3 + cos x ^ 3) / (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 5) +
+           9 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x) / (4 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5))"
+      proof -
+        have "(sin x ^ 3 + cos x ^ 3) * inverse ((sqrt (sin x ^ 3 + cos x ^ 3)) ^ 2) = 1"
+          using s3_c3_gt_zero by auto
+        then have
+          "1 / (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3) =
+           (sin x ^ 3 + cos x ^ 3) / ((sqrt (sin x ^ 3 + cos x ^ 3)) ^ 2) / (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)"
+          by (simp add: field_simps)
+        then have
+          "1 / (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3) =
+           (sin x ^ 3 + cos x ^ 3) / (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 5)"
+          by auto
+        then show ?thesis
+          by (metis (mono_tags, lifting) divide_real_def inverse_eq_divide times_divide_eq_right)
+      qed
+      also have "... =
+         ( 288 * cos x * cos x * sin x * sin x * (sin x - cos x) * (sin x - cos x) +
+           (3 * sin x + 3 * cos x + 9 * cos (3 * x) - 9 * sin (3 * x)) * (sin x ^ 3 + cos x ^ 3) * 8 +
+           144 * cos x * sin x * (sin x - cos x) * cos x * sin x * (sin x - cos x))
+         / (64 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 5)"
+        by (simp add: divide_simps)
+      also have "... =
+         ( 432 * cos x * cos x * sin x * sin x * (sin x - cos x) * (sin x - cos x) +
+           (24 * sin x + 24 * cos x + 72 * cos (3 * x) - 72 * sin (3 * x)) * (sin x ^ 3 + cos x ^ 3))
+         / (64 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 5)"
+        by (simp add: divide_simps)
+      finally show ?thesis
+        \<comment> \<open>Having matched the denominators, prove the numerators equal\<close>
+        by (simp, force intro:  disjI2 simp add: power2_eq_square power4_eq_xxxx power3_eq_cube cos_times_sin
+            sin_times_cos cos_times_cos sin_times_sin right_diff_distrib power6 power5
+            distrib_left distrib_right left_diff_distrib divide_simps)
+    qed
+    then show ?thesis
+      \<comment> \<open>Put it all together\<close>
+      by (simp add: w7b)
+  qed
+
+  let ?w8 = "?w1 * ?w7"
+  have w8: "?w8 =
+    ((exp x)/(sqrt(sin x ^ 3 + cos x ^ 3))) *\<^sub>H ba +
+    ((exp x * (3 * cos x + 5 * cos (3 * x) + 9 * sin x + sin (3 * x)))/(8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)) *\<^sub>H e1 +
+    ((exp x * (3 * cos x + 5 * cos (3 * x) + 9 * sin x + sin (3 * x)))/(8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)) *\<^sub>H e2 +
+    ((exp x * (130 - 12 * cos (2 * x) + 30 * cos (4 * x) + 12 * cos (6 * x) - 111 * sin (2 * x) + 48 * sin (4 * x) +
+               5 * sin (6 * x)))/(64 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 5)) *\<^sub>H e12"
+  proof (auto simp add: w7 w1 times_hyperdual_parts)
+    show "exp x * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) = exp x / sqrt (sin x ^ 3 + cos x ^ 3)"
+      by (simp add: divide_inverse)
+
+    have sqrt_sc3: "sqrt (sin x ^ 3 + cos x ^ 3) ^ 3 = (sin x ^ 3 + cos x ^ 3) * sqrt (sin x ^ 3 + cos x ^ 3)"
+      using s3_c3_gt_zero
+      by (simp add: power3_eq_cube)
+    then have "inverse (sqrt (sin x ^ 3 + cos x ^ 3)) -
+          (3 * cos x * sin x * (sin x - cos x)) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3) =
+          (3 * cos x + 5 * cos (3 * x) + 9 * sin x + sin (3 * x)) / (8 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3)"
+      using sqrt_pos
+      apply (simp add: right_diff_distrib' sin_times_sin cos_times_cos sin_times_cos cos_times_sin power3_eq_cube left_diff_distrib' divide_simps)
+      by (simp add: distrib_right cos_times_cos)
+    then show " exp x * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) -
+           exp x * (3 * cos x * sin x * (sin x - cos x)) / (2 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3) =
+          exp x * (3 * cos x + 5 * cos (3 * x) + 9 * sin x + sin (3 * x)) / (8 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 3)"
+      by (simp add: algebra_simps)
+
+    have "sqrt (sin x ^ 3 + cos x ^ 3) ^ 8 =  (sin x ^ 3 + cos x ^ 3) ^ 4"
+      using s3_c3_gt_zero
+      by (smt mult_2 numeral_Bit0 power2_eq_square power_even_eq real_sqrt_mult_self)
+    moreover have "sqrt (sin x ^ 3 + cos x ^ 3) ^ 5 =  (sin x ^ 3 + cos x ^ 3) ^ 2 * sqrt (sin x ^ 3 + cos x ^ 3)"
+      by (simp add: mult.assoc power2_eq_square power5)
+    ultimately
+    have "(90 + 6 * cos (4 * x) - 123 * sin (2 * x) + 9 * sin (6 * x)) / (64 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5) -
+           3 * (cos x * ((sin x * (sin x - cos x)))) / sqrt (sin x ^ 3 + cos x ^ 3) ^ 3 +
+           inverse (sqrt (sin x ^ 3 + cos x ^ 3)) =
+           (130 - 12 * cos (2 * x) + 30 * cos (4 * x) + 12 * cos (6 * x) - 111 * sin (2 * x) + 48 * sin (4 * x) + 5 * sin (6 * x)) /
+           (64 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5)"
+      using sqrt_pos
+      apply (simp add: sqrt_sc3 power2_eq_square divide_simps)
+      by (simp add: distrib_right distrib_left power3_eq_cube sin_times_sin sin_times_cos
+              cos_times_cos cos_times_sin right_diff_distrib' left_diff_distrib' divide_simps)
+    moreover have "\<forall>r a b c. (a::real) * b - c * (a * r) = a * (b - c * r)"
+      by (simp add: right_diff_distrib)
+    ultimately show "exp x * (90 + 6 * cos (4 * x) - 123 * sin (2 * x) + 9 * sin (6 * x)) / (64 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5) -
+               3 * (cos x * (exp x * (sin x * (sin x - cos x)))) / sqrt (sin x ^ 3 + cos x ^ 3) ^ 3 +
+              exp x * inverse (sqrt (sin x ^ 3 + cos x ^ 3)) =
+              exp x *
+              (130 - 12 * cos (2 * x) + 30 * cos (4 * x) + 12 * cos (6 * x) - 111 * sin (2 * x) + 48 * sin (4 * x) + 5 * sin (6 * x)) /
+              (64 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5)"
+      by (metis (mono_tags, opaque_lifting) distrib_left mult.left_commute times_divide_eq_right)
+  qed
+
+  \<comment> \<open>Check that we have indeed computed the hyperdual extension of our analytic test function \<close>
+  moreover have w8_eq_hyp_fa_test: "?w8 = (*h* fa_test) (\<beta> x)"
+    using assms
+    by (simp add: hyp_divide_inverse hyp_fa_test_def hypext_fa_test hypext_power)
+
+  \<comment> \<open>Now simply show that "extraction" of first and second derivatives are as expected\<close>
+  ultimately
+  show "First (autodiff fa_test x) =
+          (exp x * (3 * cos x + 5 * cos (3 * x) + 9 * sin x + sin (3 * x))) /
+            (8 * (sqrt (sin x ^ 3 + cos x ^ 3)) ^ 3)"
+  and
+    "Second (autodiff fa_test x) =
+        exp x * (130 - 12 * cos (2 * x) +
+                  30 * cos (4 * x) + 12 * cos (6 * x) -
+                   111 * sin (2 * x) + 48 * sin (4 * x) + 5 * sin (6 * x)) /
+       (64 * sqrt (sin x ^ 3 + cos x ^ 3) ^ 5)"
+    by (metis autodiff_sel hyperdual_comb_sel)+
+qed
+
+end
diff --git a/thys/Hyperdual/Hyperdual.thy b/thys/Hyperdual/Hyperdual.thy
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/Hyperdual.thy
@@ -0,0 +1,1454 @@
+(*  Title:   Hyperdual.thy
+    Authors: Filip Smola and Jacques D. Fleuriot, University of Edinburgh, 2019-2021
+*)
+
+theory Hyperdual
+  imports "HOL-Analysis.Analysis"
+begin
+
+section\<open>Hyperdual Numbers\<close>
+text\<open>
+  Let \<open>\<tau>\<close> be some type.
+  Second-order hyperdual numbers over \<open>\<tau>\<close> take the form \<open>a\<^sub>1 + a\<^sub>2\<epsilon>\<^sub>1 + a\<^sub>3\<epsilon>\<^sub>2 + a\<^sub>4\<epsilon>\<^sub>1\<epsilon>\<^sub>2\<close> where all
+  \<open>a\<^sub>i :: \<tau>\<close>, and \<open>\<epsilon>\<^sub>1\<close> and \<open>\<epsilon>\<^sub>2\<close> are non-zero but nilpotent infinitesimals: \<open>\<epsilon>\<^sub>1\<^sup>2 = \<epsilon>\<^sub>2\<^sup>2 = (\<epsilon>\<^sub>1\<epsilon>\<^sub>2)\<^sup>2 = 0\<close>.
+
+  We define second-order hyperdual numbers as a coinductive data type with four components: the base
+  component, two first-order hyperdual components and one second-order hyperdual component.
+\<close>
+codatatype 'a hyperdual = Hyperdual (Base: 'a) (Eps1: 'a) (Eps2: 'a) (Eps12: 'a)
+
+text\<open>Two hyperduals are equal iff all their components are equal.\<close>
+lemma hyperdual_eq_iff [iff]:
+  "x = y \<longleftrightarrow> ((Base x = Base y) \<and> (Eps1 x = Eps1 y) \<and> (Eps2 x = Eps2 y) \<and> (Eps12 x = Eps12 y))"
+  using hyperdual.expand by auto
+
+lemma hyperdual_eqI:
+  assumes "Base x = Base y"
+      and "Eps1 x = Eps1 y"
+      and "Eps2 x = Eps2 y"
+      and "Eps12 x = Eps12 y"
+    shows "x = y"
+  by (simp add: assms)
+
+text\<open>
+  The embedding from the component type to hyperduals requires the component type to have a zero
+  element.
+\<close>
+definition of_comp :: "('a :: zero) \<Rightarrow> 'a hyperdual"
+  where "of_comp a = Hyperdual a 0 0 0"
+
+lemma of_comp_simps [simp]:
+  "Base (of_comp a) = a"
+  "Eps1 (of_comp a) = 0"
+  "Eps2 (of_comp a) = 0"
+  "Eps12 (of_comp a) = 0"
+  by (simp_all add: of_comp_def)
+
+subsection \<open>Addition and Subtraction\<close>
+
+text\<open>We define hyperdual addition, subtraction and unary minus pointwise, and zero by embedding.\<close>
+(* Define addition using component addition *)
+instantiation hyperdual :: (plus) plus
+begin
+
+primcorec plus_hyperdual
+  where
+    "Base (x + y) = Base x + Base y"
+  | "Eps1 (x + y) = Eps1 x + Eps1 y"
+  | "Eps2 (x + y) = Eps2 x + Eps2 y"
+  | "Eps12 (x + y) = Eps12 x + Eps12 y"
+
+instance by standard
+end
+
+(* Define zero by embedding component zero *)
+instantiation hyperdual :: (zero) zero
+begin
+
+definition zero_hyperdual
+  where "0 = of_comp 0"
+
+instance by standard
+end
+
+lemma zero_hyperdual_simps [simp]:
+  "Base 0 = 0"
+  "Eps1 0 = 0"
+  "Eps2 0 = 0"
+  "Eps12 0 = 0"
+  "Hyperdual 0 0 0 0 = 0"
+  by (simp_all add: zero_hyperdual_def)
+
+(* Define minus using component minus *)
+instantiation hyperdual :: (uminus) uminus
+begin
+
+primcorec uminus_hyperdual
+  where
+    "Base (-x) = - Base x"
+  | "Eps1 (-x) = - Eps1 x"
+  | "Eps2 (-x) = - Eps2 x"
+  | "Eps12 (-x) = - Eps12 x"
+
+instance by standard
+end
+
+(* Define subtraction using component subtraction *)
+instantiation hyperdual :: (minus) minus
+begin
+
+primcorec minus_hyperdual
+  where
+    "Base (x - y) = Base x - Base y"
+  | "Eps1 (x - y) = Eps1 x - Eps1 y"
+  | "Eps2 (x - y) = Eps2 x - Eps2 y"
+  | "Eps12 (x - y) = Eps12 x - Eps12 y"
+
+instance by standard
+end
+
+text\<open>If the components form a commutative group under addition then so do the hyperduals.\<close>
+instance hyperdual :: (semigroup_add) semigroup_add
+  by standard (simp add: add.assoc)
+
+instance hyperdual :: (monoid_add) monoid_add
+  by standard simp_all
+
+instance hyperdual :: (ab_semigroup_add) ab_semigroup_add
+  by standard (simp_all add: add.commute)
+
+instance hyperdual :: (comm_monoid_add) comm_monoid_add
+  by standard simp
+
+instance hyperdual :: (group_add) group_add
+  by standard simp_all
+
+instance hyperdual :: (ab_group_add) ab_group_add
+  by standard simp_all
+
+lemma of_comp_add:
+  fixes a b :: "'a :: monoid_add"
+  shows "of_comp (a + b) = of_comp a + of_comp b"
+  by simp
+
+lemma
+  fixes a b :: "'a :: group_add"
+  shows of_comp_minus: "of_comp (- a) = - of_comp a"
+    and of_comp_diff: "of_comp (a - b) = of_comp a - of_comp b"
+  by simp_all
+
+subsection \<open>Multiplication and Scaling\<close>
+
+text\<open>
+  Multiplication of hyperduals is defined by distributing the expressions and using the nilpotence
+  of \<open>\<epsilon>\<^sub>1\<close> and \<open>\<epsilon>\<^sub>2\<close>, resulting in the definition used here.
+  The hyperdual one is again defined by embedding.
+\<close>
+(* Define one by embedding the component one, which also requires it to have zero *)
+instantiation hyperdual :: ("{one, zero}") one
+begin
+
+definition one_hyperdual
+  where "1 = of_comp 1"
+
+instance by standard
+end
+
+lemma one_hyperdual_simps [simp]:
+  "Base 1 = 1"
+  "Eps1 1 = 0"
+  "Eps2 1 = 0"
+  "Eps12 1 = 0"
+  "Hyperdual 1 0 0 0 = 1"
+  by (simp_all add: one_hyperdual_def)
+
+(* Define multiplication using component multiplication and addition *)
+instantiation hyperdual :: ("{times, plus}") times
+begin
+
+primcorec times_hyperdual
+  where
+    "Base (x * y) = Base x * Base y"
+  | "Eps1 (x * y) = (Base x * Eps1 y) + (Eps1 x * Base y)"
+  | "Eps2 (x * y) = (Base x * Eps2 y) + (Eps2 x * Base y)"
+  | "Eps12 (x * y) = (Base x * Eps12 y) + (Eps1 x * Eps2 y) + (Eps2 x * Eps1 y) + (Eps12 x * Base y)"
+
+instance by standard
+end
+
+text\<open>If the components form a ring then so do the hyperduals.\<close>
+instance hyperdual :: (semiring) semiring
+  by standard (simp_all add: mult.assoc distrib_left distrib_right add.assoc add.left_commute)
+
+instance hyperdual :: ("{monoid_add, mult_zero}") mult_zero
+  by standard simp_all
+
+instance hyperdual :: (ring) ring
+  by standard
+
+instance hyperdual :: (comm_ring) comm_ring
+  by standard (simp_all add: mult.commute distrib_left)
+
+instance hyperdual :: (ring_1) ring_1
+  by standard simp_all
+
+instance hyperdual :: (comm_ring_1) comm_ring_1
+  by standard simp
+
+lemma of_comp_times:
+  fixes a b :: "'a :: semiring_0"
+  shows "of_comp (a * b) = of_comp a * of_comp b"
+  by (simp add: of_comp_def times_hyperdual.code)
+
+text\<open>Hyperdual scaling is multiplying each component by a factor from the component type.\<close>
+(* Named scaleH for hyperdual like scaleR is for real *)
+primcorec scaleH :: "('a :: times) \<Rightarrow> 'a hyperdual \<Rightarrow> 'a hyperdual"  (infixr "*\<^sub>H" 75)
+  where
+    "Base (f *\<^sub>H x) = f * Base x"
+  | "Eps1 (f *\<^sub>H x) = f * Eps1 x"
+  | "Eps2 (f *\<^sub>H x) = f * Eps2 x"
+  | "Eps12 (f *\<^sub>H x) = f * Eps12 x"
+
+lemma scaleH_times:
+  fixes f :: "'a :: {monoid_add, mult_zero}"
+  shows "f *\<^sub>H x = of_comp f * x"
+  by simp
+
+lemma scaleH_add:
+  fixes a :: "'a :: semiring"
+  shows "(a + a') *\<^sub>H b = a *\<^sub>H b + a' *\<^sub>H b"
+    and "a *\<^sub>H (b + b') = a *\<^sub>H b + a *\<^sub>H b'"
+  by (simp_all add: distrib_left distrib_right)
+
+lemma scaleH_diff:
+  fixes a :: "'a :: ring"
+  shows "(a - a') *\<^sub>H b = a *\<^sub>H b - a' *\<^sub>H b"
+    and "a *\<^sub>H (b - b') = a *\<^sub>H b - a *\<^sub>H b'"
+  by (auto simp add: left_diff_distrib right_diff_distrib scaleH_times of_comp_diff)
+
+lemma scaleH_mult:
+  fixes a :: "'a :: semigroup_mult"
+  shows "(a * a') *\<^sub>H b = a *\<^sub>H a' *\<^sub>H b"
+  by (simp add: mult.assoc)
+
+lemma scaleH_one [simp]:
+  fixes b :: "('a :: monoid_mult) hyperdual"
+  shows "1 *\<^sub>H b = b"
+  by simp
+
+lemma scaleH_zero [simp]:
+  fixes b :: "('a :: {mult_zero, times}) hyperdual"
+  shows "0 *\<^sub>H b = 0"
+  by simp
+
+lemma
+  fixes b :: "('a :: ring_1) hyperdual"
+  shows scaleH_minus [simp]:"- 1 *\<^sub>H b = - b"
+    and scaleH_minus_left: "- (a *\<^sub>H b) = - a *\<^sub>H b"
+    and scaleH_minus_right: "- (a *\<^sub>H b) = a *\<^sub>H - b"
+  by simp_all
+
+text\<open>Induction rule for natural numbers that takes 0 and 1 as base cases.\<close>
+lemma nat_induct01Suc[case_names 0 1 Suc]:
+  assumes "P 0"
+      and "P 1"
+      and "\<And>n. n > 0 \<Longrightarrow> P n \<Longrightarrow> P (Suc n)"
+    shows "P n"
+  by (metis One_nat_def assms nat_induct neq0_conv)
+
+lemma hyperdual_power:
+  fixes x :: "('a :: comm_ring_1) hyperdual"
+  shows "x ^ n = Hyperdual ((Base x) ^ n)
+                           (Eps1 x * of_nat n * (Base x) ^ (n - 1))
+                           (Eps2 x * of_nat n * (Base x) ^ (n - 1))
+                           (Eps12 x * of_nat n * (Base x) ^ (n - 1) + Eps1 x * Eps2 x * of_nat n * of_nat (n - 1) * (Base x) ^ (n - 2))"
+proof (induction n rule: nat_induct01Suc)
+  case 0
+  show ?case
+    by simp
+next
+  case 1
+  show ?case
+    by simp
+next
+  case hyp: (Suc n)
+  show ?case
+  proof (simp add: hyp, intro conjI)
+    show "Base x * (Eps1 x * of_nat n * Base x ^ (n - Suc 0)) + Eps1 x * Base x ^ n = Eps1 x * (1 + of_nat n) * Base x ^ n"
+     and "Base x * (Eps2 x * of_nat n * Base x ^ (n - Suc 0)) + Eps2 x * Base x ^ n = Eps2 x * (1 + of_nat n) * Base x ^ n"
+      by (simp_all add: distrib_left distrib_right power_eq_if)
+    show
+      "2 * (Eps1 x * (Eps2 x * (of_nat n * Base x ^ (n - Suc 0)))) +
+       Base x * (Eps12 x * of_nat n * Base x ^ (n - Suc 0) + Eps1 x * Eps2 x * of_nat n * of_nat (n - Suc 0) * Base x ^ (n - 2)) +
+       Eps12 x * Base x ^ n =
+       Eps12 x * (1 + of_nat n) * Base x ^ n + Eps1 x * Eps2 x * (1 + of_nat n) * of_nat n * Base x ^ (n - Suc 0)"
+    proof -
+      have
+        "2 * (Eps1 x * (Eps2 x * (of_nat n * Base x ^ (n - Suc 0)))) +
+         Base x * (Eps12 x * of_nat n * Base x ^ (n - Suc 0) + Eps1 x * Eps2 x * of_nat n * of_nat (n - Suc 0) * Base x ^ (n - 2)) +
+         Eps12 x * Base x ^ n =
+         2 * Eps1 x * Eps2 x * of_nat n * Base x ^ (n - Suc 0) +
+         Eps12 x * of_nat (n + 1) * Base x ^ n + Eps1 x * Eps2 x * of_nat n * of_nat (n - Suc 0) * Base x ^ (n - Suc 0)"
+        by (simp add: field_simps power_eq_if)
+      also have "... = Eps12 x * of_nat (n + 1) * Base x ^ n + of_nat (n - 1 + 2) * Eps1 x * Eps2 x * of_nat n * Base x ^ (n - Suc 0)"
+        by (simp add: distrib_left mult.commute)
+      finally show ?thesis
+        by (simp add: hyp.hyps)
+    qed
+  qed
+qed
+
+lemma hyperdual_power_simps [simp]:
+  shows "Base ((x :: 'a :: comm_ring_1 hyperdual) ^ n) = Base x ^ n"
+    and "Eps1 ((x :: 'a :: comm_ring_1 hyperdual) ^ n) = Eps1 x * of_nat n * (Base x) ^ (n - 1)"
+    and "Eps2 ((x :: 'a :: comm_ring_1 hyperdual) ^ n) = Eps2 x * of_nat n * (Base x) ^ (n - 1)"
+    and "Eps12 ((x :: 'a :: comm_ring_1 hyperdual) ^ n) =
+  (Eps12 x * of_nat n * (Base x) ^ (n - 1) + Eps1 x * Eps2 x * of_nat n * of_nat (n - 1) * (Base x) ^ (n - 2))"
+  by (simp_all add: hyperdual_power)
+
+text\<open>Squaring the hyperdual one behaves as expected from the reals.\<close>
+(* Commutativity required to reorder times definition expressions, division algebra required for
+    base component's x * x = 1 \<longrightarrow> x = 1 \<or> x = -1 *)
+lemma hyperdual_square_eq_1_iff [iff]:
+  fixes x :: "('a :: {real_div_algebra, comm_ring}) hyperdual"
+  shows "x * x = 1 \<longleftrightarrow> x = 1 \<or> x = - 1"
+proof
+  assume a: "x * x = 1"
+
+  have base: "Base x * Base x = 1"
+    using a by simp
+  moreover have e1: "Eps1 x = 0"
+  proof -
+    have "Base x * Eps1 x = - (Base x * Eps1 x)"
+      using mult.commute[of "Base x"] add_eq_0_iff[of "Base x * Eps1 x"] times_hyperdual.simps(2)[of x x]
+      by (simp add: a)
+    then have "Base x * Base x * Eps1 x = - Base x * Base x * Eps1 x"
+      using mult_left_cancel[of "Base x"] base by fastforce
+    then show ?thesis
+      using base mult_right_cancel[of "Eps1 x" "Base x * Base x" "- Base x * Base x"] one_neq_neg_one
+      by auto
+  qed
+  moreover have e2: "Eps2 x = 0"
+  proof -
+    have "Base x * Eps2 x = - (Base x * Eps2 x)"
+      using a mult.commute[of "Base x" "Eps2 x"] add_eq_0_iff[of "Base x * Eps2 x"] times_hyperdual.simps(3)[of x x]
+      by simp
+    then have "Base x * Base x * Eps2 x = - Base x * Base x * Eps2 x"
+      using mult_left_cancel[of "Base x"] base by fastforce
+    then show ?thesis
+      using base mult_right_cancel[of "Eps2 x" "Base x * Base x" "- Base x * Base x"] one_neq_neg_one
+      by auto
+  qed
+  moreover have "Eps12 x = 0"
+  proof -
+    have "Base x * Eps12 x = - (Base x * Eps12 x)"
+      using a e1 e2 mult.commute[of "Base x" "Eps12 x"] add_eq_0_iff[of "Base x * Eps12 x"] times_hyperdual.simps(4)[of x x]
+      by simp
+    then have "Base x * Base x * Eps12 x = - Base x * Base x * Eps12 x"
+      using mult_left_cancel[of "Base x"] base by fastforce
+    then show ?thesis
+      using base mult_right_cancel[of "Eps12 x" "Base x * Base x" "- Base x * Base x"] one_neq_neg_one
+      by auto
+  qed
+  ultimately show "x = 1 \<or> x = - 1"
+    using square_eq_1_iff[of "Base x"] by simp
+next
+  assume "x = 1 \<or> x = - 1"
+  then show "x * x = 1"
+    by (simp, safe, simp_all)
+qed
+
+subsubsection\<open>Properties of Zero Divisors\<close>
+
+text\<open>Unlike the reals, hyperdual numbers may have non-trivial divisors of zero as we show below.\<close>
+
+text\<open>
+  First, if the components have no non-trivial zero divisors then that behaviour is preserved on the
+  base component.
+\<close>
+lemma divisors_base_zero:
+  fixes a b :: "('a :: ring_no_zero_divisors) hyperdual"
+  assumes "Base (a * b) = 0"
+  shows "Base a = 0 \<or> Base b = 0"
+  using assms by auto
+lemma hyp_base_mult_eq_0_iff [iff]:
+  fixes a b :: "('a :: ring_no_zero_divisors) hyperdual"
+  shows "Base (a * b) = 0 \<longleftrightarrow> Base a = 0 \<or> Base b = 0"
+  by simp
+
+text\<open>
+  However, the conditions are relaxed on the full hyperdual numbers.
+  This is due to some terms vanishing in the multiplication and thus not constraining the result.
+\<close>
+lemma divisors_hyperdual_zero [iff]:
+  fixes a b :: "('a :: ring_no_zero_divisors) hyperdual"
+  shows "a * b = 0 \<longleftrightarrow> (a = 0 \<or> b = 0 \<or> (Base a = 0 \<and> Base b = 0 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b))"
+proof
+  assume mult: "a * b = 0"
+  then have split: "Base a = 0 \<or> Base b = 0"
+    by simp
+  show "a = 0 \<or> b = 0 \<or> Base a = 0 \<and> Base b = 0 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b"
+  proof (cases "Base a = 0")
+    case aT: True
+    then show ?thesis
+    proof (cases "Base b = 0")
+      \<comment> \<open>@{term "Base a = 0 \<and> Base b = 0"}\<close>
+      case bT: True
+      then have "Eps12 (a * b) = Eps1 a * Eps2 b + Eps2 a * Eps1 b"
+        by (simp add: aT)
+      then show ?thesis
+        by (simp add: aT bT mult eq_neg_iff_add_eq_0)
+    next
+      \<comment> \<open>@{term "Base a = 0 \<and> Base b \<noteq> 0"}\<close>
+      case bF: False
+      then have e1: "Eps1 a = 0"
+      proof -
+        have "Eps1 (a * b) = Eps1 a * Base b"
+          by (simp add: aT)
+        then show ?thesis
+          by (simp add: bF mult)
+      qed
+      moreover have e2: "Eps2 a = 0"
+      proof -
+        have "Eps2 (a * b) = Eps2 a * Base b"
+          by (simp add: aT)
+        then show ?thesis
+          by (simp add: bF mult)
+      qed
+      moreover have "Eps12 a = 0"
+      proof -
+        have "Eps12 (a * b) = Eps1 a * Eps2 b + Eps2 a * Eps1 b"
+          by (simp add: e1 e2 mult)
+        then show ?thesis
+          by (simp add: aT bF)
+      qed
+      ultimately show ?thesis
+        by (simp add: aT)
+    qed
+  next
+    case aF: False
+    then show ?thesis
+    proof (cases "Base b = 0")
+      \<comment> \<open>@{term "Base a \<noteq> 0 \<and> Base b = 0"}\<close>
+      case bT: True
+      then have e1: "Eps1 b = 0"
+      proof -
+        have "Eps1 (a * b) = Base a * Eps1 b"
+          by (simp add: bT)
+        then show ?thesis
+          by (simp add: aF mult)
+      qed
+      moreover have e2: "Eps2 b = 0"
+      proof -
+        have "Eps2 (a * b) = Base a * Eps2 b"
+          by (simp add: bT)
+        then show ?thesis
+          by (simp add: aF mult)
+      qed
+      moreover have "Eps12 b = 0"
+      proof -
+        have "Eps12 (a * b) = Eps1 a * Eps2 b + Eps2 a * Eps1 b"
+          by (simp add: e1 e2 mult)
+        then show ?thesis
+          by (simp add: bT aF)
+      qed
+      ultimately show ?thesis
+        by (simp add: bT)
+    next
+      \<comment> \<open>@{term "Base a \<noteq> 0 \<and> Base b \<noteq> 0"}\<close>
+      case bF: False
+      then have "False"
+        using split aF by blast
+      then show ?thesis
+        by simp
+    qed
+  qed
+next
+  show "a = 0 \<or> b = 0 \<or> Base a = 0 \<and> Base b = 0 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b \<Longrightarrow> a * b = 0"
+    by (simp, auto)
+qed
+
+subsubsection\<open>Multiplication Cancellation\<close>
+
+text\<open>
+  Similarly to zero divisors, multiplication cancellation rules for hyperduals are not exactly the
+  same as those for reals.
+\<close>
+
+text\<open>
+  First, cancelling a common factor has a relaxed condition compared to reals.
+  It only requires the common factor to have base component zero, instead of requiring the whole
+  number to be zero.
+\<close>
+lemma hyp_mult_left_cancel [iff]:
+  fixes a b c :: "('a :: ring_no_zero_divisors) hyperdual"
+  assumes baseC: "Base c \<noteq> 0"
+  shows "c * a = c * b \<longleftrightarrow> a = b"
+proof
+  assume mult: "c * a = c * b"
+  show "a = b"
+  proof (simp, safe)
+    show base: "Base a = Base b"
+      using mult mult_cancel_left baseC by simp_all
+    show "Eps1 a = Eps1 b"
+     and "Eps2 a = Eps2 b"
+      using mult base mult_cancel_left baseC by simp_all
+    then show "Eps12 a = Eps12 b"
+      using mult base mult_cancel_left baseC by simp_all
+  qed
+next
+  show "a = b \<Longrightarrow> c * a = c * b"
+    by simp
+qed
+
+lemma hyp_mult_right_cancel [iff]:
+  fixes a b c :: "('a :: ring_no_zero_divisors) hyperdual"
+  assumes baseC: "Base c \<noteq> 0"
+    shows "a * c = b * c \<longleftrightarrow> a = b"
+proof
+  assume mult: "a * c = b * c"
+  show "a = b"
+  proof (simp, safe)
+    show base: "Base a = Base b"
+      using mult mult_cancel_left baseC by simp_all
+    show "Eps1 a = Eps1 b"
+     and "Eps2 a = Eps2 b"
+      using mult base mult_cancel_left baseC by simp_all
+    then show "Eps12 a = Eps12 b"
+      using mult base mult_cancel_left baseC by simp_all
+  qed
+next
+  show "a = b \<Longrightarrow> a * c = b * c"
+    by simp
+qed
+
+text\<open>
+  Next, when a factor absorbs another there are again relaxed conditions compared to reals.
+  For reals, either the absorbing factor is zero or the absorbed is the unit.
+  However, with hyperduals there are more possibilities again due to terms vanishing during the
+  multiplication.
+\<close>
+lemma hyp_mult_cancel_right1 [iff]:
+  fixes a b :: "('a :: ring_1_no_zero_divisors) hyperdual"
+  shows "a = b * a \<longleftrightarrow> a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 b * Eps2 a = - Eps2 b * Eps1 a)"
+proof
+  assume mult: "a = b * a"
+  show "a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 b * Eps2 a = - Eps2 b * Eps1 a)"
+  proof (cases "Base a = 0")
+    case aT: True
+    then show ?thesis
+    proof (cases "Base b = 1")
+      \<comment> \<open>@{term "Base a = 0 \<and> Base b = 1"}\<close>
+      case bT: True
+      then show ?thesis
+        using aT mult add_cancel_right_right add_eq_0_iff[of "Eps1 b * Eps2 a"] times_hyperdual.simps(4)[of b a]
+        by simp
+    next
+      \<comment> \<open>@{term "Base a = 0 \<and> Base b \<noteq> 1"}\<close>
+      case bF: False
+      then show ?thesis
+        using aT mult by (simp, auto)
+    qed
+  next
+    case aF: False
+    then show ?thesis
+    proof (cases "Base b = 1")
+      \<comment> \<open>@{term "Base a \<noteq> 0 \<and> Base b = 1"}\<close>
+      case bT: True
+      then show ?thesis
+        using aF mult by (simp, auto)
+    next
+      \<comment> \<open>@{term "Base a \<noteq> 0 \<and> Base b \<noteq> 1"}\<close>
+      case bF: False
+      then show ?thesis
+        using aF mult by simp
+    qed
+  qed
+next
+  have "a = 0 \<Longrightarrow> a = b * a"
+   and "b = 1 \<Longrightarrow> a = b * a"
+    by simp_all
+  moreover have "Base a = 0 \<and> Base b = 1 \<and> Eps1 b * Eps2 a = - Eps2 b * Eps1 a \<Longrightarrow> a = b * a"
+    by simp
+  ultimately show "a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 b * Eps2 a = - Eps2 b * Eps1 a) \<Longrightarrow> a = b * a"
+    by blast
+qed
+lemma hyp_mult_cancel_right2 [iff]:
+  fixes a b :: "('a :: ring_1_no_zero_divisors) hyperdual"
+  shows "b * a = a \<longleftrightarrow> a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 b * Eps2 a = - Eps2 b * Eps1 a)"
+  using hyp_mult_cancel_right1 by smt
+
+lemma hyp_mult_cancel_left1 [iff]:
+  fixes a b :: "('a :: ring_1_no_zero_divisors) hyperdual"
+  shows "a = a * b \<longleftrightarrow> a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b)"
+proof
+  assume mult: "a = a * b"
+  show "a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b)"
+  proof (cases "Base a = 0")
+    case aT: True
+    then show ?thesis
+    proof (cases "Base b = 1")
+      \<comment> \<open>@{term "Base a = 0 \<and> Base b = 1"}\<close>
+      case bT: True
+      then show ?thesis
+        using aT mult add_cancel_right_right add_eq_0_iff[of "Eps1 a * Eps2 b"] times_hyperdual.simps(4)[of a b]
+        by simp
+    next
+      \<comment> \<open>@{term "Base a = 0 \<and> Base b \<noteq> 1"}\<close>
+      case bF: False
+      then show ?thesis
+        using aT mult by (simp, auto)
+    qed
+  next
+    case aF: False
+    then show ?thesis
+    proof (cases "Base b = 1")
+      \<comment> \<open>@{term "Base a \<noteq> 0 \<and> Base b = 1"}\<close>
+      case bT: True
+      then show ?thesis
+        using aF mult by (simp, auto)
+    next
+      \<comment> \<open>@{term "Base a \<noteq> 0 \<and> Base b \<noteq> 1"}\<close>
+      case bF: False
+      then show ?thesis
+        using aF mult by simp
+    qed
+  qed
+next
+  have "a = 0 \<Longrightarrow> a = a * b"
+   and "b = 1 \<Longrightarrow> a = a * b"
+    by simp_all
+  moreover have "Base a = 0 \<and> Base b = 1 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b \<Longrightarrow> a = a * b"
+    by simp
+  ultimately show "a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b) \<Longrightarrow> a = a * b"
+    by blast
+qed
+lemma hyp_mult_cancel_left2 [iff]:
+  fixes a b :: "('a :: ring_1_no_zero_divisors) hyperdual"
+  shows "a * b = a \<longleftrightarrow> a = 0 \<or> b = 1 \<or> (Base a = 0 \<and> Base b = 1 \<and> Eps1 a * Eps2 b = - Eps2 a * Eps1 b)"
+  using hyp_mult_cancel_left1 by smt
+
+subsection\<open>Multiplicative Inverse and Division\<close>
+
+text\<open>
+  If the components form a ring with a multiplicative inverse then so do the hyperduals.
+  The hyperdual inverse of @{term a} is defined as the solution to @{term "a * x = 1"}.
+  Hyperdual division is then multiplication by divisor's inverse.
+
+  Each component of the inverse has as denominator a power of the base component.
+  Therefore this inverse is only well defined for hyperdual numbers with non-zero base components.
+\<close>
+instantiation hyperdual :: ("{inverse, ring_1}") inverse
+begin
+
+primcorec inverse_hyperdual
+  where
+    "Base (inverse a) = 1 / Base a"
+  | "Eps1 (inverse a) = - Eps1 a / (Base a)^2"
+  | "Eps2 (inverse a) = - Eps2 a / (Base a)^2"
+  | "Eps12 (inverse a) = 2 * (Eps1 a * Eps2 a / (Base a)^3) - Eps12 a / (Base a)^2"
+
+primcorec divide_hyperdual
+  where
+    "Base (divide a b) = Base a / Base b"
+  | "Eps1 (divide a b) = (Eps1 a * Base b - Base a * Eps1 b) / ((Base b)^2)"
+  | "Eps2 (divide a b) = (Eps2 a * Base b - Base a * Eps2 b) / ((Base b)^2)"
+  | "Eps12 (divide a b) = (2 * Base a * Eps1 b * Eps2 b -
+                           Base a * Base b * Eps12 b -
+                           Eps1 a * Base b * Eps2 b -
+                           Eps2 a * Base b * Eps1 b +
+                           Eps12 a * ((Base b)^2)) / ((Base b)^3)"
+instance
+  by standard
+end
+
+text\<open>
+  Because hyperduals have non-trivial zero divisors, they do not form a division ring and so we
+  can't use the @{class division_ring} type class to establish properties of hyperdual division.
+  However, if the components form a division ring as well as a commutative ring, we can prove some
+  similar facts about hyperdual division inspired by @{class division_ring}.
+\<close>
+
+text\<open>Inverse is multiplicative inverse from both sides.\<close>
+lemma
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base a \<noteq> 0"
+  shows hyp_left_inverse [simp]: "inverse a * a = 1"
+    and hyp_right_inverse [simp]: "a * inverse a = 1"
+  by (simp_all add: assms power2_eq_square power3_eq_cube field_simps)
+
+text\<open>Division is multiplication by inverse.\<close>
+lemma hyp_divide_inverse:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "a / b = a * inverse b"
+  by (cases "Base b = 0" ; simp add: field_simps power2_eq_square power3_eq_cube)
+
+text\<open>Hyperdual inverse is zero when not well defined.\<close>
+lemma zero_base_zero_inverse:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base a = 0"
+  shows "inverse a = 0"
+  by (simp add: assms)
+
+lemma zero_inverse_zero_base:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "inverse a = 0"
+  shows "Base a = 0"
+  using assms right_inverse hyp_left_inverse by force
+
+lemma hyp_inverse_zero:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "(inverse a = 0) = (Base a = 0)"
+  using zero_base_zero_inverse[of a] zero_inverse_zero_base[of a] by blast
+
+text\<open>Inverse preserves invertibility.\<close>
+lemma hyp_invertible_inverse:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "(Base a = 0) = (Base (inverse a) = 0)"
+  by (safe, simp_all add: divide_inverse)
+
+text\<open>Inverse is the only number that satisfies the defining equation.\<close>
+lemma hyp_inverse_unique:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "a * b = 1"
+  shows "b = inverse a"
+proof -
+  have "Base a \<noteq> 0"
+    using assms one_hyperdual_def of_comp_simps zero_neq_one hyp_base_mult_eq_0_iff by smt
+  then show ?thesis
+    by (metis assms hyp_right_inverse mult.left_commute mult.right_neutral)
+qed
+
+text\<open>Multiplicative inverse commutes with additive inverse.\<close>
+lemma hyp_minus_inverse_comm:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "inverse (- a) = - inverse a"
+proof (cases "Base a = 0")
+  case True
+  then show ?thesis
+    by (simp add: zero_base_zero_inverse)
+next
+  case False
+  then show ?thesis
+    by (simp, simp add: nonzero_minus_divide_right)
+qed
+
+text\<open>
+  Inverse is an involution (only) where well defined.
+  Counter-example for non-invertible is @{term "Hyperdual 0 0 0 0"} with inverse
+  @{term "Hyperdual 0 0 0 0"} which then inverts to @{term "Hyperdual 0 0 0 0"}.
+\<close>
+lemma hyp_inverse_involution:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base a \<noteq> 0"
+  shows "inverse (inverse a) = a"
+  by (metis assms hyp_inverse_unique hyp_right_inverse mult.commute)
+
+lemma inverse_inverse_neq_Ex:
+  "\<exists>a :: ('a :: {inverse, comm_ring_1, division_ring}) hyperdual . inverse (inverse a) \<noteq> a"
+proof -
+  have "\<exists>a :: 'a hyperdual . Base a = 0 \<and> a \<noteq> 0"
+    by (metis (full_types) hyperdual.sel(1) hyperdual.sel(4) zero_neq_one)
+  moreover have "\<And>a :: 'a hyperdual . (Base a = 0 \<and> a \<noteq> 0) \<Longrightarrow> (inverse (inverse a) \<noteq> a)"
+    using hyp_inverse_zero hyp_invertible_inverse by smt
+  ultimately show ?thesis
+    by blast
+qed
+
+text\<open>
+  Inverses of equal invertible numbers are equal.
+  This includes the other direction by inverse preserving invertibility and being an involution.
+
+  From a different point of view, inverse is injective on invertible numbers.
+  The other direction for is again by inverse preserving invertibility and being an involution.
+\<close>
+lemma hyp_inverse_injection:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base a \<noteq> 0"
+      and "Base b \<noteq> 0"
+    shows "(inverse a = inverse b) = (a = b)"
+  by (metis assms hyp_inverse_involution)
+
+text\<open>One is its own inverse.\<close>
+lemma hyp_inverse_1 [simp]:
+  "inverse (1 :: ('a :: {inverse, comm_ring_1, division_ring}) hyperdual) = 1"
+  using hyp_inverse_unique mult.left_neutral by metis
+
+text\<open>Inverse distributes over multiplication (even when not well defined).\<close>
+lemma hyp_inverse_mult_distrib:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "inverse (a * b) = inverse b * inverse a"
+proof (cases "Base a = 0 \<or> Base b = 0")
+  case True
+  then show ?thesis
+    by (metis hyp_base_mult_eq_0_iff mult_zero_left mult_zero_right zero_base_zero_inverse)
+next
+  case False
+  then have "a * (b * inverse b) * inverse a = 1"
+    by simp
+  then have "a * b * (inverse b * inverse a) = 1"
+    by (simp only: mult.assoc)
+  thus ?thesis
+    using hyp_inverse_unique[of "a * b" "(inverse b * inverse a)"] by simp
+qed
+
+text\<open>We derive expressions for addition and subtraction of inverses.\<close>
+lemma hyp_inverse_add:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base a \<noteq> 0"
+      and "Base b \<noteq> 0"
+    shows "inverse a + inverse b = inverse a * (a + b) * inverse b"
+  by (simp add: assms distrib_left mult.commute semiring_normalization_rules(18) add.commute)
+
+lemma hyp_inverse_diff:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes a: "Base a \<noteq> 0"
+      and b: "Base b \<noteq> 0"
+    shows "inverse a - inverse b = inverse a * (b - a) * inverse b"
+proof -
+  have x: "inverse a - inverse b = inverse b * (inverse a * b - 1)"
+    by (simp add: b mult.left_commute right_diff_distrib')
+  show ?thesis
+    by (simp add: x a mult.commute right_diff_distrib')
+qed
+
+text\<open>Division is one only when dividing by self.\<close>
+lemma hyp_divide_self:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base b \<noteq> 0"
+    shows "a / b = 1 \<longleftrightarrow> a = b"
+  by (metis assms hyp_divide_inverse hyp_inverse_unique hyp_right_inverse mult.commute)
+
+text\<open>Taking inverse is the same as division of one, even when not invertible.\<close>
+lemma hyp_inverse_divide_1 [divide_simps]:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "inverse a = 1 / a"
+  by (simp add: hyp_divide_inverse)
+
+text\<open>Division distributes over addition and subtraction.\<close>
+lemma hyp_add_divide_distrib:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "(a + b) / c = a/c + b/c"
+  by (simp add: distrib_right hyp_divide_inverse)
+
+lemma hyp_diff_divide_distrib:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "(a - b) / c = a / c - b / c"
+  by (simp add: left_diff_distrib hyp_divide_inverse)
+
+text\<open>Multiplication associates with division.\<close>
+lemma hyp_times_divide_assoc [simp]:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "a * (b / c) = (a * b) / c"
+  by (simp add: hyp_divide_inverse mult.assoc)
+
+text\<open>Additive inverse commutes with division, because it is multiplication by inverse.\<close>
+lemma hyp_divide_minus_left [simp]:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "(-a) / b = - (a / b)"
+  by (simp add: hyp_divide_inverse)
+
+lemma hyp_divide_minus_right [simp]:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "a / (-b) = - (a / b)"
+  by (simp add: hyp_divide_inverse hyp_minus_inverse_comm)
+
+text\<open>Additive inverses on both sides of division cancel out.\<close>
+lemma hyp_minus_divide_minus:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "(-a) / (-b) = a / b"
+  by simp
+
+text\<open>We can multiply both sides of equations by an invertible denominator.\<close>
+lemma hyp_denominator_eliminate [divide_simps]:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base c \<noteq> 0"
+    shows "a = b / c \<longleftrightarrow> a * c = b"
+  by (metis assms hyp_divide_self hyp_times_divide_assoc mult.commute mult.right_neutral)
+
+text\<open>We can move addition and subtraction to a common denominator in the following ways:\<close>
+lemma hyp_add_divide_eq_iff:
+  fixes x y z :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base z \<noteq> 0"
+  shows "x + y / z = (x * z + y) / z"
+  by (metis assms hyp_add_divide_distrib hyp_denominator_eliminate)
+
+text\<open>Result of division by non-invertible number is not invertible.\<close>
+lemma hyp_divide_base_zero [simp]:
+  fixes a b :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  assumes "Base b = 0"
+    shows "Base (a / b) = 0"
+  by (simp add: assms)
+
+text\<open>Division of self is 1 when invertible, 0 otherwise.\<close>
+lemma hyp_divide_self_if [simp]:
+  fixes a :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "a / a = (if Base a = 0 then 0 else 1)"
+  by (metis hyp_divide_inverse zero_base_zero_inverse hyp_divide_self mult_zero_right)
+
+text\<open>Repeated division is division by product of the denominators.\<close>
+lemma hyp_denominators_merge:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "(a / b) / c = a / (c * b)"
+  using hyp_inverse_mult_distrib[of c b]
+  by (simp add: hyp_divide_inverse mult.assoc)
+
+text\<open>Finally, we derive general simplifications for division with addition and subtraction.\<close>
+lemma hyp_add_divide_eq_if_simps [divide_simps]:
+  fixes a b z :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "a + b / z = (if Base z = 0 then a else (a * z + b) / z)"
+    and "a / z + b = (if Base z = 0 then b else (a + b * z) / z)"
+    and "- (a / z) + b = (if Base z = 0 then b else (-a + b * z) / z)"
+    and "a - b / z = (if Base z = 0 then a else (a * z - b) / z)"
+    and "a / z - b = (if Base z = 0 then -b else (a - b * z) / z)"
+    and "- (a / z) - b = (if Base z = 0 then -b else (- a - b * z) / z)"
+  by (simp_all add: algebra_simps hyp_add_divide_eq_iff hyp_divide_inverse zero_base_zero_inverse)
+
+lemma hyp_divide_eq_eq [divide_simps]:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "b / c = a \<longleftrightarrow> (if Base c \<noteq> 0 then b = a * c else a = 0)"
+  by (metis hyp_divide_inverse hyp_denominator_eliminate mult_not_zero zero_base_zero_inverse)
+
+lemma hyp_eq_divide_eq [divide_simps]:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "a = b / c \<longleftrightarrow> (if Base c \<noteq> 0 then a * c = b else a = 0)"
+  by (metis hyp_divide_eq_eq)
+
+lemma hyp_minus_divide_eq_eq [divide_simps]:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "- (b / c) = a \<longleftrightarrow> (if Base c \<noteq> 0 then - b = a * c else a = 0)"
+  by (metis hyp_divide_minus_left hyp_eq_divide_eq)
+
+lemma hyp_eq_minus_divide_eq [divide_simps]:
+  fixes a b c :: "('a :: {inverse, comm_ring_1, division_ring}) hyperdual"
+  shows "a = - (b / c) \<longleftrightarrow> (if Base c \<noteq> 0 then a * c = - b else a = 0)"
+  by (metis hyp_minus_divide_eq_eq)
+
+subsection \<open>Real Scaling, Real Vector and Real Algebra\<close>
+
+text\<open>
+  If the components can be scaled by real numbers then so can the hyperduals.
+  We define the scaling pointwise.
+\<close>
+instantiation hyperdual :: (scaleR) scaleR
+begin
+
+primcorec scaleR_hyperdual
+  where
+    "Base (f *\<^sub>R x) = f *\<^sub>R Base x"
+  | "Eps1 (f *\<^sub>R x) = f *\<^sub>R Eps1 x"
+  | "Eps2 (f *\<^sub>R x) = f *\<^sub>R Eps2 x"
+  | "Eps12 (f *\<^sub>R x) = f *\<^sub>R Eps12 x"
+
+instance
+  by standard
+end
+
+text\<open>If the components form a real vector space then so do the hyperduals.\<close>
+instance hyperdual :: (real_vector) real_vector
+  by standard (simp_all add: algebra_simps)
+
+text\<open>If the components form a real algebra then so do the hyperduals\<close>
+instance hyperdual :: (real_algebra_1) real_algebra_1
+  by standard (simp_all add: algebra_simps)
+
+text\<open>
+  If the components are reals then @{const of_real} matches our embedding @{const of_comp}, and
+  @{const scaleR} matches our scalar product @{const scaleH}.
+\<close>
+lemma "of_real = of_comp"
+  by (standard, simp add: of_real_def)
+
+lemma scaleR_eq_scale:
+  "(*\<^sub>R) = (*\<^sub>H)"
+  by (standard, standard, simp)
+
+text\<open>Hyperdual scalar product @{const scaleH} is compatible with @{const scaleR}.\<close>
+lemma scaleH_scaleR:
+  fixes a :: "'a :: real_algebra_1"
+    and b :: "'a hyperdual"
+  shows "(f *\<^sub>R a) *\<^sub>H b = f *\<^sub>R a *\<^sub>H b"
+    and "a *\<^sub>H f *\<^sub>R b = f *\<^sub>R a *\<^sub>H b"
+  by simp_all
+
+subsection\<open>Real Inner Product and Real-Normed Vector Space\<close>
+
+text\<open>
+  We now take a closer look at hyperduals as a real vector space.
+
+  If the components form a real inner product space then we can define one on the hyperduals as the
+  sum of componentwise inner products.
+  The norm is then defined as the square root of that inner product.
+  We define signum, distance, uniformity and openness similarly as they are defined for complex
+  numbers.
+\<close>
+instantiation hyperdual :: (real_inner) real_inner
+begin
+
+definition inner_hyperdual :: "'a hyperdual \<Rightarrow> 'a hyperdual \<Rightarrow> real"
+  where "x \<bullet> y = Base x \<bullet> Base y + Eps1 x \<bullet> Eps1 y + Eps2 x \<bullet> Eps2 y + Eps12 x \<bullet> Eps12 y"
+
+definition norm_hyperdual :: "'a hyperdual \<Rightarrow> real"
+  where "norm_hyperdual x = sqrt (x \<bullet> x)"
+
+definition sgn_hyperdual :: "'a hyperdual \<Rightarrow> 'a hyperdual"
+  where "sgn_hyperdual x = x /\<^sub>R norm x"
+
+definition dist_hyperdual :: "'a hyperdual \<Rightarrow> 'a hyperdual \<Rightarrow> real"
+  where "dist_hyperdual a b = norm(a - b)"
+
+definition uniformity_hyperdual :: "('a hyperdual \<times> 'a hyperdual) filter"
+  where "uniformity_hyperdual = (INF e\<in>{0 <..}. principal {(x, y). dist x y < e})"
+
+definition open_hyperdual :: "('a hyperdual) set \<Rightarrow> bool"
+  where  "open_hyperdual U \<longleftrightarrow> (\<forall>x\<in>U. eventually (\<lambda>(x', y). x' = x \<longrightarrow> y \<in> U) uniformity)"
+
+instance
+proof
+  fix x y z :: "'a hyperdual"
+  fix U :: "('a hyperdual) set"
+  fix r :: real
+
+  show "dist x y = norm (x - y)"
+    by (simp add: dist_hyperdual_def)
+  show "sgn x = x /\<^sub>R norm x"
+    by (simp add: sgn_hyperdual_def)
+  show "uniformity = (INF e\<in>{0<..}. principal {(x :: 'a hyperdual, y). dist x y < e})"
+    using uniformity_hyperdual_def by blast
+  show "open U = (\<forall>x\<in>U. \<forall>\<^sub>F (x', y) in uniformity. x' = x \<longrightarrow> y \<in> U)"
+    using open_hyperdual_def by blast
+  show  "x \<bullet> y = y \<bullet> x"
+    by (simp add: inner_hyperdual_def inner_commute)
+  show "(x + y) \<bullet> z = x \<bullet> z + y \<bullet> z"
+    and "r *\<^sub>R x \<bullet> y = r * (x \<bullet> y)"
+    and "0 \<le> x \<bullet> x"
+    and "norm x = sqrt (x \<bullet> x)"
+    by (simp_all add: inner_hyperdual_def norm_hyperdual_def algebra_simps)
+  show "(x \<bullet> x = 0) = (x = 0)"
+  proof
+    assume "x \<bullet> x = 0"
+    then have "Base x \<bullet> Base x + Eps1 x \<bullet> Eps1 x + Eps2 x \<bullet> Eps2 x + Eps12 x \<bullet> Eps12 x = 0"
+      by (simp add: inner_hyperdual_def)
+    then have "Base x = 0 \<and> Eps1 x = 0 \<and> Eps2 x = 0 \<and> Eps12 x = 0"
+      using inner_gt_zero_iff inner_ge_zero add_nonneg_nonneg
+      by smt
+    then show "x = 0"
+      by simp
+  next
+    assume "x = 0"
+    then show "x \<bullet> x = 0"
+      by (simp add: inner_hyperdual_def)
+  qed
+qed
+end
+
+text\<open>
+  We then show that with this norm hyperduals with components that form a real normed algebra do not
+  themselves form a normed algebra, by counter-example to the assumption that class adds.
+\<close>
+(* Components must be real_inner for hyperduals to have a norm *)
+(* Components must be real_normed_algebra_1 to know |1| = 1, because we need some element with known
+    norm of its square. Otherwise we would need the precise definition of the norm. *)
+lemma not_normed_algebra:
+  shows "\<not>(\<forall>x y :: ('a :: {real_normed_algebra_1, real_inner}) hyperdual . norm (x * y) \<le> norm x * norm y)"
+proof -
+  have "norm (Hyperdual (1 :: 'a) 1 1 1) = 2"
+    by (simp add: norm_hyperdual_def inner_hyperdual_def dot_square_norm)
+  moreover have "(Hyperdual (1 :: 'a) 1 1 1) * (Hyperdual 1 1 1 1) = Hyperdual 1 2 2 4"
+    by (simp add: hyperdual.expand)
+  moreover have "norm (Hyperdual (1 :: 'a) 2 2 4) > 4"
+    by (simp add: norm_hyperdual_def inner_hyperdual_def dot_square_norm)
+  ultimately have "\<exists>x y :: 'a hyperdual . norm (x * y) > norm x * norm y"
+    by (smt power2_eq_square real_sqrt_four real_sqrt_pow2)
+  then show ?thesis
+    by (simp add: not_le)
+qed
+
+subsection\<open>Euclidean Space\<close>
+
+text\<open>
+  Next we define a basis for the space, consisting of four elements one for each component with \<open>1\<close>
+  in the relevant component and \<open>0\<close> elsewhere.
+\<close>
+definition ba :: "('a :: zero_neq_one) hyperdual"
+  where "ba = Hyperdual 1 0 0 0"
+definition e1 :: "('a :: zero_neq_one) hyperdual"
+  where "e1 = Hyperdual 0 1 0 0"
+definition e2 :: "('a :: zero_neq_one) hyperdual"
+  where "e2 = Hyperdual 0 0 1 0"
+definition e12 :: "('a :: zero_neq_one) hyperdual"
+  where "e12 = Hyperdual 0 0 0 1"
+
+lemmas hyperdual_bases = ba_def e1_def e2_def e12_def
+
+text\<open>
+  Using the constructor @{const Hyperdual} is equivalent to using the linear combination with
+  coefficients the relevant arguments.
+\<close>
+lemma Hyperdual_eq:
+  fixes a b c d :: "'a :: ring_1"
+  shows "Hyperdual a b c d = a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12"
+  by (simp add: hyperdual_bases)
+
+text\<open>Projecting from the combination returns the relevant coefficient:\<close>
+lemma hyperdual_comb_sel [simp]:
+  fixes a b c d :: "'a :: ring_1"
+  shows "Base (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) = a"
+    and "Eps1 (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) = b"
+    and "Eps2 (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) = c"
+    and "Eps12 (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) = d"
+  using Hyperdual_eq hyperdual.sel by metis+
+
+text\<open>Any hyperdual number is a linear combination of these four basis elements.\<close>
+lemma hyperdual_linear_comb:
+  fixes x :: "('a :: ring_1) hyperdual"
+  obtains a b c d :: 'a where "x = a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12"
+  using hyperdual.exhaust Hyperdual_eq by metis
+
+text\<open>
+  The linear combination expressing any hyperdual number has as coefficients the projections of
+  that number onto the relevant basis element.
+\<close>
+lemma hyperdual_eq:
+  fixes x :: "('a :: ring_1) hyperdual"
+  shows "x = Base x *\<^sub>H ba + Eps1 x *\<^sub>H e1 + Eps2 x *\<^sub>H e2 + Eps12 x *\<^sub>H e12"
+  using Hyperdual_eq hyperdual.collapse by smt
+
+text\<open>Equality of hyperduals as linear combinations is equality of corresponding components.\<close>
+lemma hyperdual_eq_parts_cancel [simp]:
+  fixes a b c d :: "'a :: ring_1"
+  shows "(a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12 = a' *\<^sub>H ba + b' *\<^sub>H e1 + c' *\<^sub>H e2 + d' *\<^sub>H e12) \<equiv>
+         (a = a' \<and> b = b' \<and> c = c' \<and> d = d')"
+  by (smt Hyperdual_eq hyperdual.inject)
+
+lemma scaleH_cancel [simp]:
+  fixes a b :: "'a :: ring_1"
+  shows "(a *\<^sub>H ba = b *\<^sub>H ba) \<equiv> (a = b)"
+    and "(a *\<^sub>H e1 = b *\<^sub>H e1) \<equiv> (a = b)"
+    and "(a *\<^sub>H e2 = b *\<^sub>H e2) \<equiv> (a = b)"
+    and "(a *\<^sub>H e12 = b *\<^sub>H e12) \<equiv> (a = b)"
+  by (auto simp add: hyperdual_bases)+
+
+text\<open>We can now also show that the multiplication we use indeed has the hyperdual units nilpotent.\<close>
+(* The components must have multiplication and addition for this to be defined and have an absorbing
+   zero to prove it. *)
+lemma epsilon_squares [simp]:
+  "(e1 :: ('a :: ring_1) hyperdual) * e1 = 0"
+  "(e2 :: ('a :: ring_1) hyperdual) * e2 = 0"
+  "(e12 :: ('a :: ring_1) hyperdual) * e12 = 0"
+  by (simp_all add: hyperdual_bases)
+
+text\<open>However none of the hyperdual units is zero.\<close>
+lemma hyperdual_bases_nonzero [simp]:
+  "ba \<noteq> 0"
+  "e1 \<noteq> 0"
+  "e2 \<noteq> 0"
+  "e12 \<noteq> 0"
+  by (simp_all add: hyperdual_bases)
+
+text\<open>Hyperdual units are orthogonal.\<close>
+lemma hyperdual_bases_ortho [simp]:
+  "(ba :: ('a :: {real_inner,zero_neq_one}) hyperdual) \<bullet> e1 = 0"
+  "(ba :: ('a :: {real_inner,zero_neq_one}) hyperdual) \<bullet> e2 = 0"
+  "(ba :: ('a :: {real_inner,zero_neq_one}) hyperdual) \<bullet> e12 = 0"
+  "(e1 :: ('a :: {real_inner,zero_neq_one}) hyperdual) \<bullet> e2 = 0"
+  "(e1 :: ('a :: {real_inner,zero_neq_one}) hyperdual) \<bullet> e12 = 0"
+  "(e2 :: ('a :: {real_inner,zero_neq_one}) hyperdual) \<bullet> e12 = 0"
+  by (simp_all add: hyperdual_bases inner_hyperdual_def)
+
+text\<open>Hyperdual units of norm equal to 1.\<close>
+lemma hyperdual_bases_norm [simp]:
+  "(ba :: ('a :: {real_inner,real_normed_algebra_1}) hyperdual) \<bullet> ba = 1"
+  "(e1 :: ('a :: {real_inner,real_normed_algebra_1}) hyperdual) \<bullet> e1 = 1"
+  "(e2 :: ('a :: {real_inner,real_normed_algebra_1}) hyperdual) \<bullet> e2 = 1"
+  "(e12 :: ('a :: {real_inner,real_normed_algebra_1}) hyperdual) \<bullet> e12 = 1"
+  by (simp_all add: hyperdual_bases inner_hyperdual_def norm_eq_1[symmetric])
+
+text\<open>We can also express earlier operations in terms of the linear combination.\<close>
+lemma add_hyperdual_parts:
+  fixes a b c d :: "'a :: ring_1"
+  shows "(a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) + (a' *\<^sub>H ba + b' *\<^sub>H e1 + c' *\<^sub>H e2 + d' *\<^sub>H e12) =
+         (a + a') *\<^sub>H ba + (b + b') *\<^sub>H e1 + (c + c') *\<^sub>H e2 + (d + d') *\<^sub>H e12"
+  by (simp add: scaleH_add(1))
+
+lemma times_hyperdual_parts:
+  fixes a b c d :: "'a :: ring_1"
+  shows "(a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) * (a' *\<^sub>H ba + b' *\<^sub>H e1 + c' *\<^sub>H e2 + d' *\<^sub>H e12) =
+         (a * a') *\<^sub>H ba + (a * b' + b * a') *\<^sub>H e1 + (a * c' + c * a') *\<^sub>H e2 + (a * d' + b * c' + c * b' + d * a') *\<^sub>H e12"
+  by (simp add: hyperdual_bases)
+
+lemma inverse_hyperdual_parts:
+  fixes a b c d :: "'a :: {inverse,ring_1}"
+  shows "inverse (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) =
+         (1 / a) *\<^sub>H ba + (- b / a ^ 2) *\<^sub>H e1 + (- c / a ^ 2) *\<^sub>H e2 + (2 * (b * c / a ^ 3) - d / a ^ 2) *\<^sub>H e12"
+  by (simp add: hyperdual_bases)
+
+text\<open>
+  Next we show that hyperduals form a euclidean space with the help of the basis we defined earlier
+  and the above inner product if the component is an instance of @{class euclidean_space} and
+  @{class real_algebra_1}.
+  The basis of this space is each of the basis elements we defined scaled by each of the basis
+  elements of the component type, representing the expansion of the space for each component of the
+  hyperdual numbers.
+\<close>
+instantiation hyperdual :: ("{euclidean_space, real_algebra_1}") euclidean_space
+begin
+
+definition Basis_hyperdual :: "('a hyperdual) set"
+  where "Basis = (\<Union>i\<in>{ba,e1,e2,e12}. (\<lambda>u. u *\<^sub>H i) ` Basis)"
+
+instance
+proof
+  fix x y z :: "'a hyperdual"
+  show "(Basis :: ('a hyperdual) set) \<noteq> {}"
+    and "finite (Basis :: ('a hyperdual) set)"
+    by (simp_all add: Basis_hyperdual_def)
+  show "x \<in> Basis \<Longrightarrow> y \<in> Basis \<Longrightarrow> x \<bullet> y = (if x = y then 1 else 0)"
+    unfolding Basis_hyperdual_def inner_hyperdual_def hyperdual_bases
+    by (auto dest: inner_not_same_Basis)
+  show "(\<forall>u\<in>Basis. x \<bullet> u = 0) = (x = 0)"
+    by (auto simp add: Basis_hyperdual_def ball_Un inner_hyperdual_def hyperdual_bases euclidean_all_zero_iff)
+qed
+end
+
+subsection\<open>Bounded Linear Projections\<close>
+
+text\<open>Now we can show that each projection to a basis element is a bounded linear map.\<close>
+lemma bounded_linear_Base: "bounded_linear Base"
+proof
+  show "\<And>b1 b2. Base (b1 + b2) = Base b1 + Base b2"
+    by simp
+  show "\<And>r b. Base (r *\<^sub>R b) = r *\<^sub>R Base b"
+    by simp
+  have "\<forall>x. norm (Base x) \<le> norm x"
+    by (simp add: inner_hyperdual_def norm_eq_sqrt_inner)
+  then show "\<exists>K. \<forall>x. norm (Base x) \<le> norm x * K"
+    by (metis mult.commute mult.left_neutral)
+qed
+lemma bounded_linear_Eps1: "bounded_linear Eps1"
+proof
+  show "\<And>b1 b2. Eps1 (b1 + b2) = Eps1 b1 + Eps1 b2"
+    by simp
+  show "\<And>r b. Eps1 (r *\<^sub>R b) = r *\<^sub>R Eps1 b"
+    by simp
+  have "\<forall>x. norm (Eps1 x) \<le> norm x"
+    by (simp add: inner_hyperdual_def norm_eq_sqrt_inner)
+  then show "\<exists>K. \<forall>x. norm (Eps1 x) \<le> norm x * K"
+    by (metis mult.commute mult.left_neutral)
+qed
+lemma bounded_linear_Eps2: "bounded_linear Eps2"
+proof
+  show "\<And>b1 b2. Eps2 (b1 + b2) = Eps2 b1 + Eps2 b2"
+    by simp
+  show "\<And>r b. Eps2 (r *\<^sub>R b) = r *\<^sub>R Eps2 b"
+    by simp
+  have "\<forall>x. norm (Eps2 x) \<le> norm x"
+    by (simp add: inner_hyperdual_def norm_eq_sqrt_inner)
+  then show "\<exists>K. \<forall>x. norm (Eps2 x) \<le> norm x * K"
+    by (metis mult.commute mult.left_neutral)
+qed
+lemma bounded_linear_Eps12: "bounded_linear Eps12"
+proof
+  show "\<And>b1 b2. Eps12 (b1 + b2) = Eps12 b1 + Eps12 b2"
+    by simp
+  show "\<And>r b. Eps12 (r *\<^sub>R b) = r *\<^sub>R Eps12 b"
+    by simp
+  have "\<forall>x. norm (Eps12 x) \<le> norm x"
+    by (simp add: inner_hyperdual_def norm_eq_sqrt_inner)
+  then show "\<exists>K. \<forall>x. norm (Eps12 x) \<le> norm x * K"
+    by (metis mult.commute mult.left_neutral)
+qed
+
+text\<open>
+  This bounded linearity gives us a range of useful theorems about limits, convergence and
+  derivatives of these projections.
+\<close>
+lemmas tendsto_Base = bounded_linear.tendsto[OF bounded_linear_Base]
+lemmas tendsto_Eps1 = bounded_linear.tendsto[OF bounded_linear_Eps1]
+lemmas tendsto_Eps2 = bounded_linear.tendsto[OF bounded_linear_Eps2]
+lemmas tendsto_Eps12 = bounded_linear.tendsto[OF bounded_linear_Eps12]
+
+lemmas has_derivative_Base = bounded_linear.has_derivative[OF bounded_linear_Base]
+lemmas has_derivative_Eps1 = bounded_linear.has_derivative[OF bounded_linear_Eps1]
+lemmas has_derivative_Eps2 = bounded_linear.has_derivative[OF bounded_linear_Eps2]
+lemmas has_derivative_Eps12 = bounded_linear.has_derivative[OF bounded_linear_Eps12]
+
+subsection\<open>Convergence\<close>
+lemma inner_mult_le_mult_inner:
+  fixes a b :: "'a :: {real_inner,real_normed_algebra}"
+  shows "((a * b) \<bullet> (a * b)) \<le> (a \<bullet> a) * (b \<bullet> b)"
+  by (metis real_sqrt_le_iff norm_eq_sqrt_inner real_sqrt_mult norm_mult_ineq)
+
+lemma bounded_bilinear_scaleH:
+  "bounded_bilinear ((*\<^sub>H) :: ('a :: {real_normed_algebra_1, real_inner}) \<Rightarrow> 'a hyperdual \<Rightarrow> 'a hyperdual)"
+proof (auto simp add: bounded_bilinear_def scaleH_add scaleH_scaleR del:exI intro!:exI)
+  fix a :: 'a
+  and b :: "'a hyperdual"
+  have "norm (a *\<^sub>H b) = sqrt ((a * Base b) \<bullet> (a * Base b) + (a * Eps1 b) \<bullet> (a * Eps1 b) + (a * Eps2 b) \<bullet> (a * Eps2 b) + (a * Eps12 b) \<bullet> (a * Eps12 b))"
+    by (simp add: norm_eq_sqrt_inner inner_hyperdual_def)
+  moreover have "norm a * norm b = sqrt (a \<bullet> a * (Base b \<bullet> Base b + Eps1 b \<bullet> Eps1 b + Eps2 b \<bullet> Eps2 b + Eps12 b \<bullet> Eps12 b))"
+    by (simp add: norm_eq_sqrt_inner inner_hyperdual_def real_sqrt_mult)
+  moreover have "sqrt ((a * Base b) \<bullet> (a * Base b) + (a * Eps1 b) \<bullet> (a * Eps1 b) + (a * Eps2 b) \<bullet> (a * Eps2 b) + (a * Eps12 b) \<bullet> (a * Eps12 b)) \<le>
+                 sqrt (a \<bullet> a * (Base b \<bullet> Base b + Eps1 b \<bullet> Eps1 b + Eps2 b \<bullet> Eps2 b + Eps12 b \<bullet> Eps12 b))"
+    by (simp add: distrib_left add_mono inner_mult_le_mult_inner)
+  ultimately show "norm (a *\<^sub>H b) \<le> norm a * norm b * 1"
+    by simp
+qed
+
+lemmas tendsto_scaleH = bounded_bilinear.tendsto[OF bounded_bilinear_scaleH]
+
+text\<open>
+  We describe how limits behave for general hyperdual-valued functions.
+
+  First we prove that we can go from convergence of the four component functions to the convergence
+  of the hyperdual-valued function whose components they define.
+\<close>
+lemma tendsto_Hyperdual:
+  fixes f :: "'a \<Rightarrow> ('b :: {real_normed_algebra_1, real_inner})"
+  assumes "(f \<longlongrightarrow> a) F"
+      and "(g \<longlongrightarrow> b) F"
+      and "(h \<longlongrightarrow> c) F"
+      and "(i \<longlongrightarrow> d) F"
+    shows "((\<lambda>x. Hyperdual (f x) (g x) (h x) (i x)) \<longlongrightarrow> Hyperdual a b c d) F"
+proof -
+  have "((\<lambda>x. (f x) *\<^sub>H ba) \<longlongrightarrow> a *\<^sub>H ba) F"
+       "((\<lambda>x. (g x) *\<^sub>H e1) \<longlongrightarrow> b *\<^sub>H e1) F"
+       "((\<lambda>x. (h x) *\<^sub>H e2) \<longlongrightarrow> c *\<^sub>H e2) F"
+       "((\<lambda>x. (i x) *\<^sub>H e12) \<longlongrightarrow> d *\<^sub>H e12) F"
+    by (rule tendsto_scaleH[OF _ tendsto_const], rule assms)+
+  then have "((\<lambda>x. (f x) *\<^sub>H ba + (g x) *\<^sub>H e1 + (h x) *\<^sub>H e2 + (i x) *\<^sub>H e12) \<longlongrightarrow>
+             a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) F"
+    by (rule tendsto_add[OF tendsto_add[OF tendsto_add]] ; assumption)
+  then show ?thesis
+    by (simp add: Hyperdual_eq)
+qed
+
+text\<open>
+  Next we complete the equivalence by proving the other direction, from convergence of a
+  hyperdual-valued function to the convergence of the projected component functions.
+\<close>
+lemma tendsto_hyperdual_iff:
+  fixes f :: "'a \<Rightarrow> ('b :: {real_normed_algebra_1, real_inner}) hyperdual"
+  shows "(f \<longlongrightarrow> x) F \<longleftrightarrow>
+    ((\<lambda>x. Base (f x)) \<longlongrightarrow> Base x) F \<and>
+    ((\<lambda>x. Eps1 (f x)) \<longlongrightarrow> Eps1 x) F \<and>
+    ((\<lambda>x. Eps2 (f x)) \<longlongrightarrow> Eps2 x) F \<and>
+    ((\<lambda>x. Eps12 (f x)) \<longlongrightarrow> Eps12 x) F"
+proof safe
+  assume "(f \<longlongrightarrow> x) F"
+  then show "((\<lambda>x. Base (f x)) \<longlongrightarrow> Base x) F"
+        and "((\<lambda>x. Eps1 (f x)) \<longlongrightarrow> Eps1 x) F"
+        and "((\<lambda>x. Eps2 (f x)) \<longlongrightarrow> Eps2 x) F"
+        and "((\<lambda>x. Eps12 (f x)) \<longlongrightarrow> Eps12 x) F"
+    by (simp_all add: tendsto_Base tendsto_Eps1 tendsto_Eps2 tendsto_Eps12)
+next
+  assume "((\<lambda>x. Base (f x)) \<longlongrightarrow> Base x) F"
+     and "((\<lambda>x. Eps1 (f x)) \<longlongrightarrow> Eps1 x) F"
+     and "((\<lambda>x. Eps2 (f x)) \<longlongrightarrow> Eps2 x) F"
+     and "((\<lambda>x. Eps12 (f x)) \<longlongrightarrow> Eps12 x) F"
+  then show "(f \<longlongrightarrow> x) F"
+    using tendsto_Hyperdual[of "\<lambda>x. Base (f x)" "Base x" F "\<lambda>x. Eps1 (f x)" "Eps1 x" "\<lambda>x. Eps2 (f x)" "Eps2 x" "\<lambda>x. Eps12 (f x)" "Eps12 x"]
+    by simp
+qed
+
+subsection\<open>Derivatives\<close>
+
+text\<open>
+  We describe how derivatives of hyperdual-valued functions behave.
+  Due to hyperdual numbers not forming a normed field, the derivative relation we must use is the
+  general Fréchet derivative @{const has_derivative}.
+
+  The left to right implication of the following equivalence is easily proved by the known
+  derivative behaviour of the projections.
+  The other direction is more difficult, because we have to construct the two requirements of the
+  @{const has_derivative} relation, the limit and the bounded linearity of the derivative.
+  While the limit is simple to construct from the component functions by previous lemma, the bounded
+  linearity is more involved.
+\<close>
+(* The derivative of a hyperdual function is composed of the derivatives of its coefficient functions *)
+lemma has_derivative_hyperdual_iff:
+  fixes f :: "('a :: real_normed_vector) \<Rightarrow> ('b :: {real_normed_algebra_1, real_inner}) hyperdual"
+  shows "(f has_derivative Df) F \<longleftrightarrow>
+    ((\<lambda>x. Base (f x)) has_derivative (\<lambda>x. Base (Df x))) F \<and>
+    ((\<lambda>x. Eps1 (f x)) has_derivative (\<lambda>x. Eps1 (Df x))) F \<and>
+    ((\<lambda>x. Eps2 (f x)) has_derivative (\<lambda>x. Eps2 (Df x))) F \<and>
+    ((\<lambda>x. Eps12 (f x)) has_derivative (\<lambda>x. Eps12 (Df x))) F"
+proof safe
+  \<comment> \<open>Left to Right\<close>
+  assume assm: "(f has_derivative Df) F"
+  show "((\<lambda>x. Base (f x)) has_derivative (\<lambda>x. Base (Df x))) F"
+    using assm has_derivative_Base by blast
+  show "((\<lambda>x. Eps1 (f x)) has_derivative (\<lambda>x. Eps1 (Df x))) F"
+    using assm has_derivative_Eps1 by blast
+  show "((\<lambda>x. Eps2 (f x)) has_derivative (\<lambda>x. Eps2 (Df x))) F"
+    using assm has_derivative_Eps2 by blast
+  show "((\<lambda>x. Eps12 (f x)) has_derivative (\<lambda>x. Eps12 (Df x))) F"
+    using assm has_derivative_Eps12 by blast
+next
+  \<comment> \<open>Right to Left\<close>
+  assume assm:
+    "((\<lambda>x. Base (f x)) has_derivative (\<lambda>x. Base (Df x))) F"
+    "((\<lambda>x. Eps1 (f x)) has_derivative (\<lambda>x. Eps1 (Df x))) F"
+    "((\<lambda>x. Eps2 (f x)) has_derivative (\<lambda>x. Eps2 (Df x))) F"
+    "((\<lambda>x. Eps12 (f x)) has_derivative (\<lambda>x. Eps12 (Df x))) F"
+
+  \<comment> \<open>First prove the limit from component function limits\<close>
+  have "((\<lambda>y. Base (((f y - f (Lim F (\<lambda>x. x))) - Df (y - Lim F (\<lambda>x. x))) /\<^sub>R norm (y - Lim F (\<lambda>x. x)))) \<longlongrightarrow> Base 0) F"
+    using assm has_derivative_def[of "(\<lambda>x. Base (f x))" "(\<lambda>x. Base (Df x))" F] by simp
+  moreover have "((\<lambda>y. Eps1 (((f y - f (Lim F (\<lambda>x. x))) - Df (y - Lim F (\<lambda>x. x))) /\<^sub>R norm (y - Lim F (\<lambda>x. x)))) \<longlongrightarrow> Base 0) F"
+    using assm has_derivative_def[of "(\<lambda>x. Eps1 (f x))" "(\<lambda>x. Eps1 (Df x))" F] by simp
+  moreover have "((\<lambda>y. Eps2 (((f y - f (Lim F (\<lambda>x. x))) - Df (y - Lim F (\<lambda>x. x))) /\<^sub>R norm (y - Lim F (\<lambda>x. x)))) \<longlongrightarrow> Base 0) F"
+    using assm has_derivative_def[of "(\<lambda>x. Eps2 (f x))" "(\<lambda>x. Eps2 (Df x))" F] by simp
+  moreover have "((\<lambda>y. Eps12 (((f y - f (Lim F (\<lambda>x. x))) - Df (y - Lim F (\<lambda>x. x))) /\<^sub>R norm (y - Lim F (\<lambda>x. x)))) \<longlongrightarrow> Base 0) F"
+    using assm has_derivative_def[of "(\<lambda>x. Eps12 (f x))" "(\<lambda>x. Eps12 (Df x))" F] by simp
+  ultimately have "((\<lambda>y. ((f y - f (Lim F (\<lambda>x. x))) - Df (y - Lim F (\<lambda>x. x))) /\<^sub>R norm (y - Lim F (\<lambda>x. x))) \<longlongrightarrow> 0) F"
+    by (simp add: tendsto_hyperdual_iff)
+
+  \<comment> \<open>Next prove bounded linearity of the composed derivative by proving each of that class'
+      assumptions from bounded linearity of the component derivatives\<close>
+  moreover have "bounded_linear Df"
+  proof
+    have bl:
+      "bounded_linear (\<lambda>x. Base (Df x))"
+      "bounded_linear (\<lambda>x. Eps1 (Df x))"
+      "bounded_linear (\<lambda>x. Eps2 (Df x))"
+      "bounded_linear (\<lambda>x. Eps12 (Df x))"
+      using assm has_derivative_def by blast+
+    then have "linear (\<lambda>x. Base (Df x))"
+     and "linear (\<lambda>x. Eps1 (Df x))"
+     and "linear (\<lambda>x. Eps2 (Df x))"
+     and "linear (\<lambda>x. Eps12 (Df x))"
+      using bounded_linear.linear by blast+
+    then show "\<And>x y. Df (x + y) = Df x + Df y"
+          and "\<And>x y. Df (x *\<^sub>R y) = x *\<^sub>R Df y"
+      using plus_hyperdual.code scaleR_hyperdual.code by (simp_all add: linear_iff)
+
+    show "\<exists>K. \<forall>x. norm (Df x) \<le> norm x * K"
+    proof -
+      obtain k_re k_eps1 k_eps2 k_eps12
+       where "\<forall>x. (norm (Base (Df x)))^2 \<le> (norm x * k_re)^2"
+         and "\<forall>x. (norm (Eps1 (Df x)))^2 \<le> (norm x * k_eps1)^2"
+         and "\<forall>x. (norm (Eps2 (Df x)))^2 \<le> (norm x * k_eps2)^2"
+         and "\<forall>x. (norm (Eps12 (Df x)))^2 \<le> (norm x * k_eps12)^2"
+        using bl bounded_linear.bounded norm_ge_zero power_mono by metis
+      moreover have "\<forall>x. (norm (Df x))^2 = (norm (Base (Df x)))^2 + (norm (Eps1 (Df x)))^2 + (norm (Eps2 (Df x)))^2 + (norm (Eps12 (Df x)))^2"
+        using inner_hyperdual_def power2_norm_eq_inner by metis
+      ultimately have "\<forall>x. (norm (Df x))^2 \<le> (norm x * k_re)^2 + (norm x * k_eps1)^2 + (norm x * k_eps2)^2 + (norm x * k_eps12)^2"
+        by smt
+      then have "\<forall>x. (norm (Df x))^2 \<le> (norm x)^2 * (k_re^2 + k_eps1^2 + k_eps2^2 + k_eps12^2)"
+        by (simp add: distrib_left power_mult_distrib)
+      then have final: "\<forall>x. norm (Df x) \<le> norm x * sqrt(k_re^2 + k_eps1^2 + k_eps2^2 + k_eps12^2)"
+        using real_le_rsqrt real_sqrt_mult real_sqrt_pow2 by fastforce
+      then show "\<exists>K. \<forall>x. norm (Df x) \<le> norm x * K"
+        by blast
+    qed
+  qed
+  \<comment> \<open>Finally put the two together to finish the proof\<close>
+  ultimately show "(f has_derivative Df) F"
+    by (simp add: has_derivative_def)
+qed
+
+text\<open>Stop automatically unfolding hyperduals into components outside this theory:\<close>
+lemmas [iff del] = hyperdual_eq_iff
+
+end
diff --git a/thys/Hyperdual/HyperdualFunctionExtension.thy b/thys/Hyperdual/HyperdualFunctionExtension.thy
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/HyperdualFunctionExtension.thy
@@ -0,0 +1,647 @@
+(*  Title:   HyperdualFunctionExtension.thy
+    Authors: Jacques D. Fleuriot and Filip Smola, University of Edinburgh, 2021
+*)
+
+section \<open>Hyperdual Extension of Functions\<close>
+
+theory HyperdualFunctionExtension
+  imports Hyperdual TwiceFieldDifferentiable
+begin
+
+text\<open>The following is an important fact in the derivation of the hyperdual extension.\<close>
+lemma
+    fixes x :: "('a :: comm_ring_1) hyperdual" and n :: nat
+  assumes "Base x = 0"
+    shows "x ^ (n + 3) = 0"
+proof (induct n)
+  case 0
+  then show ?case
+    using assms hyperdual_power[of x 3] by simp
+next
+  case (Suc n)
+  then show ?case
+    using assms power_Suc[of x "n + 3"] mult_zero_right add_Suc by simp
+qed
+
+text\<open>We define the extension of a function to the hyperdual numbers.\<close>
+primcorec hypext :: "(('a :: real_normed_field) \<Rightarrow> 'a) \<Rightarrow> 'a hyperdual \<Rightarrow> 'a hyperdual" (\<open>*h* _\<close> [80] 80)
+  where
+    "Base ((*h* f) x) = f (Base x)"
+  | "Eps1 ((*h* f) x) = Eps1 x * deriv f (Base x)"
+  | "Eps2 ((*h* f) x) = Eps2 x * deriv f (Base x)"
+  | "Eps12 ((*h* f) x) = Eps12 x * deriv f (Base x) + Eps1 x * Eps2 x * deriv (deriv f) (Base x)"
+
+text\<open>This has the expected behaviour when expressed in terms of the units.\<close>
+lemma hypext_Hyperdual_eq:
+  "(*h* f) (Hyperdual a b c d) =
+     Hyperdual (f a) (b * deriv f a) (c * deriv f a) (d * deriv f a + b * c * deriv (deriv f) a)"
+  by (simp add: hypext.code)
+
+lemma hypext_Hyperdual_eq_parts:
+  "(*h* f) (Hyperdual a b c d) =
+      f a *\<^sub>H ba + (b * deriv f a) *\<^sub>H e1 + (c * deriv f a) *\<^sub>H e2 +
+         (d * deriv f a + b * c * deriv (deriv f) a) *\<^sub>H e12 "
+  by (metis Hyperdual_eq hypext_Hyperdual_eq)
+
+text\<open>
+  The extension can be used to extract the function value, and first and second derivatives at x
+  when applied to @{term "x *\<^sub>H re + e1 + e2 + 0 *\<^sub>H e12"}, which we denote by @{term "\<beta> x"}.
+\<close>
+definition hyperdualx :: "('a :: real_normed_field) \<Rightarrow> 'a hyperdual" ("\<beta>")
+  where "\<beta> x = (Hyperdual x 1 1 0)"
+
+lemma hyperdualx_sel [simp]:
+  shows "Base (\<beta> x) = x"
+    and "Eps1 (\<beta> x) = 1"
+    and "Eps2 (\<beta> x) = 1"
+    and "Eps12 (\<beta> x) = 0"
+  by (simp_all add: hyperdualx_def)
+
+lemma hypext_extract_eq:
+  "(*h* f) (\<beta> x) = f x *\<^sub>H ba + deriv f x *\<^sub>H e1 + deriv f x *\<^sub>H e2 + deriv (deriv f) x *\<^sub>H e12"
+  by (simp add: hypext_Hyperdual_eq_parts hyperdualx_def)
+
+lemma Base_hypext:
+  "Base ((*h* f) (\<beta> x)) = f x"
+  by (simp add: hyperdualx_def)
+
+lemma Eps1_hypext:
+  "Eps1 ((*h* f) (\<beta> x)) = deriv f x"
+  by (simp add: hyperdualx_def)
+
+lemma Eps2_hypext:
+  "Eps2 ((*h* f) (\<beta> x)) = deriv f x"
+  by (simp add: hyperdualx_def)
+
+lemma Eps12_hypext:
+  "Eps12 ((*h* f) (\<beta> x)) = deriv (deriv f) x"
+  by (simp add: hyperdualx_def)
+
+subsubsection\<open>Convenience Interface\<close>
+
+text\<open>Define a datatype to hold the function value, and the first and second derivative values.\<close>
+datatype ('a :: real_normed_field) derivs = Derivs (Value: 'a) (First: 'a) (Second: 'a)
+
+text\<open>
+  Then we convert a hyperdual number to derivative values by extracting the base component, one of
+  the first-order components, and the second-order component.
+\<close>
+fun hyperdual_to_derivs :: "('a :: real_normed_field) hyperdual \<Rightarrow> 'a derivs"
+  where "hyperdual_to_derivs x = Derivs (Base x) (Eps1 x) (Eps12 x)"
+
+text\<open>
+  Finally we define way of converting any compatible function into one that yields the value and the
+  derivatives.
+\<close>
+fun autodiff :: "('a :: real_normed_field \<Rightarrow> 'a) \<Rightarrow> 'a \<Rightarrow> 'a derivs"
+  where "autodiff f = (\<lambda>x. hyperdual_to_derivs ((*h* f) (\<beta> x)))"
+
+lemma autodiff_sel:
+  "Value (autodiff f x) = Base ((*h* f) (\<beta> x))"
+  "First (autodiff f x) = Eps1 ((*h* f) (\<beta> x))"
+  "Second (autodiff f x) = Eps12 ((*h* f) (\<beta> x))"
+  by simp_all
+
+text\<open>The result contains the expected values.\<close>
+lemma autodiff_extract_value:
+  "Value (autodiff f x) = f x"
+  by (simp del: hypext.simps add: Base_hypext)
+
+lemma autodiff_extract_first:
+  "First (autodiff f x) = deriv f x"
+  by (simp del: hypext.simps add: Eps1_hypext)
+
+lemma autodiff_extract_second:
+  "Second (autodiff f x) = deriv (deriv f) x"
+  by (simp del: hypext.simps add: Eps12_hypext)
+
+text\<open>
+  The derivative components of the result are actual derivatives if the function is sufficiently
+  differentiable on that argument.
+\<close>
+lemma autodiff_first_derivative:
+  assumes "f field_differentiable (at x)"
+  shows "(f has_field_derivative First (autodiff f x)) (at x)"
+  by (simp add: autodiff_extract_first DERIV_deriv_iff_field_differentiable assms)
+
+lemma autodiff_second_derivative:
+  assumes "f twice_field_differentiable_at x"
+  shows "((deriv f) has_field_derivative Second (autodiff f x)) (at x)"
+  by (simp add: autodiff_extract_second DERIV_deriv_iff_field_differentiable assms deriv_field_differentiable_at)
+
+subsubsection\<open>Composition\<close>
+
+text\<open>Composition of hyperdual extensions is the hyperdual extension of composition:\<close>
+lemma hypext_compose:
+  assumes "f twice_field_differentiable_at (Base x)"
+      and "g twice_field_differentiable_at (f (Base x))"
+    shows "(*h* (\<lambda>x. g (f x))) x = (*h* g) ((*h* f) x)"
+proof (simp add: hyperdual_eq_iff, intro conjI disjI2)
+  show goal1: "deriv (\<lambda>x. g (f x)) (Base x) = deriv f (Base x) * deriv g (f (Base x))"
+  proof -
+    have "deriv (\<lambda>x. g (f x)) (Base x) = deriv (g \<circ> f) (Base x)"
+      by (simp add: comp_def)
+    also have "... = deriv g (f (Base x)) * deriv f (Base x)"
+      using assms by (simp add: deriv_chain once_field_differentiable_at)
+    finally show ?thesis
+      by (simp add: mult.commute deriv_chain)
+  qed
+  then show "deriv (\<lambda>x. g (f x)) (Base x) = deriv f (Base x) * deriv g (f (Base x))" .
+
+  have first_diff: "(\<lambda>x. deriv g (f x)) field_differentiable at (Base x)"
+    by (metis DERIV_chain2 assms deriv_field_differentiable_at field_differentiable_def once_field_differentiable_at)
+
+  have "deriv (deriv g \<circ> f) (Base x) = deriv (deriv g) (f (Base x)) * deriv f (Base x)"
+    using deriv_chain assms once_field_differentiable_at deriv_field_differentiable_at
+    by blast
+  then have deriv_deriv_comp: "deriv (\<lambda>x. deriv g (f x)) (Base x) = deriv (deriv g) (f (Base x)) * deriv f (Base x)"
+    by (simp add: comp_def)
+
+  have "deriv (deriv (\<lambda>x. g (f x))) (Base x) = deriv ((\<lambda>x. deriv f x * deriv g (f x))) (Base x)"
+    using assms eventually_deriv_compose'[of f "Base x" g]
+    by (simp add: mult.commute deriv_cong_ev)
+  also have "... = deriv f (Base x) * deriv (\<lambda>x. deriv g (f x)) (Base x) + deriv (deriv f) (Base x) * deriv g (f (Base x))"
+    using assms(1) first_diff by (simp add: deriv_field_differentiable_at)
+  also have "... = deriv f (Base x) * deriv (deriv g) (f (Base x)) * deriv f (Base x) + deriv (deriv f) (Base x) * deriv g (f (Base x))"
+    using deriv_deriv_comp by simp
+  finally show "Eps12 x * deriv (\<lambda>x. g (f x)) (Base x) + Eps1 x * Eps2 x * deriv (deriv (\<lambda>x. g (f x))) (Base x) =
+    (Eps12 x * deriv f (Base x) + Eps1 x * Eps2 x * deriv (deriv f) (Base x)) * deriv g (f (Base x)) +
+    Eps1 x * deriv f (Base x) * (Eps2 x * deriv f (Base x)) * deriv (deriv g) (f (Base x))"
+    by (simp add: goal1 field_simps)
+qed
+
+subsection\<open>Concrete Instances\<close>
+
+subsubsection\<open>Constant\<close>
+
+text\<open>Component embedding is an extension of the constant function.\<close>
+lemma hypext_const [simp]:
+  "(*h* (\<lambda>x. a)) x = of_comp a"
+  by (simp add: of_comp_def hyperdual_eq_iff)
+
+lemma "autodiff (\<lambda>x. a) = (\<lambda>x. Derivs a 0 0)"
+  by simp
+
+subsubsection\<open>Identity\<close>
+
+text\<open>Identity is an extension of the component identity.\<close>
+lemma hypext_ident:
+  "(*h* (\<lambda>x. x)) x = x"
+  by (simp add: hyperdual_eq_iff)
+
+subsubsection\<open>Component Scalar Multiplication\<close>
+
+text\<open>Component scaling is an extension of component constant multiplication:\<close>
+lemma hypext_scaleH:
+  "(*h* (\<lambda>x. k * x)) x = k *\<^sub>H x"
+  by (simp add: hyperdual_eq_iff)
+
+lemma hypext_fun_scaleH:
+  assumes "f twice_field_differentiable_at (Base x)"
+  shows "(*h* (\<lambda>x. k * f x)) x = k *\<^sub>H (*h* f) x"
+  using assms by (simp add: hypext_compose hypext_scaleH)
+
+text\<open>Unary minus is just an instance of constant multiplication:\<close>
+lemma hypext_uminus:
+  "(*h* uminus) x = - x"
+  using hypext_scaleH[of "-1" x] by simp
+
+subsubsection\<open>Real Scalar Multiplication\<close>
+
+text\<open>Real scaling is an extension of component real scaling:\<close>
+lemma hypext_scaleR:
+  "(*h* (\<lambda>x. k *\<^sub>R x)) x = k *\<^sub>R x"
+  by (auto simp add: hyperdual_eq_iff)
+
+lemma hypext_fun_scaleR:
+  assumes "f twice_field_differentiable_at (Base x)"
+  shows "(*h* (\<lambda>x. k *\<^sub>R f x)) x = k *\<^sub>R (*h* f) x"
+  using assms by (simp add: hypext_compose hypext_scaleR)
+
+subsubsection\<open>Addition\<close>
+
+text\<open>Addition of hyperdual extensions is a hyperdual extension of addition of functions.\<close>
+lemma hypext_fun_add:
+  assumes "f twice_field_differentiable_at (Base x)"
+      and "g twice_field_differentiable_at (Base x)"
+    shows "(*h* (\<lambda>x. f x + g x)) x = (*h* f) x + (*h* g) x"
+proof (simp add: hyperdual_eq_iff distrib_left[symmetric], intro conjI disjI2)
+  show goal1: "deriv (\<lambda>x. f x + g x) (Base x) = deriv f (Base x) + deriv g (Base x)"
+    by (simp add: assms once_field_differentiable_at distrib_left)
+  then show "deriv (\<lambda>x. f x + g x) (Base x) = deriv f (Base x) + deriv g (Base x)" .
+
+  have "deriv (deriv (\<lambda>x. f x + g x)) (Base x) = deriv (\<lambda>w. deriv f w + deriv g w) (Base x)"
+    by (simp add: assms deriv_cong_ev eventually_deriv_add)
+  moreover have "Eps12 x * deriv f (Base x) + Eps12 x * deriv g (Base x) = Eps12 x * deriv (\<lambda>x. f x + g x) (Base x)"
+    by (metis distrib_left goal1)
+  ultimately show "Eps12 x * deriv (\<lambda>x. f x + g x) (Base x) +
+    Eps1 x * Eps2 x * deriv (deriv (\<lambda>x. f x + g x)) (Base x) =
+    Eps12 x * deriv f (Base x) + Eps1 x * Eps2 x * deriv (deriv f) (Base x) +
+    (Eps12 x * deriv g (Base x) + Eps1 x * Eps2 x * deriv (deriv g) (Base x))"
+    using deriv_add[OF deriv_field_differentiable_at deriv_field_differentiable_at, OF assms]
+    by (simp add: distrib_left add.left_commute)
+qed
+
+lemma hypext_cadd [simp]:
+  "(*h* (\<lambda>x. x + a)) x = x + of_comp a"
+  by (auto simp add: hyperdual_eq_iff of_comp_def)
+
+lemma hypext_fun_cadd:
+  assumes "f twice_field_differentiable_at (Base x)"
+  shows "(*h* (\<lambda>x. f x + a)) x = (*h* f) x + of_comp a"
+  using assms hypext_compose[of f x "\<lambda>x. x + a"] by simp
+
+subsubsection\<open>Component Linear Function\<close>
+
+text\<open>Hyperdual linear function is an extension of the component linear function:\<close>
+lemma hypext_linear:
+  "(*h* (\<lambda>x. k * x + a)) x = k *\<^sub>H x + of_comp a"
+  using hypext_fun_add[of "(*) k" x "\<lambda>x. a"]
+  by (simp add: hypext_scaleH)
+
+lemma hypext_fun_linear:
+  assumes "f twice_field_differentiable_at (Base x)"
+    shows "(*h* (\<lambda>x. k * f x + a)) x = k *\<^sub>H (*h* f) x + of_comp a"
+  using assms hypext_compose[of f x "\<lambda>x. k * x + a"] by (simp add: hypext_linear)
+
+subsubsection\<open>Real Linear Function\<close>
+
+text\<open>We have the same for real scaling instead of component multiplication:\<close>
+lemma hypext_linearR:
+  "(*h* (\<lambda>x. k *\<^sub>R x + a)) x = k *\<^sub>R x + of_comp a"
+  using hypext_fun_add[of "(*\<^sub>R) k" x "\<lambda>x. a"]
+  by (simp add: hypext_scaleR)
+
+lemma hypext_fun_linearR:
+  assumes "f twice_field_differentiable_at (Base x)"
+    shows "(*h* (\<lambda>x. k *\<^sub>R f x + a)) x = k *\<^sub>R (*h* f) x + of_comp a"
+  using assms hypext_compose[of f x "\<lambda>x. k *\<^sub>R x + a"] by (simp add: hypext_linearR)
+
+subsubsection\<open>Multiplication\<close>
+
+text\<open>Extension of multiplication is multiplication of the functions' extensions.\<close>
+lemma hypext_fun_mult:
+  assumes "f twice_field_differentiable_at (Base x)"
+      and "g twice_field_differentiable_at (Base x)"
+    shows "(*h* (\<lambda>z. f z * g z)) x = (*h* f) x * (*h* g) x"
+proof (simp add: hyperdual_eq_iff distrib_left[symmetric], intro conjI)
+  show "Eps1 x * deriv (\<lambda>z. f z * g z) (Base x) =
+        f (Base x) * (Eps1 x * deriv g (Base x)) + Eps1 x * deriv f (Base x) * g (Base x)"
+   and "Eps2 x * deriv (\<lambda>z. f z * g z) (Base x) =
+        f (Base x) * (Eps2 x * deriv g (Base x)) + Eps2 x * deriv f (Base x) * g (Base x)"
+    using assms by (simp_all add: once_field_differentiable_at distrib_left)
+
+  have
+    "deriv (deriv (\<lambda>z. f z * g z)) (Base x) =
+     f (Base x) * deriv (deriv g) (Base x) + 2 * deriv f (Base x) * deriv g (Base x) + deriv (deriv f) (Base x) * g (Base x)"
+  proof -
+    have "deriv (deriv (\<lambda>z. f z * g z)) (Base x) = deriv (\<lambda>z. f z * deriv g z + deriv f z * g z) (Base x)"
+      using assms by (simp add: eventually_deriv_mult deriv_cong_ev)
+    also have "... = (\<lambda>z. f z * deriv (deriv g) z + deriv f z * deriv g z + deriv f z * deriv g z + deriv (deriv f) z * g z) (Base x)"
+      by (simp add: assms deriv_field_differentiable_at field_differentiable_mult once_field_differentiable_at)
+    finally show ?thesis
+      by simp
+  qed
+  then show
+    "Eps12 x * deriv (\<lambda>z. f z * g z) (Base x) + Eps1 x * Eps2 x * deriv (deriv (\<lambda>z. f z * g z)) (Base x) =
+     2 * (Eps1 x * (Eps2 x * (deriv f (Base x) * deriv g (Base x)))) +
+     f (Base x) * (Eps12 x * deriv g (Base x) + Eps1 x * Eps2 x * deriv (deriv g) (Base x)) +
+     (Eps12 x * deriv f (Base x) + Eps1 x * Eps2 x * deriv (deriv f) (Base x)) * g (Base x)"
+     using assms by (simp add: once_field_differentiable_at field_simps)
+qed
+
+subsubsection\<open>Sine and Cosine\<close>
+
+text\<open>The extended sin and cos at an arbitrary hyperdual.\<close>
+
+lemma hypext_sin_Hyperdual:
+  "(*h* sin) (Hyperdual a b c d) = sin a *\<^sub>H ba + (b *cos a) *\<^sub>H e1 + (c * cos a) *\<^sub>H e2 + (d * cos a - b * c * sin a) *\<^sub>H e12 "
+  by (simp add: hypext_Hyperdual_eq_parts)
+
+lemma hypext_cos_Hyperdual:
+  "(*h* cos) (Hyperdual a b c d) = cos a *\<^sub>H ba - (b * sin a) *\<^sub>H e1 - (c * sin a) *\<^sub>H e2 - (d * sin a + b * c * cos a) *\<^sub>H e12 "
+proof -
+  have "of_comp (- (d * sin a) - b * c * cos a) * e12 = - (of_comp (d * sin a + b * c * cos a) * e12)"
+    by (metis add_uminus_conv_diff minus_add_distrib mult_minus_left of_comp_minus)
+  then show ?thesis
+    by (simp add: hypext_Hyperdual_eq_parts of_comp_minus scaleH_times)
+qed
+
+lemma Eps1_hypext_sin [simp]:
+  "Eps1 ((*h* sin) x) = Eps1 x * cos (Base x)"
+  by simp
+
+lemma Eps2_hypext_sin [simp]:
+  "Eps2 ((*h* sin) x) = Eps2 x * cos (Base x)"
+  by simp
+
+lemma Eps12_hypext_sin [simp]:
+  "Eps12 ((*h* sin) x) = Eps12 x * cos (Base x) - Eps1 x * Eps2 x * sin (Base x)"
+  by simp
+
+lemma hypext_sin_e1 [simp]:
+  "(*h* sin) (x * e1) = e1 * x"
+  by (simp add: e1_def hyperdual_eq_iff one_hyperdual_def)
+
+lemma hypext_sin_e2 [simp]:
+  "(*h* sin) (x * e2) = e2 * x"
+  by (simp add: e2_def hyperdual_eq_iff one_hyperdual_def)
+
+lemma hypext_sin_e12 [simp]:
+  "(*h* sin) (x * e12) = e12 * x"
+  by (simp add: e12_def hyperdual_eq_iff one_hyperdual_def)
+
+lemma hypext_cos_e1 [simp]:
+  "(*h* cos) (x * e1) = 1"
+  by (simp add: e1_def hyperdual_eq_iff one_hyperdual_def)
+
+lemma hypext_cos_e2 [simp]:
+  "(*h* cos) (x * e2) = 1"
+  by (simp add: e2_def hyperdual_eq_iff one_hyperdual_def)
+
+lemma hypext_cos_e12 [simp]:
+  "(*h* cos) (x * e12) = 1"
+  by (simp add: e12_def hyperdual_eq_iff one_hyperdual_def)
+
+text\<open>The extended sin and cos at @{term "\<beta> x"}.\<close>
+
+lemma hypext_sin_extract:
+  "(*h* sin) (\<beta> x) = sin x *\<^sub>H ba + cos x *\<^sub>H e1 + cos x *\<^sub>H e2 - sin x *\<^sub>H e12"
+  by (simp add: hypext_sin_Hyperdual of_comp_minus scaleH_times hyperdualx_def)
+
+lemma hypext_cos_extract:
+  "(*h* cos) (\<beta> x) = cos x *\<^sub>H ba - sin x *\<^sub>H e1 - sin x *\<^sub>H e2 - cos x *\<^sub>H e12"
+  by (simp add: hypext_cos_Hyperdual hyperdualx_def)
+
+text\<open>Extracting the extended sin components at @{term "\<beta> x"}.\<close>
+
+lemma Base_hypext_sin_extract [simp]:
+  "Base ((*h* sin) (\<beta> x)) = sin x"
+  by (rule Base_hypext)
+
+lemma Eps2_hypext_sin_extract [simp]:
+  "Eps2 ((*h* sin) (\<beta> x)) = cos x"
+  using Eps2_hypext[of sin] by simp
+
+lemma Eps12_hypext_sin_extract [simp]:
+  "Eps12 ((*h* sin) (\<beta> x)) = - sin x"
+  using Eps12_hypext[of sin] by simp
+
+text\<open>Extracting the extended cos components at @{term "\<beta> x"}.\<close>
+
+lemma Base_hypext_cos_extract [simp]:
+  "Base ((*h* cos) (\<beta> x)) = cos x"
+  by (rule Base_hypext)
+
+lemma Eps2_hypext_cos_extract [simp]:
+  "Eps2 ((*h* cos) (\<beta> x)) = - sin x"
+  using Eps2_hypext[of cos] by simp
+
+lemma Eps12_hypext_cos_extract [simp]:
+  "Eps12 ((*h* cos) (\<beta> x)) = - cos x"
+  using Eps12_hypext[of cos] by simp
+
+text\<open>We get one of the key trigonometric properties for the extensions of sin and cos.\<close>
+
+lemma "((*h* sin) x)\<^sup>2 + ((*h* cos) x)\<^sup>2 = 1"
+  by (simp add: hyperdual_eq_iff one_hyperdual_def power2_eq_square field_simps)
+
+(* example *)
+lemma "(*h* sin) x + (*h* cos) x = (*h* (\<lambda>x. sin x + cos x)) x"
+  by (simp add: hypext_fun_add)
+
+subsubsection\<open>Exponential\<close>
+
+text\<open>The exponential function extension behaves as expected.\<close>
+
+lemma hypext_exp_Hyperdual:
+  "(*h* exp) (Hyperdual a b c d) =
+       exp a *\<^sub>H ba + (b * exp a) *\<^sub>H e1 + (c * exp a) *\<^sub>H e2 + (d * exp a + b * c * exp a) *\<^sub>H e12"
+  by (simp add: hypext_Hyperdual_eq_parts)
+
+lemma hypext_exp_extract:
+  "(*h* exp) (\<beta> x) = exp x *\<^sub>H ba + exp x *\<^sub>H e1 + exp x *\<^sub>H e2 + exp x *\<^sub>H e12"
+  by (simp add: hypext_extract_eq)
+
+lemma hypext_exp_e1 [simp]:
+  "(*h* exp) (x * e1) = 1 + e1 * x"
+  by (simp add: e1_def hyperdual_eq_iff)
+
+lemma hypext_exp_e2 [simp]:
+  "(*h* exp) (x * e2) = 1 + e2 * x"
+  by (simp add: e2_def hyperdual_eq_iff)
+
+lemma hypext_exp_e12 [simp]:
+  "(*h* exp) (x * e12) = 1 + e12 * x"
+  by (simp add: e12_def hyperdual_eq_iff)
+
+text\<open>Extracting the parts for the exponential function extension.\<close>
+
+lemma Eps1_hypext_exp_extract [simp]:
+  "Eps1 ((*h* exp) (\<beta> x)) = exp x"
+  using Eps1_hypext[of exp] by simp
+
+lemma Eps2_hypext_exp_extract [simp]:
+  "Eps2 ((*h* exp) (\<beta> x)) = exp x"
+  using Eps2_hypext[of exp] by simp
+
+lemma Eps12_hypext_exp_extract [simp]:
+  "Eps12 ((*h* exp) (\<beta> x)) = exp x"
+  using Eps12_hypext[of exp] by simp
+
+subsubsection\<open>Square Root\<close>
+text\<open>Square root function extension.\<close>
+
+lemma hypext_sqrt_Hyperdual_Hyperdual:
+  assumes "a > 0"
+  shows "(*h* sqrt) (Hyperdual a b c d) =
+         Hyperdual (sqrt a) (b * inverse (sqrt a) / 2) (c * inverse (sqrt a) / 2)
+           (d * inverse (sqrt a) / 2 - b * c * inverse (sqrt a ^ 3) / 4)"
+  by (simp add: assms hypext_Hyperdual_eq)
+
+lemma hypext_sqrt_Hyperdual:
+      "a > 0 \<Longrightarrow> (*h* sqrt) (Hyperdual a b c d) =
+       sqrt a *\<^sub>H ba + (b * inverse (sqrt a) / 2) *\<^sub>H e1 + (c * inverse (sqrt a) / 2) *\<^sub>H e2 +
+       (d * inverse (sqrt a) / 2 - b * c * inverse (sqrt a ^ 3) / 4) *\<^sub>H e12"
+  by (auto simp add: hypext_Hyperdual_eq_parts)
+
+lemma hypext_sqrt_extract:
+  "x > 0 \<Longrightarrow> (*h* sqrt) (\<beta> x) = sqrt x *\<^sub>H ba + (inverse (sqrt x) / 2) *\<^sub>H e1 +
+        (inverse (sqrt x) / 2) *\<^sub>H e2 - (inverse (sqrt x ^ 3) / 4) *\<^sub>H e12"
+  by (simp add: hypext_sqrt_Hyperdual hyperdualx_def of_comp_minus scaleH_times)
+
+text\<open>Extracting the parts for the square root extension.\<close>
+
+lemma Eps1_hypext_sqrt_extract [simp]:
+  "x > 0 \<Longrightarrow> Eps1 ((*h* sqrt) (\<beta> x)) = inverse (sqrt x) / 2"
+  using Eps1_hypext[of sqrt] by simp
+
+lemma Eps2_hypext_sqrt_extract [simp]:
+  "x > 0 \<Longrightarrow> Eps2 ((*h* sqrt) (\<beta> x)) = inverse (sqrt x) / 2"
+  using Eps2_hypext[of sqrt] by simp
+
+lemma Eps12_hypext_sqrt_extract [simp]:
+  "x > 0 \<Longrightarrow> Eps12 ((*h* sqrt) (\<beta> x)) = - (inverse (sqrt x ^ 3) / 4)"
+  using Eps12_hypext[of sqrt] by simp
+
+(* example *)
+lemma "Base x > 0 \<Longrightarrow> (*h* sin) x + (*h* sqrt) x = (*h* (\<lambda>x. sin x + sqrt x)) x"
+  by (simp add: hypext_fun_add)
+
+subsubsection\<open>Natural Power\<close>
+
+lemma hypext_power:
+  "(*h* (\<lambda>x. x ^ n)) x = x ^ n"
+  by (simp add: hyperdual_eq_iff hyperdual_power)
+
+lemma hypext_fun_power:
+  assumes "f twice_field_differentiable_at (Base x)"
+    shows "(*h* (\<lambda>x. (f x) ^ n)) x = ((*h* f) x) ^ n"
+  using assms hypext_compose[of f x "\<lambda>x. x ^ n"] by (simp add: hypext_power)
+
+lemma hypext_power_Hyperdual:
+  "(*h* (\<lambda>x. x ^ n)) (Hyperdual a b c d) =
+        a ^ n *\<^sub>H ba + (of_nat n * b * a ^ (n - 1)) *\<^sub>H e1 + (of_nat n * c * a ^ (n - 1)) *\<^sub>H e2 +
+        (d * (of_nat n * a ^ (n - 1)) + b * c * (of_nat n * of_nat (n - 1) * a ^ (n - 2))) *\<^sub>H e12"
+  by (simp add: hypext_Hyperdual_eq_parts algebra_simps)
+
+lemma hypext_power_Hyperdual_parts:
+  "(*h* (\<lambda>x. x ^ n)) (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) =
+        a ^ n *\<^sub>H ba + (of_nat n * b * a ^ (n - 1)) *\<^sub>H e1 + (of_nat n * c * a ^ (n - 1)) *\<^sub>H e2 +
+        (d * (of_nat n * a ^ (n - 1)) + b * c * (of_nat n * of_nat (n - 1) * a ^ (n - 2))) *\<^sub>H e12"
+  by (simp add: Hyperdual_eq [symmetric] hypext_power_Hyperdual)
+
+lemma hypext_power_extract:
+  "(*h* (\<lambda>x. x ^ n)) (\<beta> x) =
+      x ^ n *\<^sub>H ba + (of_nat n * x ^ (n - 1)) *\<^sub>H e1 + (of_nat n * x ^ (n - 1)) *\<^sub>H e2 +
+      (of_nat n * of_nat (n - 1) * x ^ (n - 2)) *\<^sub>H e12"
+  by (simp add: hypext_extract_eq)
+
+lemma Eps1_hypext_power [simp]:
+  "Eps1 ((*h* (\<lambda>x. x ^ n)) x) = of_nat n * Eps1 x * (Base x) ^ (n - 1)"
+  by simp
+
+lemma Eps2_hypext_power [simp]:
+  "Eps2 ((*h* (\<lambda>x. x ^ n)) x) = of_nat n * Eps2 x * (Base x) ^ (n - 1)"
+  by simp
+
+lemma Eps12_hypext_power [simp]:
+  "Eps12 ((*h* (\<lambda>x. x ^ n)) x) =
+   Eps12 x * (of_nat n * Base x ^ (n - 1)) + Eps1 x * Eps2 x * (of_nat n * of_nat (n - 1) * Base x ^ (n - 2))"
+  by simp
+
+subsubsection\<open>Inverse\<close>
+
+lemma hypext_inverse:
+  assumes "Base x \<noteq> 0"
+  shows "(*h* inverse) x = inverse x"
+  using assms by (simp add: hyperdual_eq_iff inverse_eq_divide)
+
+lemma hypext_fun_inverse:
+  assumes "f twice_field_differentiable_at (Base x)"
+      and "f (Base x) \<noteq> 0"
+    shows "(*h* (\<lambda>x. inverse (f x))) x = inverse ((*h* f) x)"
+  using assms hypext_compose[of f x inverse] by (simp add: hypext_inverse)
+
+lemma hypext_inverse_Hyperdual:
+  "a \<noteq> 0 \<Longrightarrow>
+    (*h* inverse) (Hyperdual a b c d) =
+    Hyperdual (inverse a) (- (b / a\<^sup>2)) (- (c / a\<^sup>2)) (2 * b * c / (a ^ 3) - d / a\<^sup>2)"
+  by (simp add: hypext_Hyperdual_eq divide_inverse)
+
+lemma hypext_inverse_Hyperdual_parts:
+  "a \<noteq> 0 \<Longrightarrow>
+    (*h* inverse) (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) =
+    inverse a *\<^sub>H ba + - (b / a\<^sup>2) *\<^sub>H e1 + - (c / a\<^sup>2) *\<^sub>H e2 + (2 * b * c / a ^ 3 - d / a\<^sup>2) *\<^sub>H e12"
+  by (metis Hyperdual_eq hypext_inverse_Hyperdual)
+
+lemma inverse_Hyperdual_parts:
+  "(a::'a::real_normed_field) \<noteq> 0 \<Longrightarrow>
+    inverse (a *\<^sub>H ba + b *\<^sub>H e1 + c *\<^sub>H e2 + d *\<^sub>H e12) =
+    inverse a *\<^sub>H ba + - (b / a\<^sup>2) *\<^sub>H e1 + - (c / a\<^sup>2) *\<^sub>H e2 + (2 * b * c / a ^ 3 - d / a\<^sup>2) *\<^sub>H e12"
+  by (metis Hyperdual_eq hyperdual.sel(1) hypext_inverse hypext_inverse_Hyperdual_parts)
+
+lemma hypext_inverse_extract:
+  "x \<noteq> 0 \<Longrightarrow> (*h* inverse) (\<beta> x) = inverse x *\<^sub>H ba - (1 / x\<^sup>2) *\<^sub>H e1 - (1 / x\<^sup>2) *\<^sub>H e2 + (2 / x ^ 3) *\<^sub>H e12"
+  by (simp add: hypext_extract_eq divide_inverse of_comp_minus scaleH_times)
+
+lemma inverse_extract:
+  "x \<noteq> 0 \<Longrightarrow> inverse (\<beta> x) = inverse x *\<^sub>H ba - (1 / x\<^sup>2) *\<^sub>H e1 - (1 / x\<^sup>2) *\<^sub>H e2 + (2 / x ^ 3) *\<^sub>H e12"
+  by (metis hyperdual.sel(1) hyperdualx_def hypext_inverse hypext_inverse_extract)
+
+lemma Eps1_hypext_inverse [simp]:
+  "Base x \<noteq> 0 \<Longrightarrow> Eps1 ((*h* inverse) x) = - Eps1 x * (1 / (Base x)\<^sup>2)"
+  by simp
+lemma Eps1_inverse [simp]:
+  "Base (x::'a::real_normed_field hyperdual) \<noteq> 0 \<Longrightarrow> Eps1 (inverse x) = - Eps1 x * (1 / (Base x)\<^sup>2)"
+  by simp
+
+lemma Eps2_hypext_inverse [simp]:
+  "Base (x::'a::real_normed_field hyperdual) \<noteq> 0 \<Longrightarrow> Eps2 (inverse x) = - Eps2 x * (1 / (Base x)\<^sup>2)"
+  by simp
+
+lemma Eps12_hypext_inverse [simp]:
+  "Base  (x::'a::real_normed_field hyperdual) \<noteq> 0
+   \<Longrightarrow> Eps12 (inverse x) = Eps1 x * Eps2 x * (2/ (Base x ^ 3)) - Eps12 x / (Base x)\<^sup>2"
+  by simp
+
+subsubsection\<open>Division\<close>
+
+lemma hypext_fun_divide:
+  assumes "f twice_field_differentiable_at (Base x)"
+      and "g twice_field_differentiable_at (Base x)"
+      and "g (Base x) \<noteq> 0"
+    shows "(*h* (\<lambda>x. f x / g x)) x = (*h* f) x / (*h* g) x"
+proof -
+  have "(\<lambda>x. inverse (g x)) twice_field_differentiable_at Base x"
+    by (simp add: assms(2) assms(3) twice_field_differentiable_at_compose)
+  moreover have "(*h* f) x * (*h* (\<lambda>x. inverse (g x))) x = (*h* f) x * inverse ((*h* g) x)"
+    by (simp add: assms(2) assms(3) hypext_fun_inverse)
+  ultimately have "(*h* (\<lambda>x. f x * inverse (g x))) x = (*h* f) x * inverse ((*h* g) x)"
+    by (simp add: assms(1) hypext_fun_mult)
+  then show ?thesis
+    by (simp add: divide_inverse hyp_divide_inverse)
+qed
+
+subsubsection\<open>Polynomial\<close>
+
+lemma hypext_polyn:
+  fixes coef :: "nat \<Rightarrow> 'a :: {real_normed_field}"
+    and n :: nat
+  shows "(*h* (\<lambda>x. \<Sum>i<n. coef i * x^i)) x = (\<Sum>i<n. (coef i) *\<^sub>H (x^i))"
+proof (induction n)
+  case 0
+  then show ?case
+    by (simp add: zero_hyperdual_def)
+next
+  case hyp: (Suc n)
+
+  have "(\<lambda>x. \<Sum>i<Suc n. coef i * x ^ i) = (\<lambda>x. (\<Sum>i<n. coef i * x ^ i) + coef n * x ^ n)"
+   and "(\<lambda>x. \<Sum>i<Suc n. coef i *\<^sub>H x ^ i) = (\<lambda>x. (\<Sum>i<n. coef i *\<^sub>H x ^ i) + coef n *\<^sub>H x ^ n)"
+    by (simp_all add: field_simps)
+
+  then show ?case
+  proof (simp)
+    have "(\<lambda>x. coef n * x ^ n) twice_field_differentiable_at Base x"
+      using twice_field_differentiable_at_compose[of "\<lambda>x. x ^ n" "Base x" "(*) (coef n)"]
+      by simp
+    then have "(*h* (\<lambda>x. (\<Sum>i<n. coef i * x ^ i) + coef n * x ^ n)) x =
+          (*h* (\<lambda>x. (\<Sum>i<n. coef i * x ^ i))) x + (*h* (\<lambda>x. coef n * x ^ n)) x"
+      by (simp add: hypext_fun_add)
+    moreover have "(*h* (\<lambda>x. coef n * x ^ n)) x = coef n *\<^sub>H x ^ n"
+      by (simp add: hypext_fun_scaleH hypext_power)
+    ultimately have "(*h* (\<lambda>x. (\<Sum>i<n. coef i * x ^ i) + coef n * x ^ n)) x = (\<Sum>i<n. coef i *\<^sub>H x ^ i) + coef n *\<^sub>H x ^ n"
+      using hyp by simp
+    then show "(*h* (\<lambda>x. (\<Sum>i<n. coef i * x ^ i) + coef n * x ^ n)) x = (\<Sum>i<n. coef i *\<^sub>H x ^ i) + coef n *\<^sub>H x ^ n"
+      by simp
+  qed
+qed
+
+lemma hypext_fun_polyn:
+    fixes coef :: "nat \<Rightarrow> 'a :: {real_normed_field}"
+      and n :: nat
+  assumes "f twice_field_differentiable_at (Base x)"
+    shows "(*h* (\<lambda>x. \<Sum>i<n. coef i * (f x)^i)) x = (\<Sum>i<n. (coef i) *\<^sub>H (((*h* f) x)^i))"
+  using assms hypext_compose[of f x "\<lambda>x. (\<Sum>i<n. coef i * x^i)"] by (simp add: hypext_polyn)
+
+end
diff --git a/thys/Hyperdual/LogisticFunction.thy b/thys/Hyperdual/LogisticFunction.thy
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/LogisticFunction.thy
@@ -0,0 +1,94 @@
+(*  Title:       LogisticFunction.thy
+    Author:      Filip Smola, 2019-2021
+*)
+
+theory LogisticFunction
+  imports HyperdualFunctionExtension
+begin
+
+subsection\<open>Logistic Function\<close>
+
+text\<open>Define the standard logistic function and its hyperdual variant:\<close>
+definition logistic :: "real \<Rightarrow> real"
+  where "logistic x = inverse (1 + exp (-x))"
+definition hyp_logistic :: "real hyperdual \<Rightarrow> real hyperdual"
+  where "hyp_logistic x = inverse (1 + (*h* exp) (-x))"
+
+text\<open>Hyperdual extension of the logistic function is its hyperdual variant:\<close>
+lemma hypext_logistic:
+  "(*h* logistic) x = hyp_logistic x"
+proof -
+  have "(*h* (\<lambda>x. exp (- x) + 1)) x = (*h* exp) (- x) + of_comp 1"
+    by (simp add: hypext_compose hypext_uminus hypext_fun_cadd twice_field_differentiable_at_compose)
+  then have "(*h* (\<lambda>x. 1 + exp (- x))) x = 1 + (*h* exp) (- x)"
+    by (simp add: one_hyperdual_def add.commute)
+  moreover have "1 + exp (- Base x) \<noteq> 0"
+    by (metis exp_ge_zero add_eq_0_iff neg_0_le_iff_le not_one_le_zero)
+  moreover have "(\<lambda>x. 1 + exp (- x)) twice_field_differentiable_at Base x"
+  proof -
+    have "(\<lambda>x. exp (- x)) twice_field_differentiable_at Base x"
+      by (simp add: twice_field_differentiable_at_compose)
+    then have "(\<lambda>x. exp (- x) + 1) twice_field_differentiable_at Base x"
+      using twice_field_differentiable_at_compose[of "\<lambda>x. exp (- x)" "Base x" "\<lambda>x. x + 1"]
+      by simp
+    then show ?thesis
+      by (simp add: add.commute)
+  qed
+  ultimately have "(*h* (\<lambda>x. inverse (1 + exp (- x)))) x = inverse (1 + (*h* exp) (- x))"
+    by (simp add: hypext_fun_inverse)
+  then show ?thesis
+    unfolding logistic_def hyp_logistic_def .
+qed
+
+text\<open>From properties of autodiff we know it gives us the derivative:\<close>
+lemma "Eps1 (hyp_logistic (\<beta> x)) = deriv logistic x"
+  by (metis Eps1_hypext hypext_logistic)
+text\<open>which is equal to the known derivative of the standard logistic function:\<close>
+lemma "First (autodiff logistic x) = exp (- x) / (1 + exp (- x)) ^ 2"
+  (* Move to hyperdual variant: *)
+  apply (simp only: autodiff.simps hyperdual_to_derivs.simps derivs.sel hypext_logistic)
+  (* Unfold extensions of functions that have a hyperdual variant (all except exp): *)
+  apply (simp only: hyp_logistic_def inverse_hyperdual.code hyperdual.sel)
+  (* Finish by expanding the extension of exp and hyperdual computations: *)
+  apply (simp add: hyperdualx_def hypext_exp_Hyperdual hyperdual_bases)
+  done
+
+text\<open>Similarly we can get the second derivative:\<close>
+lemma "Second (autodiff logistic x) = deriv (deriv logistic) x"
+  by (rule autodiff_extract_second)
+text\<open>and derive its value:\<close>
+lemma "Second (autodiff logistic x) = ((exp (- x) - 1) * exp (- x)) / ((1 + exp (- x)) ^ 3)"
+  (* Move to hyperdual variant: *)
+  apply (simp only: autodiff.simps hyperdual_to_derivs.simps derivs.sel hypext_logistic)
+  (* Unfold extensions of functions that have a hyperdual variant (all except exp): *)
+  apply (simp only: hyp_logistic_def inverse_hyperdual.code hyperdual.sel)
+  (* Finish by expanding the extension of exp and hyperdual computations: *)
+  apply (simp add: hyperdualx_def hypext_exp_Hyperdual hyperdual_bases)
+  (* Simplify the resulting expression: *)
+proof -
+  have
+    "2 * (exp (- x) * exp (- x)) / (1 + exp (- x)) ^ 3 - exp (- x) / (1 + exp (- x)) ^ 2 =
+     (2 * exp (- x) / (1 + exp (- x)) ^ 3 - 1 / (1 + exp (- x)) ^ 2) * exp (- x)"
+    by (simp add: field_simps)
+  also have "... = (2 * exp (- x) / (1 + exp (- x)) ^ 3 - (1 + exp (- x)) / (1 + exp (- x)) ^ 3) * exp (- x)"
+  proof -
+    have "inverse ((1 + exp (- x)) ^ 2) = inverse (1 + exp (- x)) ^ 2"
+      by (simp add: power_inverse)
+    also have "... = (1 + exp (- x)) * inverse (1 + exp (- x)) * inverse (1 + exp (- x)) ^ 2"
+      by (simp add: inverse_eq_divide)
+    also have "... = (1 + exp (- x)) * inverse (1 + exp (- x)) ^ 3"
+      by (simp add: power2_eq_square power3_eq_cube)
+    finally have "inverse ((1 + exp (- x)) ^ 2) = (1 + exp (- x)) * inverse ((1 + exp (- x)) ^ 3)"
+      by (simp add: power_inverse)
+    then show ?thesis
+      by (simp add: inverse_eq_divide)
+  qed
+  also have "... = (2 * exp (- x) - (1 + exp (- x))) / (1 + exp (- x)) ^ 3 * exp (- x)"
+    by (metis diff_divide_distrib)
+  finally show
+    "2 * (exp (- x) * exp (- x)) / (1 + exp (- x)) ^ 3 - exp (- x) / (1 + exp (- x))\<^sup>2 =
+     (exp (- x) - 1) * exp (- x) / (1 + exp (- x)) ^ 3"
+    by (simp add: field_simps)
+qed
+
+end
\ No newline at end of file
diff --git a/thys/Hyperdual/ROOT b/thys/Hyperdual/ROOT
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/ROOT
@@ -0,0 +1,15 @@
+chapter AFP
+
+session "Hyperdual" (AFP) = "HOL-Analysis" +
+  options [timeout = 600]
+  sessions
+    "HOL-Decision_Procs"
+  theories
+    Hyperdual
+    TwiceFieldDifferentiable
+    HyperdualFunctionExtension
+    LogisticFunction
+    AnalyticTestFunction
+  document_files
+    "root.tex"
+    "root.bib"
diff --git a/thys/Hyperdual/TwiceFieldDifferentiable.thy b/thys/Hyperdual/TwiceFieldDifferentiable.thy
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/TwiceFieldDifferentiable.thy
@@ -0,0 +1,741 @@
+(*  Title:   TwiceFieldDifferentiable.thy
+    Authors: Jacques Fleuriot and Filip Smola, University of Edinburgh, 2020
+*)
+
+section \<open>Twice Field Differentiable\<close>
+
+theory TwiceFieldDifferentiable
+  imports "HOL-Analysis.Analysis"
+begin
+
+subsection\<open>Differentiability on a Set\<close>
+
+text\<open>A function is differentiable on a set iff it is differentiable at any point within that set.\<close>
+definition field_differentiable_on :: "('a \<Rightarrow> 'a::real_normed_field) \<Rightarrow> 'a set \<Rightarrow> bool"
+  (infix "field'_differentiable'_on" 50)
+  where "f field_differentiable_on s \<equiv> \<forall>x\<in>s. f field_differentiable (at x within s)"
+
+text\<open>This is preserved for subsets.\<close>
+lemma field_differentiable_on_subset:
+  assumes "f field_differentiable_on S"
+      and "T \<subseteq> S"
+    shows "f field_differentiable_on T"
+  by (meson assms field_differentiable_on_def field_differentiable_within_subset in_mono)
+
+subsection\<open>Twice Differentiability\<close>
+text\<open>
+  Informally, a function is twice differentiable at x iff it is differentiable on some neighbourhood
+  of x and its derivative is differentiable at x.
+\<close>
+definition twice_field_differentiable_at :: "['a \<Rightarrow> 'a::real_normed_field, 'a ] \<Rightarrow> bool"
+  (infixr "(twice'_field'_differentiable'_at)" 50)
+  where "f twice_field_differentiable_at x \<equiv>
+           \<exists>S. f field_differentiable_on S \<and> x \<in> interior S \<and> (deriv f) field_differentiable (at x)"
+
+lemma once_field_differentiable_at:
+  "f twice_field_differentiable_at x \<Longrightarrow> f field_differentiable (at x)"
+  by (metis at_within_interior field_differentiable_on_def interior_subset subsetD twice_field_differentiable_at_def)
+
+lemma deriv_field_differentiable_at:
+  "f twice_field_differentiable_at x \<Longrightarrow> deriv f field_differentiable (at x)"
+  using twice_field_differentiable_at_def by blast
+
+text\<open>
+  For a composition of two functions twice differentiable at x, the chain rule eventually holds on
+  some neighbourhood of x.
+\<close>
+lemma eventually_deriv_compose:
+  assumes "\<exists>S. f field_differentiable_on S \<and> x \<in> interior S"
+      and "g twice_field_differentiable_at (f x)"
+    shows "\<forall>\<^sub>F x in nhds x. deriv (\<lambda>x. g (f x)) x = deriv g (f x) * deriv f x"
+proof -
+  obtain S S'
+   where Df_on_S:  "f field_differentiable_on S" and x_int_S: "x \<in> interior S"
+     and Dg_on_S': "g field_differentiable_on S'" and fx_int_S': "f x \<in> interior S'"
+    using assms twice_field_differentiable_at_def by blast
+
+  let ?T = "{x \<in> interior S. f x \<in> interior S'}"
+
+  have "continuous_on (interior S) f"
+    by (meson Df_on_S continuous_on_eq_continuous_within continuous_on_subset field_differentiable_imp_continuous_at
+         field_differentiable_on_def interior_subset)
+  then have "open (interior S \<inter> {x. f x \<in> interior S'})"
+    by (metis continuous_open_preimage open_interior vimage_def)
+  then have x_int_T: "x \<in> interior ?T"
+    by (metis (no_types) Collect_conj_eq Collect_mem_eq Int_Collect fx_int_S' interior_eq x_int_S)
+  moreover have  Dg_on_fT: "g field_differentiable_on f`?T"
+   by (metis (no_types, lifting) Dg_on_S' field_differentiable_on_subset image_Collect_subsetI interior_subset)
+  moreover have Df_on_T: "f field_differentiable_on ?T"
+    using  field_differentiable_on_subset Df_on_S
+    by (metis Collect_subset interior_subset)
+  moreover have "\<forall>x \<in> interior ?T. deriv (\<lambda>x. g (f x)) x = deriv g (f x) * deriv f x"
+  proof
+    fix x
+    assume x_int_T: "x \<in> interior ?T"
+    have "f field_differentiable at x"
+      by (metis (no_types, lifting) Df_on_T at_within_interior field_differentiable_on_def
+          interior_subset subsetD x_int_T)
+    moreover have "g field_differentiable at (f x)"
+      by (metis (no_types, lifting) Dg_on_S' at_within_interior field_differentiable_on_def
+          interior_subset mem_Collect_eq subsetD x_int_T)
+    ultimately have "deriv (g \<circ> f) x = deriv g (f x) * deriv f x"
+      using deriv_chain[of f x g] by simp
+    then show "deriv (\<lambda>x. g (f x)) x = deriv g (f x) * deriv f x"
+      by (simp add: comp_def)
+  qed
+  ultimately show ?thesis
+    using eventually_nhds by blast
+qed
+
+lemma eventually_deriv_compose':
+  assumes "f twice_field_differentiable_at x"
+      and "g twice_field_differentiable_at (f x)"
+    shows "\<forall>\<^sub>F x in nhds x. deriv (\<lambda>x. g (f x)) x = deriv g (f x) * deriv f x"
+  using assms eventually_deriv_compose twice_field_differentiable_at_def by blast
+
+text\<open>Composition of twice differentiable functions is twice differentiable.\<close>
+lemma twice_field_differentiable_at_compose:
+  assumes "f twice_field_differentiable_at x"
+      and "g twice_field_differentiable_at (f x)"
+    shows "(\<lambda>x. g (f x)) twice_field_differentiable_at x"
+proof -
+  obtain S S'
+   where Df_on_S:  "f field_differentiable_on S" and x_int_S: "x \<in> interior S"
+     and Dg_on_S': "g field_differentiable_on S'" and fx_int_S': "f x \<in> interior S'"
+    using assms twice_field_differentiable_at_def by blast
+
+  let ?T = "{x \<in> interior S. f x \<in> interior S'}"
+
+  have "continuous_on (interior S) f"
+    by (meson Df_on_S continuous_on_eq_continuous_within continuous_on_subset field_differentiable_imp_continuous_at
+         field_differentiable_on_def interior_subset)
+  then have "open (interior S \<inter> {x. f x \<in> interior S'})"
+    by (metis continuous_open_preimage open_interior vimage_def)
+  then have x_int_T: "x \<in> interior ?T"
+    by (metis (no_types) Collect_conj_eq Collect_mem_eq Int_Collect fx_int_S' interior_eq x_int_S)
+
+  have  Dg_on_fT: "g field_differentiable_on f`?T"
+    by (metis (no_types, lifting) Dg_on_S' field_differentiable_on_subset image_Collect_subsetI interior_subset)
+
+  have Df_on_T: "f field_differentiable_on ?T"
+    using  field_differentiable_on_subset Df_on_S
+    by (metis Collect_subset interior_subset)
+
+  have "(\<lambda>x. g (f x)) field_differentiable_on ?T"
+    unfolding field_differentiable_on_def
+  proof
+    fix x assume x_int: "x \<in> {x \<in> interior S. f x \<in> interior S'}"
+    have "f field_differentiable at x"
+      by (metis Df_on_S at_within_interior field_differentiable_on_def interior_subset mem_Collect_eq subsetD x_int)
+    moreover have "g field_differentiable at (f x)"
+      by (metis Dg_on_S' at_within_interior field_differentiable_on_def interior_subset mem_Collect_eq subsetD x_int)
+    ultimately have "(g \<circ> f) field_differentiable at x"
+      by (simp add: field_differentiable_compose)
+    then have "(\<lambda>x. g (f x)) field_differentiable at x"
+      by (simp add: comp_def)
+    then show "(\<lambda>x. g (f x)) field_differentiable at x within {x \<in> interior S. f x \<in> interior S'}"
+      using field_differentiable_at_within by blast
+  qed
+  moreover have "deriv (\<lambda>x. g (f x)) field_differentiable at x"
+  proof -
+    have "(\<lambda>x. deriv g (f x)) field_differentiable at x"
+      by (metis DERIV_chain2 assms deriv_field_differentiable_at field_differentiable_def once_field_differentiable_at)
+    then have "(\<lambda>x. deriv g (f x) * deriv f x) field_differentiable at x"
+      using assms field_differentiable_mult[of "\<lambda>x. deriv g (f x)"]
+      by (simp add: deriv_field_differentiable_at)
+    moreover have "deriv (deriv (\<lambda>x. g (f x))) x = deriv (\<lambda>x. deriv g (f x) * deriv f x) x"
+      using assms Df_on_S x_int_S deriv_cong_ev eventually_deriv_compose by fastforce
+    ultimately show ?thesis
+      using assms eventually_deriv_compose DERIV_deriv_iff_field_differentiable Df_on_S x_int_S
+            DERIV_cong_ev[of x x "deriv (\<lambda>x. g (f x))" "\<lambda>x. deriv g (f x) * deriv f x"]
+      by blast
+  qed
+  ultimately show ?thesis
+    using twice_field_differentiable_at_def x_int_T by blast
+qed
+
+subsubsection\<open>Constant\<close>
+lemma twice_field_differentiable_at_const [simp, intro]:
+  "(\<lambda>x. a) twice_field_differentiable_at x"
+  by (auto intro: exI [of _ UNIV] simp add: twice_field_differentiable_at_def field_differentiable_on_def)
+
+subsubsection\<open>Identity\<close>
+lemma twice_field_differentiable_at_ident [simp, intro]:
+  "(\<lambda>x. x) twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. (\<lambda>x. x) field_differentiable at x"
+   and "deriv ((\<lambda>x. x)) field_differentiable at x"
+    by simp_all
+  then show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by fastforce
+qed
+
+subsubsection\<open>Constant Multiplication\<close>
+lemma twice_field_differentiable_at_cmult [simp, intro]:
+"(*) k twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. (*) k field_differentiable at x"
+    by simp
+  moreover have "deriv ((*) k) field_differentiable at x"
+    by simp
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by fastforce
+qed
+
+lemma twice_field_differentiable_at_uminus [simp, intro]:
+  "uminus twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. uminus field_differentiable at x"
+    by (simp add: field_differentiable_minus)
+  moreover have "deriv uminus field_differentiable at x"
+    by simp
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by auto
+qed
+
+lemma twice_field_differentiable_at_uminus_fun [intro]:
+  assumes "f twice_field_differentiable_at x"
+    shows "(\<lambda>x. - f x) twice_field_differentiable_at x"
+  by (simp add: assms twice_field_differentiable_at_compose)
+
+subsubsection\<open>Real Scaling\<close>
+
+lemma deriv_scaleR_right_id [simp]:
+  "(deriv ((*\<^sub>R) k)) = (\<lambda>z. k *\<^sub>R 1)"
+  using DERIV_imp_deriv has_field_derivative_scaleR_right DERIV_ident by blast
+
+lemma deriv_deriv_scaleR_right_id [simp]:
+  "deriv (deriv ((*\<^sub>R) k)) = (\<lambda>z. 0)"
+  by simp
+
+lemma deriv_scaleR_right:
+  "f field_differentiable (at z) \<Longrightarrow> deriv (\<lambda>x. k *\<^sub>R f x) z = k *\<^sub>R deriv f z"
+  by (simp add: DERIV_imp_deriv field_differentiable_derivI has_field_derivative_scaleR_right)
+
+lemma field_differentiable_scaleR_right [intro]:
+  "f field_differentiable F \<Longrightarrow> (\<lambda>x. c *\<^sub>R f x) field_differentiable F"
+  using field_differentiable_def has_field_derivative_scaleR_right by blast
+
+lemma has_field_derivative_scaleR_deriv_right:
+  assumes "f twice_field_differentiable_at z"
+  shows "((\<lambda>x. k *\<^sub>R deriv f x) has_field_derivative k *\<^sub>R deriv (deriv f) z) (at z)"
+  by (simp add: DERIV_deriv_iff_field_differentiable assms deriv_field_differentiable_at has_field_derivative_scaleR_right)
+
+lemma deriv_scaleR_deriv_right:
+  assumes "f twice_field_differentiable_at z"
+  shows "deriv (\<lambda>x. k *\<^sub>R deriv f x) z = k *\<^sub>R deriv (deriv f) z"
+  using assms deriv_scaleR_right twice_field_differentiable_at_def by blast
+
+lemma twice_field_differentiable_at_scaleR [simp, intro]:
+  "(*\<^sub>R) k twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. (*\<^sub>R) k field_differentiable at x"
+    by (simp add: field_differentiable_scaleR_right)
+  moreover have "deriv ((*\<^sub>R) k) field_differentiable at x"
+    by simp
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by auto
+qed
+
+lemma twice_field_differentiable_at_scaleR_fun [simp, intro]:
+  assumes "f twice_field_differentiable_at x"
+  shows "(\<lambda>x. k *\<^sub>R f x) twice_field_differentiable_at x"
+  by (simp add: assms twice_field_differentiable_at_compose)
+
+subsubsection\<open>Addition\<close>
+
+lemma eventually_deriv_add:
+  assumes "f twice_field_differentiable_at x"
+      and "g twice_field_differentiable_at x"
+    shows "\<forall>\<^sub>F x in nhds x. deriv (\<lambda>x. f x + g x) x = deriv f x + deriv g x"
+proof -
+  obtain S where Df_on_S: "f field_differentiable_on S" and x_int_S: "x \<in> interior S"
+    using assms twice_field_differentiable_at_def by blast
+  obtain S' where Dg_on_S: "g field_differentiable_on S'" and x_int_S': "x \<in> interior S'"
+    using assms twice_field_differentiable_at_def by blast
+  have "x \<in> interior (S \<inter> S')"
+     by (simp add: x_int_S x_int_S')
+  moreover have Df_on_SS': "f field_differentiable_on (S \<inter> S')"
+     by (meson Df_on_S IntD1 field_differentiable_on_def field_differentiable_within_subset inf_sup_ord(1))
+  moreover have Dg_on_SS': "g field_differentiable_on (S \<inter> S')"
+     by (meson Dg_on_S IntD2 field_differentiable_on_def field_differentiable_within_subset inf_le2)
+  moreover have "open (interior (S \<inter> S'))"
+    by blast
+  moreover have "\<forall>x\<in> interior (S \<inter> S'). deriv (\<lambda>x. f x + g x) x = deriv f x + deriv g x"
+    by (metis (full_types) Df_on_SS' Dg_on_SS' at_within_interior deriv_add field_differentiable_on_def in_mono interior_subset)
+  ultimately show ?thesis
+    using eventually_nhds by blast
+qed
+
+lemma twice_field_differentiable_at_add [intro]:
+  assumes "f twice_field_differentiable_at x"
+      and "g twice_field_differentiable_at x"
+    shows "(\<lambda>x. f x + g x) twice_field_differentiable_at x"
+proof -
+  obtain S S'
+   where Df_on_S:  "f field_differentiable_on S" and x_int_S: "x \<in> interior S"
+     and Dg_on_S': "g field_differentiable_on S'" and x_int_S': " x \<in> interior S'"
+    using assms twice_field_differentiable_at_def by blast
+
+  let ?T = "interior (S \<inter> S')"
+
+  have x_int_T: "x \<in> interior ?T"
+     by (simp add: x_int_S x_int_S')
+
+  have Df_on_T: "f field_differentiable_on ?T"
+    by (meson Df_on_S field_differentiable_on_subset inf_sup_ord(1) interior_subset)
+  have  Dg_on_fT: "g field_differentiable_on ?T"
+    by (meson Dg_on_S' field_differentiable_on_subset interior_subset le_infE)
+
+   have "(\<lambda>x. f x + g x) field_differentiable_on ?T"
+     unfolding field_differentiable_on_def
+  proof
+    fix x assume x_in_T: "x \<in> ?T"
+    have "f field_differentiable at x"
+      by (metis x_in_T Df_on_T at_within_open field_differentiable_on_def open_interior)
+    moreover have "g field_differentiable at x"
+      by (metis x_in_T Dg_on_fT at_within_open field_differentiable_on_def open_interior)
+    ultimately have "(\<lambda>x. f x + g x) field_differentiable at x"
+      by (simp add: field_differentiable_add)
+    then show "(\<lambda>x. f x + g x) field_differentiable at x within ?T"
+      using field_differentiable_at_within by blast
+  qed
+  moreover have "deriv (\<lambda>x. f x + g x) field_differentiable at x"
+  proof -
+    have "deriv (\<lambda>x. f x + g x) x = deriv f x + deriv g x"
+      by (simp add: assms once_field_differentiable_at)
+    moreover have "(\<lambda>x. deriv f x + deriv g x) field_differentiable at x"
+      by (simp add: field_differentiable_add assms deriv_field_differentiable_at)
+    ultimately show ?thesis
+      using assms DERIV_deriv_iff_field_differentiable
+            DERIV_cong_ev[of x x "deriv (\<lambda>x. f x + g x)" "\<lambda>x. deriv f x + deriv g x"]
+      by (simp add: eventually_deriv_add field_differentiable_def)
+  qed
+  ultimately show ?thesis
+    using twice_field_differentiable_at_def x_int_T by blast
+qed
+
+lemma deriv_add_id_const [simp]:
+  "deriv (\<lambda>x. x + a) = (\<lambda>z. 1)"
+  using ext trans[OF deriv_add] by force
+
+lemma deriv_deriv_add_id_const [simp]:
+  "deriv (deriv (\<lambda>x. x + a)) z = 0"
+  by simp
+
+lemma twice_field_differentiable_at_cadd [simp]:
+  "(\<lambda>x. x + a) twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. (\<lambda>x. x + a) field_differentiable at x"
+    by (simp add: field_differentiable_add)
+  moreover have "deriv ((\<lambda>x. x + a)) field_differentiable at x"
+    by (simp add: ext)
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by auto
+qed
+
+subsubsection\<open>Linear Function\<close>
+lemma twice_field_differentiable_at_linear [simp, intro]:
+  "(\<lambda>x. k * x + a) twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. (\<lambda>x. k * x + a) field_differentiable at x"
+    by (simp add: field_differentiable_add)
+  moreover have "deriv ((\<lambda>x. k * x + a)) field_differentiable at x"
+  proof -
+    have "deriv ((\<lambda>x. k * x + a)) = (\<lambda>x. k)"
+      by (simp add: ext)
+    then show ?thesis
+      by simp
+  qed
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by auto
+qed
+
+lemma twice_field_differentiable_at_linearR [simp, intro]:
+  "(\<lambda>x. k *\<^sub>R x + a) twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. (\<lambda>x. k *\<^sub>R x + a) field_differentiable at x"
+    by (simp add: field_differentiable_scaleR_right field_differentiable_add)
+  moreover have "deriv ((\<lambda>x. k *\<^sub>R x + a)) field_differentiable at x"
+  proof -
+    have "deriv ((\<lambda>x. k *\<^sub>R x + a)) = (\<lambda>x. k *\<^sub>R 1)"
+      by (simp add: ext once_field_differentiable_at)
+    then show ?thesis
+      by simp
+  qed
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by auto
+qed
+
+subsubsection\<open>Multiplication\<close>
+
+lemma eventually_deriv_mult:
+  assumes "f twice_field_differentiable_at x"
+      and "g twice_field_differentiable_at x"
+    shows "\<forall>\<^sub>F x in nhds x. deriv (\<lambda>x. f x * g x) x = f x * deriv g x + deriv f x * g x"
+proof -
+  obtain S and S'
+    where "f field_differentiable_on S"   and in_S:  "x \<in> interior S"
+      and "g field_differentiable_on S'"  and in_S': "x \<in> interior S'"
+    using assms twice_field_differentiable_at_def by blast
+  then have Df_on_SS': "f field_differentiable_on (S \<inter> S')"
+        and Dg_on_SS': "g field_differentiable_on (S \<inter> S')"
+    using field_differentiable_on_subset by blast+
+
+  have "\<forall>x\<in> interior (S \<inter> S'). deriv (\<lambda>x. f x * g x) x = f x * deriv g x + deriv f x * g x"
+  proof
+    fix x assume "x \<in> interior (S \<inter> S')"
+    then have "f field_differentiable (at x)"
+          and "g field_differentiable (at x)"
+      using Df_on_SS' Dg_on_SS' field_differentiable_on_def at_within_interior interior_subset subsetD by metis+
+    then show "deriv (\<lambda>x. f x * g x) x = f x * deriv g x + deriv f x * g x"
+      by simp
+  qed
+  moreover have "x \<in> interior (S \<inter> S')"
+    by (simp add: in_S in_S')
+  moreover have "open (interior (S \<inter> S'))"
+    by blast
+  ultimately show ?thesis
+    using eventually_nhds by blast
+qed
+
+lemma twice_field_differentiable_at_mult [intro]:
+  assumes "f twice_field_differentiable_at x"
+      and "g twice_field_differentiable_at x"
+    shows "(\<lambda>x. f x * g x) twice_field_differentiable_at x"
+proof -
+  obtain S S'
+   where Df_on_S:  "f field_differentiable_on S" and x_int_S: "x \<in> interior S"
+     and Dg_on_S': "g field_differentiable_on S'" and x_int_S': " x \<in> interior S'"
+    using assms twice_field_differentiable_at_def by blast
+
+  let ?T = "interior (S \<inter> S')"
+
+  have x_int_T: "x \<in> interior ?T"
+     by (simp add: x_int_S x_int_S')
+
+  have Df_on_T: "f field_differentiable_on ?T"
+    by (meson Df_on_S field_differentiable_on_subset inf_sup_ord(1) interior_subset)
+  have  Dg_on_fT: "g field_differentiable_on ?T"
+    by (meson Dg_on_S' field_differentiable_on_subset interior_subset le_infE)
+
+   have "(\<lambda>x. f x * g x) field_differentiable_on ?T"
+     unfolding field_differentiable_on_def
+  proof
+    fix x assume x_in_T: "x \<in> ?T"
+    have "f field_differentiable at x"
+      by (metis x_in_T Df_on_T at_within_open field_differentiable_on_def open_interior)
+    moreover have "g field_differentiable at x"
+      by (metis x_in_T Dg_on_fT at_within_open field_differentiable_on_def open_interior)
+    ultimately have "(\<lambda>x. f x * g x) field_differentiable at x"
+      by (simp add: field_differentiable_mult)
+    then show "(\<lambda>x. f x * g x) field_differentiable at x within ?T"
+      using field_differentiable_at_within by blast
+  qed
+  moreover have "deriv (\<lambda>x. f x * g x) field_differentiable at x"
+  proof -
+    have "deriv (\<lambda>x. f x * g x) x = f x * deriv g x + deriv f x * g x"
+      by (simp add: assms once_field_differentiable_at)
+    moreover have "(\<lambda>x. f x * deriv g x + deriv f x * g x) field_differentiable at x"
+      by (rule field_differentiable_add, simp_all add: field_differentiable_mult assms once_field_differentiable_at deriv_field_differentiable_at)
+    ultimately show ?thesis
+      using assms DERIV_deriv_iff_field_differentiable
+            DERIV_cong_ev[of x x "deriv (\<lambda>x. f x * g x)" "\<lambda>x. f x * deriv g x + deriv f x * g x"]
+      by (simp add: eventually_deriv_mult field_differentiable_def)
+  qed
+  ultimately show ?thesis
+    using twice_field_differentiable_at_def x_int_T by blast
+qed
+
+subsubsection\<open>Sine and Cosine\<close>
+
+lemma deriv_sin [simp]: "deriv sin a = cos a"
+  by (simp add: DERIV_imp_deriv)
+
+lemma deriv_sinf [simp]: "deriv sin = (\<lambda>x. cos x)"
+  by auto
+
+lemma deriv_cos [simp]: "deriv cos a = - sin a"
+  by (simp add: DERIV_imp_deriv)
+
+lemma deriv_cosf [simp]: "deriv cos = (\<lambda>x. - sin x)"
+  by auto
+
+lemma deriv_sin_minus [simp]:
+  "deriv (\<lambda>x. - sin x) a = - deriv (\<lambda>x. sin x) a"
+  by (simp add: DERIV_imp_deriv Deriv.field_differentiable_minus)
+
+lemma twice_field_differentiable_at_sin [simp, intro]:
+  "sin twice_field_differentiable_at x"
+  by (auto intro!: exI [of _ UNIV] simp add: field_differentiable_at_sin
+      field_differentiable_on_def twice_field_differentiable_at_def field_differentiable_at_cos)
+
+lemma twice_field_differentiable_at_sin_fun [intro]:
+  assumes "f twice_field_differentiable_at x"
+  shows   "(\<lambda>x. sin (f x)) twice_field_differentiable_at x"
+  by (simp add: assms twice_field_differentiable_at_compose)
+
+lemma twice_field_differentiable_at_cos [simp, intro]:
+  "cos twice_field_differentiable_at x"
+  by (auto intro!: exI [of _ UNIV] simp add:  field_differentiable_within_sin field_differentiable_minus
+      field_differentiable_on_def twice_field_differentiable_at_def field_differentiable_at_cos)
+
+lemma twice_field_differentiable_at_cos_fun [intro]:
+  assumes "f twice_field_differentiable_at x"
+  shows   "(\<lambda>x. cos (f x)) twice_field_differentiable_at x"
+  by (simp add: assms twice_field_differentiable_at_compose)
+
+subsubsection\<open>Exponential\<close>
+
+lemma deriv_exp [simp]: "deriv exp x = exp x"
+  using DERIV_exp DERIV_imp_deriv by blast
+
+lemma deriv_expf [simp]: "deriv exp = exp"
+  by (simp add: ext)
+
+lemma deriv_deriv_exp [simp]: "deriv (deriv exp) x = exp x"
+  by simp
+
+lemma twice_field_differentiable_at_exp [simp, intro]:
+  "exp twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. exp field_differentiable at x"
+   and "deriv exp field_differentiable at x"
+    by (simp_all add: field_differentiable_within_exp)
+  then show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by auto
+qed
+
+lemma twice_field_differentiable_at_exp_fun [simp, intro]:
+  assumes "f twice_field_differentiable_at x"
+  shows "(\<lambda>x. exp (f x)) twice_field_differentiable_at x"
+  by (simp add: assms twice_field_differentiable_at_compose)
+
+subsubsection\<open>Square Root\<close>
+
+lemma deriv_real_sqrt [simp]: "x > 0 \<Longrightarrow> deriv sqrt x = inverse (sqrt x) / 2"
+  using DERIV_imp_deriv DERIV_real_sqrt by blast
+
+lemma has_real_derivative_inverse_sqrt:
+  assumes "x > 0"
+  shows "((\<lambda>x. inverse (sqrt x) / 2) has_real_derivative - (inverse (sqrt x ^ 3) / 4)) (at x)"
+proof -
+  have inv_sqrt_mult: "(inverse (sqrt x)/2) * (sqrt x * 2) = 1"
+    using assms by simp
+  have inv_sqrt_mult2: "(- inverse ((sqrt x)^3)/2)* x * (sqrt x * 2) = -1"
+    using assms by (simp add: field_simps power3_eq_cube)
+  then show ?thesis
+    using assms by (safe intro!: DERIV_imp_deriv derivative_eq_intros)
+        (auto intro: derivative_eq_intros inv_sqrt_mult [THEN ssubst] inv_sqrt_mult2 [THEN ssubst]
+              simp add: divide_simps power3_eq_cube)
+qed
+
+lemma deriv_deriv_real_sqrt':
+  assumes "x > 0"
+  shows "deriv (\<lambda>x. inverse (sqrt x) / 2) x = - inverse ((sqrt x)^3)/4"
+  by (simp add: DERIV_imp_deriv assms has_real_derivative_inverse_sqrt)
+
+lemma has_real_derivative_deriv_sqrt:
+  assumes "x > 0"
+  shows "(deriv sqrt has_real_derivative - inverse (sqrt x ^ 3) / 4) (at x)"
+proof -
+  have "((\<lambda>x. inverse (sqrt x) / 2) has_real_derivative - inverse (sqrt x ^ 3) / 4) (at x)"
+    using assms has_real_derivative_inverse_sqrt by auto
+  moreover
+  {fix xa :: real
+   assume "xa \<in> {0<..}"
+   then have "inverse (sqrt xa) / 2 = deriv sqrt xa"
+     by simp
+  }
+  ultimately show ?thesis
+    using has_field_derivative_transform_within_open [where S="{0 <..}" and f="(\<lambda>x. inverse (sqrt x) / 2)"]
+    by (meson assms greaterThan_iff open_greaterThan)
+qed
+
+lemma deriv_deriv_real_sqrt [simp]:
+  assumes "x > 0"
+  shows "deriv(deriv sqrt) x = - inverse ((sqrt x)^3)/4"
+  using DERIV_imp_deriv assms has_real_derivative_deriv_sqrt by blast
+
+lemma twice_field_differentiable_at_sqrt [simp, intro]:
+  assumes "x > 0"
+    shows "sqrt twice_field_differentiable_at x"
+proof -
+  have "sqrt field_differentiable_on {0<..}"
+    by (metis DERIV_real_sqrt at_within_open field_differentiable_def field_differentiable_on_def
+        greaterThan_iff open_greaterThan)
+  moreover have "x \<in> interior {0<..}"
+    by (metis assms greaterThan_iff interior_interior interior_real_atLeast)
+  moreover have "deriv sqrt field_differentiable at x"
+    using assms field_differentiable_def has_real_derivative_deriv_sqrt by blast
+  ultimately show ?thesis
+    using twice_field_differentiable_at_def by blast
+qed
+
+lemma twice_field_differentiable_at_sqrt_fun [intro]:
+  assumes "f twice_field_differentiable_at x"
+    and "f x > 0"
+  shows "(\<lambda>x. sqrt (f x)) twice_field_differentiable_at x"
+  by (simp add: assms(1) assms(2) twice_field_differentiable_at_compose)
+
+subsubsection\<open>Natural Power\<close>
+
+lemma field_differentiable_power [simp]:
+  "(\<lambda>x. x ^ n) field_differentiable at x"
+  using DERIV_power DERIV_ident field_differentiable_def
+  by blast
+
+lemma deriv_power_fun [simp]:
+  assumes "f field_differentiable at x"
+    shows "deriv (\<lambda>x. f x ^ n) x = of_nat n * deriv f x * f x ^ (n - 1)"
+  using DERIV_power[of f "deriv f x"]
+  by (simp add:  DERIV_imp_deriv assms field_differentiable_derivI mult.assoc [symmetric])
+
+lemma deriv_power [simp]:
+  "deriv (\<lambda>x. x ^ n) x = of_nat n * x ^ (n - 1)"
+  using DERIV_power[of "\<lambda>x. x" 1] DERIV_imp_deriv by force
+
+lemma deriv_deriv_power [simp]:
+  "deriv (deriv (\<lambda>x. x ^ n)) x = of_nat n * of_nat (n - Suc 0) * x ^ (n - 2)"
+proof -
+  have "(\<lambda>x. x ^ (n - 1)) field_differentiable at x"
+    by simp
+  then have "deriv (\<lambda>x. of_nat n * x ^ (n - 1)) x = of_nat n * of_nat (n - Suc 0) * x ^ (n - 2)"
+    by (simp add: diff_diff_add mult.assoc numeral_2_eq_2)
+  then show ?thesis
+    by (simp add: ext[OF deriv_power])
+qed
+
+lemma twice_field_differentiable_at_power [simp, intro]:
+  "(\<lambda>x. x ^ n) twice_field_differentiable_at x"
+proof -
+  have "\<forall>x\<in>UNIV. (\<lambda>x. x ^ n) field_differentiable at x"
+    by simp
+  moreover have "deriv ((\<lambda>x. x ^ n)) field_differentiable at x"
+  proof -
+    have "deriv ((\<lambda>x. x ^ n)) = (\<lambda>x. of_nat n * x ^ (n - 1))"
+      by (simp add: ext)
+    then show ?thesis
+      using field_differentiable_mult[of "\<lambda>x. of_nat n" x UNIV "\<lambda>x. x ^ (n - 1)"]
+      by (simp add: field_differentiable_caratheodory_at)
+  qed
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by force
+qed
+
+lemma twice_field_differentiable_at_power_fun [intro]:
+  assumes "f twice_field_differentiable_at x"
+    shows "(\<lambda>x. f x ^ n) twice_field_differentiable_at x"
+  by (blast intro: assms twice_field_differentiable_at_compose [OF _ twice_field_differentiable_at_power])
+
+subsubsection\<open>Inverse\<close>
+
+lemma eventually_deriv_inverse:
+  assumes "x \<noteq> 0"
+    shows "\<forall>\<^sub>F x in nhds x. deriv inverse x = - 1 / (x ^ 2)"
+proof -
+  obtain T where open_T: "open T" and "\<forall>z\<in>T. z \<noteq> 0" and x_in_T: "x \<in> T"
+    using assms t1_space by blast
+
+  then have "\<forall>x \<in> T. deriv inverse x = - 1 / (x ^ 2)"
+    by simp
+  then show ?thesis
+    using eventually_nhds open_T x_in_T by blast
+qed
+
+lemma deriv_deriv_inverse [simp]:
+  assumes "x \<noteq> 0"
+  shows "deriv (deriv inverse) x = 2 * inverse (x ^ 3)"
+proof -
+  have "deriv (\<lambda>x. inverse (x ^ 2)) x = - (of_nat 2 * x) / ((x ^ 2) ^ 2)"
+    using assms by simp
+  moreover have "(\<lambda>x. inverse (x ^ 2)) field_differentiable at x"
+    using assms by (simp add: field_differentiable_inverse)
+  ultimately have "deriv (\<lambda>x. - (inverse (x ^ 2))) x = of_nat 2 * x / (x ^ 4)"
+    using deriv_chain[of "\<lambda>x. inverse (x ^ 2)" x]
+    by (simp add: comp_def field_differentiable_minus field_simps)
+  then have "deriv (\<lambda>x. - 1 / (x ^ 2)) x = 2 * inverse (x ^ 3)"
+    by (simp add: power4_eq_xxxx power3_eq_cube field_simps)
+  then show ?thesis
+    using assms eventually_deriv_inverse deriv_cong_ev by fastforce
+qed
+
+lemma twice_field_differentiable_at_inverse [simp, intro]:
+  assumes "x \<noteq> 0"
+  shows "inverse twice_field_differentiable_at x"
+proof -
+  obtain T where zero_T: "0 \<notin> T" and x_in_T: "x \<in> T" and open_T: "open T"
+    using assms t1_space by blast
+  then have "T \<subseteq> {z. z \<noteq> 0}"
+    by blast
+  then have "\<forall>x\<in>T. inverse field_differentiable at x within T"
+    using DERIV_inverse field_differentiable_def by blast
+  moreover have "deriv inverse field_differentiable at x"
+  proof -
+    have "(\<lambda>x. - inverse (x ^ 2)) field_differentiable at x"
+      using assms by (simp add: field_differentiable_inverse field_differentiable_minus)
+    then have "(\<lambda>x. - 1 / (x ^ 2)) field_differentiable at x"
+      by (simp add: inverse_eq_divide)
+    then show ?thesis
+      using eventually_deriv_inverse[OF assms]
+      by (simp add: DERIV_cong_ev field_differentiable_def)
+  qed
+  moreover have "x \<in> interior T"
+    by (simp add: x_in_T open_T interior_open)
+  ultimately show ?thesis
+    unfolding twice_field_differentiable_at_def field_differentiable_on_def
+    by blast
+qed
+
+lemma twice_field_differentiable_at_inverse_fun [simp, intro]:
+  assumes "f twice_field_differentiable_at x"
+          "f x \<noteq> 0"
+  shows "(\<lambda>x. inverse (f x)) twice_field_differentiable_at x"
+  by (simp add: assms twice_field_differentiable_at_compose)
+
+lemma twice_field_differentiable_at_divide [intro]:
+  assumes "f twice_field_differentiable_at x"
+      and "g twice_field_differentiable_at x"
+      and "g x \<noteq> 0"
+    shows "(\<lambda>x. f x / g x) twice_field_differentiable_at x"
+  by (simp add: assms divide_inverse twice_field_differentiable_at_mult)
+
+subsubsection\<open>Polynomial\<close>
+
+lemma twice_field_differentiable_at_polyn [simp, intro]:
+  fixes coef :: "nat \<Rightarrow> 'a :: {real_normed_field}"
+    and n :: nat
+  shows "(\<lambda>x. \<Sum>i<n. coef i * x ^ i) twice_field_differentiable_at x"
+proof (induction n)
+  case 0
+  then show ?case
+    by simp
+next
+  case hyp: (Suc n)
+  show ?case
+  proof (simp, rule twice_field_differentiable_at_add)
+    show "(\<lambda>x. \<Sum>i<n. coef i * x ^ i) twice_field_differentiable_at x"
+      by (rule hyp)
+    show "(\<lambda>x. coef n * x ^ n) twice_field_differentiable_at x"
+      using twice_field_differentiable_at_compose[of "\<lambda>x. x ^ n" x "(*) (coef n)"]
+      by simp
+  qed
+qed
+
+lemma twice_field_differentiable_at_polyn_fun [simp]:
+  fixes coef :: "nat \<Rightarrow> 'a :: {real_normed_field}"
+    and n :: nat
+  assumes "f twice_field_differentiable_at x"
+  shows "(\<lambda>x. \<Sum>i<n. coef i * f x ^ i) twice_field_differentiable_at x"
+  by (blast intro: assms twice_field_differentiable_at_compose [OF _ twice_field_differentiable_at_polyn])
+
+end
diff --git a/thys/Hyperdual/document/root.bib b/thys/Hyperdual/document/root.bib
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/document/root.bib
@@ -0,0 +1,6 @@
+@InProceedings{fike_alonso-2011,
+  author={{Fike}, J.~A. and {Alonso}, J.~J.},
+  title={The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations},
+  booktitle={AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting},
+  year=2011
+}
diff --git a/thys/Hyperdual/document/root.tex b/thys/Hyperdual/document/root.tex
new file mode 100644
--- /dev/null
+++ b/thys/Hyperdual/document/root.tex
@@ -0,0 +1,70 @@
+\documentclass[11pt,a4paper]{article}
+\usepackage{isabelle,isabellesym}
+
+% further packages required for unusual symbols (see also
+% isabellesym.sty), use only when needed
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+  %for \<leadsto>, \<box>, \<diamond>, \<sqsupset>, \<mho>, \<Join>,
+  %\<lhd>, \<lesssim>, \<greatersim>, \<lessapprox>, \<greaterapprox>,
+  %\<triangleq>, \<yen>, \<lozenge>
+
+%\usepackage{eurosym}
+  %for \<euro>
+
+%\usepackage[only,bigsqcap]{stmaryrd}
+  %for \<Sqinter>
+
+%\usepackage{eufrak}
+  %for \<AA> ... \<ZZ>, \<aa> ... \<zz> (also included in amssymb)
+
+%\usepackage{textcomp}
+  %for \<onequarter>, \<onehalf>, \<threequarters>, \<degree>, \<cent>,
+  %\<currency>
+
+% this should be the last package used
+\usepackage{pdfsetup}
+
+% urls in roman style, theory text in math-similar italics
+\urlstyle{rm}
+\isabellestyle{it}
+
+% for uniform font size
+%\renewcommand{\isastyle}{\isastyleminor}
+
+
+\begin{document}
+
+\title{Hyperdual Numbers and Forward Differentiation}
+\author{Filip Smola, Jacques Fleuriot}
+\maketitle
+
+\begin{abstract}
+  Hyperdual numbers are ones with a real component and a number of infinitesimal components, usually written as $a_0 + a_1 \cdot \epsilon_1 + a_2 \cdot \epsilon_2 + a_3 \cdot \epsilon_1\epsilon_2$.
+  They have been proposed by Fike and Alonso~\cite{fike_alonso-2011} in an approach to automatic differentiation.
+
+  In this entry we formalise hyperdual numbers and their application to forward differentiation.
+  We show them to be an instance of multiple algebraic structures and then, along with facts about twice-differentiability, we define what we call the hyperdual extensions of functions on real-normed fields.
+  This extension formally represents the proposed way that the first and second derivatives of a function can be automatically calculated.
+  We demonstrate it on the standard logistic function $f(x) = \frac{1}{1 + e^{-x}}$ and also reproduce the example analytic function $f(x) = \frac{e^x}{\sqrt{sin(x)^3 + cos(x)^3}}$ used for demonstration by Fike and Alonso.
+\end{abstract}
+
+\tableofcontents
+
+% sane default for proof documents
+\parindent 0pt\parskip 0.5ex
+
+% generated text of all theories
+\input{session}
+
+% optional bibliography
+\bibliographystyle{abbrv}
+\bibliography{root}
+
+\end{document}
+
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: t
+%%% End:
diff --git a/thys/ROOTS b/thys/ROOTS
--- a/thys/ROOTS
+++ b/thys/ROOTS
@@ -1,651 +1,652 @@
 ADS_Functor
 AI_Planning_Languages_Semantics
 AODV
 AVL-Trees
 AWN
 Abortable_Linearizable_Modules
 Abs_Int_ITP2012
 Abstract-Hoare-Logics
 Abstract-Rewriting
 Abstract_Completeness
 Abstract_Soundness
 Adaptive_State_Counting
 Affine_Arithmetic
 Aggregation_Algebras
 Akra_Bazzi
 Algebraic_Numbers
 Algebraic_VCs
 Allen_Calculus
 Amicable_Numbers
 Amortized_Complexity
 AnselmGod
 Applicative_Lifting
 Approximation_Algorithms
 Architectural_Design_Patterns
 Aristotles_Assertoric_Syllogistic
 Arith_Prog_Rel_Primes
 ArrowImpossibilityGS
 Attack_Trees
 Auto2_HOL
 Auto2_Imperative_HOL
 AutoFocus-Stream
 Automated_Stateful_Protocol_Verification
 Automatic_Refinement
 AxiomaticCategoryTheory
 BDD
 BD_Security_Compositional
 BNF_CC
 BNF_Operations
 BTree
 Banach_Steinhaus
 Belief_Revision
 Bell_Numbers_Spivey
 BenOr_Kozen_Reif
 Berlekamp_Zassenhaus
 Bernoulli
 Bertrands_Postulate
 Bicategory
 BinarySearchTree
 Binding_Syntax_Theory
 Binomial-Heaps
 Binomial-Queues
 BirdKMP
 Blue_Eyes
 Bondy
 Boolean_Expression_Checkers
 Bounded_Deducibility_Security
 Buchi_Complementation
 Budan_Fourier
 Buffons_Needle
 Buildings
 BytecodeLogicJmlTypes
 C2KA_DistributedSystems
 CAVA_Automata
 CAVA_LTL_Modelchecker
 CCS
 CISC-Kernel
 CRDT
 CSP_RefTK
 CYK
 CZH_Elementary_Categories
 CZH_Foundations
 CZH_Universal_Constructions
 CakeML
 CakeML_Codegen
 Call_Arity
 Card_Equiv_Relations
 Card_Multisets
 Card_Number_Partitions
 Card_Partitions
 Cartan_FP
 Case_Labeling
 Catalan_Numbers
 Category
 Category2
 Category3
 Cauchy
 Cayley_Hamilton
 Certification_Monads
 Chandy_Lamport
 Chord_Segments
 Circus
 Clean
 ClockSynchInst
 Closest_Pair_Points
 CoCon
 CoSMeDis
 CoSMed
 CofGroups
 Coinductive
 Coinductive_Languages
 Collections
 Combinatorics_Words
 Combinatorics_Words_Graph_Lemma
 Combinatorics_Words_Lyndon
 Comparison_Sort_Lower_Bound
 Compiling-Exceptions-Correctly
 Complete_Non_Orders
 Completeness
 Complex_Bounded_Operators
 Complex_Geometry
 Complx
 ComponentDependencies
 ConcurrentGC
 ConcurrentIMP
 Concurrent_Ref_Alg
 Concurrent_Revisions
 Conditional_Simplification
 Conditional_Transfer_Rule
 Consensus_Refined
 Constructive_Cryptography
 Constructive_Cryptography_CM
 Constructor_Funs
 Containers
 CoreC++
 Core_DOM
 Core_SC_DOM
 Correctness_Algebras
 Count_Complex_Roots
 CryptHOL
 CryptoBasedCompositionalProperties
 Cubic_Quartic_Equations
 DFS_Framework
 DOM_Components
 DPT-SAT-Solver
 DataRefinementIBP
 Datatype_Order_Generator
 Decl_Sem_Fun_PL
 Decreasing-Diagrams
 Decreasing-Diagrams-II
 Deep_Learning
 Delta_System_Lemma
 Density_Compiler
 Dependent_SIFUM_Refinement
 Dependent_SIFUM_Type_Systems
 Depth-First-Search
 Derangements
 Deriving
 Descartes_Sign_Rule
 Design_Theory
 Dict_Construction
 Differential_Dynamic_Logic
 Differential_Game_Logic
 Dijkstra_Shortest_Path
 Diophantine_Eqns_Lin_Hom
 Dirichlet_L
 Dirichlet_Series
 DiscretePricing
 Discrete_Summation
 DiskPaxos
 Dominance_CHK
 DynamicArchitectures
 Dynamic_Tables
 E_Transcendental
 Echelon_Form
 EdmondsKarp_Maxflow
 Efficient-Mergesort
 Elliptic_Curves_Group_Law
 Encodability_Process_Calculi
 Epistemic_Logic
 Ergodic_Theory
 Error_Function
 Euler_MacLaurin
 Euler_Partition
 Example-Submission
 Extended_Finite_State_Machine_Inference
 Extended_Finite_State_Machines
 FFT
 FLP
 FOL-Fitting
 FOL_Axiomatic
 FOL_Harrison
 FOL_Seq_Calc1
 Factor_Algebraic_Polynomial
 Factored_Transition_System_Bounding
 Falling_Factorial_Sum
 Farkas
 FeatherweightJava
 Featherweight_OCL
 Fermat3_4
 FileRefinement
 FinFun
 Finger-Trees
 Finite-Map-Extras
 Finite_Automata_HF
 Finitely_Generated_Abelian_Groups
 First_Order_Terms
 First_Welfare_Theorem
 Fishburn_Impossibility
 Fisher_Yates
 Flow_Networks
 Floyd_Warshall
 Flyspeck-Tame
 FocusStreamsCaseStudies
 Forcing
 Formal_Puiseux_Series
 Formal_SSA
 Formula_Derivatives
 Foundation_of_geometry
 Fourier
 Free-Boolean-Algebra
 Free-Groups
 Fresh_Identifiers
 FunWithFunctions
 FunWithTilings
 Functional-Automata
 Functional_Ordered_Resolution_Prover
 Furstenberg_Topology
 GPU_Kernel_PL
 Gabow_SCC
 GaleStewart_Games
 Gale_Shapley
 Game_Based_Crypto
 Gauss-Jordan-Elim-Fun
 Gauss_Jordan
 Gauss_Sums
 Gaussian_Integers
 GenClock
 General-Triangle
 Generalized_Counting_Sort
 Generic_Deriving
 Generic_Join
 GewirthPGCProof
 Girth_Chromatic
 GoedelGod
 Goedel_HFSet_Semantic
 Goedel_HFSet_Semanticless
 Goedel_Incompleteness
 Goodstein_Lambda
 GraphMarkingIBP
 Graph_Saturation
 Graph_Theory
 Green
 Groebner_Bases
 Groebner_Macaulay
 Gromov_Hyperbolicity
 Grothendieck_Schemes
 Group-Ring-Module
 HOL-CSP
 HOLCF-Prelude
 HRB-Slicing
 Hahn_Jordan_Decomposition
 Heard_Of
 Hello_World
 HereditarilyFinite
 Hermite
 Hermite_Lindemann
 Hidden_Markov_Models
 Higher_Order_Terms
 Hoare_Time
 Hood_Melville_Queue
 HotelKeyCards
 Huffman
 Hybrid_Logic
 Hybrid_Multi_Lane_Spatial_Logic
 Hybrid_Systems_VCs
 HyperCTL
+Hyperdual
 IEEE_Floating_Point
 IFC_Tracking
 IMAP-CRDT
 IMO2019
 IMP2
 IMP2_Binary_Heap
 IMP_Compiler
 IP_Addresses
 Imperative_Insertion_Sort
 Impossible_Geometry
 Incompleteness
 Incredible_Proof_Machine
 Inductive_Confidentiality
 Inductive_Inference
 InfPathElimination
 InformationFlowSlicing
 InformationFlowSlicing_Inter
 Integration
 Interpreter_Optimizations
 Interval_Arithmetic_Word32
 Intro_Dest_Elim
 Iptables_Semantics
 Irrational_Series_Erdos_Straus
 Irrationality_J_Hancl
 IsaGeoCoq
 Isabelle_C
 Isabelle_Marries_Dirac
 Isabelle_Meta_Model
 Jacobson_Basic_Algebra
 Jinja
 JinjaDCI
 JinjaThreads
 JiveDataStoreModel
 Jordan_Hoelder
 Jordan_Normal_Form
 KAD
 KAT_and_DRA
 KBPs
 KD_Tree
 Key_Agreement_Strong_Adversaries
 Kleene_Algebra
 Knights_Tour
 Knot_Theory
 Knuth_Bendix_Order
 Knuth_Morris_Pratt
 Koenigsberg_Friendship
 Kruskal
 Kuratowski_Closure_Complement
 LLL_Basis_Reduction
 LLL_Factorization
 LOFT
 LTL
 LTL_Master_Theorem
 LTL_Normal_Form
 LTL_to_DRA
 LTL_to_GBA
 Lam-ml-Normalization
 LambdaAuth
 LambdaMu
 Lambda_Free_EPO
 Lambda_Free_KBOs
 Lambda_Free_RPOs
 Lambert_W
 Landau_Symbols
 Laplace_Transform
 Latin_Square
 LatticeProperties
 Launchbury
 Laws_of_Large_Numbers
 Lazy-Lists-II
 Lazy_Case
 Lehmer
 Lifting_Definition_Option
 Lifting_the_Exponent
 LightweightJava
 LinearQuantifierElim
 Linear_Inequalities
 Linear_Programming
 Linear_Recurrences
 Liouville_Numbers
 List-Index
 List-Infinite
 List_Interleaving
 List_Inversions
 List_Update
 LocalLexing
 Localization_Ring
 Locally-Nameless-Sigma
 Logging_Independent_Anonymity
 Lowe_Ontological_Argument
 Lower_Semicontinuous
 Lp
 Lucas_Theorem
 MDP-Algorithms
 MDP-Rewards
 MFMC_Countable
 MFODL_Monitor_Optimized
 MFOTL_Monitor
 MSO_Regex_Equivalence
 Markov_Models
 Marriage
 Mason_Stothers
 Matrices_for_ODEs
 Matrix
 Matrix_Tensor
 Matroids
 Max-Card-Matching
 Median_Of_Medians_Selection
 Menger
 Mereology
 Mersenne_Primes
 Metalogic_ProofChecker
 MiniML
 MiniSail
 Minimal_SSA
 Minkowskis_Theorem
 Minsky_Machines
 Modal_Logics_for_NTS
 Modular_Assembly_Kit_Security
 Modular_arithmetic_LLL_and_HNF_algorithms
 Monad_Memo_DP
 Monad_Normalisation
 MonoBoolTranAlgebra
 MonoidalCategory
 Monomorphic_Monad
 MuchAdoAboutTwo
 Multi_Party_Computation
 Multirelations
 Myhill-Nerode
 Name_Carrying_Type_Inference
 Nash_Williams
 Nat-Interval-Logic
 Native_Word
 Nested_Multisets_Ordinals
 Network_Security_Policy_Verification
 Neumann_Morgenstern_Utility
 No_FTL_observers
 Nominal2
 Noninterference_CSP
 Noninterference_Concurrent_Composition
 Noninterference_Generic_Unwinding
 Noninterference_Inductive_Unwinding
 Noninterference_Ipurge_Unwinding
 Noninterference_Sequential_Composition
 NormByEval
 Nullstellensatz
 Octonions
 OpSets
 Open_Induction
 Optics
 Optimal_BST
 Orbit_Stabiliser
 Order_Lattice_Props
 Ordered_Resolution_Prover
 Ordinal
 Ordinal_Partitions
 Ordinals_and_Cardinals
 Ordinary_Differential_Equations
 PAC_Checker
 PAL
 PCF
 PLM
 POPLmark-deBruijn
 PSemigroupsConvolution
 Padic_Ints
 Pairing_Heap
 Paraconsistency
 Parity_Game
 Partial_Function_MR
 Partial_Order_Reduction
 Password_Authentication_Protocol
 Pell
 Perfect-Number-Thm
 Perron_Frobenius
 Physical_Quantities
 Pi_Calculus
 Pi_Transcendental
 Planarity_Certificates
 Poincare_Bendixson
 Poincare_Disc
 Polynomial_Factorization
 Polynomial_Interpolation
 Polynomials
 Pop_Refinement
 Posix-Lexing
 Possibilistic_Noninterference
 Power_Sum_Polynomials
 Pratt_Certificate
 Presburger-Automata
 Prim_Dijkstra_Simple
 Prime_Distribution_Elementary
 Prime_Harmonic_Series
 Prime_Number_Theorem
 Priority_Queue_Braun
 Priority_Search_Trees
 Probabilistic_Noninterference
 Probabilistic_Prime_Tests
 Probabilistic_System_Zoo
 Probabilistic_Timed_Automata
 Probabilistic_While
 Program-Conflict-Analysis
 Progress_Tracking
 Projective_Geometry
 Projective_Measurements
 Promela
 Proof_Strategy_Language
 PropResPI
 Propositional_Proof_Systems
 Prpu_Maxflow
 PseudoHoops
 Psi_Calculi
 Ptolemys_Theorem
 Public_Announcement_Logic
 QHLProver
 QR_Decomposition
 Quantales
 Quaternions
 Quick_Sort_Cost
 RIPEMD-160-SPARK
 ROBDD
 RSAPSS
 Ramsey-Infinite
 Random_BSTs
 Random_Graph_Subgraph_Threshold
 Randomised_BSTs
 Randomised_Social_Choice
 Rank_Nullity_Theorem
 Real_Impl
 Real_Power
 Recursion-Addition
 Recursion-Theory-I
 Refine_Imperative_HOL
 Refine_Monadic
 RefinementReactive
 Regex_Equivalence
 Registers
 Regression_Test_Selection
 Regular-Sets
 Regular_Algebras
 Regular_Tree_Relations
 Relation_Algebra
 Relational-Incorrectness-Logic
 Relational_Disjoint_Set_Forests
 Relational_Forests
 Relational_Method
 Relational_Minimum_Spanning_Trees
 Relational_Paths
 Rep_Fin_Groups
 Residuated_Lattices
 Resolution_FOL
 Rewriting_Z
 Ribbon_Proofs
 Robbins-Conjecture
 Robinson_Arithmetic
 Root_Balanced_Tree
 Roth_Arithmetic_Progressions
 Routing
 Roy_Floyd_Warshall
 SATSolverVerification
 SC_DOM_Components
 SDS_Impossibility
 SIFPL
 SIFUM_Type_Systems
 SPARCv8
 Safe_Distance
 Safe_OCL
 Saturation_Framework
 Saturation_Framework_Extensions
 Schutz_Spacetime
 Secondary_Sylow
 Security_Protocol_Refinement
 Selection_Heap_Sort
 SenSocialChoice
 Separata
 Separation_Algebra
 Separation_Logic_Imperative_HOL
 SequentInvertibility
 Shadow_DOM
 Shadow_SC_DOM
 Shivers-CFA
 ShortestPath
 Show
 Sigma_Commit_Crypto
 Signature_Groebner
 Simpl
 Simple_Firewall
 Simplex
 Simplicial_complexes_and_boolean_functions
 SimplifiedOntologicalArgument
 Skew_Heap
 Skip_Lists
 Slicing
 Sliding_Window_Algorithm
 Smith_Normal_Form
 Smooth_Manifolds
 Sort_Encodings
 Source_Coding_Theorem
 SpecCheck
 Special_Function_Bounds
 Splay_Tree
 Sqrt_Babylonian
 Stable_Matching
 Statecharts
 Stateful_Protocol_Composition_and_Typing
 Stellar_Quorums
 Stern_Brocot
 Stewart_Apollonius
 Stirling_Formula
 Stochastic_Matrices
 Stone_Algebras
 Stone_Kleene_Relation_Algebras
 Stone_Relation_Algebras
 Store_Buffer_Reduction
 Stream-Fusion
 Stream_Fusion_Code
 Strong_Security
 Sturm_Sequences
 Sturm_Tarski
 Stuttering_Equivalence
 Subresultants
 Subset_Boolean_Algebras
 SumSquares
 Sunflowers
 SuperCalc
 Surprise_Paradox
 Symmetric_Polynomials
 Syntax_Independent_Logic
 Szemeredi_Regularity
 Szpilrajn
 TESL_Language
 TLA
 Tail_Recursive_Functions
 Tarskis_Geometry
 Taylor_Models
 Three_Circles
 Timed_Automata
 Topological_Semantics
 Topology
 TortoiseHare
 Transcendence_Series_Hancl_Rucki
 Transformer_Semantics
 Transition_Systems_and_Automata
 Transitive-Closure
 Transitive-Closure-II
 Treaps
 Tree-Automata
 Tree_Decomposition
 Triangle
 Trie
 Twelvefold_Way
 Tycon
 Types_Tableaus_and_Goedels_God
 Types_To_Sets_Extension
 UPF
 UPF_Firewall
 UTP
 Universal_Turing_Machine
 UpDown_Scheme
 Valuation
 Van_Emde_Boas_Trees
 Van_der_Waerden
 VectorSpace
 VeriComp
 Verified-Prover
 Verified_SAT_Based_AI_Planning
 VerifyThis2018
 VerifyThis2019
 Vickrey_Clarke_Groves
 Virtual_Substitution
 VolpanoSmith
 WHATandWHERE_Security
 WOOT_Strong_Eventual_Consistency
 WebAssembly
 Weight_Balanced_Trees
 Weighted_Path_Order
 Well_Quasi_Orders
 Winding_Number_Eval
 Word_Lib
 WorkerWrapper
 X86_Semantics
 XML
 ZFC_in_HOL
 Zeta_3_Irrational
 Zeta_Function
 pGCL