
(FPCore (x y z t a b) :precision binary64 (/ (+ x (/ (* y z) t)) (+ (+ a 1.0) (/ (* y b) t))))
double code(double x, double y, double z, double t, double a, double b) {
return (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t));
}
real(8) function code(x, y, z, t, a, b)
real(8), intent (in) :: x
real(8), intent (in) :: y
real(8), intent (in) :: z
real(8), intent (in) :: t
real(8), intent (in) :: a
real(8), intent (in) :: b
code = (x + ((y * z) / t)) / ((a + 1.0d0) + ((y * b) / t))
end function
public static double code(double x, double y, double z, double t, double a, double b) {
return (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t));
}
def code(x, y, z, t, a, b): return (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t))
function code(x, y, z, t, a, b) return Float64(Float64(x + Float64(Float64(y * z) / t)) / Float64(Float64(a + 1.0) + Float64(Float64(y * b) / t))) end
function tmp = code(x, y, z, t, a, b) tmp = (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t)); end
code[x_, y_, z_, t_, a_, b_] := N[(N[(x + N[(N[(y * z), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision] / N[(N[(a + 1.0), $MachinePrecision] + N[(N[(y * b), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{x + \frac{y \cdot z}{t}}{\left(a + 1\right) + \frac{y \cdot b}{t}}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 15 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (x y z t a b) :precision binary64 (/ (+ x (/ (* y z) t)) (+ (+ a 1.0) (/ (* y b) t))))
double code(double x, double y, double z, double t, double a, double b) {
return (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t));
}
real(8) function code(x, y, z, t, a, b)
real(8), intent (in) :: x
real(8), intent (in) :: y
real(8), intent (in) :: z
real(8), intent (in) :: t
real(8), intent (in) :: a
real(8), intent (in) :: b
code = (x + ((y * z) / t)) / ((a + 1.0d0) + ((y * b) / t))
end function
public static double code(double x, double y, double z, double t, double a, double b) {
return (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t));
}
def code(x, y, z, t, a, b): return (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t))
function code(x, y, z, t, a, b) return Float64(Float64(x + Float64(Float64(y * z) / t)) / Float64(Float64(a + 1.0) + Float64(Float64(y * b) / t))) end
function tmp = code(x, y, z, t, a, b) tmp = (x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t)); end
code[x_, y_, z_, t_, a_, b_] := N[(N[(x + N[(N[(y * z), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision] / N[(N[(a + 1.0), $MachinePrecision] + N[(N[(y * b), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{x + \frac{y \cdot z}{t}}{\left(a + 1\right) + \frac{y \cdot b}{t}}
\end{array}
(FPCore (x y z t a b) :precision binary64 (if (<= (/ (+ x (/ (* y z) t)) (+ (+ a 1.0) (/ (* y b) t))) INFINITY) (fma z (/ y (fma (* y b) 1.0 (fma t a t))) (/ x (fma y (/ b t) (+ a 1.0)))) (/ (fma t (/ x y) z) b)))
double code(double x, double y, double z, double t, double a, double b) {
double tmp;
if (((x + ((y * z) / t)) / ((a + 1.0) + ((y * b) / t))) <= ((double) INFINITY)) {
tmp = fma(z, (y / fma((y * b), 1.0, fma(t, a, t))), (x / fma(y, (b / t), (a + 1.0))));
} else {
tmp = fma(t, (x / y), z) / b;
}
return tmp;
}
function code(x, y, z, t, a, b) tmp = 0.0 if (Float64(Float64(x + Float64(Float64(y * z) / t)) / Float64(Float64(a + 1.0) + Float64(Float64(y * b) / t))) <= Inf) tmp = fma(z, Float64(y / fma(Float64(y * b), 1.0, fma(t, a, t))), Float64(x / fma(y, Float64(b / t), Float64(a + 1.0)))); else tmp = Float64(fma(t, Float64(x / y), z) / b); end return tmp end
code[x_, y_, z_, t_, a_, b_] := If[LessEqual[N[(N[(x + N[(N[(y * z), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision] / N[(N[(a + 1.0), $MachinePrecision] + N[(N[(y * b), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], Infinity], N[(z * N[(y / N[(N[(y * b), $MachinePrecision] * 1.0 + N[(t * a + t), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(x / N[(y * N[(b / t), $MachinePrecision] + N[(a + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(t * N[(x / y), $MachinePrecision] + z), $MachinePrecision] / b), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;\frac{x + \frac{y \cdot z}{t}}{\left(a + 1\right) + \frac{y \cdot b}{t}} \leq \infty:\\
\;\;\;\;\mathsf{fma}\left(z, \frac{y}{\mathsf{fma}\left(y \cdot b, 1, \mathsf{fma}\left(t, a, t\right)\right)}, \frac{x}{\mathsf{fma}\left(y, \frac{b}{t}, a + 1\right)}\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{\mathsf{fma}\left(t, \frac{x}{y}, z\right)}{b}\\
\end{array}
\end{array}
if (/.f64 (+.f64 x (/.f64 (*.f64 y z) t)) (+.f64 (+.f64 a #s(literal 1 binary64)) (/.f64 (*.f64 y b) t))) < +inf.0Initial program 82.6%
Taylor expanded in x around 0
+-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
lower-/.f64N/A
lower-+.f64N/A
Applied rewrites87.7%
Applied rewrites94.1%
if +inf.0 < (/.f64 (+.f64 x (/.f64 (*.f64 y z) t)) (+.f64 (+.f64 a #s(literal 1 binary64)) (/.f64 (*.f64 y b) t))) Initial program 0.0%
Taylor expanded in x around 0
+-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
lower-/.f64N/A
lower-+.f64N/A
Applied rewrites9.6%
Taylor expanded in b around inf
Applied rewrites99.9%
(FPCore (x y z t a b)
:precision binary64
(let* ((t_1 (+ x (/ (* y z) t))) (t_2 (/ t_1 (+ (+ a 1.0) (/ (* y b) t)))))
(if (<= t_2 (- INFINITY))
(fma z (/ y (fma (* y b) 1.0 (fma t a t))) (/ x a))
(if (<= t_2 5e+303)
(/ t_1 (fma b (/ y t) (+ a 1.0)))
(if (<= t_2 INFINITY)
(fma (/ z (fma y (/ b t) (+ a 1.0))) (/ y t) (/ x a))
(/ (fma t (/ x y) z) b))))))
double code(double x, double y, double z, double t, double a, double b) {
double t_1 = x + ((y * z) / t);
double t_2 = t_1 / ((a + 1.0) + ((y * b) / t));
double tmp;
if (t_2 <= -((double) INFINITY)) {
tmp = fma(z, (y / fma((y * b), 1.0, fma(t, a, t))), (x / a));
} else if (t_2 <= 5e+303) {
tmp = t_1 / fma(b, (y / t), (a + 1.0));
} else if (t_2 <= ((double) INFINITY)) {
tmp = fma((z / fma(y, (b / t), (a + 1.0))), (y / t), (x / a));
} else {
tmp = fma(t, (x / y), z) / b;
}
return tmp;
}
function code(x, y, z, t, a, b) t_1 = Float64(x + Float64(Float64(y * z) / t)) t_2 = Float64(t_1 / Float64(Float64(a + 1.0) + Float64(Float64(y * b) / t))) tmp = 0.0 if (t_2 <= Float64(-Inf)) tmp = fma(z, Float64(y / fma(Float64(y * b), 1.0, fma(t, a, t))), Float64(x / a)); elseif (t_2 <= 5e+303) tmp = Float64(t_1 / fma(b, Float64(y / t), Float64(a + 1.0))); elseif (t_2 <= Inf) tmp = fma(Float64(z / fma(y, Float64(b / t), Float64(a + 1.0))), Float64(y / t), Float64(x / a)); else tmp = Float64(fma(t, Float64(x / y), z) / b); end return tmp end
code[x_, y_, z_, t_, a_, b_] := Block[{t$95$1 = N[(x + N[(N[(y * z), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(t$95$1 / N[(N[(a + 1.0), $MachinePrecision] + N[(N[(y * b), $MachinePrecision] / t), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$2, (-Infinity)], N[(z * N[(y / N[(N[(y * b), $MachinePrecision] * 1.0 + N[(t * a + t), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(x / a), $MachinePrecision]), $MachinePrecision], If[LessEqual[t$95$2, 5e+303], N[(t$95$1 / N[(b * N[(y / t), $MachinePrecision] + N[(a + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[t$95$2, Infinity], N[(N[(z / N[(y * N[(b / t), $MachinePrecision] + N[(a + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(y / t), $MachinePrecision] + N[(x / a), $MachinePrecision]), $MachinePrecision], N[(N[(t * N[(x / y), $MachinePrecision] + z), $MachinePrecision] / b), $MachinePrecision]]]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_1 := x + \frac{y \cdot z}{t}\\
t_2 := \frac{t\_1}{\left(a + 1\right) + \frac{y \cdot b}{t}}\\
\mathbf{if}\;t\_2 \leq -\infty:\\
\;\;\;\;\mathsf{fma}\left(z, \frac{y}{\mathsf{fma}\left(y \cdot b, 1, \mathsf{fma}\left(t, a, t\right)\right)}, \frac{x}{a}\right)\\
\mathbf{elif}\;t\_2 \leq 5 \cdot 10^{+303}:\\
\;\;\;\;\frac{t\_1}{\mathsf{fma}\left(b, \frac{y}{t}, a + 1\right)}\\
\mathbf{elif}\;t\_2 \leq \infty:\\
\;\;\;\;\mathsf{fma}\left(\frac{z}{\mathsf{fma}\left(y, \frac{b}{t}, a + 1\right)}, \frac{y}{t}, \frac{x}{a}\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{\mathsf{fma}\left(t, \frac{x}{y}, z\right)}{b}\\
\end{array}
\end{array}
if (/.f64 (+.f64 x (/.f64 (*.f64 y z) t)) (+.f64 (+.f64 a #s(literal 1 binary64)) (/.f64 (*.f64 y b) t))) < -inf.0Initial program 33.2%
Taylor expanded in x around 0
+-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
lower-/.f64N/A
lower-+.f64N/A
Applied rewrites86.7%
Applied rewrites90.2%
Taylor expanded in a around inf
Applied rewrites76.6%
if -inf.0 < (/.f64 (+.f64 x (/.f64 (*.f64 y z) t)) (+.f64 (+.f64 a #s(literal 1 binary64)) (/.f64 (*.f64 y b) t))) < 4.9999999999999997e303Initial program 90.6%
lift-+.f64N/A
+-commutativeN/A
lift-/.f64N/A
lift-*.f64N/A
*-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f6491.9
Applied rewrites91.9%
if 4.9999999999999997e303 < (/.f64 (+.f64 x (/.f64 (*.f64 y z) t)) (+.f64 (+.f64 a #s(literal 1 binary64)) (/.f64 (*.f64 y b) t))) < +inf.0Initial program 33.6%
Taylor expanded in x around 0
+-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
lower-/.f64N/A
lower-+.f64N/A
Applied rewrites88.7%
Applied rewrites85.2%
Taylor expanded in a around inf
Applied rewrites65.0%
if +inf.0 < (/.f64 (+.f64 x (/.f64 (*.f64 y z) t)) (+.f64 (+.f64 a #s(literal 1 binary64)) (/.f64 (*.f64 y b) t))) Initial program 0.0%
Taylor expanded in x around 0
+-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
lower-/.f64N/A
lower-+.f64N/A
Applied rewrites8.4%
Taylor expanded in b around inf
Applied rewrites97.5%
(FPCore (x y z t a b)
:precision binary64
(let* ((t_1
(* 1.0 (* (+ x (* (/ y t) z)) (/ 1.0 (+ (+ a 1.0) (* (/ y t) b)))))))
(if (< t -1.3659085366310088e-271)
t_1
(if (< t 3.036967103737246e-130) (/ z b) t_1))))
double code(double x, double y, double z, double t, double a, double b) {
double t_1 = 1.0 * ((x + ((y / t) * z)) * (1.0 / ((a + 1.0) + ((y / t) * b))));
double tmp;
if (t < -1.3659085366310088e-271) {
tmp = t_1;
} else if (t < 3.036967103737246e-130) {
tmp = z / b;
} else {
tmp = t_1;
}
return tmp;
}
real(8) function code(x, y, z, t, a, b)
real(8), intent (in) :: x
real(8), intent (in) :: y
real(8), intent (in) :: z
real(8), intent (in) :: t
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8) :: t_1
real(8) :: tmp
t_1 = 1.0d0 * ((x + ((y / t) * z)) * (1.0d0 / ((a + 1.0d0) + ((y / t) * b))))
if (t < (-1.3659085366310088d-271)) then
tmp = t_1
else if (t < 3.036967103737246d-130) then
tmp = z / b
else
tmp = t_1
end if
code = tmp
end function
public static double code(double x, double y, double z, double t, double a, double b) {
double t_1 = 1.0 * ((x + ((y / t) * z)) * (1.0 / ((a + 1.0) + ((y / t) * b))));
double tmp;
if (t < -1.3659085366310088e-271) {
tmp = t_1;
} else if (t < 3.036967103737246e-130) {
tmp = z / b;
} else {
tmp = t_1;
}
return tmp;
}
def code(x, y, z, t, a, b): t_1 = 1.0 * ((x + ((y / t) * z)) * (1.0 / ((a + 1.0) + ((y / t) * b)))) tmp = 0 if t < -1.3659085366310088e-271: tmp = t_1 elif t < 3.036967103737246e-130: tmp = z / b else: tmp = t_1 return tmp
function code(x, y, z, t, a, b) t_1 = Float64(1.0 * Float64(Float64(x + Float64(Float64(y / t) * z)) * Float64(1.0 / Float64(Float64(a + 1.0) + Float64(Float64(y / t) * b))))) tmp = 0.0 if (t < -1.3659085366310088e-271) tmp = t_1; elseif (t < 3.036967103737246e-130) tmp = Float64(z / b); else tmp = t_1; end return tmp end
function tmp_2 = code(x, y, z, t, a, b) t_1 = 1.0 * ((x + ((y / t) * z)) * (1.0 / ((a + 1.0) + ((y / t) * b)))); tmp = 0.0; if (t < -1.3659085366310088e-271) tmp = t_1; elseif (t < 3.036967103737246e-130) tmp = z / b; else tmp = t_1; end tmp_2 = tmp; end
code[x_, y_, z_, t_, a_, b_] := Block[{t$95$1 = N[(1.0 * N[(N[(x + N[(N[(y / t), $MachinePrecision] * z), $MachinePrecision]), $MachinePrecision] * N[(1.0 / N[(N[(a + 1.0), $MachinePrecision] + N[(N[(y / t), $MachinePrecision] * b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[Less[t, -1.3659085366310088e-271], t$95$1, If[Less[t, 3.036967103737246e-130], N[(z / b), $MachinePrecision], t$95$1]]]
\begin{array}{l}
\\
\begin{array}{l}
t_1 := 1 \cdot \left(\left(x + \frac{y}{t} \cdot z\right) \cdot \frac{1}{\left(a + 1\right) + \frac{y}{t} \cdot b}\right)\\
\mathbf{if}\;t < -1.3659085366310088 \cdot 10^{-271}:\\
\;\;\;\;t\_1\\
\mathbf{elif}\;t < 3.036967103737246 \cdot 10^{-130}:\\
\;\;\;\;\frac{z}{b}\\
\mathbf{else}:\\
\;\;\;\;t\_1\\
\end{array}
\end{array}
herbie shell --seed 2024228
(FPCore (x y z t a b)
:name "Diagrams.Solve.Tridiagonal:solveCyclicTriDiagonal from diagrams-solve-0.1, B"
:precision binary64
:alt
(! :herbie-platform default (if (< t -1707385670788761/12500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000) (* 1 (* (+ x (* (/ y t) z)) (/ 1 (+ (+ a 1) (* (/ y t) b))))) (if (< t 1518483551868623/5000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000) (/ z b) (* 1 (* (+ x (* (/ y t) z)) (/ 1 (+ (+ a 1) (* (/ y t) b))))))))
(/ (+ x (/ (* y z) t)) (+ (+ a 1.0) (/ (* y b) t))))