diff --git a/crates/argmin-testfunctions/proptest-regressions/picheny.txt b/crates/argmin-testfunctions/proptest-regressions/picheny.txt
new file mode 100644
index 000000000..03e681b0d
--- /dev/null
+++ b/crates/argmin-testfunctions/proptest-regressions/picheny.txt
@@ -0,0 +1,11 @@
+# Seeds for failure cases proptest has generated in the past. It is
+# automatically read and these particular cases re-run before any
+# novel cases are generated.
+#
+# It is recommended to check this file in to source control so that
+# everyone who runs the test benefits from these saved cases.
+cc 5aa3ba0e4323a8dc0dd12ac1dbc1c38b376037a0c0bb90ccc29c39d7eaf085b0 # shrinks to a = 0.0, b = 0.0
+cc eee01c0df03e5b5723d2865397ab5db4834e8bd8e0e97a6bc9cc8788d26edbdf # shrinks to a = 0.8734667101799871, b = 0.5043381978510337
+cc 99447ad3249d3b6b719ee686fbd85095e8488bfd0e433c56571eadf8539bff21 # shrinks to a = 0.5555316208817483, b = 0.15793555431430475
+cc 3b77493a4ce1cd91948096939c65ca5859b1b7f8b2714eb94bd442561ca373fe
+cc d8e9c4c61c5dd81da980ed3275484ddf7f558553f059e0513d90ffe4dfe5dddb # shrinks to a = 0.9727914197428091, b = 0.21845982572353503
diff --git a/crates/argmin-testfunctions/src/picheny.rs b/crates/argmin-testfunctions/src/picheny.rs
index f817ccc43..1a2e78630 100644
--- a/crates/argmin-testfunctions/src/picheny.rs
+++ b/crates/argmin-testfunctions/src/picheny.rs
@@ -12,7 +12,7 @@
 //! Defined as
 //!
 //! `f(x_1, x_2) = (1/2.427) * log([1 + (\bar{x}_1 + \bar{x}_2 + 1)^2 * (19 - 14*\bar{x}_2 +
-//!                3*\bar{x}_1^2 - 14*\bar{x}_2 6*\bar{x}_1*\bar{x}_2 + 3*\bar{x}_2^2)]
+//!                3*\bar{x}_1^2 - 14*\bar{x}_2 + 6*\bar{x}_1*\bar{x}_2 + 3*\bar{x}_2^2)]
 //!                * [30 + (2*\bar{x}_1 - 3*\bar{x}_2)^2(18 - 32 * \bar{x}_1 + 12* \bar{x}_1^2 +
 //!                48 * \bar{x}_2 - 36 * \bar{x}_1 * \bar{x}_2 + 27 * \bar{x}_2^2) ] - 8.693)`
 //!
@@ -20,6 +20,13 @@
 //!
 //! The global minimum is at `f(x_1, x_2) = f(0.5, 0.25) = 3.3851993182036826`.
 
+//  (1/2.427) * log([1 + ((4*x_1-2) + (4*x_2-2) + 1)^2 * (19 - 14*(4*x_1-2) +   3*(4*x_2-2)^2 - 14*(4*x_2-2) + 6*(4*x_1-2)*(4*x_2-2) + 3*(4*x_2-2)^2)]    * [30 + (2*(4*x_1-2) - 3*(4*x_2-2))^2*(18 - 32 * (4*x_1-2) + 12* (4*x_1-2)^2 +   48 * (4*x_2-2) - 36 * (4*x_1-2) * (4*x_2-2) + 27 * (4*x_2-2)^2) ] - 8.693)
+//  (1/2.427) * log([1 + (x_1 + x_2 + 1)^2 * (19 - 14*x_1 +   3*x_2^2 - 14*x_2 + 6*x_1*x_2 + 3*x_2^2)]    * [30 + (2*x_1 - 3*x_2)^2*(18 - 32 * x_1 + 12* x_1^2 +   48 * x_2 - 36 * x_1 * x_2 + 27 * x_2^2) ] - 8.693)
+//  (1/2.427) * log([1 + (x_1 + x_2 + 1)^2 * (19 - 14*x_1 +   3*x_1^2 - 14*x_2 + 6*x_1*x_2 + 3*x_2^2)]    * [30 + (2*x_1 - 3*x_2)^2*(18 - 32 * x_1 + 12* x_1^2 +   48 * x_2 - 36 * x_1 * x_2 + 27 * x_2^2) ]) - 8.693
+//  (1/2.427) * log_10([1 + (x_1 + x_2 + 1)^2 * (19 - 14*x_1 +   3*x_1^2 - 14*x_2 + 6*x_1*x_2 + 3*x_2^2)]    * [30 + (2*x_1 - 3*x_2)^2*(18 - 32 * x_1 + 12* x_1^2 +   48 * x_2 - 36 * x_1 * x_2 + 27 * x_2^2) ]) - 8.693
+//  (1/2.427) * log_10([1 + ((4*x_1-2) + (4*x_2-2) + 1)^2 * (19 - 14*(4*x_1-2) +   3*(4*x_1-2)^2 - 14*(4*x_2-2) + 6*(4*x_1-2)*(4*x_2-2) + 3*(4*x_2-2)^2)]    * [30 + (2*(4*x_1-2) - 3*(4*x_2-2))^2*(18 - 32 * (4*x_1-2) + 12* (4*x_1-2)^2 +   48 * (4*x_2-2) - 36 * (4*x_1-2) * (4*x_2-2) + 27 * (4*x_2-2)^2) ]) - 8.693
+//  (1/2.427) * (log_10([1 + ((4*x_1-2) + (4*x_2-2) + 1)^2 * (19 - 14*(4*x_1-2) +   3*(4*x_1-2)^2 - 14*(4*x_2-2) + 6*(4*x_1-2)*(4*x_2-2) + 3*(4*x_2-2)^2)]    * [30 + (2*(4*x_1-2) - 3*(4*x_2-2))^2*(18 - 32 * (4*x_1-2) + 12* (4*x_1-2)^2 +   48 * (4*x_2-2) - 36 * (4*x_1-2) * (4*x_2-2) + 27 * (4*x_2-2)^2) ]) - 8.693)
+
 use num::{Float, FromPrimitive};
 
 /// Picheny test function
@@ -68,11 +75,760 @@ where
         .log10()
             - T::from_f64(8.693).unwrap())
 }
+/// Derivative of Picheny test function
+pub fn picheny_derivative<T>(param: &[T; 2]) -> [T; 2]
+where
+    T: Float + FromPrimitive,
+{
+    let n1 = T::from_f64(1.0).unwrap();
+    let n2 = T::from_f64(2.0).unwrap();
+    let n3 = T::from_f64(3.0).unwrap();
+    let n4 = T::from_f64(4.0).unwrap();
+    let n6 = T::from_f64(6.0).unwrap();
+    let n8 = T::from_f64(8.0).unwrap();
+    let n10 = T::from_f64(10.0).unwrap();
+    let n12 = T::from_f64(12.0).unwrap();
+    let n14 = T::from_f64(14.0).unwrap();
+    let n16 = T::from_f64(16.0).unwrap();
+    let n18 = T::from_f64(18.0).unwrap();
+    let n19 = T::from_f64(19.0).unwrap();
+    let n24 = T::from_f64(24.0).unwrap();
+    let n27 = T::from_f64(27.0).unwrap();
+    let n30 = T::from_f64(30.0).unwrap();
+    let n32 = T::from_f64(32.0).unwrap();
+    let n36 = T::from_f64(36.0).unwrap();
+    let n48 = T::from_f64(48.0).unwrap();
+    let n56 = T::from_f64(56.0).unwrap();
+    let n96 = T::from_f64(96.0).unwrap();
+    let n128 = T::from_f64(128.0).unwrap();
+    let n144 = T::from_f64(144.0).unwrap();
+    let n192 = T::from_f64(192.0).unwrap();
+    let n216 = T::from_f64(216.0).unwrap();
+    let factor = T::from_f64(1.0 / 2.427).unwrap();
+
+    let [x1, x2] = *param;
+
+    let a = (factor
+        * ((n8
+            * (n4 * x1 + n4 * x2 - n3)
+            * (n3 * (n4 * x1 - n2).powi(2) + n6 * (n4 * x2 - n2) * (n4 * x1 - n2)
+                - n14 * (n4 * x1 - n2)
+                + n3 * (n4 * x2 - n2).powi(2)
+                - n14 * (n4 * x2 - n2)
+                + n19)
+            + (n4 * x1 + n4 * x2 - n3).powi(2)
+                * (n24 * (n4 * x1 - n2) + n24 * (n4 * x2 - n2) - n56))
+            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
+                * (n12 * (n4 * x1 - n2).powi(2)
+                    - n36 * (n4 * x2 - n2) * (n4 * x1 - n2)
+                    - n32 * (n4 * x1 - n2)
+                    + n27 * (n4 * x2 - n2).powi(2)
+                    + n48 * (n4 * x2 - n2)
+                    + n18)
+                + n30)
+            + ((n4 * x1 + n4 * x2 - n3).powi(2)
+                * (n3 * (n4 * x1 - n2).powi(2) + n6 * (n4 * x2 - n2) * (n4 * x1 - n2)
+                    - n14 * (n4 * x1 - n2)
+                    + n3 * (n4 * x2 - n2).powi(2)
+                    - n14 * (n4 * x2 - n2)
+                    + n19)
+                + n1)
+                * (n16
+                    * (n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2))
+                    * (n12 * (n4 * x1 - n2).powi(2)
+                        - n36 * (n4 * x2 - n2) * (n4 * x1 - n2)
+                        - n32 * (n4 * x1 - n2)
+                        + n27 * (n4 * x2 - n2).powi(2)
+                        + n48 * (n4 * x2 - n2)
+                        + n18)
+                    + (n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
+                        * (n96 * (n4 * x1 - n2) - n144 * (n4 * x2 - n2) - n128))))
+        / (n10.ln()
+            * ((n4 * x1 + n4 * x2 - n3).powi(2)
+                * (n3 * (n4 * x1 - n2).powi(2) + n6 * (n4 * x2 - n2) * (n4 * x1 - n2)
+                    - n14 * (n4 * x1 - n2)
+                    + n3 * (n4 * x2 - n2).powi(2)
+                    - n14 * (n4 * x2 - n2)
+                    + n19)
+                + n1)
+            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
+                * (n12 * (n4 * x1 - n2).powi(2)
+                    - n36 * (n4 * x2 - n2) * (n4 * x1 - n2)
+                    - n32 * (n4 * x1 - n2)
+                    + n27 * (n4 * x2 - n2).powi(2)
+                    + n48 * (n4 * x2 - n2)
+                    + n18)
+                + n30));
+
+    let b = (factor
+        * ((n8
+            * (n4 * x2 + n4 * x1 - n3)
+            * (n3 * (n4 * x2 - n2).powi(2) + n6 * (n4 * x1 - n2) * (n4 * x2 - n2)
+                - n14 * (n4 * x2 - n2)
+                + n3 * (n4 * x1 - n2).powi(2)
+                - n14 * (n4 * x1 - n2)
+                + n19)
+            + (n4 * x2 + n4 * x1 - n3).powi(2)
+                * (n24 * (n4 * x2 - n2) + n24 * (n4 * x1 - n2) - n56))
+            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
+                * (n27 * (n4 * x2 - n2).powi(2) - n36 * (n4 * x1 - n2) * (n4 * x2 - n2)
+                    + n48 * (n4 * x2 - n2)
+                    + n12 * (n4 * x1 - n2).powi(2)
+                    - n32 * (n4 * x1 - n2)
+                    + n18)
+                + n30)
+            + ((n4 * x2 + n4 * x1 - n3).powi(2)
+                * (n3 * (n4 * x2 - n2).powi(2) + n6 * (n4 * x1 - n2) * (n4 * x2 - n2)
+                    - n14 * (n4 * x2 - n2)
+                    + n3 * (n4 * x1 - n2).powi(2)
+                    - n14 * (n4 * x1 - n2)
+                    + n19)
+                + n1)
+                * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
+                    * (n216 * (n4 * x2 - n2) - n144 * (n4 * x1 - n2) + n192)
+                    - n24
+                        * (n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2))
+                        * (n27 * (n4 * x2 - n2).powi(2) - n36 * (n4 * x1 - n2) * (n4 * x2 - n2)
+                            + n48 * (n4 * x2 - n2)
+                            + n12 * (n4 * x1 - n2).powi(2)
+                            - n32 * (n4 * x1 - n2)
+                            + n18))))
+        / (n10.ln()
+            * ((n4 * x2 + n4 * x1 - n3).powi(2)
+                * (n3 * (n4 * x2 - n2).powi(2) + n6 * (n4 * x1 - n2) * (n4 * x2 - n2)
+                    - n14 * (n4 * x2 - n2)
+                    + n3 * (n4 * x1 - n2).powi(2)
+                    - n14 * (n4 * x1 - n2)
+                    + n19)
+                + n1)
+            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
+                * (n27 * (n4 * x2 - n2).powi(2) - n36 * (n4 * x1 - n2) * (n4 * x2 - n2)
+                    + n48 * (n4 * x2 - n2)
+                    + n12 * (n4 * x1 - n2).powi(2)
+                    - n32 * (n4 * x1 - n2)
+                    + n18)
+                + n30));
+
+    [a, b]
+}
+
+/// Hessian of Picheny test function
+pub fn picheny_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
+where
+    T: Float + FromPrimitive,
+{
+    // oh dear...
+    let n2 = T::from_f64(2.0).unwrap();
+    let n3 = T::from_f64(3.0).unwrap();
+    let n4 = T::from_f64(4.0).unwrap();
+    let n5 = T::from_f64(5.0).unwrap();
+    let n6 = T::from_f64(6.0).unwrap();
+    let n10 = T::from_f64(10.0).unwrap();
+    let n11 = T::from_f64(11.0).unwrap();
+    let n96 = T::from_f64(96.0).unwrap();
+    let n144 = T::from_f64(144.0).unwrap();
+    let n192 = T::from_f64(192.0).unwrap();
+    let n277 = T::from_f64(277.0).unwrap();
+    let n432 = T::from_f64(432.0).unwrap();
+    let n576 = T::from_f64(576.0).unwrap();
+    let n768 = T::from_f64(768.0).unwrap();
+    let n864 = T::from_f64(864.0).unwrap();
+    let n896 = T::from_f64(896.0).unwrap();
+    let n1080 = T::from_f64(1080.0).unwrap();
+    let n1152 = T::from_f64(1152.0).unwrap();
+    let n1167 = T::from_f64(1167.0).unwrap();
+    let n1512 = T::from_f64(1512.0).unwrap();
+    let n1603 = T::from_f64(1603.0).unwrap();
+    let n2048 = T::from_f64(2048.0).unwrap();
+    let n2304 = T::from_f64(2304.0).unwrap();
+    let n2688 = T::from_f64(2688.0).unwrap();
+    let n3024 = T::from_f64(3024.0).unwrap();
+    let n5376 = T::from_f64(5376.0).unwrap();
+    let n5594 = T::from_f64(5594.0).unwrap();
+    let n5874 = T::from_f64(5874.0).unwrap();
+    let n6144 = T::from_f64(6144.0).unwrap();
+    let n6912 = T::from_f64(6912.0).unwrap();
+    let n9216 = T::from_f64(9216.0).unwrap();
+    let n13824 = T::from_f64(13824.0).unwrap();
+    let n20736 = T::from_f64(20736.0).unwrap();
+    let n21546 = T::from_f64(21546.0).unwrap();
+    let n24576 = T::from_f64(24576.0).unwrap();
+    let n27648 = T::from_f64(27648.0).unwrap();
+    let n31104 = T::from_f64(31104.0).unwrap();
+    let n36864 = T::from_f64(36864.0).unwrap();
+    let n77456 = T::from_f64(77456.0).unwrap();
+    let n82944 = T::from_f64(82944.0).unwrap();
+    let n85088 = T::from_f64(85088.0).unwrap();
+    let n98304 = T::from_f64(98304.0).unwrap();
+    let n118688 = T::from_f64(118688.0).unwrap();
+    let n124416 = T::from_f64(124416.0).unwrap();
+    let n127904 = T::from_f64(127904.0).unwrap();
+    let n155024 = T::from_f64(155024.0).unwrap();
+    let n163936 = T::from_f64(163936.0).unwrap();
+    let n165888 = T::from_f64(165888.0).unwrap();
+    let n248832 = T::from_f64(248832.0).unwrap();
+    let n255264 = T::from_f64(255264.0).unwrap();
+    let n255808 = T::from_f64(255808.0).unwrap();
+    let n327872 = T::from_f64(327872.0).unwrap();
+    let n331776 = T::from_f64(331776.0).unwrap();
+    let n393216 = T::from_f64(393216.0).unwrap();
+    let n442368 = T::from_f64(442368.0).unwrap();
+    let n497664 = T::from_f64(497664.0).unwrap();
+    let n516096 = T::from_f64(516096.0).unwrap();
+    let n679936 = T::from_f64(679936.0).unwrap();
+    let n688128 = T::from_f64(688128.0).unwrap();
+    let n705280 = T::from_f64(705280.0).unwrap();
+    let n995328 = T::from_f64(995328.0).unwrap();
+    let n1057536 = T::from_f64(1057536.0).unwrap();
+    let n1057920 = T::from_f64(1057920.0).unwrap();
+    let n1136640 = T::from_f64(1136640.0).unwrap();
+    let n1253376 = T::from_f64(1253376.0).unwrap();
+    let n1290240 = T::from_f64(1290240.0).unwrap();
+    let n1298688 = T::from_f64(1298688.0).unwrap();
+    let n1327104 = T::from_f64(1327104.0).unwrap();
+    let n1435136 = T::from_f64(1435136.0).unwrap();
+    let n1441280 = T::from_f64(1441280.0).unwrap();
+    let n1490944 = T::from_f64(1490944.0).unwrap();
+    let n1515520 = T::from_f64(1515520.0).unwrap();
+    let n1556736 = T::from_f64(1556736.0).unwrap();
+    let n1781760 = T::from_f64(1781760.0).unwrap();
+    let n1854464 = T::from_f64(1854464.0).unwrap();
+    let n1880064 = T::from_f64(1880064.0).unwrap();
+    let n1990656 = T::from_f64(1990656.0).unwrap();
+    let n2088960 = T::from_f64(2088960.0).unwrap();
+    let n2161920 = T::from_f64(2161920.0).unwrap();
+    let n2322432 = T::from_f64(2322432.0).unwrap();
+    let n2363520 = T::from_f64(2363520.0).unwrap();
+    let n2520320 = T::from_f64(2520320.0).unwrap();
+    let n2550112 = T::from_f64(2550112.0).unwrap();
+    let n2580480 = T::from_f64(2580480.0).unwrap();
+    let n2654208 = T::from_f64(2654208.0).unwrap();
+    let n2752512 = T::from_f64(2752512.0).unwrap();
+    let n2985984 = T::from_f64(2985984.0).unwrap();
+    let n3113472 = T::from_f64(3113472.0).unwrap();
+    let n3133440 = T::from_f64(3133440.0).unwrap();
+    let n3896064 = T::from_f64(3896064.0).unwrap();
+    let n4079616 = T::from_f64(4079616.0).unwrap();
+    let n4231680 = T::from_f64(4231680.0).unwrap();
+    let n4323840 = T::from_f64(4323840.0).unwrap();
+    let n4523520 = T::from_f64(4523520.0).unwrap();
+    let n4546560 = T::from_f64(4546560.0).unwrap();
+    let n5040640 = T::from_f64(5040640.0).unwrap();
+    let n5287680 = T::from_f64(5287680.0).unwrap();
+    let n5492736 = T::from_f64(5492736.0).unwrap();
+    let n5971968 = T::from_f64(5971968.0).unwrap();
+    let n6266880 = T::from_f64(6266880.0).unwrap();
+    let n6785280 = T::from_f64(6785280.0).unwrap();
+    let n6850560 = T::from_f64(6850560.0).unwrap();
+    let n7127040 = T::from_f64(7127040.0).unwrap();
+    let n7175680 = T::from_f64(7175680.0).unwrap();
+    let n7399680 = T::from_f64(7399680.0).unwrap();
+    let n7728640 = T::from_f64(7728640.0).unwrap();
+    let n8515584 = T::from_f64(8515584.0).unwrap();
+    let n9089280 = T::from_f64(9089280.0).unwrap();
+    let n9134080 = T::from_f64(9134080.0).unwrap();
+    let n9400320 = T::from_f64(9400320.0).unwrap();
+    let n9454080 = T::from_f64(9454080.0).unwrap();
+    let n10081280 = T::from_f64(10081280.0).unwrap();
+    let n13284864 = T::from_f64(13284864.0).unwrap();
+    let n13570560 = T::from_f64(13570560.0).unwrap();
+    let n13731840 = T::from_f64(13731840.0).unwrap();
+    let n13934592 = T::from_f64(13934592.0).unwrap();
+    let n14799360 = T::from_f64(14799360.0).unwrap();
+    let n23185920 = T::from_f64(23185920.0).unwrap();
+    let n27463680 = T::from_f64(27463680.0).unwrap();
+    let n29598720 = T::from_f64(29598720.0).unwrap();
+    let n27402240 = T::from_f64(27402240.0).unwrap();
+
+    let factor = T::from_f64(9.88875154511743).unwrap();
+
+    let [x_1, x_2] = *param;
+
+    let offdiag = -(factor
+        * (n768 * x_2.powi(3) + x_1 * (n2304 * x_2.powi(2) - n5376 * x_2 + n3024)
+            - n2688 * x_2.powi(2)
+            + x_1.powi(2) * (n2304 * x_2 - n2688)
+            + n3024 * x_2
+            + n768 * x_1.powi(3)
+            - n1080)
+        * (n331776 * x_2.powi(7) - n497664 * x_2.powi(6)
+            + x_1
+                * (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
+                    + n1441280 * x_2.powi(3)
+                    + n1556736 * x_2.powi(2)
+                    - n255808 * x_2
+                    + n5594)
+            - n1057536 * x_2.powi(5)
+            + x_1.powi(2)
+                * (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
+                    - n6785280 * x_2.powi(2)
+                    - n255264 * x_2
+                    + n77456)
+            + x_1.powi(3)
+                * (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
+                    + n5040640 * x_2.powi(2)
+                    + n4231680 * x_2
+                    - n118688)
+            + n2363520 * x_2.powi(4)
+            + x_1.powi(4)
+                * (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2) - n7175680 * x_2 - n705280)
+            - n1298688 * x_2.powi(3)
+            + n163936 * x_2.powi(2)
+            + x_1.powi(5) * (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464)
+            + n5874 * x_2
+            + x_1.powi(6) * (-n688128 * x_2 - n1490944)
+            + n393216 * x_1.powi(7)
+            - n1603))
+        / (n10.ln()
+            * (n192 * x_2.powi(4)
+                + x_1 * (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080)
+                - n896 * x_2.powi(3)
+                + x_1.powi(2) * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512)
+                + n1512 * x_2.powi(2)
+                + x_1.powi(3) * (n768 * x_2 - n896)
+                - n1080 * x_2
+                + n192 * x_1.powi(4)
+                + n277)
+                .powi(2)
+            * (n31104 * x_2.powi(4) - n6912 * x_2.powi(3)
+                + x_1 * (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96)
+                + x_1.powi(2) * (n82944 * x_2.powi(2) - n9216 * x_2 - n192)
+                - n432 * x_2.powi(2)
+                + n144 * x_2
+                + x_1.powi(3) * (n2048 - n36864 * x_2)
+                + n6144 * x_1.powi(4)
+                + n11))
+        - (factor
+            * (n124416 * x_2.powi(3) - n20736 * x_2.powi(2)
+                + x_1 * (-n248832 * x_2.powi(2) + n27648 * x_2 + n576)
+                + x_1.powi(2) * (n165888 * x_2 - n9216)
+                - n864 * x_2
+                - n36864 * x_1.powi(3)
+                + n144)
+            * (n331776 * x_2.powi(7) - n497664 * x_2.powi(6)
+                + x_1
+                    * (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
+                        + n1441280 * x_2.powi(3)
+                        + n1556736 * x_2.powi(2)
+                        - n255808 * x_2
+                        + n5594)
+                - n1057536 * x_2.powi(5)
+                + x_1.powi(2)
+                    * (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
+                        - n6785280 * x_2.powi(2)
+                        - n255264 * x_2
+                        + n77456)
+                + x_1.powi(3)
+                    * (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
+                        + n5040640 * x_2.powi(2)
+                        + n4231680 * x_2
+                        - n118688)
+                + n2363520 * x_2.powi(4)
+                + x_1.powi(4)
+                    * (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2)
+                        - n7175680 * x_2
+                        - n705280)
+                - n1298688 * x_2.powi(3)
+                + n163936 * x_2.powi(2)
+                + x_1.powi(5) * (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464)
+                + n5874 * x_2
+                + x_1.powi(6) * (-n688128 * x_2 - n1490944)
+                + n393216 * x_1.powi(7)
+                - n1603))
+            / (n10.ln()
+                * (n192 * x_2.powi(4)
+                    + x_1 * (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080)
+                    - n896 * x_2.powi(3)
+                    + x_1.powi(2) * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512)
+                    + n1512 * x_2.powi(2)
+                    + x_1.powi(3) * (n768 * x_2 - n896)
+                    - n1080 * x_2
+                    + n192 * x_1.powi(4)
+                    + n277)
+                * (n31104 * x_2.powi(4) - n6912 * x_2.powi(3)
+                    + x_1 * (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96)
+                    + x_1.powi(2) * (n82944 * x_2.powi(2) - n9216 * x_2 - n192)
+                    - n432 * x_2.powi(2)
+                    + n144 * x_2
+                    + x_1.powi(3) * (n2048 - n36864 * x_2)
+                    + n6144 * x_1.powi(4)
+                    + n11)
+                    .powi(2))
+        + (factor
+            * (n2322432 * x_2.powi(6) - n2985984 * x_2.powi(5)
+                + x_1
+                    * (-n5971968 * x_2.powi(5) + n27463680 * x_2.powi(4)
+                        - n29598720 * x_2.powi(3)
+                        + n4323840 * x_2.powi(2)
+                        + n3113472 * x_2
+                        - n255808)
+                - n5287680 * x_2.powi(4)
+                + x_1.powi(2)
+                    * (-n9400320 * x_2.powi(4)
+                        + n4546560 * x_2.powi(3)
+                        + n23185920 * x_2.powi(2)
+                        - n13570560 * x_2
+                        - n255264)
+                + x_1.powi(3)
+                    * (n7127040 * x_2.powi(3) - n27402240 * x_2.powi(2)
+                        + n10081280 * x_2
+                        + n4231680)
+                + n9454080 * x_2.powi(3)
+                + x_1.powi(4) * (n6266880 * x_2.powi(2) + n2580480 * x_2 - n7175680)
+                - n3896064 * x_2.powi(2)
+                + n327872 * x_2
+                + x_1.powi(5) * (n4079616 - n2654208 * x_2)
+                - n688128 * x_1.powi(6)
+                + n5874))
+            / (n10.ln()
+                * (n192 * x_2.powi(4)
+                    + x_1 * (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080)
+                    - n896 * x_2.powi(3)
+                    + x_1.powi(2) * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512)
+                    + n1512 * x_2.powi(2)
+                    + x_1.powi(3) * (n768 * x_2 - n896)
+                    - n1080 * x_2
+                    + n192 * x_1.powi(4)
+                    + n277)
+                * (n31104 * x_2.powi(4) - n6912 * x_2.powi(3)
+                    + x_1 * (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96)
+                    + x_1.powi(2) * (n82944 * x_2.powi(2) - n9216 * x_2 - n192)
+                    - n432 * x_2.powi(2)
+                    + n144 * x_2
+                    + x_1.powi(3) * (n2048 - n36864 * x_2)
+                    + n6144 * x_1.powi(4)
+                    + n11));
+
+    let b = -(factor
+        * (n768 * x_2.powi(3)
+            + n3 * (n768 * x_1 - n896) * x_2.powi(2)
+            + n2 * (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2
+            + n768 * x_1.powi(3)
+            - n2688 * x_1.powi(2)
+            + n3024 * x_1
+            - n1080)
+        * (n1990656 * x_2.powi(7)
+            + (n2322432 * x_1 - n8515584) * x_2.powi(6)
+            + (-n2985984 * x_1.powi(2) - n2985984 * x_1 + n13284864) * x_2.powi(5)
+            + (-n3133440 * x_1.powi(3) + n13731840 * x_1.powi(2) - n5287680 * x_1 - n9089280)
+                * x_2.powi(4)
+            + (n1781760 * x_1.powi(4) + n1515520 * x_1.powi(3) - n14799360 * x_1.powi(2)
+                + n9454080 * x_1
+                + n2550112)
+                * x_2.powi(3)
+            + (n1253376 * x_1.powi(5) - n6850560 * x_1.powi(4)
+                + n7728640 * x_1.powi(3)
+                + n2161920 * x_1.powi(2)
+                - n3896064 * x_1
+                - n155024)
+                * x_2.powi(2)
+            + (-n442368 * x_1.powi(6) + n516096 * x_1.powi(5) + n2520320 * x_1.powi(4)
+                - n4523520 * x_1.powi(3)
+                + n1556736 * x_1.powi(2)
+                + n327872 * x_1
+                - n21546)
+                * x_2
+            - n98304 * x_1.powi(7)
+            + n679936 * x_1.powi(6)
+            - n1435136 * x_1.powi(5)
+            + n1057920 * x_1.powi(4)
+            - n85088 * x_1.powi(3)
+            - n127904 * x_1.powi(2)
+            + n5874 * x_1
+            + n1167))
+        / (n10.ln()
+            * (n192 * x_2.powi(4)
+                + (n768 * x_1 - n896) * x_2.powi(3)
+                + (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2.powi(2)
+                + (n768 * x_1.powi(3) - n2688 * x_1.powi(2) + n3024 * x_1 - n1080) * x_2
+                + n192 * x_1.powi(4)
+                - n896 * x_1.powi(3)
+                + n1512 * x_1.powi(2)
+                - n1080 * x_1
+                + n277)
+                .powi(2)
+            * (n31104 * x_2.powi(4)
+                + (-n82944 * x_1 - n6912) * x_2.powi(3)
+                + (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2.powi(2)
+                + (-n36864 * x_1.powi(3) - n9216 * x_1.powi(2) + n576 * x_1 + n144) * x_2
+                + n6144 * x_1.powi(4)
+                + n2048 * x_1.powi(3)
+                - n192 * x_1.powi(2)
+                - n96 * x_1
+                + n11))
+        - (factor
+            * (n124416 * x_2.powi(3)
+                + n3 * (-n82944 * x_1 - n6912) * x_2.powi(2)
+                + n2 * (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2
+                - n36864 * x_1.powi(3)
+                - n9216 * x_1.powi(2)
+                + n576 * x_1
+                + n144)
+            * (n1990656 * x_2.powi(7)
+                + (n2322432 * x_1 - n8515584) * x_2.powi(6)
+                + (-n2985984 * x_1.powi(2) - n2985984 * x_1 + n13284864) * x_2.powi(5)
+                + (-n3133440 * x_1.powi(3) + n13731840 * x_1.powi(2) - n5287680 * x_1 - n9089280)
+                    * x_2.powi(4)
+                + (n1781760 * x_1.powi(4) + n1515520 * x_1.powi(3) - n14799360 * x_1.powi(2)
+                    + n9454080 * x_1
+                    + n2550112)
+                    * x_2.powi(3)
+                + (n1253376 * x_1.powi(5) - n6850560 * x_1.powi(4)
+                    + n7728640 * x_1.powi(3)
+                    + n2161920 * x_1.powi(2)
+                    - n3896064 * x_1
+                    - n155024)
+                    * x_2.powi(2)
+                + (-n442368 * x_1.powi(6) + n516096 * x_1.powi(5) + n2520320 * x_1.powi(4)
+                    - n4523520 * x_1.powi(3)
+                    + n1556736 * x_1.powi(2)
+                    + n327872 * x_1
+                    - n21546)
+                    * x_2
+                - n98304 * x_1.powi(7)
+                + n679936 * x_1.powi(6)
+                - n1435136 * x_1.powi(5)
+                + n1057920 * x_1.powi(4)
+                - n85088 * x_1.powi(3)
+                - n127904 * x_1.powi(2)
+                + n5874 * x_1
+                + n1167))
+            / (n10.ln()
+                * (n192 * x_2.powi(4)
+                    + (n768 * x_1 - n896) * x_2.powi(3)
+                    + (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2.powi(2)
+                    + (n768 * x_1.powi(3) - n2688 * x_1.powi(2) + n3024 * x_1 - n1080) * x_2
+                    + n192 * x_1.powi(4)
+                    - n896 * x_1.powi(3)
+                    + n1512 * x_1.powi(2)
+                    - n1080 * x_1
+                    + n277)
+                * (n31104 * x_2.powi(4)
+                    + (-n82944 * x_1 - n6912) * x_2.powi(3)
+                    + (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2.powi(2)
+                    + (-n36864 * x_1.powi(3) - n9216 * x_1.powi(2) + n576 * x_1 + n144) * x_2
+                    + n6144 * x_1.powi(4)
+                    + n2048 * x_1.powi(3)
+                    - n192 * x_1.powi(2)
+                    - n96 * x_1
+                    + n11)
+                    .powi(2))
+        + (factor
+            * (n13934592 * x_2.powi(6)
+                + n6 * (n2322432 * x_1 - n8515584) * x_2.powi(5)
+                + n5 * (-n2985984 * x_1.powi(2) - n2985984 * x_1 + n13284864) * x_2.powi(4)
+                + n4 * (-n3133440 * x_1.powi(3) + n13731840 * x_1.powi(2)
+                    - n5287680 * x_1
+                    - n9089280)
+                    * x_2.powi(3)
+                + n3 * (n1781760 * x_1.powi(4) + n1515520 * x_1.powi(3)
+                    - n14799360 * x_1.powi(2)
+                    + n9454080 * x_1
+                    + n2550112)
+                    * x_2.powi(2)
+                + n2 * (n1253376 * x_1.powi(5) - n6850560 * x_1.powi(4)
+                    + n7728640 * x_1.powi(3)
+                    + n2161920 * x_1.powi(2)
+                    - n3896064 * x_1
+                    - n155024)
+                    * x_2
+                - n442368 * x_1.powi(6)
+                + n516096 * x_1.powi(5)
+                + n2520320 * x_1.powi(4)
+                - n4523520 * x_1.powi(3)
+                + n1556736 * x_1.powi(2)
+                + n327872 * x_1
+                - n21546))
+            / (n10.ln()
+                * (n192 * x_2.powi(4)
+                    + (n768 * x_1 - n896) * x_2.powi(3)
+                    + (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2.powi(2)
+                    + (n768 * x_1.powi(3) - n2688 * x_1.powi(2) + n3024 * x_1 - n1080) * x_2
+                    + n192 * x_1.powi(4)
+                    - n896 * x_1.powi(3)
+                    + n1512 * x_1.powi(2)
+                    - n1080 * x_1
+                    + n277)
+                * (n31104 * x_2.powi(4)
+                    + (-n82944 * x_1 - n6912) * x_2.powi(3)
+                    + (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2.powi(2)
+                    + (-n36864 * x_1.powi(3) - n9216 * x_1.powi(2) + n576 * x_1 + n144) * x_2
+                    + n6144 * x_1.powi(4)
+                    + n2048 * x_1.powi(3)
+                    - n192 * x_1.powi(2)
+                    - n96 * x_1
+                    + n11));
+
+    let a = -(factor
+        * (n768 * x_1.powi(3)
+            + n3 * (n768 * x_2 - n896) * x_1.powi(2)
+            + n2 * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1
+            + n768 * x_2.powi(3)
+            - n2688 * x_2.powi(2)
+            + n3024 * x_2
+            - n1080)
+        * (n393216 * x_1.powi(7)
+            + (-n688128 * x_2 - n1490944) * x_1.powi(6)
+            + (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464) * x_1.powi(5)
+            + (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2) - n7175680 * x_2 - n705280)
+                * x_1.powi(4)
+            + (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
+                + n5040640 * x_2.powi(2)
+                + n4231680 * x_2
+                - n118688)
+                * x_1.powi(3)
+            + (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
+                - n6785280 * x_2.powi(2)
+                - n255264 * x_2
+                + n77456)
+                * x_1.powi(2)
+            + (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
+                + n1441280 * x_2.powi(3)
+                + n1556736 * x_2.powi(2)
+                - n255808 * x_2
+                + n5594)
+                * x_1
+            + n331776 * x_2.powi(7)
+            - n497664 * x_2.powi(6)
+            - n1057536 * x_2.powi(5)
+            + n2363520 * x_2.powi(4)
+            - n1298688 * x_2.powi(3)
+            + n163936 * x_2.powi(2)
+            + n5874 * x_2
+            - n1603))
+        / (n10.ln()
+            * (n192 * x_1.powi(4)
+                + (n768 * x_2 - n896) * x_1.powi(3)
+                + (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1.powi(2)
+                + (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080) * x_1
+                + n192 * x_2.powi(4)
+                - n896 * x_2.powi(3)
+                + n1512 * x_2.powi(2)
+                - n1080 * x_2
+                + n277)
+                .powi(2)
+            * (n6144 * x_1.powi(4)
+                + (n2048 - n36864 * x_2) * x_1.powi(3)
+                + (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1.powi(2)
+                + (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96) * x_1
+                + n31104 * x_2.powi(4)
+                - n6912 * x_2.powi(3)
+                - n432 * x_2.powi(2)
+                + n144 * x_2
+                + n11))
+        - (factor
+            * (n24576 * x_1.powi(3)
+                + n3 * (n2048 - n36864 * x_2) * x_1.powi(2)
+                + n2 * (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1
+                - n82944 * x_2.powi(3)
+                + n13824 * x_2.powi(2)
+                + n576 * x_2
+                - n96)
+            * (n393216 * x_1.powi(7)
+                + (-n688128 * x_2 - n1490944) * x_1.powi(6)
+                + (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464) * x_1.powi(5)
+                + (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2) - n7175680 * x_2 - n705280)
+                    * x_1.powi(4)
+                + (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
+                    + n5040640 * x_2.powi(2)
+                    + n4231680 * x_2
+                    - n118688)
+                    * x_1.powi(3)
+                + (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
+                    - n6785280 * x_2.powi(2)
+                    - n255264 * x_2
+                    + n77456)
+                    * x_1.powi(2)
+                + (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
+                    + n1441280 * x_2.powi(3)
+                    + n1556736 * x_2.powi(2)
+                    - n255808 * x_2
+                    + n5594)
+                    * x_1
+                + n331776 * x_2.powi(7)
+                - n497664 * x_2.powi(6)
+                - n1057536 * x_2.powi(5)
+                + n2363520 * x_2.powi(4)
+                - n1298688 * x_2.powi(3)
+                + n163936 * x_2.powi(2)
+                + n5874 * x_2
+                - n1603))
+            / (n10.ln()
+                * (n192 * x_1.powi(4)
+                    + (n768 * x_2 - n896) * x_1.powi(3)
+                    + (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1.powi(2)
+                    + (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080) * x_1
+                    + n192 * x_2.powi(4)
+                    - n896 * x_2.powi(3)
+                    + n1512 * x_2.powi(2)
+                    - n1080 * x_2
+                    + n277)
+                * (n6144 * x_1.powi(4)
+                    + (n2048 - n36864 * x_2) * x_1.powi(3)
+                    + (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1.powi(2)
+                    + (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96) * x_1
+                    + n31104 * x_2.powi(4)
+                    - n6912 * x_2.powi(3)
+                    - n432 * x_2.powi(2)
+                    + n144 * x_2
+                    + n11)
+                    .powi(2))
+        + (factor
+            * (n2752512 * x_1.powi(6)
+                + n6 * (-n688128 * x_2 - n1490944) * x_1.powi(5)
+                + n5 * (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464) * x_1.powi(4)
+                + n4 * (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2)
+                    - n7175680 * x_2
+                    - n705280)
+                    * x_1.powi(3)
+                + n3 * (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
+                    + n5040640 * x_2.powi(2)
+                    + n4231680 * x_2
+                    - n118688)
+                    * x_1.powi(2)
+                + n2 * (-n1880064 * x_2.powi(5)
+                    + n1136640 * x_2.powi(4)
+                    + n7728640 * x_2.powi(3)
+                    - n6785280 * x_2.powi(2)
+                    - n255264 * x_2
+                    + n77456)
+                    * x_1
+                - n995328 * x_2.powi(6)
+                + n5492736 * x_2.powi(5)
+                - n7399680 * x_2.powi(4)
+                + n1441280 * x_2.powi(3)
+                + n1556736 * x_2.powi(2)
+                - n255808 * x_2
+                + n5594))
+            / (n10.ln()
+                * (n192 * x_1.powi(4)
+                    + (n768 * x_2 - n896) * x_1.powi(3)
+                    + (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1.powi(2)
+                    + (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080) * x_1
+                    + n192 * x_2.powi(4)
+                    - n896 * x_2.powi(3)
+                    + n1512 * x_2.powi(2)
+                    - n1080 * x_2
+                    + n277)
+                * (n6144 * x_1.powi(4)
+                    + (n2048 - n36864 * x_2) * x_1.powi(3)
+                    + (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1.powi(2)
+                    + (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96) * x_1
+                    + n31104 * x_2.powi(4)
+                    - n6912 * x_2.powi(3)
+                    - n432 * x_2.powi(2)
+                    + n144 * x_2
+                    + n11));
+
+    [[a, offdiag], [offdiag, b]]
+}
 
 #[cfg(test)]
 mod tests {
     use super::*;
     use approx::assert_relative_eq;
+    use finitediff::FiniteDiff;
+    use proptest::prelude::*;
     use std::{f32, f64};
 
     #[test]
@@ -87,5 +843,56 @@ mod tests {
             -3.3851993182036826,
             epsilon = f64::EPSILON
         );
+
+        let deriv = picheny_derivative(&[0.5_f64, 0.25_f64]);
+        // println!("1: {deriv:?}");
+        for i in 0..2 {
+            assert_relative_eq!(deriv[i], 0.0, epsilon = f64::EPSILON);
+        }
+    }
+
+    proptest! {
+        #[test]
+        fn test_picheny_derivative(a in 0.0..1.0, b in 0.0..1.0) {
+            let param = [a, b];
+            let derivative = picheny_derivative(&param);
+            let derivative_fd = Vec::from(param).central_diff(&|x| picheny(&[x[0], x[1]]));
+            // println!("1: {derivative:?} at {a}/{b}");
+            // println!("2: {derivative_fd:?} at {a}/{b}");
+            for i in 0..derivative.len() {
+                assert_relative_eq!(
+                    derivative[i],
+                    derivative_fd[i],
+                    epsilon = f64::EPSILON,
+                    max_relative = 1e-3
+                );
+            }
+        }
+    }
+
+    proptest! {
+        #[test]
+        fn test_picheny_hessian_finitediff(a in 0.0..1.0, b in 0.0..1.0) {
+            let param = [a, b];
+            let hessian = picheny_hessian(&param);
+            let hessian_fd =
+                Vec::from(param).central_hessian(&|x| picheny_derivative(&[x[0], x[1]]).to_vec());
+            let n = hessian.len();
+            // println!("1: {hessian:?} at {a}/{b}");
+            // println!("2: {hessian_fd:?} at {a}/{b}");
+            for i in 0..n {
+                assert_eq!(hessian[i].len(), n);
+                for j in 0..n {
+                    if hessian_fd[i][j].is_finite() {
+                        assert_relative_eq!(
+                            hessian[i][j],
+                            hessian_fd[i][j],
+                            epsilon = f64::EPSILON,
+                            max_relative = 1e-4
+                        );
+                    }
+                }
+            }
+        }
     }
 }