Diagrams.TwoD.Apollonian:initialConfig from diagrams-contrib-1.3.0.5, A

Percentage Accurate: 69.0% → 95.2%
Time: 5.1s
Alternatives: 8
Speedup: 0.6×

Specification

?
\[\begin{array}{l} \\ \frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \end{array} \]
(FPCore (x y z)
 :precision binary64
 (/ (- (+ (* x x) (* y y)) (* z z)) (* y 2.0)))
double code(double x, double y, double z) {
	return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
module fmin_fmax_functions
    implicit none
    private
    public fmax
    public fmin

    interface fmax
        module procedure fmax88
        module procedure fmax44
        module procedure fmax84
        module procedure fmax48
    end interface
    interface fmin
        module procedure fmin88
        module procedure fmin44
        module procedure fmin84
        module procedure fmin48
    end interface
contains
    real(8) function fmax88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(4) function fmax44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(8) function fmax84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmax48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
    end function
    real(8) function fmin88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(4) function fmin44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(8) function fmin84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmin48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
    end function
end module

real(8) function code(x, y, z)
use fmin_fmax_functions
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    real(8), intent (in) :: z
    code = (((x * x) + (y * y)) - (z * z)) / (y * 2.0d0)
end function
public static double code(double x, double y, double z) {
	return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
def code(x, y, z):
	return (((x * x) + (y * y)) - (z * z)) / (y * 2.0)
function code(x, y, z)
	return Float64(Float64(Float64(Float64(x * x) + Float64(y * y)) - Float64(z * z)) / Float64(y * 2.0))
end
function tmp = code(x, y, z)
	tmp = (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
end
code[x_, y_, z_] := N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y * y), $MachinePrecision]), $MachinePrecision] - N[(z * z), $MachinePrecision]), $MachinePrecision] / N[(y * 2.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}

\\
\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2}
\end{array}

Sampling outcomes in binary64 precision:

Local Percentage Accuracy vs ?

The average percentage accuracy by input value. Horizontal axis shows value of an input variable; the variable is choosen in the title. Vertical axis is accuracy; higher is better. Red represent the original program, while blue represents Herbie's suggestion. These can be toggled with buttons below the plot. The line is an average while dots represent individual samples.

Accuracy vs Speed?

Herbie found 8 alternatives:

AlternativeAccuracySpeedup
The accuracy (vertical axis) and speed (horizontal axis) of each alternatives. Up and to the right is better. The red square shows the initial program, and each blue circle shows an alternative.The line shows the best available speed-accuracy tradeoffs.

Initial Program: 69.0% accurate, 1.0× speedup?

\[\begin{array}{l} \\ \frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \end{array} \]
(FPCore (x y z)
 :precision binary64
 (/ (- (+ (* x x) (* y y)) (* z z)) (* y 2.0)))
double code(double x, double y, double z) {
	return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
module fmin_fmax_functions
    implicit none
    private
    public fmax
    public fmin

    interface fmax
        module procedure fmax88
        module procedure fmax44
        module procedure fmax84
        module procedure fmax48
    end interface
    interface fmin
        module procedure fmin88
        module procedure fmin44
        module procedure fmin84
        module procedure fmin48
    end interface
contains
    real(8) function fmax88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(4) function fmax44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(8) function fmax84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmax48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
    end function
    real(8) function fmin88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(4) function fmin44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(8) function fmin84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmin48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
    end function
end module

real(8) function code(x, y, z)
use fmin_fmax_functions
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    real(8), intent (in) :: z
    code = (((x * x) + (y * y)) - (z * z)) / (y * 2.0d0)
end function
public static double code(double x, double y, double z) {
	return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
def code(x, y, z):
	return (((x * x) + (y * y)) - (z * z)) / (y * 2.0)
function code(x, y, z)
	return Float64(Float64(Float64(Float64(x * x) + Float64(y * y)) - Float64(z * z)) / Float64(y * 2.0))
end
function tmp = code(x, y, z)
	tmp = (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
end
code[x_, y_, z_] := N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y * y), $MachinePrecision]), $MachinePrecision] - N[(z * z), $MachinePrecision]), $MachinePrecision] / N[(y * 2.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}

\\
\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2}
\end{array}

Alternative 1: 95.2% accurate, 0.3× speedup?

\[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ \begin{array}{l} t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\ y\_s \cdot \begin{array}{l} \mathbf{if}\;t\_0 \leq 0:\\ \;\;\;\;\left(y\_m - \frac{z\_m}{y\_m} \cdot z\_m\right) \cdot 0.5\\ \mathbf{elif}\;t\_0 \leq \infty:\\ \;\;\;\;\mathsf{fma}\left(\frac{x\_m}{y\_m}, x\_m, y\_m\right) \cdot 0.5\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\frac{0.5}{y\_m \cdot y\_m} \cdot \left(x\_m - z\_m\right), x\_m, 0.5\right) \cdot y\_m\\ \end{array} \end{array} \end{array} \]
x_m = (fabs.f64 x)
z_m = (fabs.f64 z)
y\_m = (fabs.f64 y)
y\_s = (copysign.f64 #s(literal 1 binary64) y)
(FPCore (y_s x_m y_m z_m)
 :precision binary64
 (let* ((t_0 (/ (- (+ (* x_m x_m) (* y_m y_m)) (* z_m z_m)) (* y_m 2.0))))
   (*
    y_s
    (if (<= t_0 0.0)
      (* (- y_m (* (/ z_m y_m) z_m)) 0.5)
      (if (<= t_0 INFINITY)
        (* (fma (/ x_m y_m) x_m y_m) 0.5)
        (* (fma (* (/ 0.5 (* y_m y_m)) (- x_m z_m)) x_m 0.5) y_m))))))
x_m = fabs(x);
z_m = fabs(z);
y\_m = fabs(y);
y\_s = copysign(1.0, y);
double code(double y_s, double x_m, double y_m, double z_m) {
	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
	double tmp;
	if (t_0 <= 0.0) {
		tmp = (y_m - ((z_m / y_m) * z_m)) * 0.5;
	} else if (t_0 <= ((double) INFINITY)) {
		tmp = fma((x_m / y_m), x_m, y_m) * 0.5;
	} else {
		tmp = fma(((0.5 / (y_m * y_m)) * (x_m - z_m)), x_m, 0.5) * y_m;
	}
	return y_s * tmp;
}
x_m = abs(x)
z_m = abs(z)
y\_m = abs(y)
y\_s = copysign(1.0, y)
function code(y_s, x_m, y_m, z_m)
	t_0 = Float64(Float64(Float64(Float64(x_m * x_m) + Float64(y_m * y_m)) - Float64(z_m * z_m)) / Float64(y_m * 2.0))
	tmp = 0.0
	if (t_0 <= 0.0)
		tmp = Float64(Float64(y_m - Float64(Float64(z_m / y_m) * z_m)) * 0.5);
	elseif (t_0 <= Inf)
		tmp = Float64(fma(Float64(x_m / y_m), x_m, y_m) * 0.5);
	else
		tmp = Float64(fma(Float64(Float64(0.5 / Float64(y_m * y_m)) * Float64(x_m - z_m)), x_m, 0.5) * y_m);
	end
	return Float64(y_s * tmp)
end
x_m = N[Abs[x], $MachinePrecision]
z_m = N[Abs[z], $MachinePrecision]
y\_m = N[Abs[y], $MachinePrecision]
y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := Block[{t$95$0 = N[(N[(N[(N[(x$95$m * x$95$m), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z$95$m * z$95$m), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision]}, N[(y$95$s * If[LessEqual[t$95$0, 0.0], N[(N[(y$95$m - N[(N[(z$95$m / y$95$m), $MachinePrecision] * z$95$m), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision], If[LessEqual[t$95$0, Infinity], N[(N[(N[(x$95$m / y$95$m), $MachinePrecision] * x$95$m + y$95$m), $MachinePrecision] * 0.5), $MachinePrecision], N[(N[(N[(N[(0.5 / N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] * N[(x$95$m - z$95$m), $MachinePrecision]), $MachinePrecision] * x$95$m + 0.5), $MachinePrecision] * y$95$m), $MachinePrecision]]]), $MachinePrecision]]
\begin{array}{l}
x_m = \left|x\right|
\\
z_m = \left|z\right|
\\
y\_m = \left|y\right|
\\
y\_s = \mathsf{copysign}\left(1, y\right)

\\
\begin{array}{l}
t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\
y\_s \cdot \begin{array}{l}
\mathbf{if}\;t\_0 \leq 0:\\
\;\;\;\;\left(y\_m - \frac{z\_m}{y\_m} \cdot z\_m\right) \cdot 0.5\\

\mathbf{elif}\;t\_0 \leq \infty:\\
\;\;\;\;\mathsf{fma}\left(\frac{x\_m}{y\_m}, x\_m, y\_m\right) \cdot 0.5\\

\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(\frac{0.5}{y\_m \cdot y\_m} \cdot \left(x\_m - z\_m\right), x\_m, 0.5\right) \cdot y\_m\\


\end{array}
\end{array}
\end{array}
Derivation
  1. Split input into 3 regimes
  2. if (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < 0.0

    1. Initial program 81.1%

      \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
    2. Add Preprocessing
    3. Taylor expanded in x around 0

      \[\leadsto \color{blue}{\frac{1}{2} \cdot \frac{{y}^{2} - {z}^{2}}{y}} \]
    4. Step-by-step derivation
      1. Applied rewrites65.3%

        \[\leadsto \color{blue}{\left(y - \frac{z \cdot z}{y}\right) \cdot 0.5} \]
      2. Step-by-step derivation
        1. Applied rewrites67.6%

          \[\leadsto \left(y - \frac{z}{y} \cdot z\right) \cdot 0.5 \]

        if 0.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < +inf.0

        1. Initial program 76.9%

          \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
        2. Add Preprocessing
        3. Taylor expanded in z around 0

          \[\leadsto \color{blue}{\frac{1}{2} \cdot \frac{{x}^{2} + {y}^{2}}{y}} \]
        4. Step-by-step derivation
          1. Applied rewrites73.2%

            \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{x}{y}, x, y\right) \cdot 0.5} \]

          if +inf.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64)))

          1. Initial program 0.0%

            \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
          2. Add Preprocessing
          3. Taylor expanded in y around inf

            \[\leadsto \color{blue}{y \cdot \left(\frac{1}{2} + \frac{1}{2} \cdot \frac{{x}^{2} - {z}^{2}}{{y}^{2}}\right)} \]
          4. Applied rewrites41.0%

            \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{\frac{0.5}{y}}{y}, \left(x - z\right) \cdot \left(z + x\right), 0.5\right) \cdot y} \]
          5. Step-by-step derivation
            1. Applied rewrites37.9%

              \[\leadsto \mathsf{fma}\left(\frac{0.5}{y \cdot y}, \left(x - z\right) \cdot \left(z + x\right), 0.5\right) \cdot y \]
            2. Taylor expanded in x around inf

              \[\leadsto \mathsf{fma}\left(\frac{\frac{1}{2}}{y \cdot y}, \left(x - z\right) \cdot x, \frac{1}{2}\right) \cdot y \]
            3. Step-by-step derivation
              1. Applied rewrites62.9%

                \[\leadsto \mathsf{fma}\left(\frac{0.5}{y \cdot y}, \left(x - z\right) \cdot x, 0.5\right) \cdot y \]
              2. Step-by-step derivation
                1. Applied rewrites70.9%

                  \[\leadsto \mathsf{fma}\left(\frac{0.5}{y \cdot y} \cdot \left(x - z\right), x, 0.5\right) \cdot y \]
              3. Recombined 3 regimes into one program.
              4. Add Preprocessing

              Alternative 2: 69.7% accurate, 0.3× speedup?

              \[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ \begin{array}{l} t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\ y\_s \cdot \begin{array}{l} \mathbf{if}\;t\_0 \leq 0:\\ \;\;\;\;\left(z\_m \cdot \frac{z\_m}{y\_m}\right) \cdot -0.5\\ \mathbf{elif}\;t\_0 \leq 5 \cdot 10^{+153} \lor \neg \left(t\_0 \leq \infty\right):\\ \;\;\;\;0.5 \cdot y\_m\\ \mathbf{else}:\\ \;\;\;\;\left(x\_m \cdot \frac{x\_m}{y\_m}\right) \cdot 0.5\\ \end{array} \end{array} \end{array} \]
              x_m = (fabs.f64 x)
              z_m = (fabs.f64 z)
              y\_m = (fabs.f64 y)
              y\_s = (copysign.f64 #s(literal 1 binary64) y)
              (FPCore (y_s x_m y_m z_m)
               :precision binary64
               (let* ((t_0 (/ (- (+ (* x_m x_m) (* y_m y_m)) (* z_m z_m)) (* y_m 2.0))))
                 (*
                  y_s
                  (if (<= t_0 0.0)
                    (* (* z_m (/ z_m y_m)) -0.5)
                    (if (or (<= t_0 5e+153) (not (<= t_0 INFINITY)))
                      (* 0.5 y_m)
                      (* (* x_m (/ x_m y_m)) 0.5))))))
              x_m = fabs(x);
              z_m = fabs(z);
              y\_m = fabs(y);
              y\_s = copysign(1.0, y);
              double code(double y_s, double x_m, double y_m, double z_m) {
              	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
              	double tmp;
              	if (t_0 <= 0.0) {
              		tmp = (z_m * (z_m / y_m)) * -0.5;
              	} else if ((t_0 <= 5e+153) || !(t_0 <= ((double) INFINITY))) {
              		tmp = 0.5 * y_m;
              	} else {
              		tmp = (x_m * (x_m / y_m)) * 0.5;
              	}
              	return y_s * tmp;
              }
              
              x_m = Math.abs(x);
              z_m = Math.abs(z);
              y\_m = Math.abs(y);
              y\_s = Math.copySign(1.0, y);
              public static double code(double y_s, double x_m, double y_m, double z_m) {
              	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
              	double tmp;
              	if (t_0 <= 0.0) {
              		tmp = (z_m * (z_m / y_m)) * -0.5;
              	} else if ((t_0 <= 5e+153) || !(t_0 <= Double.POSITIVE_INFINITY)) {
              		tmp = 0.5 * y_m;
              	} else {
              		tmp = (x_m * (x_m / y_m)) * 0.5;
              	}
              	return y_s * tmp;
              }
              
              x_m = math.fabs(x)
              z_m = math.fabs(z)
              y\_m = math.fabs(y)
              y\_s = math.copysign(1.0, y)
              def code(y_s, x_m, y_m, z_m):
              	t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0)
              	tmp = 0
              	if t_0 <= 0.0:
              		tmp = (z_m * (z_m / y_m)) * -0.5
              	elif (t_0 <= 5e+153) or not (t_0 <= math.inf):
              		tmp = 0.5 * y_m
              	else:
              		tmp = (x_m * (x_m / y_m)) * 0.5
              	return y_s * tmp
              
              x_m = abs(x)
              z_m = abs(z)
              y\_m = abs(y)
              y\_s = copysign(1.0, y)
              function code(y_s, x_m, y_m, z_m)
              	t_0 = Float64(Float64(Float64(Float64(x_m * x_m) + Float64(y_m * y_m)) - Float64(z_m * z_m)) / Float64(y_m * 2.0))
              	tmp = 0.0
              	if (t_0 <= 0.0)
              		tmp = Float64(Float64(z_m * Float64(z_m / y_m)) * -0.5);
              	elseif ((t_0 <= 5e+153) || !(t_0 <= Inf))
              		tmp = Float64(0.5 * y_m);
              	else
              		tmp = Float64(Float64(x_m * Float64(x_m / y_m)) * 0.5);
              	end
              	return Float64(y_s * tmp)
              end
              
              x_m = abs(x);
              z_m = abs(z);
              y\_m = abs(y);
              y\_s = sign(y) * abs(1.0);
              function tmp_2 = code(y_s, x_m, y_m, z_m)
              	t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
              	tmp = 0.0;
              	if (t_0 <= 0.0)
              		tmp = (z_m * (z_m / y_m)) * -0.5;
              	elseif ((t_0 <= 5e+153) || ~((t_0 <= Inf)))
              		tmp = 0.5 * y_m;
              	else
              		tmp = (x_m * (x_m / y_m)) * 0.5;
              	end
              	tmp_2 = y_s * tmp;
              end
              
              x_m = N[Abs[x], $MachinePrecision]
              z_m = N[Abs[z], $MachinePrecision]
              y\_m = N[Abs[y], $MachinePrecision]
              y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
              code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := Block[{t$95$0 = N[(N[(N[(N[(x$95$m * x$95$m), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z$95$m * z$95$m), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision]}, N[(y$95$s * If[LessEqual[t$95$0, 0.0], N[(N[(z$95$m * N[(z$95$m / y$95$m), $MachinePrecision]), $MachinePrecision] * -0.5), $MachinePrecision], If[Or[LessEqual[t$95$0, 5e+153], N[Not[LessEqual[t$95$0, Infinity]], $MachinePrecision]], N[(0.5 * y$95$m), $MachinePrecision], N[(N[(x$95$m * N[(x$95$m / y$95$m), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision]]]), $MachinePrecision]]
              
              \begin{array}{l}
              x_m = \left|x\right|
              \\
              z_m = \left|z\right|
              \\
              y\_m = \left|y\right|
              \\
              y\_s = \mathsf{copysign}\left(1, y\right)
              
              \\
              \begin{array}{l}
              t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\
              y\_s \cdot \begin{array}{l}
              \mathbf{if}\;t\_0 \leq 0:\\
              \;\;\;\;\left(z\_m \cdot \frac{z\_m}{y\_m}\right) \cdot -0.5\\
              
              \mathbf{elif}\;t\_0 \leq 5 \cdot 10^{+153} \lor \neg \left(t\_0 \leq \infty\right):\\
              \;\;\;\;0.5 \cdot y\_m\\
              
              \mathbf{else}:\\
              \;\;\;\;\left(x\_m \cdot \frac{x\_m}{y\_m}\right) \cdot 0.5\\
              
              
              \end{array}
              \end{array}
              \end{array}
              
              Derivation
              1. Split input into 3 regimes
              2. if (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < 0.0

                1. Initial program 81.1%

                  \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                2. Add Preprocessing
                3. Taylor expanded in y around inf

                  \[\leadsto \color{blue}{y \cdot \left(\frac{1}{2} + \frac{1}{2} \cdot \frac{{x}^{2} - {z}^{2}}{{y}^{2}}\right)} \]
                4. Applied rewrites84.0%

                  \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{\frac{0.5}{y}}{y}, \left(x - z\right) \cdot \left(z + x\right), 0.5\right) \cdot y} \]
                5. Taylor expanded in z around inf

                  \[\leadsto \color{blue}{\frac{-1}{2} \cdot \frac{{z}^{2}}{y}} \]
                6. Step-by-step derivation
                  1. Applied rewrites29.1%

                    \[\leadsto \color{blue}{\frac{z \cdot z}{y} \cdot -0.5} \]
                  2. Step-by-step derivation
                    1. Applied rewrites31.3%

                      \[\leadsto \left(z \cdot \frac{z}{y}\right) \cdot -0.5 \]

                    if 0.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < 5.00000000000000018e153 or +inf.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64)))

                    1. Initial program 50.6%

                      \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                    2. Add Preprocessing
                    3. Taylor expanded in y around inf

                      \[\leadsto \color{blue}{\frac{1}{2} \cdot y} \]
                    4. Step-by-step derivation
                      1. Applied rewrites48.6%

                        \[\leadsto \color{blue}{0.5 \cdot y} \]

                      if 5.00000000000000018e153 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < +inf.0

                      1. Initial program 66.3%

                        \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                      2. Add Preprocessing
                      3. Taylor expanded in y around inf

                        \[\leadsto \color{blue}{y \cdot \left(\frac{1}{2} + \frac{1}{2} \cdot \frac{{x}^{2} - {z}^{2}}{{y}^{2}}\right)} \]
                      4. Applied rewrites68.0%

                        \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{\frac{0.5}{y}}{y}, \left(x - z\right) \cdot \left(z + x\right), 0.5\right) \cdot y} \]
                      5. Taylor expanded in x around inf

                        \[\leadsto \color{blue}{\frac{1}{2} \cdot \frac{{x}^{2}}{y}} \]
                      6. Step-by-step derivation
                        1. Applied rewrites40.7%

                          \[\leadsto \color{blue}{\left(x \cdot \frac{x}{y}\right) \cdot 0.5} \]
                      7. Recombined 3 regimes into one program.
                      8. Final simplification38.3%

                        \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq 0:\\ \;\;\;\;\left(z \cdot \frac{z}{y}\right) \cdot -0.5\\ \mathbf{elif}\;\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq 5 \cdot 10^{+153} \lor \neg \left(\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq \infty\right):\\ \;\;\;\;0.5 \cdot y\\ \mathbf{else}:\\ \;\;\;\;\left(x \cdot \frac{x}{y}\right) \cdot 0.5\\ \end{array} \]
                      9. Add Preprocessing

                      Alternative 3: 68.6% accurate, 0.3× speedup?

                      \[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ \begin{array}{l} t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\ y\_s \cdot \begin{array}{l} \mathbf{if}\;t\_0 \leq -5 \cdot 10^{-103}:\\ \;\;\;\;-0.5 \cdot \frac{z\_m \cdot z\_m}{y\_m}\\ \mathbf{elif}\;t\_0 \leq 5 \cdot 10^{+153} \lor \neg \left(t\_0 \leq \infty\right):\\ \;\;\;\;0.5 \cdot y\_m\\ \mathbf{else}:\\ \;\;\;\;\left(x\_m \cdot \frac{x\_m}{y\_m}\right) \cdot 0.5\\ \end{array} \end{array} \end{array} \]
                      x_m = (fabs.f64 x)
                      z_m = (fabs.f64 z)
                      y\_m = (fabs.f64 y)
                      y\_s = (copysign.f64 #s(literal 1 binary64) y)
                      (FPCore (y_s x_m y_m z_m)
                       :precision binary64
                       (let* ((t_0 (/ (- (+ (* x_m x_m) (* y_m y_m)) (* z_m z_m)) (* y_m 2.0))))
                         (*
                          y_s
                          (if (<= t_0 -5e-103)
                            (* -0.5 (/ (* z_m z_m) y_m))
                            (if (or (<= t_0 5e+153) (not (<= t_0 INFINITY)))
                              (* 0.5 y_m)
                              (* (* x_m (/ x_m y_m)) 0.5))))))
                      x_m = fabs(x);
                      z_m = fabs(z);
                      y\_m = fabs(y);
                      y\_s = copysign(1.0, y);
                      double code(double y_s, double x_m, double y_m, double z_m) {
                      	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
                      	double tmp;
                      	if (t_0 <= -5e-103) {
                      		tmp = -0.5 * ((z_m * z_m) / y_m);
                      	} else if ((t_0 <= 5e+153) || !(t_0 <= ((double) INFINITY))) {
                      		tmp = 0.5 * y_m;
                      	} else {
                      		tmp = (x_m * (x_m / y_m)) * 0.5;
                      	}
                      	return y_s * tmp;
                      }
                      
                      x_m = Math.abs(x);
                      z_m = Math.abs(z);
                      y\_m = Math.abs(y);
                      y\_s = Math.copySign(1.0, y);
                      public static double code(double y_s, double x_m, double y_m, double z_m) {
                      	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
                      	double tmp;
                      	if (t_0 <= -5e-103) {
                      		tmp = -0.5 * ((z_m * z_m) / y_m);
                      	} else if ((t_0 <= 5e+153) || !(t_0 <= Double.POSITIVE_INFINITY)) {
                      		tmp = 0.5 * y_m;
                      	} else {
                      		tmp = (x_m * (x_m / y_m)) * 0.5;
                      	}
                      	return y_s * tmp;
                      }
                      
                      x_m = math.fabs(x)
                      z_m = math.fabs(z)
                      y\_m = math.fabs(y)
                      y\_s = math.copysign(1.0, y)
                      def code(y_s, x_m, y_m, z_m):
                      	t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0)
                      	tmp = 0
                      	if t_0 <= -5e-103:
                      		tmp = -0.5 * ((z_m * z_m) / y_m)
                      	elif (t_0 <= 5e+153) or not (t_0 <= math.inf):
                      		tmp = 0.5 * y_m
                      	else:
                      		tmp = (x_m * (x_m / y_m)) * 0.5
                      	return y_s * tmp
                      
                      x_m = abs(x)
                      z_m = abs(z)
                      y\_m = abs(y)
                      y\_s = copysign(1.0, y)
                      function code(y_s, x_m, y_m, z_m)
                      	t_0 = Float64(Float64(Float64(Float64(x_m * x_m) + Float64(y_m * y_m)) - Float64(z_m * z_m)) / Float64(y_m * 2.0))
                      	tmp = 0.0
                      	if (t_0 <= -5e-103)
                      		tmp = Float64(-0.5 * Float64(Float64(z_m * z_m) / y_m));
                      	elseif ((t_0 <= 5e+153) || !(t_0 <= Inf))
                      		tmp = Float64(0.5 * y_m);
                      	else
                      		tmp = Float64(Float64(x_m * Float64(x_m / y_m)) * 0.5);
                      	end
                      	return Float64(y_s * tmp)
                      end
                      
                      x_m = abs(x);
                      z_m = abs(z);
                      y\_m = abs(y);
                      y\_s = sign(y) * abs(1.0);
                      function tmp_2 = code(y_s, x_m, y_m, z_m)
                      	t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
                      	tmp = 0.0;
                      	if (t_0 <= -5e-103)
                      		tmp = -0.5 * ((z_m * z_m) / y_m);
                      	elseif ((t_0 <= 5e+153) || ~((t_0 <= Inf)))
                      		tmp = 0.5 * y_m;
                      	else
                      		tmp = (x_m * (x_m / y_m)) * 0.5;
                      	end
                      	tmp_2 = y_s * tmp;
                      end
                      
                      x_m = N[Abs[x], $MachinePrecision]
                      z_m = N[Abs[z], $MachinePrecision]
                      y\_m = N[Abs[y], $MachinePrecision]
                      y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
                      code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := Block[{t$95$0 = N[(N[(N[(N[(x$95$m * x$95$m), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z$95$m * z$95$m), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision]}, N[(y$95$s * If[LessEqual[t$95$0, -5e-103], N[(-0.5 * N[(N[(z$95$m * z$95$m), $MachinePrecision] / y$95$m), $MachinePrecision]), $MachinePrecision], If[Or[LessEqual[t$95$0, 5e+153], N[Not[LessEqual[t$95$0, Infinity]], $MachinePrecision]], N[(0.5 * y$95$m), $MachinePrecision], N[(N[(x$95$m * N[(x$95$m / y$95$m), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision]]]), $MachinePrecision]]
                      
                      \begin{array}{l}
                      x_m = \left|x\right|
                      \\
                      z_m = \left|z\right|
                      \\
                      y\_m = \left|y\right|
                      \\
                      y\_s = \mathsf{copysign}\left(1, y\right)
                      
                      \\
                      \begin{array}{l}
                      t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\
                      y\_s \cdot \begin{array}{l}
                      \mathbf{if}\;t\_0 \leq -5 \cdot 10^{-103}:\\
                      \;\;\;\;-0.5 \cdot \frac{z\_m \cdot z\_m}{y\_m}\\
                      
                      \mathbf{elif}\;t\_0 \leq 5 \cdot 10^{+153} \lor \neg \left(t\_0 \leq \infty\right):\\
                      \;\;\;\;0.5 \cdot y\_m\\
                      
                      \mathbf{else}:\\
                      \;\;\;\;\left(x\_m \cdot \frac{x\_m}{y\_m}\right) \cdot 0.5\\
                      
                      
                      \end{array}
                      \end{array}
                      \end{array}
                      
                      Derivation
                      1. Split input into 3 regimes
                      2. if (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < -4.99999999999999966e-103

                        1. Initial program 83.0%

                          \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                        2. Add Preprocessing
                        3. Taylor expanded in z around inf

                          \[\leadsto \color{blue}{\frac{-1}{2} \cdot \frac{{z}^{2}}{y}} \]
                        4. Step-by-step derivation
                          1. Applied rewrites30.1%

                            \[\leadsto \color{blue}{-0.5 \cdot \frac{z \cdot z}{y}} \]

                          if -4.99999999999999966e-103 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < 5.00000000000000018e153 or +inf.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64)))

                          1. Initial program 49.5%

                            \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                          2. Add Preprocessing
                          3. Taylor expanded in y around inf

                            \[\leadsto \color{blue}{\frac{1}{2} \cdot y} \]
                          4. Step-by-step derivation
                            1. Applied rewrites48.2%

                              \[\leadsto \color{blue}{0.5 \cdot y} \]

                            if 5.00000000000000018e153 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < +inf.0

                            1. Initial program 66.3%

                              \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                            2. Add Preprocessing
                            3. Taylor expanded in y around inf

                              \[\leadsto \color{blue}{y \cdot \left(\frac{1}{2} + \frac{1}{2} \cdot \frac{{x}^{2} - {z}^{2}}{{y}^{2}}\right)} \]
                            4. Applied rewrites68.0%

                              \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{\frac{0.5}{y}}{y}, \left(x - z\right) \cdot \left(z + x\right), 0.5\right) \cdot y} \]
                            5. Taylor expanded in x around inf

                              \[\leadsto \color{blue}{\frac{1}{2} \cdot \frac{{x}^{2}}{y}} \]
                            6. Step-by-step derivation
                              1. Applied rewrites40.7%

                                \[\leadsto \color{blue}{\left(x \cdot \frac{x}{y}\right) \cdot 0.5} \]
                            7. Recombined 3 regimes into one program.
                            8. Final simplification38.0%

                              \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq -5 \cdot 10^{-103}:\\ \;\;\;\;-0.5 \cdot \frac{z \cdot z}{y}\\ \mathbf{elif}\;\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq 5 \cdot 10^{+153} \lor \neg \left(\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq \infty\right):\\ \;\;\;\;0.5 \cdot y\\ \mathbf{else}:\\ \;\;\;\;\left(x \cdot \frac{x}{y}\right) \cdot 0.5\\ \end{array} \]
                            9. Add Preprocessing

                            Alternative 4: 66.9% accurate, 0.3× speedup?

                            \[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ \begin{array}{l} t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\ y\_s \cdot \begin{array}{l} \mathbf{if}\;t\_0 \leq -5 \cdot 10^{-103}:\\ \;\;\;\;-0.5 \cdot \frac{z\_m \cdot z\_m}{y\_m}\\ \mathbf{elif}\;t\_0 \leq 5 \cdot 10^{+153} \lor \neg \left(t\_0 \leq \infty\right):\\ \;\;\;\;0.5 \cdot y\_m\\ \mathbf{else}:\\ \;\;\;\;\frac{x\_m \cdot x\_m}{y\_m + y\_m}\\ \end{array} \end{array} \end{array} \]
                            x_m = (fabs.f64 x)
                            z_m = (fabs.f64 z)
                            y\_m = (fabs.f64 y)
                            y\_s = (copysign.f64 #s(literal 1 binary64) y)
                            (FPCore (y_s x_m y_m z_m)
                             :precision binary64
                             (let* ((t_0 (/ (- (+ (* x_m x_m) (* y_m y_m)) (* z_m z_m)) (* y_m 2.0))))
                               (*
                                y_s
                                (if (<= t_0 -5e-103)
                                  (* -0.5 (/ (* z_m z_m) y_m))
                                  (if (or (<= t_0 5e+153) (not (<= t_0 INFINITY)))
                                    (* 0.5 y_m)
                                    (/ (* x_m x_m) (+ y_m y_m)))))))
                            x_m = fabs(x);
                            z_m = fabs(z);
                            y\_m = fabs(y);
                            y\_s = copysign(1.0, y);
                            double code(double y_s, double x_m, double y_m, double z_m) {
                            	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
                            	double tmp;
                            	if (t_0 <= -5e-103) {
                            		tmp = -0.5 * ((z_m * z_m) / y_m);
                            	} else if ((t_0 <= 5e+153) || !(t_0 <= ((double) INFINITY))) {
                            		tmp = 0.5 * y_m;
                            	} else {
                            		tmp = (x_m * x_m) / (y_m + y_m);
                            	}
                            	return y_s * tmp;
                            }
                            
                            x_m = Math.abs(x);
                            z_m = Math.abs(z);
                            y\_m = Math.abs(y);
                            y\_s = Math.copySign(1.0, y);
                            public static double code(double y_s, double x_m, double y_m, double z_m) {
                            	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
                            	double tmp;
                            	if (t_0 <= -5e-103) {
                            		tmp = -0.5 * ((z_m * z_m) / y_m);
                            	} else if ((t_0 <= 5e+153) || !(t_0 <= Double.POSITIVE_INFINITY)) {
                            		tmp = 0.5 * y_m;
                            	} else {
                            		tmp = (x_m * x_m) / (y_m + y_m);
                            	}
                            	return y_s * tmp;
                            }
                            
                            x_m = math.fabs(x)
                            z_m = math.fabs(z)
                            y\_m = math.fabs(y)
                            y\_s = math.copysign(1.0, y)
                            def code(y_s, x_m, y_m, z_m):
                            	t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0)
                            	tmp = 0
                            	if t_0 <= -5e-103:
                            		tmp = -0.5 * ((z_m * z_m) / y_m)
                            	elif (t_0 <= 5e+153) or not (t_0 <= math.inf):
                            		tmp = 0.5 * y_m
                            	else:
                            		tmp = (x_m * x_m) / (y_m + y_m)
                            	return y_s * tmp
                            
                            x_m = abs(x)
                            z_m = abs(z)
                            y\_m = abs(y)
                            y\_s = copysign(1.0, y)
                            function code(y_s, x_m, y_m, z_m)
                            	t_0 = Float64(Float64(Float64(Float64(x_m * x_m) + Float64(y_m * y_m)) - Float64(z_m * z_m)) / Float64(y_m * 2.0))
                            	tmp = 0.0
                            	if (t_0 <= -5e-103)
                            		tmp = Float64(-0.5 * Float64(Float64(z_m * z_m) / y_m));
                            	elseif ((t_0 <= 5e+153) || !(t_0 <= Inf))
                            		tmp = Float64(0.5 * y_m);
                            	else
                            		tmp = Float64(Float64(x_m * x_m) / Float64(y_m + y_m));
                            	end
                            	return Float64(y_s * tmp)
                            end
                            
                            x_m = abs(x);
                            z_m = abs(z);
                            y\_m = abs(y);
                            y\_s = sign(y) * abs(1.0);
                            function tmp_2 = code(y_s, x_m, y_m, z_m)
                            	t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
                            	tmp = 0.0;
                            	if (t_0 <= -5e-103)
                            		tmp = -0.5 * ((z_m * z_m) / y_m);
                            	elseif ((t_0 <= 5e+153) || ~((t_0 <= Inf)))
                            		tmp = 0.5 * y_m;
                            	else
                            		tmp = (x_m * x_m) / (y_m + y_m);
                            	end
                            	tmp_2 = y_s * tmp;
                            end
                            
                            x_m = N[Abs[x], $MachinePrecision]
                            z_m = N[Abs[z], $MachinePrecision]
                            y\_m = N[Abs[y], $MachinePrecision]
                            y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
                            code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := Block[{t$95$0 = N[(N[(N[(N[(x$95$m * x$95$m), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z$95$m * z$95$m), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision]}, N[(y$95$s * If[LessEqual[t$95$0, -5e-103], N[(-0.5 * N[(N[(z$95$m * z$95$m), $MachinePrecision] / y$95$m), $MachinePrecision]), $MachinePrecision], If[Or[LessEqual[t$95$0, 5e+153], N[Not[LessEqual[t$95$0, Infinity]], $MachinePrecision]], N[(0.5 * y$95$m), $MachinePrecision], N[(N[(x$95$m * x$95$m), $MachinePrecision] / N[(y$95$m + y$95$m), $MachinePrecision]), $MachinePrecision]]]), $MachinePrecision]]
                            
                            \begin{array}{l}
                            x_m = \left|x\right|
                            \\
                            z_m = \left|z\right|
                            \\
                            y\_m = \left|y\right|
                            \\
                            y\_s = \mathsf{copysign}\left(1, y\right)
                            
                            \\
                            \begin{array}{l}
                            t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\
                            y\_s \cdot \begin{array}{l}
                            \mathbf{if}\;t\_0 \leq -5 \cdot 10^{-103}:\\
                            \;\;\;\;-0.5 \cdot \frac{z\_m \cdot z\_m}{y\_m}\\
                            
                            \mathbf{elif}\;t\_0 \leq 5 \cdot 10^{+153} \lor \neg \left(t\_0 \leq \infty\right):\\
                            \;\;\;\;0.5 \cdot y\_m\\
                            
                            \mathbf{else}:\\
                            \;\;\;\;\frac{x\_m \cdot x\_m}{y\_m + y\_m}\\
                            
                            
                            \end{array}
                            \end{array}
                            \end{array}
                            
                            Derivation
                            1. Split input into 3 regimes
                            2. if (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < -4.99999999999999966e-103

                              1. Initial program 83.0%

                                \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                              2. Add Preprocessing
                              3. Taylor expanded in z around inf

                                \[\leadsto \color{blue}{\frac{-1}{2} \cdot \frac{{z}^{2}}{y}} \]
                              4. Step-by-step derivation
                                1. Applied rewrites30.1%

                                  \[\leadsto \color{blue}{-0.5 \cdot \frac{z \cdot z}{y}} \]

                                if -4.99999999999999966e-103 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < 5.00000000000000018e153 or +inf.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64)))

                                1. Initial program 49.5%

                                  \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                2. Add Preprocessing
                                3. Taylor expanded in y around inf

                                  \[\leadsto \color{blue}{\frac{1}{2} \cdot y} \]
                                4. Step-by-step derivation
                                  1. Applied rewrites48.2%

                                    \[\leadsto \color{blue}{0.5 \cdot y} \]

                                  if 5.00000000000000018e153 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < +inf.0

                                  1. Initial program 66.3%

                                    \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                  2. Add Preprocessing
                                  3. Taylor expanded in x around inf

                                    \[\leadsto \frac{\color{blue}{{x}^{2}}}{y \cdot 2} \]
                                  4. Step-by-step derivation
                                    1. Applied rewrites34.3%

                                      \[\leadsto \frac{\color{blue}{x \cdot x}}{y \cdot 2} \]
                                    2. Step-by-step derivation
                                      1. lift-*.f64N/A

                                        \[\leadsto \frac{x \cdot x}{\color{blue}{y \cdot 2}} \]
                                      2. *-commutativeN/A

                                        \[\leadsto \frac{x \cdot x}{\color{blue}{2 \cdot y}} \]
                                      3. count-2-revN/A

                                        \[\leadsto \frac{x \cdot x}{\color{blue}{y + y}} \]
                                      4. lower-+.f6434.3

                                        \[\leadsto \frac{x \cdot x}{\color{blue}{y + y}} \]
                                    3. Applied rewrites34.3%

                                      \[\leadsto \frac{x \cdot x}{\color{blue}{y + y}} \]
                                  5. Recombined 3 regimes into one program.
                                  6. Final simplification36.2%

                                    \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq -5 \cdot 10^{-103}:\\ \;\;\;\;-0.5 \cdot \frac{z \cdot z}{y}\\ \mathbf{elif}\;\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq 5 \cdot 10^{+153} \lor \neg \left(\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq \infty\right):\\ \;\;\;\;0.5 \cdot y\\ \mathbf{else}:\\ \;\;\;\;\frac{x \cdot x}{y + y}\\ \end{array} \]
                                  7. Add Preprocessing

                                  Alternative 5: 96.3% accurate, 0.3× speedup?

                                  \[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ \begin{array}{l} t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\ y\_s \cdot \begin{array}{l} \mathbf{if}\;t\_0 \leq 0 \lor \neg \left(t\_0 \leq \infty\right):\\ \;\;\;\;\left(y\_m - \frac{z\_m}{y\_m} \cdot z\_m\right) \cdot 0.5\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\frac{x\_m}{y\_m}, x\_m, y\_m\right) \cdot 0.5\\ \end{array} \end{array} \end{array} \]
                                  x_m = (fabs.f64 x)
                                  z_m = (fabs.f64 z)
                                  y\_m = (fabs.f64 y)
                                  y\_s = (copysign.f64 #s(literal 1 binary64) y)
                                  (FPCore (y_s x_m y_m z_m)
                                   :precision binary64
                                   (let* ((t_0 (/ (- (+ (* x_m x_m) (* y_m y_m)) (* z_m z_m)) (* y_m 2.0))))
                                     (*
                                      y_s
                                      (if (or (<= t_0 0.0) (not (<= t_0 INFINITY)))
                                        (* (- y_m (* (/ z_m y_m) z_m)) 0.5)
                                        (* (fma (/ x_m y_m) x_m y_m) 0.5)))))
                                  x_m = fabs(x);
                                  z_m = fabs(z);
                                  y\_m = fabs(y);
                                  y\_s = copysign(1.0, y);
                                  double code(double y_s, double x_m, double y_m, double z_m) {
                                  	double t_0 = (((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0);
                                  	double tmp;
                                  	if ((t_0 <= 0.0) || !(t_0 <= ((double) INFINITY))) {
                                  		tmp = (y_m - ((z_m / y_m) * z_m)) * 0.5;
                                  	} else {
                                  		tmp = fma((x_m / y_m), x_m, y_m) * 0.5;
                                  	}
                                  	return y_s * tmp;
                                  }
                                  
                                  x_m = abs(x)
                                  z_m = abs(z)
                                  y\_m = abs(y)
                                  y\_s = copysign(1.0, y)
                                  function code(y_s, x_m, y_m, z_m)
                                  	t_0 = Float64(Float64(Float64(Float64(x_m * x_m) + Float64(y_m * y_m)) - Float64(z_m * z_m)) / Float64(y_m * 2.0))
                                  	tmp = 0.0
                                  	if ((t_0 <= 0.0) || !(t_0 <= Inf))
                                  		tmp = Float64(Float64(y_m - Float64(Float64(z_m / y_m) * z_m)) * 0.5);
                                  	else
                                  		tmp = Float64(fma(Float64(x_m / y_m), x_m, y_m) * 0.5);
                                  	end
                                  	return Float64(y_s * tmp)
                                  end
                                  
                                  x_m = N[Abs[x], $MachinePrecision]
                                  z_m = N[Abs[z], $MachinePrecision]
                                  y\_m = N[Abs[y], $MachinePrecision]
                                  y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
                                  code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := Block[{t$95$0 = N[(N[(N[(N[(x$95$m * x$95$m), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z$95$m * z$95$m), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision]}, N[(y$95$s * If[Or[LessEqual[t$95$0, 0.0], N[Not[LessEqual[t$95$0, Infinity]], $MachinePrecision]], N[(N[(y$95$m - N[(N[(z$95$m / y$95$m), $MachinePrecision] * z$95$m), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision], N[(N[(N[(x$95$m / y$95$m), $MachinePrecision] * x$95$m + y$95$m), $MachinePrecision] * 0.5), $MachinePrecision]]), $MachinePrecision]]
                                  
                                  \begin{array}{l}
                                  x_m = \left|x\right|
                                  \\
                                  z_m = \left|z\right|
                                  \\
                                  y\_m = \left|y\right|
                                  \\
                                  y\_s = \mathsf{copysign}\left(1, y\right)
                                  
                                  \\
                                  \begin{array}{l}
                                  t_0 := \frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2}\\
                                  y\_s \cdot \begin{array}{l}
                                  \mathbf{if}\;t\_0 \leq 0 \lor \neg \left(t\_0 \leq \infty\right):\\
                                  \;\;\;\;\left(y\_m - \frac{z\_m}{y\_m} \cdot z\_m\right) \cdot 0.5\\
                                  
                                  \mathbf{else}:\\
                                  \;\;\;\;\mathsf{fma}\left(\frac{x\_m}{y\_m}, x\_m, y\_m\right) \cdot 0.5\\
                                  
                                  
                                  \end{array}
                                  \end{array}
                                  \end{array}
                                  
                                  Derivation
                                  1. Split input into 2 regimes
                                  2. if (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < 0.0 or +inf.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64)))

                                    1. Initial program 64.0%

                                      \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                    2. Add Preprocessing
                                    3. Taylor expanded in x around 0

                                      \[\leadsto \color{blue}{\frac{1}{2} \cdot \frac{{y}^{2} - {z}^{2}}{y}} \]
                                    4. Step-by-step derivation
                                      1. Applied rewrites59.4%

                                        \[\leadsto \color{blue}{\left(y - \frac{z \cdot z}{y}\right) \cdot 0.5} \]
                                      2. Step-by-step derivation
                                        1. Applied rewrites69.3%

                                          \[\leadsto \left(y - \frac{z}{y} \cdot z\right) \cdot 0.5 \]

                                        if 0.0 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < +inf.0

                                        1. Initial program 76.9%

                                          \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                        2. Add Preprocessing
                                        3. Taylor expanded in z around 0

                                          \[\leadsto \color{blue}{\frac{1}{2} \cdot \frac{{x}^{2} + {y}^{2}}{y}} \]
                                        4. Step-by-step derivation
                                          1. Applied rewrites73.2%

                                            \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{x}{y}, x, y\right) \cdot 0.5} \]
                                        5. Recombined 2 regimes into one program.
                                        6. Final simplification70.9%

                                          \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq 0 \lor \neg \left(\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \leq \infty\right):\\ \;\;\;\;\left(y - \frac{z}{y} \cdot z\right) \cdot 0.5\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\frac{x}{y}, x, y\right) \cdot 0.5\\ \end{array} \]
                                        7. Add Preprocessing

                                        Alternative 6: 93.3% accurate, 0.6× speedup?

                                        \[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ y\_s \cdot \begin{array}{l} \mathbf{if}\;\frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2} \leq -5 \cdot 10^{-103}:\\ \;\;\;\;\left(z\_m \cdot \frac{z\_m}{y\_m}\right) \cdot -0.5\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\frac{x\_m}{y\_m}, x\_m, y\_m\right) \cdot 0.5\\ \end{array} \end{array} \]
                                        x_m = (fabs.f64 x)
                                        z_m = (fabs.f64 z)
                                        y\_m = (fabs.f64 y)
                                        y\_s = (copysign.f64 #s(literal 1 binary64) y)
                                        (FPCore (y_s x_m y_m z_m)
                                         :precision binary64
                                         (*
                                          y_s
                                          (if (<= (/ (- (+ (* x_m x_m) (* y_m y_m)) (* z_m z_m)) (* y_m 2.0)) -5e-103)
                                            (* (* z_m (/ z_m y_m)) -0.5)
                                            (* (fma (/ x_m y_m) x_m y_m) 0.5))))
                                        x_m = fabs(x);
                                        z_m = fabs(z);
                                        y\_m = fabs(y);
                                        y\_s = copysign(1.0, y);
                                        double code(double y_s, double x_m, double y_m, double z_m) {
                                        	double tmp;
                                        	if (((((x_m * x_m) + (y_m * y_m)) - (z_m * z_m)) / (y_m * 2.0)) <= -5e-103) {
                                        		tmp = (z_m * (z_m / y_m)) * -0.5;
                                        	} else {
                                        		tmp = fma((x_m / y_m), x_m, y_m) * 0.5;
                                        	}
                                        	return y_s * tmp;
                                        }
                                        
                                        x_m = abs(x)
                                        z_m = abs(z)
                                        y\_m = abs(y)
                                        y\_s = copysign(1.0, y)
                                        function code(y_s, x_m, y_m, z_m)
                                        	tmp = 0.0
                                        	if (Float64(Float64(Float64(Float64(x_m * x_m) + Float64(y_m * y_m)) - Float64(z_m * z_m)) / Float64(y_m * 2.0)) <= -5e-103)
                                        		tmp = Float64(Float64(z_m * Float64(z_m / y_m)) * -0.5);
                                        	else
                                        		tmp = Float64(fma(Float64(x_m / y_m), x_m, y_m) * 0.5);
                                        	end
                                        	return Float64(y_s * tmp)
                                        end
                                        
                                        x_m = N[Abs[x], $MachinePrecision]
                                        z_m = N[Abs[z], $MachinePrecision]
                                        y\_m = N[Abs[y], $MachinePrecision]
                                        y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
                                        code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := N[(y$95$s * If[LessEqual[N[(N[(N[(N[(x$95$m * x$95$m), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z$95$m * z$95$m), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision], -5e-103], N[(N[(z$95$m * N[(z$95$m / y$95$m), $MachinePrecision]), $MachinePrecision] * -0.5), $MachinePrecision], N[(N[(N[(x$95$m / y$95$m), $MachinePrecision] * x$95$m + y$95$m), $MachinePrecision] * 0.5), $MachinePrecision]]), $MachinePrecision]
                                        
                                        \begin{array}{l}
                                        x_m = \left|x\right|
                                        \\
                                        z_m = \left|z\right|
                                        \\
                                        y\_m = \left|y\right|
                                        \\
                                        y\_s = \mathsf{copysign}\left(1, y\right)
                                        
                                        \\
                                        y\_s \cdot \begin{array}{l}
                                        \mathbf{if}\;\frac{\left(x\_m \cdot x\_m + y\_m \cdot y\_m\right) - z\_m \cdot z\_m}{y\_m \cdot 2} \leq -5 \cdot 10^{-103}:\\
                                        \;\;\;\;\left(z\_m \cdot \frac{z\_m}{y\_m}\right) \cdot -0.5\\
                                        
                                        \mathbf{else}:\\
                                        \;\;\;\;\mathsf{fma}\left(\frac{x\_m}{y\_m}, x\_m, y\_m\right) \cdot 0.5\\
                                        
                                        
                                        \end{array}
                                        \end{array}
                                        
                                        Derivation
                                        1. Split input into 2 regimes
                                        2. if (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < -4.99999999999999966e-103

                                          1. Initial program 83.0%

                                            \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                          2. Add Preprocessing
                                          3. Taylor expanded in y around inf

                                            \[\leadsto \color{blue}{y \cdot \left(\frac{1}{2} + \frac{1}{2} \cdot \frac{{x}^{2} - {z}^{2}}{{y}^{2}}\right)} \]
                                          4. Applied rewrites86.8%

                                            \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{\frac{0.5}{y}}{y}, \left(x - z\right) \cdot \left(z + x\right), 0.5\right) \cdot y} \]
                                          5. Taylor expanded in z around inf

                                            \[\leadsto \color{blue}{\frac{-1}{2} \cdot \frac{{z}^{2}}{y}} \]
                                          6. Step-by-step derivation
                                            1. Applied rewrites30.1%

                                              \[\leadsto \color{blue}{\frac{z \cdot z}{y} \cdot -0.5} \]
                                            2. Step-by-step derivation
                                              1. Applied rewrites31.7%

                                                \[\leadsto \left(z \cdot \frac{z}{y}\right) \cdot -0.5 \]

                                              if -4.99999999999999966e-103 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64)))

                                              1. Initial program 58.0%

                                                \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                              2. Add Preprocessing
                                              3. Taylor expanded in z around 0

                                                \[\leadsto \color{blue}{\frac{1}{2} \cdot \frac{{x}^{2} + {y}^{2}}{y}} \]
                                              4. Step-by-step derivation
                                                1. Applied rewrites70.9%

                                                  \[\leadsto \color{blue}{\mathsf{fma}\left(\frac{x}{y}, x, y\right) \cdot 0.5} \]
                                              5. Recombined 2 regimes into one program.
                                              6. Add Preprocessing

                                              Alternative 7: 50.1% accurate, 1.5× speedup?

                                              \[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ y\_s \cdot \begin{array}{l} \mathbf{if}\;x\_m \leq 260000000000:\\ \;\;\;\;0.5 \cdot y\_m\\ \mathbf{else}:\\ \;\;\;\;\frac{x\_m \cdot x\_m}{y\_m + y\_m}\\ \end{array} \end{array} \]
                                              x_m = (fabs.f64 x)
                                              z_m = (fabs.f64 z)
                                              y\_m = (fabs.f64 y)
                                              y\_s = (copysign.f64 #s(literal 1 binary64) y)
                                              (FPCore (y_s x_m y_m z_m)
                                               :precision binary64
                                               (* y_s (if (<= x_m 260000000000.0) (* 0.5 y_m) (/ (* x_m x_m) (+ y_m y_m)))))
                                              x_m = fabs(x);
                                              z_m = fabs(z);
                                              y\_m = fabs(y);
                                              y\_s = copysign(1.0, y);
                                              double code(double y_s, double x_m, double y_m, double z_m) {
                                              	double tmp;
                                              	if (x_m <= 260000000000.0) {
                                              		tmp = 0.5 * y_m;
                                              	} else {
                                              		tmp = (x_m * x_m) / (y_m + y_m);
                                              	}
                                              	return y_s * tmp;
                                              }
                                              
                                              x_m =     private
                                              z_m =     private
                                              y\_m =     private
                                              y\_s =     private
                                              module fmin_fmax_functions
                                                  implicit none
                                                  private
                                                  public fmax
                                                  public fmin
                                              
                                                  interface fmax
                                                      module procedure fmax88
                                                      module procedure fmax44
                                                      module procedure fmax84
                                                      module procedure fmax48
                                                  end interface
                                                  interface fmin
                                                      module procedure fmin88
                                                      module procedure fmin44
                                                      module procedure fmin84
                                                      module procedure fmin48
                                                  end interface
                                              contains
                                                  real(8) function fmax88(x, y) result (res)
                                                      real(8), intent (in) :: x
                                                      real(8), intent (in) :: y
                                                      res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                  end function
                                                  real(4) function fmax44(x, y) result (res)
                                                      real(4), intent (in) :: x
                                                      real(4), intent (in) :: y
                                                      res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                  end function
                                                  real(8) function fmax84(x, y) result(res)
                                                      real(8), intent (in) :: x
                                                      real(4), intent (in) :: y
                                                      res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
                                                  end function
                                                  real(8) function fmax48(x, y) result(res)
                                                      real(4), intent (in) :: x
                                                      real(8), intent (in) :: y
                                                      res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
                                                  end function
                                                  real(8) function fmin88(x, y) result (res)
                                                      real(8), intent (in) :: x
                                                      real(8), intent (in) :: y
                                                      res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                  end function
                                                  real(4) function fmin44(x, y) result (res)
                                                      real(4), intent (in) :: x
                                                      real(4), intent (in) :: y
                                                      res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                  end function
                                                  real(8) function fmin84(x, y) result(res)
                                                      real(8), intent (in) :: x
                                                      real(4), intent (in) :: y
                                                      res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
                                                  end function
                                                  real(8) function fmin48(x, y) result(res)
                                                      real(4), intent (in) :: x
                                                      real(8), intent (in) :: y
                                                      res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
                                                  end function
                                              end module
                                              
                                              real(8) function code(y_s, x_m, y_m, z_m)
                                              use fmin_fmax_functions
                                                  real(8), intent (in) :: y_s
                                                  real(8), intent (in) :: x_m
                                                  real(8), intent (in) :: y_m
                                                  real(8), intent (in) :: z_m
                                                  real(8) :: tmp
                                                  if (x_m <= 260000000000.0d0) then
                                                      tmp = 0.5d0 * y_m
                                                  else
                                                      tmp = (x_m * x_m) / (y_m + y_m)
                                                  end if
                                                  code = y_s * tmp
                                              end function
                                              
                                              x_m = Math.abs(x);
                                              z_m = Math.abs(z);
                                              y\_m = Math.abs(y);
                                              y\_s = Math.copySign(1.0, y);
                                              public static double code(double y_s, double x_m, double y_m, double z_m) {
                                              	double tmp;
                                              	if (x_m <= 260000000000.0) {
                                              		tmp = 0.5 * y_m;
                                              	} else {
                                              		tmp = (x_m * x_m) / (y_m + y_m);
                                              	}
                                              	return y_s * tmp;
                                              }
                                              
                                              x_m = math.fabs(x)
                                              z_m = math.fabs(z)
                                              y\_m = math.fabs(y)
                                              y\_s = math.copysign(1.0, y)
                                              def code(y_s, x_m, y_m, z_m):
                                              	tmp = 0
                                              	if x_m <= 260000000000.0:
                                              		tmp = 0.5 * y_m
                                              	else:
                                              		tmp = (x_m * x_m) / (y_m + y_m)
                                              	return y_s * tmp
                                              
                                              x_m = abs(x)
                                              z_m = abs(z)
                                              y\_m = abs(y)
                                              y\_s = copysign(1.0, y)
                                              function code(y_s, x_m, y_m, z_m)
                                              	tmp = 0.0
                                              	if (x_m <= 260000000000.0)
                                              		tmp = Float64(0.5 * y_m);
                                              	else
                                              		tmp = Float64(Float64(x_m * x_m) / Float64(y_m + y_m));
                                              	end
                                              	return Float64(y_s * tmp)
                                              end
                                              
                                              x_m = abs(x);
                                              z_m = abs(z);
                                              y\_m = abs(y);
                                              y\_s = sign(y) * abs(1.0);
                                              function tmp_2 = code(y_s, x_m, y_m, z_m)
                                              	tmp = 0.0;
                                              	if (x_m <= 260000000000.0)
                                              		tmp = 0.5 * y_m;
                                              	else
                                              		tmp = (x_m * x_m) / (y_m + y_m);
                                              	end
                                              	tmp_2 = y_s * tmp;
                                              end
                                              
                                              x_m = N[Abs[x], $MachinePrecision]
                                              z_m = N[Abs[z], $MachinePrecision]
                                              y\_m = N[Abs[y], $MachinePrecision]
                                              y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
                                              code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := N[(y$95$s * If[LessEqual[x$95$m, 260000000000.0], N[(0.5 * y$95$m), $MachinePrecision], N[(N[(x$95$m * x$95$m), $MachinePrecision] / N[(y$95$m + y$95$m), $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]
                                              
                                              \begin{array}{l}
                                              x_m = \left|x\right|
                                              \\
                                              z_m = \left|z\right|
                                              \\
                                              y\_m = \left|y\right|
                                              \\
                                              y\_s = \mathsf{copysign}\left(1, y\right)
                                              
                                              \\
                                              y\_s \cdot \begin{array}{l}
                                              \mathbf{if}\;x\_m \leq 260000000000:\\
                                              \;\;\;\;0.5 \cdot y\_m\\
                                              
                                              \mathbf{else}:\\
                                              \;\;\;\;\frac{x\_m \cdot x\_m}{y\_m + y\_m}\\
                                              
                                              
                                              \end{array}
                                              \end{array}
                                              
                                              Derivation
                                              1. Split input into 2 regimes
                                              2. if x < 2.6e11

                                                1. Initial program 70.9%

                                                  \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                                2. Add Preprocessing
                                                3. Taylor expanded in y around inf

                                                  \[\leadsto \color{blue}{\frac{1}{2} \cdot y} \]
                                                4. Step-by-step derivation
                                                  1. Applied rewrites42.4%

                                                    \[\leadsto \color{blue}{0.5 \cdot y} \]

                                                  if 2.6e11 < x

                                                  1. Initial program 63.2%

                                                    \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                                  2. Add Preprocessing
                                                  3. Taylor expanded in x around inf

                                                    \[\leadsto \frac{\color{blue}{{x}^{2}}}{y \cdot 2} \]
                                                  4. Step-by-step derivation
                                                    1. Applied rewrites54.6%

                                                      \[\leadsto \frac{\color{blue}{x \cdot x}}{y \cdot 2} \]
                                                    2. Step-by-step derivation
                                                      1. lift-*.f64N/A

                                                        \[\leadsto \frac{x \cdot x}{\color{blue}{y \cdot 2}} \]
                                                      2. *-commutativeN/A

                                                        \[\leadsto \frac{x \cdot x}{\color{blue}{2 \cdot y}} \]
                                                      3. count-2-revN/A

                                                        \[\leadsto \frac{x \cdot x}{\color{blue}{y + y}} \]
                                                      4. lower-+.f6454.6

                                                        \[\leadsto \frac{x \cdot x}{\color{blue}{y + y}} \]
                                                    3. Applied rewrites54.6%

                                                      \[\leadsto \frac{x \cdot x}{\color{blue}{y + y}} \]
                                                  5. Recombined 2 regimes into one program.
                                                  6. Add Preprocessing

                                                  Alternative 8: 34.4% accurate, 6.3× speedup?

                                                  \[\begin{array}{l} x_m = \left|x\right| \\ z_m = \left|z\right| \\ y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ y\_s \cdot \left(0.5 \cdot y\_m\right) \end{array} \]
                                                  x_m = (fabs.f64 x)
                                                  z_m = (fabs.f64 z)
                                                  y\_m = (fabs.f64 y)
                                                  y\_s = (copysign.f64 #s(literal 1 binary64) y)
                                                  (FPCore (y_s x_m y_m z_m) :precision binary64 (* y_s (* 0.5 y_m)))
                                                  x_m = fabs(x);
                                                  z_m = fabs(z);
                                                  y\_m = fabs(y);
                                                  y\_s = copysign(1.0, y);
                                                  double code(double y_s, double x_m, double y_m, double z_m) {
                                                  	return y_s * (0.5 * y_m);
                                                  }
                                                  
                                                  x_m =     private
                                                  z_m =     private
                                                  y\_m =     private
                                                  y\_s =     private
                                                  module fmin_fmax_functions
                                                      implicit none
                                                      private
                                                      public fmax
                                                      public fmin
                                                  
                                                      interface fmax
                                                          module procedure fmax88
                                                          module procedure fmax44
                                                          module procedure fmax84
                                                          module procedure fmax48
                                                      end interface
                                                      interface fmin
                                                          module procedure fmin88
                                                          module procedure fmin44
                                                          module procedure fmin84
                                                          module procedure fmin48
                                                      end interface
                                                  contains
                                                      real(8) function fmax88(x, y) result (res)
                                                          real(8), intent (in) :: x
                                                          real(8), intent (in) :: y
                                                          res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                      end function
                                                      real(4) function fmax44(x, y) result (res)
                                                          real(4), intent (in) :: x
                                                          real(4), intent (in) :: y
                                                          res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                      end function
                                                      real(8) function fmax84(x, y) result(res)
                                                          real(8), intent (in) :: x
                                                          real(4), intent (in) :: y
                                                          res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
                                                      end function
                                                      real(8) function fmax48(x, y) result(res)
                                                          real(4), intent (in) :: x
                                                          real(8), intent (in) :: y
                                                          res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
                                                      end function
                                                      real(8) function fmin88(x, y) result (res)
                                                          real(8), intent (in) :: x
                                                          real(8), intent (in) :: y
                                                          res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                      end function
                                                      real(4) function fmin44(x, y) result (res)
                                                          real(4), intent (in) :: x
                                                          real(4), intent (in) :: y
                                                          res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                      end function
                                                      real(8) function fmin84(x, y) result(res)
                                                          real(8), intent (in) :: x
                                                          real(4), intent (in) :: y
                                                          res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
                                                      end function
                                                      real(8) function fmin48(x, y) result(res)
                                                          real(4), intent (in) :: x
                                                          real(8), intent (in) :: y
                                                          res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
                                                      end function
                                                  end module
                                                  
                                                  real(8) function code(y_s, x_m, y_m, z_m)
                                                  use fmin_fmax_functions
                                                      real(8), intent (in) :: y_s
                                                      real(8), intent (in) :: x_m
                                                      real(8), intent (in) :: y_m
                                                      real(8), intent (in) :: z_m
                                                      code = y_s * (0.5d0 * y_m)
                                                  end function
                                                  
                                                  x_m = Math.abs(x);
                                                  z_m = Math.abs(z);
                                                  y\_m = Math.abs(y);
                                                  y\_s = Math.copySign(1.0, y);
                                                  public static double code(double y_s, double x_m, double y_m, double z_m) {
                                                  	return y_s * (0.5 * y_m);
                                                  }
                                                  
                                                  x_m = math.fabs(x)
                                                  z_m = math.fabs(z)
                                                  y\_m = math.fabs(y)
                                                  y\_s = math.copysign(1.0, y)
                                                  def code(y_s, x_m, y_m, z_m):
                                                  	return y_s * (0.5 * y_m)
                                                  
                                                  x_m = abs(x)
                                                  z_m = abs(z)
                                                  y\_m = abs(y)
                                                  y\_s = copysign(1.0, y)
                                                  function code(y_s, x_m, y_m, z_m)
                                                  	return Float64(y_s * Float64(0.5 * y_m))
                                                  end
                                                  
                                                  x_m = abs(x);
                                                  z_m = abs(z);
                                                  y\_m = abs(y);
                                                  y\_s = sign(y) * abs(1.0);
                                                  function tmp = code(y_s, x_m, y_m, z_m)
                                                  	tmp = y_s * (0.5 * y_m);
                                                  end
                                                  
                                                  x_m = N[Abs[x], $MachinePrecision]
                                                  z_m = N[Abs[z], $MachinePrecision]
                                                  y\_m = N[Abs[y], $MachinePrecision]
                                                  y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
                                                  code[y$95$s_, x$95$m_, y$95$m_, z$95$m_] := N[(y$95$s * N[(0.5 * y$95$m), $MachinePrecision]), $MachinePrecision]
                                                  
                                                  \begin{array}{l}
                                                  x_m = \left|x\right|
                                                  \\
                                                  z_m = \left|z\right|
                                                  \\
                                                  y\_m = \left|y\right|
                                                  \\
                                                  y\_s = \mathsf{copysign}\left(1, y\right)
                                                  
                                                  \\
                                                  y\_s \cdot \left(0.5 \cdot y\_m\right)
                                                  \end{array}
                                                  
                                                  Derivation
                                                  1. Initial program 69.2%

                                                    \[\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2} \]
                                                  2. Add Preprocessing
                                                  3. Taylor expanded in y around inf

                                                    \[\leadsto \color{blue}{\frac{1}{2} \cdot y} \]
                                                  4. Step-by-step derivation
                                                    1. Applied rewrites37.4%

                                                      \[\leadsto \color{blue}{0.5 \cdot y} \]
                                                    2. Add Preprocessing

                                                    Developer Target 1: 99.9% accurate, 1.1× speedup?

                                                    \[\begin{array}{l} \\ y \cdot 0.5 - \left(\frac{0.5}{y} \cdot \left(z + x\right)\right) \cdot \left(z - x\right) \end{array} \]
                                                    (FPCore (x y z)
                                                     :precision binary64
                                                     (- (* y 0.5) (* (* (/ 0.5 y) (+ z x)) (- z x))))
                                                    double code(double x, double y, double z) {
                                                    	return (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x));
                                                    }
                                                    
                                                    module fmin_fmax_functions
                                                        implicit none
                                                        private
                                                        public fmax
                                                        public fmin
                                                    
                                                        interface fmax
                                                            module procedure fmax88
                                                            module procedure fmax44
                                                            module procedure fmax84
                                                            module procedure fmax48
                                                        end interface
                                                        interface fmin
                                                            module procedure fmin88
                                                            module procedure fmin44
                                                            module procedure fmin84
                                                            module procedure fmin48
                                                        end interface
                                                    contains
                                                        real(8) function fmax88(x, y) result (res)
                                                            real(8), intent (in) :: x
                                                            real(8), intent (in) :: y
                                                            res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                        end function
                                                        real(4) function fmax44(x, y) result (res)
                                                            real(4), intent (in) :: x
                                                            real(4), intent (in) :: y
                                                            res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                        end function
                                                        real(8) function fmax84(x, y) result(res)
                                                            real(8), intent (in) :: x
                                                            real(4), intent (in) :: y
                                                            res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
                                                        end function
                                                        real(8) function fmax48(x, y) result(res)
                                                            real(4), intent (in) :: x
                                                            real(8), intent (in) :: y
                                                            res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
                                                        end function
                                                        real(8) function fmin88(x, y) result (res)
                                                            real(8), intent (in) :: x
                                                            real(8), intent (in) :: y
                                                            res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                        end function
                                                        real(4) function fmin44(x, y) result (res)
                                                            real(4), intent (in) :: x
                                                            real(4), intent (in) :: y
                                                            res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                        end function
                                                        real(8) function fmin84(x, y) result(res)
                                                            real(8), intent (in) :: x
                                                            real(4), intent (in) :: y
                                                            res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
                                                        end function
                                                        real(8) function fmin48(x, y) result(res)
                                                            real(4), intent (in) :: x
                                                            real(8), intent (in) :: y
                                                            res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
                                                        end function
                                                    end module
                                                    
                                                    real(8) function code(x, y, z)
                                                    use fmin_fmax_functions
                                                        real(8), intent (in) :: x
                                                        real(8), intent (in) :: y
                                                        real(8), intent (in) :: z
                                                        code = (y * 0.5d0) - (((0.5d0 / y) * (z + x)) * (z - x))
                                                    end function
                                                    
                                                    public static double code(double x, double y, double z) {
                                                    	return (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x));
                                                    }
                                                    
                                                    def code(x, y, z):
                                                    	return (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x))
                                                    
                                                    function code(x, y, z)
                                                    	return Float64(Float64(y * 0.5) - Float64(Float64(Float64(0.5 / y) * Float64(z + x)) * Float64(z - x)))
                                                    end
                                                    
                                                    function tmp = code(x, y, z)
                                                    	tmp = (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x));
                                                    end
                                                    
                                                    code[x_, y_, z_] := N[(N[(y * 0.5), $MachinePrecision] - N[(N[(N[(0.5 / y), $MachinePrecision] * N[(z + x), $MachinePrecision]), $MachinePrecision] * N[(z - x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
                                                    
                                                    \begin{array}{l}
                                                    
                                                    \\
                                                    y \cdot 0.5 - \left(\frac{0.5}{y} \cdot \left(z + x\right)\right) \cdot \left(z - x\right)
                                                    \end{array}
                                                    

                                                    Reproduce

                                                    ?
                                                    herbie shell --seed 2025026 
                                                    (FPCore (x y z)
                                                      :name "Diagrams.TwoD.Apollonian:initialConfig from diagrams-contrib-1.3.0.5, A"
                                                      :precision binary64
                                                    
                                                      :alt
                                                      (! :herbie-platform default (- (* y 1/2) (* (* (/ 1/2 y) (+ z x)) (- z x))))
                                                    
                                                      (/ (- (+ (* x x) (* y y)) (* z z)) (* y 2.0)))