
(FPCore (a b eps) :precision binary64 (/ (* eps (- (exp (* (+ a b) eps)) 1.0)) (* (- (exp (* a eps)) 1.0) (- (exp (* b eps)) 1.0))))
double code(double a, double b, double eps) {
return (eps * (exp(((a + b) * eps)) - 1.0)) / ((exp((a * eps)) - 1.0) * (exp((b * eps)) - 1.0));
}
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
code = (eps * (exp(((a + b) * eps)) - 1.0d0)) / ((exp((a * eps)) - 1.0d0) * (exp((b * eps)) - 1.0d0))
end function
public static double code(double a, double b, double eps) {
return (eps * (Math.exp(((a + b) * eps)) - 1.0)) / ((Math.exp((a * eps)) - 1.0) * (Math.exp((b * eps)) - 1.0));
}
def code(a, b, eps): return (eps * (math.exp(((a + b) * eps)) - 1.0)) / ((math.exp((a * eps)) - 1.0) * (math.exp((b * eps)) - 1.0))
function code(a, b, eps) return Float64(Float64(eps * Float64(exp(Float64(Float64(a + b) * eps)) - 1.0)) / Float64(Float64(exp(Float64(a * eps)) - 1.0) * Float64(exp(Float64(b * eps)) - 1.0))) end
function tmp = code(a, b, eps) tmp = (eps * (exp(((a + b) * eps)) - 1.0)) / ((exp((a * eps)) - 1.0) * (exp((b * eps)) - 1.0)); end
code[a_, b_, eps_] := N[(N[(eps * N[(N[Exp[N[(N[(a + b), $MachinePrecision] * eps), $MachinePrecision]], $MachinePrecision] - 1.0), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Exp[N[(a * eps), $MachinePrecision]], $MachinePrecision] - 1.0), $MachinePrecision] * N[(N[Exp[N[(b * eps), $MachinePrecision]], $MachinePrecision] - 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\varepsilon \cdot \left(e^{\left(a + b\right) \cdot \varepsilon} - 1\right)}{\left(e^{a \cdot \varepsilon} - 1\right) \cdot \left(e^{b \cdot \varepsilon} - 1\right)}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 6 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b eps) :precision binary64 (/ (* eps (- (exp (* (+ a b) eps)) 1.0)) (* (- (exp (* a eps)) 1.0) (- (exp (* b eps)) 1.0))))
double code(double a, double b, double eps) {
return (eps * (exp(((a + b) * eps)) - 1.0)) / ((exp((a * eps)) - 1.0) * (exp((b * eps)) - 1.0));
}
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
code = (eps * (exp(((a + b) * eps)) - 1.0d0)) / ((exp((a * eps)) - 1.0d0) * (exp((b * eps)) - 1.0d0))
end function
public static double code(double a, double b, double eps) {
return (eps * (Math.exp(((a + b) * eps)) - 1.0)) / ((Math.exp((a * eps)) - 1.0) * (Math.exp((b * eps)) - 1.0));
}
def code(a, b, eps): return (eps * (math.exp(((a + b) * eps)) - 1.0)) / ((math.exp((a * eps)) - 1.0) * (math.exp((b * eps)) - 1.0))
function code(a, b, eps) return Float64(Float64(eps * Float64(exp(Float64(Float64(a + b) * eps)) - 1.0)) / Float64(Float64(exp(Float64(a * eps)) - 1.0) * Float64(exp(Float64(b * eps)) - 1.0))) end
function tmp = code(a, b, eps) tmp = (eps * (exp(((a + b) * eps)) - 1.0)) / ((exp((a * eps)) - 1.0) * (exp((b * eps)) - 1.0)); end
code[a_, b_, eps_] := N[(N[(eps * N[(N[Exp[N[(N[(a + b), $MachinePrecision] * eps), $MachinePrecision]], $MachinePrecision] - 1.0), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Exp[N[(a * eps), $MachinePrecision]], $MachinePrecision] - 1.0), $MachinePrecision] * N[(N[Exp[N[(b * eps), $MachinePrecision]], $MachinePrecision] - 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\varepsilon \cdot \left(e^{\left(a + b\right) \cdot \varepsilon} - 1\right)}{\left(e^{a \cdot \varepsilon} - 1\right) \cdot \left(e^{b \cdot \varepsilon} - 1\right)}
\end{array}
NOTE: a and b should be sorted in increasing order before calling this function. (FPCore (a b eps) :precision binary64 (if (<= a 1.35e+121) (+ (/ 1.0 b) (/ 1.0 a)) (/ eps (expm1 (* a eps)))))
assert(a < b);
double code(double a, double b, double eps) {
double tmp;
if (a <= 1.35e+121) {
tmp = (1.0 / b) + (1.0 / a);
} else {
tmp = eps / expm1((a * eps));
}
return tmp;
}
assert a < b;
public static double code(double a, double b, double eps) {
double tmp;
if (a <= 1.35e+121) {
tmp = (1.0 / b) + (1.0 / a);
} else {
tmp = eps / Math.expm1((a * eps));
}
return tmp;
}
[a, b] = sort([a, b]) def code(a, b, eps): tmp = 0 if a <= 1.35e+121: tmp = (1.0 / b) + (1.0 / a) else: tmp = eps / math.expm1((a * eps)) return tmp
a, b = sort([a, b]) function code(a, b, eps) tmp = 0.0 if (a <= 1.35e+121) tmp = Float64(Float64(1.0 / b) + Float64(1.0 / a)); else tmp = Float64(eps / expm1(Float64(a * eps))); end return tmp end
NOTE: a and b should be sorted in increasing order before calling this function. code[a_, b_, eps_] := If[LessEqual[a, 1.35e+121], N[(N[(1.0 / b), $MachinePrecision] + N[(1.0 / a), $MachinePrecision]), $MachinePrecision], N[(eps / N[(Exp[N[(a * eps), $MachinePrecision]] - 1), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\begin{array}{l}
\mathbf{if}\;a \leq 1.35 \cdot 10^{+121}:\\
\;\;\;\;\frac{1}{b} + \frac{1}{a}\\
\mathbf{else}:\\
\;\;\;\;\frac{\varepsilon}{\mathsf{expm1}\left(a \cdot \varepsilon\right)}\\
\end{array}
\end{array}
if a < 1.3500000000000001e121Initial program 1.1%
associate-*l/1.1%
*-commutative1.1%
expm1-def2.9%
*-commutative2.9%
expm1-def13.0%
*-commutative13.0%
expm1-def33.8%
*-commutative33.8%
Simplified33.8%
Taylor expanded in eps around 0 79.8%
Taylor expanded in a around 0 99.5%
if 1.3500000000000001e121 < a Initial program 30.2%
times-frac30.2%
expm1-def34.4%
*-commutative34.4%
expm1-def35.2%
*-commutative35.2%
expm1-def71.0%
*-commutative71.0%
Simplified71.0%
Taylor expanded in eps around 0 33.4%
Taylor expanded in b around inf 22.6%
expm1-def26.5%
*-commutative26.5%
Simplified26.5%
Final simplification92.1%
NOTE: a and b should be sorted in increasing order before calling this function. (FPCore (a b eps) :precision binary64 (if (<= b 2.3e-46) (/ 1.0 b) (+ (/ 1.0 a) (* eps -0.5))))
assert(a < b);
double code(double a, double b, double eps) {
double tmp;
if (b <= 2.3e-46) {
tmp = 1.0 / b;
} else {
tmp = (1.0 / a) + (eps * -0.5);
}
return tmp;
}
NOTE: a and b should be sorted in increasing order before calling this function.
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
real(8) :: tmp
if (b <= 2.3d-46) then
tmp = 1.0d0 / b
else
tmp = (1.0d0 / a) + (eps * (-0.5d0))
end if
code = tmp
end function
assert a < b;
public static double code(double a, double b, double eps) {
double tmp;
if (b <= 2.3e-46) {
tmp = 1.0 / b;
} else {
tmp = (1.0 / a) + (eps * -0.5);
}
return tmp;
}
[a, b] = sort([a, b]) def code(a, b, eps): tmp = 0 if b <= 2.3e-46: tmp = 1.0 / b else: tmp = (1.0 / a) + (eps * -0.5) return tmp
a, b = sort([a, b]) function code(a, b, eps) tmp = 0.0 if (b <= 2.3e-46) tmp = Float64(1.0 / b); else tmp = Float64(Float64(1.0 / a) + Float64(eps * -0.5)); end return tmp end
a, b = num2cell(sort([a, b])){:}
function tmp_2 = code(a, b, eps)
tmp = 0.0;
if (b <= 2.3e-46)
tmp = 1.0 / b;
else
tmp = (1.0 / a) + (eps * -0.5);
end
tmp_2 = tmp;
end
NOTE: a and b should be sorted in increasing order before calling this function. code[a_, b_, eps_] := If[LessEqual[b, 2.3e-46], N[(1.0 / b), $MachinePrecision], N[(N[(1.0 / a), $MachinePrecision] + N[(eps * -0.5), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\begin{array}{l}
\mathbf{if}\;b \leq 2.3 \cdot 10^{-46}:\\
\;\;\;\;\frac{1}{b}\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{a} + \varepsilon \cdot -0.5\\
\end{array}
\end{array}
if b < 2.2999999999999999e-46Initial program 2.1%
associate-*l/2.1%
*-commutative2.1%
expm1-def3.9%
*-commutative3.9%
expm1-def11.1%
*-commutative11.1%
expm1-def32.4%
*-commutative32.4%
Simplified32.4%
Taylor expanded in b around 0 53.8%
if 2.2999999999999999e-46 < b Initial program 10.5%
times-frac10.5%
expm1-def27.0%
*-commutative27.0%
expm1-def28.8%
*-commutative28.8%
expm1-def68.2%
*-commutative68.2%
Simplified68.2%
Taylor expanded in eps around 0 58.2%
Taylor expanded in b around inf 10.6%
expm1-def47.5%
*-commutative47.5%
Simplified47.5%
Taylor expanded in eps around 0 72.5%
Final simplification58.2%
NOTE: a and b should be sorted in increasing order before calling this function. (FPCore (a b eps) :precision binary64 (+ (/ 1.0 b) (/ 1.0 a)))
assert(a < b);
double code(double a, double b, double eps) {
return (1.0 / b) + (1.0 / a);
}
NOTE: a and b should be sorted in increasing order before calling this function.
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
code = (1.0d0 / b) + (1.0d0 / a)
end function
assert a < b;
public static double code(double a, double b, double eps) {
return (1.0 / b) + (1.0 / a);
}
[a, b] = sort([a, b]) def code(a, b, eps): return (1.0 / b) + (1.0 / a)
a, b = sort([a, b]) function code(a, b, eps) return Float64(Float64(1.0 / b) + Float64(1.0 / a)) end
a, b = num2cell(sort([a, b])){:}
function tmp = code(a, b, eps)
tmp = (1.0 / b) + (1.0 / a);
end
NOTE: a and b should be sorted in increasing order before calling this function. code[a_, b_, eps_] := N[(N[(1.0 / b), $MachinePrecision] + N[(1.0 / a), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\frac{1}{b} + \frac{1}{a}
\end{array}
Initial program 4.1%
associate-*l/4.1%
*-commutative4.1%
expm1-def5.8%
*-commutative5.8%
expm1-def15.2%
*-commutative15.2%
expm1-def37.8%
*-commutative37.8%
Simplified37.8%
Taylor expanded in eps around 0 78.4%
Taylor expanded in a around 0 97.0%
Final simplification97.0%
NOTE: a and b should be sorted in increasing order before calling this function. (FPCore (a b eps) :precision binary64 (if (<= b 1.4e-44) (/ 1.0 b) (/ 1.0 a)))
assert(a < b);
double code(double a, double b, double eps) {
double tmp;
if (b <= 1.4e-44) {
tmp = 1.0 / b;
} else {
tmp = 1.0 / a;
}
return tmp;
}
NOTE: a and b should be sorted in increasing order before calling this function.
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
real(8) :: tmp
if (b <= 1.4d-44) then
tmp = 1.0d0 / b
else
tmp = 1.0d0 / a
end if
code = tmp
end function
assert a < b;
public static double code(double a, double b, double eps) {
double tmp;
if (b <= 1.4e-44) {
tmp = 1.0 / b;
} else {
tmp = 1.0 / a;
}
return tmp;
}
[a, b] = sort([a, b]) def code(a, b, eps): tmp = 0 if b <= 1.4e-44: tmp = 1.0 / b else: tmp = 1.0 / a return tmp
a, b = sort([a, b]) function code(a, b, eps) tmp = 0.0 if (b <= 1.4e-44) tmp = Float64(1.0 / b); else tmp = Float64(1.0 / a); end return tmp end
a, b = num2cell(sort([a, b])){:}
function tmp_2 = code(a, b, eps)
tmp = 0.0;
if (b <= 1.4e-44)
tmp = 1.0 / b;
else
tmp = 1.0 / a;
end
tmp_2 = tmp;
end
NOTE: a and b should be sorted in increasing order before calling this function. code[a_, b_, eps_] := If[LessEqual[b, 1.4e-44], N[(1.0 / b), $MachinePrecision], N[(1.0 / a), $MachinePrecision]]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\begin{array}{l}
\mathbf{if}\;b \leq 1.4 \cdot 10^{-44}:\\
\;\;\;\;\frac{1}{b}\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{a}\\
\end{array}
\end{array}
if b < 1.4e-44Initial program 2.1%
associate-*l/2.1%
*-commutative2.1%
expm1-def3.9%
*-commutative3.9%
expm1-def11.0%
*-commutative11.0%
expm1-def32.6%
*-commutative32.6%
Simplified32.6%
Taylor expanded in b around 0 54.3%
if 1.4e-44 < b Initial program 10.8%
associate-*l/10.8%
*-commutative10.8%
expm1-def12.3%
*-commutative12.3%
expm1-def29.7%
*-commutative29.7%
expm1-def55.7%
*-commutative55.7%
Simplified55.7%
Taylor expanded in a around 0 73.7%
Final simplification58.7%
NOTE: a and b should be sorted in increasing order before calling this function. (FPCore (a b eps) :precision binary64 (* eps -0.5))
assert(a < b);
double code(double a, double b, double eps) {
return eps * -0.5;
}
NOTE: a and b should be sorted in increasing order before calling this function.
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
code = eps * (-0.5d0)
end function
assert a < b;
public static double code(double a, double b, double eps) {
return eps * -0.5;
}
[a, b] = sort([a, b]) def code(a, b, eps): return eps * -0.5
a, b = sort([a, b]) function code(a, b, eps) return Float64(eps * -0.5) end
a, b = num2cell(sort([a, b])){:}
function tmp = code(a, b, eps)
tmp = eps * -0.5;
end
NOTE: a and b should be sorted in increasing order before calling this function. code[a_, b_, eps_] := N[(eps * -0.5), $MachinePrecision]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\varepsilon \cdot -0.5
\end{array}
Initial program 4.1%
times-frac4.1%
expm1-def13.7%
*-commutative13.7%
expm1-def15.2%
*-commutative15.2%
expm1-def53.0%
*-commutative53.0%
Simplified53.0%
Taylor expanded in eps around 0 55.3%
Taylor expanded in b around inf 4.7%
expm1-def30.2%
*-commutative30.2%
Simplified30.2%
Taylor expanded in eps around 0 51.3%
Taylor expanded in eps around inf 3.0%
Final simplification3.0%
NOTE: a and b should be sorted in increasing order before calling this function. (FPCore (a b eps) :precision binary64 (/ 1.0 a))
assert(a < b);
double code(double a, double b, double eps) {
return 1.0 / a;
}
NOTE: a and b should be sorted in increasing order before calling this function.
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
code = 1.0d0 / a
end function
assert a < b;
public static double code(double a, double b, double eps) {
return 1.0 / a;
}
[a, b] = sort([a, b]) def code(a, b, eps): return 1.0 / a
a, b = sort([a, b]) function code(a, b, eps) return Float64(1.0 / a) end
a, b = num2cell(sort([a, b])){:}
function tmp = code(a, b, eps)
tmp = 1.0 / a;
end
NOTE: a and b should be sorted in increasing order before calling this function. code[a_, b_, eps_] := N[(1.0 / a), $MachinePrecision]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\frac{1}{a}
\end{array}
Initial program 4.1%
associate-*l/4.1%
*-commutative4.1%
expm1-def5.8%
*-commutative5.8%
expm1-def15.2%
*-commutative15.2%
expm1-def37.8%
*-commutative37.8%
Simplified37.8%
Taylor expanded in a around 0 51.0%
Final simplification51.0%
(FPCore (a b eps) :precision binary64 (/ (+ a b) (* a b)))
double code(double a, double b, double eps) {
return (a + b) / (a * b);
}
real(8) function code(a, b, eps)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: eps
code = (a + b) / (a * b)
end function
public static double code(double a, double b, double eps) {
return (a + b) / (a * b);
}
def code(a, b, eps): return (a + b) / (a * b)
function code(a, b, eps) return Float64(Float64(a + b) / Float64(a * b)) end
function tmp = code(a, b, eps) tmp = (a + b) / (a * b); end
code[a_, b_, eps_] := N[(N[(a + b), $MachinePrecision] / N[(a * b), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{a + b}{a \cdot b}
\end{array}
herbie shell --seed 2023256
(FPCore (a b eps)
:name "expq3 (problem 3.4.2)"
:precision binary64
:pre (and (< -1.0 eps) (< eps 1.0))
:herbie-target
(/ (+ a b) (* a b))
(/ (* eps (- (exp (* (+ a b) eps)) 1.0)) (* (- (exp (* a eps)) 1.0) (- (exp (* b eps)) 1.0))))