Statistics.Distribution.Beta:$cvariance from math-functions-0.1.5.2

Percentage Accurate: 82.9% → 98.2%
Time: 10.0s
Alternatives: 13
Speedup: 0.9×

Specification

?
\[\begin{array}{l} \\ \frac{x \cdot y}{\left(z \cdot z\right) \cdot \left(z + 1\right)} \end{array} \]
(FPCore (x y z) :precision binary64 (/ (* x y) (* (* z z) (+ z 1.0))))
double code(double x, double y, double z) {
	return (x * y) / ((z * z) * (z + 1.0));
}
real(8) function code(x, y, z)
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    real(8), intent (in) :: z
    code = (x * y) / ((z * z) * (z + 1.0d0))
end function
public static double code(double x, double y, double z) {
	return (x * y) / ((z * z) * (z + 1.0));
}
def code(x, y, z):
	return (x * y) / ((z * z) * (z + 1.0))
function code(x, y, z)
	return Float64(Float64(x * y) / Float64(Float64(z * z) * Float64(z + 1.0)))
end
function tmp = code(x, y, z)
	tmp = (x * y) / ((z * z) * (z + 1.0));
end
code[x_, y_, z_] := N[(N[(x * y), $MachinePrecision] / N[(N[(z * z), $MachinePrecision] * N[(z + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}

\\
\frac{x \cdot y}{\left(z \cdot z\right) \cdot \left(z + 1\right)}
\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 13 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: 82.9% accurate, 1.0× speedup?

\[\begin{array}{l} \\ \frac{x \cdot y}{\left(z \cdot z\right) \cdot \left(z + 1\right)} \end{array} \]
(FPCore (x y z) :precision binary64 (/ (* x y) (* (* z z) (+ z 1.0))))
double code(double x, double y, double z) {
	return (x * y) / ((z * z) * (z + 1.0));
}
real(8) function code(x, y, z)
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    real(8), intent (in) :: z
    code = (x * y) / ((z * z) * (z + 1.0d0))
end function
public static double code(double x, double y, double z) {
	return (x * y) / ((z * z) * (z + 1.0));
}
def code(x, y, z):
	return (x * y) / ((z * z) * (z + 1.0))
function code(x, y, z)
	return Float64(Float64(x * y) / Float64(Float64(z * z) * Float64(z + 1.0)))
end
function tmp = code(x, y, z)
	tmp = (x * y) / ((z * z) * (z + 1.0));
end
code[x_, y_, z_] := N[(N[(x * y), $MachinePrecision] / N[(N[(z * z), $MachinePrecision] * N[(z + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}

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

Alternative 1: 98.2% accurate, 0.5× speedup?

\[\begin{array}{l} y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ x\_m = \left|x\right| \\ x\_s = \mathsf{copysign}\left(1, x\right) \\ [x_m, y_m, z] = \mathsf{sort}([x_m, y_m, z])\\ \\ x\_s \cdot \left(y\_s \cdot \begin{array}{l} \mathbf{if}\;y\_m \cdot x\_m \leq 5 \cdot 10^{-312}:\\ \;\;\;\;\frac{\frac{x\_m}{z}}{\frac{z}{y\_m}}\\ \mathbf{elif}\;y\_m \cdot x\_m \leq 10^{+177}:\\ \;\;\;\;\frac{\frac{y\_m \cdot x\_m}{z}}{\mathsf{fma}\left(z, z, z\right)}\\ \mathbf{else}:\\ \;\;\;\;\frac{x\_m}{z \cdot \left(z \cdot \frac{z + 1}{y\_m}\right)}\\ \end{array}\right) \end{array} \]
y\_m = (fabs.f64 y)
y\_s = (copysign.f64 #s(literal 1 binary64) y)
x\_m = (fabs.f64 x)
x\_s = (copysign.f64 #s(literal 1 binary64) x)
NOTE: x_m, y_m, and z should be sorted in increasing order before calling this function.
(FPCore (x_s y_s x_m y_m z)
 :precision binary64
 (*
  x_s
  (*
   y_s
   (if (<= (* y_m x_m) 5e-312)
     (/ (/ x_m z) (/ z y_m))
     (if (<= (* y_m x_m) 1e+177)
       (/ (/ (* y_m x_m) z) (fma z z z))
       (/ x_m (* z (* z (/ (+ z 1.0) y_m)))))))))
y\_m = fabs(y);
y\_s = copysign(1.0, y);
x\_m = fabs(x);
x\_s = copysign(1.0, x);
assert(x_m < y_m && y_m < z);
double code(double x_s, double y_s, double x_m, double y_m, double z) {
	double tmp;
	if ((y_m * x_m) <= 5e-312) {
		tmp = (x_m / z) / (z / y_m);
	} else if ((y_m * x_m) <= 1e+177) {
		tmp = ((y_m * x_m) / z) / fma(z, z, z);
	} else {
		tmp = x_m / (z * (z * ((z + 1.0) / y_m)));
	}
	return x_s * (y_s * tmp);
}
y\_m = abs(y)
y\_s = copysign(1.0, y)
x\_m = abs(x)
x\_s = copysign(1.0, x)
x_m, y_m, z = sort([x_m, y_m, z])
function code(x_s, y_s, x_m, y_m, z)
	tmp = 0.0
	if (Float64(y_m * x_m) <= 5e-312)
		tmp = Float64(Float64(x_m / z) / Float64(z / y_m));
	elseif (Float64(y_m * x_m) <= 1e+177)
		tmp = Float64(Float64(Float64(y_m * x_m) / z) / fma(z, z, z));
	else
		tmp = Float64(x_m / Float64(z * Float64(z * Float64(Float64(z + 1.0) / y_m))));
	end
	return Float64(x_s * Float64(y_s * tmp))
end
y\_m = N[Abs[y], $MachinePrecision]
y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
x\_m = N[Abs[x], $MachinePrecision]
x\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[x]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: x_m, y_m, and z should be sorted in increasing order before calling this function.
code[x$95$s_, y$95$s_, x$95$m_, y$95$m_, z_] := N[(x$95$s * N[(y$95$s * If[LessEqual[N[(y$95$m * x$95$m), $MachinePrecision], 5e-312], N[(N[(x$95$m / z), $MachinePrecision] / N[(z / y$95$m), $MachinePrecision]), $MachinePrecision], If[LessEqual[N[(y$95$m * x$95$m), $MachinePrecision], 1e+177], N[(N[(N[(y$95$m * x$95$m), $MachinePrecision] / z), $MachinePrecision] / N[(z * z + z), $MachinePrecision]), $MachinePrecision], N[(x$95$m / N[(z * N[(z * N[(N[(z + 1.0), $MachinePrecision] / y$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
y\_m = \left|y\right|
\\
y\_s = \mathsf{copysign}\left(1, y\right)
\\
x\_m = \left|x\right|
\\
x\_s = \mathsf{copysign}\left(1, x\right)
\\
[x_m, y_m, z] = \mathsf{sort}([x_m, y_m, z])\\
\\
x\_s \cdot \left(y\_s \cdot \begin{array}{l}
\mathbf{if}\;y\_m \cdot x\_m \leq 5 \cdot 10^{-312}:\\
\;\;\;\;\frac{\frac{x\_m}{z}}{\frac{z}{y\_m}}\\

\mathbf{elif}\;y\_m \cdot x\_m \leq 10^{+177}:\\
\;\;\;\;\frac{\frac{y\_m \cdot x\_m}{z}}{\mathsf{fma}\left(z, z, z\right)}\\

\mathbf{else}:\\
\;\;\;\;\frac{x\_m}{z \cdot \left(z \cdot \frac{z + 1}{y\_m}\right)}\\


\end{array}\right)
\end{array}
Derivation
  1. Split input into 3 regimes
  2. if (*.f64 x y) < 5.0000000000022e-312

    1. Initial program 68.2%

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

      \[\leadsto \color{blue}{\frac{x \cdot y}{{z}^{2}}} \]
    4. Step-by-step derivation
      1. associate-/l*N/A

        \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{2}}} \]
      2. lower-*.f64N/A

        \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{2}}} \]
      3. lower-/.f64N/A

        \[\leadsto x \cdot \color{blue}{\frac{y}{{z}^{2}}} \]
      4. unpow2N/A

        \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot z}} \]
      5. lower-*.f6482.5

        \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot z}} \]
    5. Applied rewrites82.5%

      \[\leadsto \color{blue}{x \cdot \frac{y}{z \cdot z}} \]
    6. Step-by-step derivation
      1. Applied rewrites99.8%

        \[\leadsto \frac{\frac{x}{z}}{\color{blue}{\frac{z}{y}}} \]

      if 5.0000000000022e-312 < (*.f64 x y) < 1e177

      1. Initial program 92.5%

        \[\frac{x \cdot y}{\left(z \cdot z\right) \cdot \left(z + 1\right)} \]
      2. Add Preprocessing
      3. Step-by-step derivation
        1. lift-/.f64N/A

          \[\leadsto \color{blue}{\frac{x \cdot y}{\left(z \cdot z\right) \cdot \left(z + 1\right)}} \]
        2. lift-*.f64N/A

          \[\leadsto \frac{x \cdot y}{\color{blue}{\left(z \cdot z\right) \cdot \left(z + 1\right)}} \]
        3. lift-*.f64N/A

          \[\leadsto \frac{x \cdot y}{\color{blue}{\left(z \cdot z\right)} \cdot \left(z + 1\right)} \]
        4. associate-*l*N/A

          \[\leadsto \frac{x \cdot y}{\color{blue}{z \cdot \left(z \cdot \left(z + 1\right)\right)}} \]
        5. associate-/r*N/A

          \[\leadsto \color{blue}{\frac{\frac{x \cdot y}{z}}{z \cdot \left(z + 1\right)}} \]
        6. lower-/.f64N/A

          \[\leadsto \color{blue}{\frac{\frac{x \cdot y}{z}}{z \cdot \left(z + 1\right)}} \]
        7. lower-/.f64N/A

          \[\leadsto \frac{\color{blue}{\frac{x \cdot y}{z}}}{z \cdot \left(z + 1\right)} \]
        8. *-commutativeN/A

          \[\leadsto \frac{\frac{x \cdot y}{z}}{\color{blue}{\left(z + 1\right) \cdot z}} \]
        9. lift-+.f64N/A

          \[\leadsto \frac{\frac{x \cdot y}{z}}{\color{blue}{\left(z + 1\right)} \cdot z} \]
        10. distribute-lft1-inN/A

          \[\leadsto \frac{\frac{x \cdot y}{z}}{\color{blue}{z \cdot z + z}} \]
        11. lower-fma.f6499.7

          \[\leadsto \frac{\frac{x \cdot y}{z}}{\color{blue}{\mathsf{fma}\left(z, z, z\right)}} \]
      4. Applied rewrites99.7%

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

      if 1e177 < (*.f64 x y)

      1. Initial program 66.7%

        \[\frac{x \cdot y}{\left(z \cdot z\right) \cdot \left(z + 1\right)} \]
      2. Add Preprocessing
      3. Step-by-step derivation
        1. lift-/.f64N/A

          \[\leadsto \color{blue}{\frac{x \cdot y}{\left(z \cdot z\right) \cdot \left(z + 1\right)}} \]
        2. lift-*.f64N/A

          \[\leadsto \frac{\color{blue}{x \cdot y}}{\left(z \cdot z\right) \cdot \left(z + 1\right)} \]
        3. lift-*.f64N/A

          \[\leadsto \frac{x \cdot y}{\color{blue}{\left(z \cdot z\right) \cdot \left(z + 1\right)}} \]
        4. times-fracN/A

          \[\leadsto \color{blue}{\frac{x}{z \cdot z} \cdot \frac{y}{z + 1}} \]
        5. *-commutativeN/A

          \[\leadsto \color{blue}{\frac{y}{z + 1} \cdot \frac{x}{z \cdot z}} \]
        6. lift-*.f64N/A

          \[\leadsto \frac{y}{z + 1} \cdot \frac{x}{\color{blue}{z \cdot z}} \]
        7. associate-/r*N/A

          \[\leadsto \frac{y}{z + 1} \cdot \color{blue}{\frac{\frac{x}{z}}{z}} \]
        8. clear-numN/A

          \[\leadsto \color{blue}{\frac{1}{\frac{z + 1}{y}}} \cdot \frac{\frac{x}{z}}{z} \]
        9. frac-timesN/A

          \[\leadsto \color{blue}{\frac{1 \cdot \frac{x}{z}}{\frac{z + 1}{y} \cdot z}} \]
        10. metadata-evalN/A

          \[\leadsto \frac{\color{blue}{{1}^{-1}} \cdot \frac{x}{z}}{\frac{z + 1}{y} \cdot z} \]
        11. clear-numN/A

          \[\leadsto \frac{{1}^{-1} \cdot \color{blue}{\frac{1}{\frac{z}{x}}}}{\frac{z + 1}{y} \cdot z} \]
        12. inv-powN/A

          \[\leadsto \frac{{1}^{-1} \cdot \color{blue}{{\left(\frac{z}{x}\right)}^{-1}}}{\frac{z + 1}{y} \cdot z} \]
        13. unpow-prod-downN/A

          \[\leadsto \frac{\color{blue}{{\left(1 \cdot \frac{z}{x}\right)}^{-1}}}{\frac{z + 1}{y} \cdot z} \]
        14. associate-/l*N/A

          \[\leadsto \frac{{\color{blue}{\left(\frac{1 \cdot z}{x}\right)}}^{-1}}{\frac{z + 1}{y} \cdot z} \]
        15. *-lft-identityN/A

          \[\leadsto \frac{{\left(\frac{\color{blue}{z}}{x}\right)}^{-1}}{\frac{z + 1}{y} \cdot z} \]
        16. inv-powN/A

          \[\leadsto \frac{\color{blue}{\frac{1}{\frac{z}{x}}}}{\frac{z + 1}{y} \cdot z} \]
        17. clear-numN/A

          \[\leadsto \frac{\color{blue}{\frac{x}{z}}}{\frac{z + 1}{y} \cdot z} \]
        18. lower-/.f64N/A

          \[\leadsto \color{blue}{\frac{\frac{x}{z}}{\frac{z + 1}{y} \cdot z}} \]
        19. lower-/.f64N/A

          \[\leadsto \frac{\color{blue}{\frac{x}{z}}}{\frac{z + 1}{y} \cdot z} \]
        20. lower-*.f64N/A

          \[\leadsto \frac{\frac{x}{z}}{\color{blue}{\frac{z + 1}{y} \cdot z}} \]
        21. lower-/.f6499.8

          \[\leadsto \frac{\frac{x}{z}}{\color{blue}{\frac{z + 1}{y}} \cdot z} \]
      4. Applied rewrites99.8%

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

          \[\leadsto \color{blue}{\frac{\frac{x}{z}}{\frac{z + 1}{y} \cdot z}} \]
        2. lift-/.f64N/A

          \[\leadsto \frac{\color{blue}{\frac{x}{z}}}{\frac{z + 1}{y} \cdot z} \]
        3. associate-/l/N/A

          \[\leadsto \color{blue}{\frac{x}{\left(\frac{z + 1}{y} \cdot z\right) \cdot z}} \]
        4. frac-2negN/A

          \[\leadsto \color{blue}{\frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\left(\frac{z + 1}{y} \cdot z\right) \cdot z\right)}} \]
        5. lift-*.f64N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\color{blue}{\left(\frac{z + 1}{y} \cdot z\right)} \cdot z\right)} \]
        6. lift-+.f64N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\left(\frac{\color{blue}{z + 1}}{y} \cdot z\right) \cdot z\right)} \]
        7. lift-/.f64N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\left(\color{blue}{\frac{z + 1}{y}} \cdot z\right) \cdot z\right)} \]
        8. associate-*l/N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\color{blue}{\frac{\left(z + 1\right) \cdot z}{y}} \cdot z\right)} \]
        9. distribute-lft1-inN/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\frac{\color{blue}{z \cdot z + z}}{y} \cdot z\right)} \]
        10. lift-fma.f64N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\frac{\color{blue}{\mathsf{fma}\left(z, z, z\right)}}{y} \cdot z\right)} \]
        11. associate-*l/N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\color{blue}{\frac{\mathsf{fma}\left(z, z, z\right) \cdot z}{y}}\right)} \]
        12. *-commutativeN/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\frac{\color{blue}{z \cdot \mathsf{fma}\left(z, z, z\right)}}{y}\right)} \]
        13. lift-*.f64N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\frac{\color{blue}{z \cdot \mathsf{fma}\left(z, z, z\right)}}{y}\right)} \]
        14. lower-/.f64N/A

          \[\leadsto \color{blue}{\frac{\mathsf{neg}\left(x\right)}{\mathsf{neg}\left(\frac{z \cdot \mathsf{fma}\left(z, z, z\right)}{y}\right)}} \]
        15. lower-neg.f64N/A

          \[\leadsto \frac{\color{blue}{\mathsf{neg}\left(x\right)}}{\mathsf{neg}\left(\frac{z \cdot \mathsf{fma}\left(z, z, z\right)}{y}\right)} \]
        16. distribute-neg-fracN/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\color{blue}{\frac{\mathsf{neg}\left(z \cdot \mathsf{fma}\left(z, z, z\right)\right)}{y}}} \]
        17. lift-*.f64N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\frac{\mathsf{neg}\left(\color{blue}{z \cdot \mathsf{fma}\left(z, z, z\right)}\right)}{y}} \]
        18. distribute-lft-neg-inN/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\frac{\color{blue}{\left(\mathsf{neg}\left(z\right)\right) \cdot \mathsf{fma}\left(z, z, z\right)}}{y}} \]
        19. associate-*r/N/A

          \[\leadsto \frac{\mathsf{neg}\left(x\right)}{\color{blue}{\left(\mathsf{neg}\left(z\right)\right) \cdot \frac{\mathsf{fma}\left(z, z, z\right)}{y}}} \]
      6. Applied rewrites98.2%

        \[\leadsto \color{blue}{\frac{-x}{z \cdot \left(z \cdot \frac{-1 - z}{y}\right)}} \]
    7. Recombined 3 regimes into one program.
    8. Final simplification99.4%

      \[\leadsto \begin{array}{l} \mathbf{if}\;y \cdot x \leq 5 \cdot 10^{-312}:\\ \;\;\;\;\frac{\frac{x}{z}}{\frac{z}{y}}\\ \mathbf{elif}\;y \cdot x \leq 10^{+177}:\\ \;\;\;\;\frac{\frac{y \cdot x}{z}}{\mathsf{fma}\left(z, z, z\right)}\\ \mathbf{else}:\\ \;\;\;\;\frac{x}{z \cdot \left(z \cdot \frac{z + 1}{y}\right)}\\ \end{array} \]
    9. Add Preprocessing

    Alternative 2: 91.6% accurate, 0.4× speedup?

    \[\begin{array}{l} y\_m = \left|y\right| \\ y\_s = \mathsf{copysign}\left(1, y\right) \\ x\_m = \left|x\right| \\ x\_s = \mathsf{copysign}\left(1, x\right) \\ [x_m, y_m, z] = \mathsf{sort}([x_m, y_m, z])\\ \\ \begin{array}{l} t_0 := x\_m \cdot \frac{y\_m}{z \cdot \left(z \cdot z\right)}\\ t_1 := \left(z + 1\right) \cdot \left(z \cdot z\right)\\ x\_s \cdot \left(y\_s \cdot \begin{array}{l} \mathbf{if}\;t\_1 \leq -5000:\\ \;\;\;\;t\_0\\ \mathbf{elif}\;t\_1 \leq 0:\\ \;\;\;\;\frac{x\_m}{z} \cdot \frac{y\_m}{z}\\ \mathbf{elif}\;t\_1 \leq 5 \cdot 10^{-18}:\\ \;\;\;\;y\_m \cdot \frac{x\_m}{z \cdot z}\\ \mathbf{else}:\\ \;\;\;\;t\_0\\ \end{array}\right) \end{array} \end{array} \]
    y\_m = (fabs.f64 y)
    y\_s = (copysign.f64 #s(literal 1 binary64) y)
    x\_m = (fabs.f64 x)
    x\_s = (copysign.f64 #s(literal 1 binary64) x)
    NOTE: x_m, y_m, and z should be sorted in increasing order before calling this function.
    (FPCore (x_s y_s x_m y_m z)
     :precision binary64
     (let* ((t_0 (* x_m (/ y_m (* z (* z z))))) (t_1 (* (+ z 1.0) (* z z))))
       (*
        x_s
        (*
         y_s
         (if (<= t_1 -5000.0)
           t_0
           (if (<= t_1 0.0)
             (* (/ x_m z) (/ y_m z))
             (if (<= t_1 5e-18) (* y_m (/ x_m (* z z))) t_0)))))))
    y\_m = fabs(y);
    y\_s = copysign(1.0, y);
    x\_m = fabs(x);
    x\_s = copysign(1.0, x);
    assert(x_m < y_m && y_m < z);
    double code(double x_s, double y_s, double x_m, double y_m, double z) {
    	double t_0 = x_m * (y_m / (z * (z * z)));
    	double t_1 = (z + 1.0) * (z * z);
    	double tmp;
    	if (t_1 <= -5000.0) {
    		tmp = t_0;
    	} else if (t_1 <= 0.0) {
    		tmp = (x_m / z) * (y_m / z);
    	} else if (t_1 <= 5e-18) {
    		tmp = y_m * (x_m / (z * z));
    	} else {
    		tmp = t_0;
    	}
    	return x_s * (y_s * tmp);
    }
    
    y\_m = abs(y)
    y\_s = copysign(1.0d0, y)
    x\_m = abs(x)
    x\_s = copysign(1.0d0, x)
    NOTE: x_m, y_m, and z should be sorted in increasing order before calling this function.
    real(8) function code(x_s, y_s, x_m, y_m, z)
        real(8), intent (in) :: x_s
        real(8), intent (in) :: y_s
        real(8), intent (in) :: x_m
        real(8), intent (in) :: y_m
        real(8), intent (in) :: z
        real(8) :: t_0
        real(8) :: t_1
        real(8) :: tmp
        t_0 = x_m * (y_m / (z * (z * z)))
        t_1 = (z + 1.0d0) * (z * z)
        if (t_1 <= (-5000.0d0)) then
            tmp = t_0
        else if (t_1 <= 0.0d0) then
            tmp = (x_m / z) * (y_m / z)
        else if (t_1 <= 5d-18) then
            tmp = y_m * (x_m / (z * z))
        else
            tmp = t_0
        end if
        code = x_s * (y_s * tmp)
    end function
    
    y\_m = Math.abs(y);
    y\_s = Math.copySign(1.0, y);
    x\_m = Math.abs(x);
    x\_s = Math.copySign(1.0, x);
    assert x_m < y_m && y_m < z;
    public static double code(double x_s, double y_s, double x_m, double y_m, double z) {
    	double t_0 = x_m * (y_m / (z * (z * z)));
    	double t_1 = (z + 1.0) * (z * z);
    	double tmp;
    	if (t_1 <= -5000.0) {
    		tmp = t_0;
    	} else if (t_1 <= 0.0) {
    		tmp = (x_m / z) * (y_m / z);
    	} else if (t_1 <= 5e-18) {
    		tmp = y_m * (x_m / (z * z));
    	} else {
    		tmp = t_0;
    	}
    	return x_s * (y_s * tmp);
    }
    
    y\_m = math.fabs(y)
    y\_s = math.copysign(1.0, y)
    x\_m = math.fabs(x)
    x\_s = math.copysign(1.0, x)
    [x_m, y_m, z] = sort([x_m, y_m, z])
    def code(x_s, y_s, x_m, y_m, z):
    	t_0 = x_m * (y_m / (z * (z * z)))
    	t_1 = (z + 1.0) * (z * z)
    	tmp = 0
    	if t_1 <= -5000.0:
    		tmp = t_0
    	elif t_1 <= 0.0:
    		tmp = (x_m / z) * (y_m / z)
    	elif t_1 <= 5e-18:
    		tmp = y_m * (x_m / (z * z))
    	else:
    		tmp = t_0
    	return x_s * (y_s * tmp)
    
    y\_m = abs(y)
    y\_s = copysign(1.0, y)
    x\_m = abs(x)
    x\_s = copysign(1.0, x)
    x_m, y_m, z = sort([x_m, y_m, z])
    function code(x_s, y_s, x_m, y_m, z)
    	t_0 = Float64(x_m * Float64(y_m / Float64(z * Float64(z * z))))
    	t_1 = Float64(Float64(z + 1.0) * Float64(z * z))
    	tmp = 0.0
    	if (t_1 <= -5000.0)
    		tmp = t_0;
    	elseif (t_1 <= 0.0)
    		tmp = Float64(Float64(x_m / z) * Float64(y_m / z));
    	elseif (t_1 <= 5e-18)
    		tmp = Float64(y_m * Float64(x_m / Float64(z * z)));
    	else
    		tmp = t_0;
    	end
    	return Float64(x_s * Float64(y_s * tmp))
    end
    
    y\_m = abs(y);
    y\_s = sign(y) * abs(1.0);
    x\_m = abs(x);
    x\_s = sign(x) * abs(1.0);
    x_m, y_m, z = num2cell(sort([x_m, y_m, z])){:}
    function tmp_2 = code(x_s, y_s, x_m, y_m, z)
    	t_0 = x_m * (y_m / (z * (z * z)));
    	t_1 = (z + 1.0) * (z * z);
    	tmp = 0.0;
    	if (t_1 <= -5000.0)
    		tmp = t_0;
    	elseif (t_1 <= 0.0)
    		tmp = (x_m / z) * (y_m / z);
    	elseif (t_1 <= 5e-18)
    		tmp = y_m * (x_m / (z * z));
    	else
    		tmp = t_0;
    	end
    	tmp_2 = x_s * (y_s * tmp);
    end
    
    y\_m = N[Abs[y], $MachinePrecision]
    y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
    x\_m = N[Abs[x], $MachinePrecision]
    x\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[x]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
    NOTE: x_m, y_m, and z should be sorted in increasing order before calling this function.
    code[x$95$s_, y$95$s_, x$95$m_, y$95$m_, z_] := Block[{t$95$0 = N[(x$95$m * N[(y$95$m / N[(z * N[(z * z), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(z + 1.0), $MachinePrecision] * N[(z * z), $MachinePrecision]), $MachinePrecision]}, N[(x$95$s * N[(y$95$s * If[LessEqual[t$95$1, -5000.0], t$95$0, If[LessEqual[t$95$1, 0.0], N[(N[(x$95$m / z), $MachinePrecision] * N[(y$95$m / z), $MachinePrecision]), $MachinePrecision], If[LessEqual[t$95$1, 5e-18], N[(y$95$m * N[(x$95$m / N[(z * z), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], t$95$0]]]), $MachinePrecision]), $MachinePrecision]]]
    
    \begin{array}{l}
    y\_m = \left|y\right|
    \\
    y\_s = \mathsf{copysign}\left(1, y\right)
    \\
    x\_m = \left|x\right|
    \\
    x\_s = \mathsf{copysign}\left(1, x\right)
    \\
    [x_m, y_m, z] = \mathsf{sort}([x_m, y_m, z])\\
    \\
    \begin{array}{l}
    t_0 := x\_m \cdot \frac{y\_m}{z \cdot \left(z \cdot z\right)}\\
    t_1 := \left(z + 1\right) \cdot \left(z \cdot z\right)\\
    x\_s \cdot \left(y\_s \cdot \begin{array}{l}
    \mathbf{if}\;t\_1 \leq -5000:\\
    \;\;\;\;t\_0\\
    
    \mathbf{elif}\;t\_1 \leq 0:\\
    \;\;\;\;\frac{x\_m}{z} \cdot \frac{y\_m}{z}\\
    
    \mathbf{elif}\;t\_1 \leq 5 \cdot 10^{-18}:\\
    \;\;\;\;y\_m \cdot \frac{x\_m}{z \cdot z}\\
    
    \mathbf{else}:\\
    \;\;\;\;t\_0\\
    
    
    \end{array}\right)
    \end{array}
    \end{array}
    
    Derivation
    1. Split input into 3 regimes
    2. if (*.f64 (*.f64 z z) (+.f64 z #s(literal 1 binary64))) < -5e3 or 5.00000000000000036e-18 < (*.f64 (*.f64 z z) (+.f64 z #s(literal 1 binary64)))

      1. Initial program 83.6%

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

        \[\leadsto \color{blue}{\frac{x \cdot y}{{z}^{3}}} \]
      4. Step-by-step derivation
        1. associate-/l*N/A

          \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{3}}} \]
        2. lower-*.f64N/A

          \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{3}}} \]
        3. lower-/.f64N/A

          \[\leadsto x \cdot \color{blue}{\frac{y}{{z}^{3}}} \]
        4. cube-multN/A

          \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot \left(z \cdot z\right)}} \]
        5. unpow2N/A

          \[\leadsto x \cdot \frac{y}{z \cdot \color{blue}{{z}^{2}}} \]
        6. lower-*.f64N/A

          \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot {z}^{2}}} \]
        7. unpow2N/A

          \[\leadsto x \cdot \frac{y}{z \cdot \color{blue}{\left(z \cdot z\right)}} \]
        8. lower-*.f6485.7

          \[\leadsto x \cdot \frac{y}{z \cdot \color{blue}{\left(z \cdot z\right)}} \]
      5. Applied rewrites85.7%

        \[\leadsto \color{blue}{x \cdot \frac{y}{z \cdot \left(z \cdot z\right)}} \]

      if -5e3 < (*.f64 (*.f64 z z) (+.f64 z #s(literal 1 binary64))) < 0.0

      1. Initial program 72.9%

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

        \[\leadsto \color{blue}{\frac{x \cdot y}{{z}^{2}}} \]
      4. Step-by-step derivation
        1. associate-/l*N/A

          \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{2}}} \]
        2. lower-*.f64N/A

          \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{2}}} \]
        3. lower-/.f64N/A

          \[\leadsto x \cdot \color{blue}{\frac{y}{{z}^{2}}} \]
        4. unpow2N/A

          \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot z}} \]
        5. lower-*.f6472.9

          \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot z}} \]
      5. Applied rewrites72.9%

        \[\leadsto \color{blue}{x \cdot \frac{y}{z \cdot z}} \]
      6. Step-by-step derivation
        1. Applied rewrites97.8%

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

        if 0.0 < (*.f64 (*.f64 z z) (+.f64 z #s(literal 1 binary64))) < 5.00000000000000036e-18

        1. Initial program 91.6%

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

          \[\leadsto \color{blue}{\frac{x \cdot y}{{z}^{2}}} \]
        4. Step-by-step derivation
          1. associate-/l*N/A

            \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{2}}} \]
          2. lower-*.f64N/A

            \[\leadsto \color{blue}{x \cdot \frac{y}{{z}^{2}}} \]
          3. lower-/.f64N/A

            \[\leadsto x \cdot \color{blue}{\frac{y}{{z}^{2}}} \]
          4. unpow2N/A

            \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot z}} \]
          5. lower-*.f6481.2

            \[\leadsto x \cdot \frac{y}{\color{blue}{z \cdot z}} \]
        5. Applied rewrites81.2%

          \[\leadsto \color{blue}{x \cdot \frac{y}{z \cdot z}} \]
        6. Step-by-step derivation
          1. Applied rewrites97.9%

            \[\leadsto \frac{x}{z \cdot z} \cdot \color{blue}{y} \]
        7. Recombined 3 regimes into one program.
        8. Final simplification91.6%

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

        Developer Target 1: 96.8% accurate, 0.6× speedup?

        \[\begin{array}{l} \\ \begin{array}{l} \mathbf{if}\;z < 249.6182814532307:\\ \;\;\;\;\frac{y \cdot \frac{x}{z}}{z + z \cdot z}\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{\frac{y}{z}}{1 + z} \cdot x}{z}\\ \end{array} \end{array} \]
        (FPCore (x y z)
         :precision binary64
         (if (< z 249.6182814532307)
           (/ (* y (/ x z)) (+ z (* z z)))
           (/ (* (/ (/ y z) (+ 1.0 z)) x) z)))
        double code(double x, double y, double z) {
        	double tmp;
        	if (z < 249.6182814532307) {
        		tmp = (y * (x / z)) / (z + (z * z));
        	} else {
        		tmp = (((y / z) / (1.0 + z)) * x) / z;
        	}
        	return tmp;
        }
        
        real(8) function code(x, y, z)
            real(8), intent (in) :: x
            real(8), intent (in) :: y
            real(8), intent (in) :: z
            real(8) :: tmp
            if (z < 249.6182814532307d0) then
                tmp = (y * (x / z)) / (z + (z * z))
            else
                tmp = (((y / z) / (1.0d0 + z)) * x) / z
            end if
            code = tmp
        end function
        
        public static double code(double x, double y, double z) {
        	double tmp;
        	if (z < 249.6182814532307) {
        		tmp = (y * (x / z)) / (z + (z * z));
        	} else {
        		tmp = (((y / z) / (1.0 + z)) * x) / z;
        	}
        	return tmp;
        }
        
        def code(x, y, z):
        	tmp = 0
        	if z < 249.6182814532307:
        		tmp = (y * (x / z)) / (z + (z * z))
        	else:
        		tmp = (((y / z) / (1.0 + z)) * x) / z
        	return tmp
        
        function code(x, y, z)
        	tmp = 0.0
        	if (z < 249.6182814532307)
        		tmp = Float64(Float64(y * Float64(x / z)) / Float64(z + Float64(z * z)));
        	else
        		tmp = Float64(Float64(Float64(Float64(y / z) / Float64(1.0 + z)) * x) / z);
        	end
        	return tmp
        end
        
        function tmp_2 = code(x, y, z)
        	tmp = 0.0;
        	if (z < 249.6182814532307)
        		tmp = (y * (x / z)) / (z + (z * z));
        	else
        		tmp = (((y / z) / (1.0 + z)) * x) / z;
        	end
        	tmp_2 = tmp;
        end
        
        code[x_, y_, z_] := If[Less[z, 249.6182814532307], N[(N[(y * N[(x / z), $MachinePrecision]), $MachinePrecision] / N[(z + N[(z * z), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(N[(y / z), $MachinePrecision] / N[(1.0 + z), $MachinePrecision]), $MachinePrecision] * x), $MachinePrecision] / z), $MachinePrecision]]
        
        \begin{array}{l}
        
        \\
        \begin{array}{l}
        \mathbf{if}\;z < 249.6182814532307:\\
        \;\;\;\;\frac{y \cdot \frac{x}{z}}{z + z \cdot z}\\
        
        \mathbf{else}:\\
        \;\;\;\;\frac{\frac{\frac{y}{z}}{1 + z} \cdot x}{z}\\
        
        
        \end{array}
        \end{array}
        

        Reproduce

        ?
        herbie shell --seed 2024228 
        (FPCore (x y z)
          :name "Statistics.Distribution.Beta:$cvariance from math-functions-0.1.5.2"
          :precision binary64
        
          :alt
          (! :herbie-platform default (if (< z 2496182814532307/10000000000000) (/ (* y (/ x z)) (+ z (* z z))) (/ (* (/ (/ y z) (+ 1 z)) x) z)))
        
          (/ (* x y) (* (* z z) (+ z 1.0))))