
(FPCore (wj x) :precision binary64 (let* ((t_0 (* wj (exp wj)))) (- wj (/ (- t_0 x) (+ (exp wj) t_0)))))
double code(double wj, double x) {
double t_0 = wj * exp(wj);
return wj - ((t_0 - x) / (exp(wj) + t_0));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
real(8) :: t_0
t_0 = wj * exp(wj)
code = wj - ((t_0 - x) / (exp(wj) + t_0))
end function
public static double code(double wj, double x) {
double t_0 = wj * Math.exp(wj);
return wj - ((t_0 - x) / (Math.exp(wj) + t_0));
}
def code(wj, x): t_0 = wj * math.exp(wj) return wj - ((t_0 - x) / (math.exp(wj) + t_0))
function code(wj, x) t_0 = Float64(wj * exp(wj)) return Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0))) end
function tmp = code(wj, x) t_0 = wj * exp(wj); tmp = wj - ((t_0 - x) / (exp(wj) + t_0)); end
code[wj_, x_] := Block[{t$95$0 = N[(wj * N[Exp[wj], $MachinePrecision]), $MachinePrecision]}, N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := wj \cdot e^{wj}\\
wj - \frac{t_0 - x}{e^{wj} + t_0}
\end{array}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 7 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (wj x) :precision binary64 (let* ((t_0 (* wj (exp wj)))) (- wj (/ (- t_0 x) (+ (exp wj) t_0)))))
double code(double wj, double x) {
double t_0 = wj * exp(wj);
return wj - ((t_0 - x) / (exp(wj) + t_0));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
real(8) :: t_0
t_0 = wj * exp(wj)
code = wj - ((t_0 - x) / (exp(wj) + t_0))
end function
public static double code(double wj, double x) {
double t_0 = wj * Math.exp(wj);
return wj - ((t_0 - x) / (Math.exp(wj) + t_0));
}
def code(wj, x): t_0 = wj * math.exp(wj) return wj - ((t_0 - x) / (math.exp(wj) + t_0))
function code(wj, x) t_0 = Float64(wj * exp(wj)) return Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0))) end
function tmp = code(wj, x) t_0 = wj * exp(wj); tmp = wj - ((t_0 - x) / (exp(wj) + t_0)); end
code[wj_, x_] := Block[{t$95$0 = N[(wj * N[Exp[wj], $MachinePrecision]), $MachinePrecision]}, N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := wj \cdot e^{wj}\\
wj - \frac{t_0 - x}{e^{wj} + t_0}
\end{array}
\end{array}
(FPCore (wj x) :precision binary64 (if (or (<= wj -7.5e-9) (not (<= wj 8e-11))) (fma (- (/ x (exp wj)) wj) (/ 1.0 (+ wj 1.0)) wj) (+ (+ x (* -2.0 (* wj x))) (* wj wj))))
double code(double wj, double x) {
double tmp;
if ((wj <= -7.5e-9) || !(wj <= 8e-11)) {
tmp = fma(((x / exp(wj)) - wj), (1.0 / (wj + 1.0)), wj);
} else {
tmp = (x + (-2.0 * (wj * x))) + (wj * wj);
}
return tmp;
}
function code(wj, x) tmp = 0.0 if ((wj <= -7.5e-9) || !(wj <= 8e-11)) tmp = fma(Float64(Float64(x / exp(wj)) - wj), Float64(1.0 / Float64(wj + 1.0)), wj); else tmp = Float64(Float64(x + Float64(-2.0 * Float64(wj * x))) + Float64(wj * wj)); end return tmp end
code[wj_, x_] := If[Or[LessEqual[wj, -7.5e-9], N[Not[LessEqual[wj, 8e-11]], $MachinePrecision]], N[(N[(N[(x / N[Exp[wj], $MachinePrecision]), $MachinePrecision] - wj), $MachinePrecision] * N[(1.0 / N[(wj + 1.0), $MachinePrecision]), $MachinePrecision] + wj), $MachinePrecision], N[(N[(x + N[(-2.0 * N[(wj * x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(wj * wj), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;wj \leq -7.5 \cdot 10^{-9} \lor \neg \left(wj \leq 8 \cdot 10^{-11}\right):\\
\;\;\;\;\mathsf{fma}\left(\frac{x}{e^{wj}} - wj, \frac{1}{wj + 1}, wj\right)\\
\mathbf{else}:\\
\;\;\;\;\left(x + -2 \cdot \left(wj \cdot x\right)\right) + wj \cdot wj\\
\end{array}
\end{array}
if wj < -7.49999999999999933e-9 or 7.99999999999999952e-11 < wj Initial program 60.5%
sub-neg60.5%
div-sub60.5%
sub-neg60.5%
+-commutative60.5%
distribute-neg-in60.5%
remove-double-neg60.5%
sub-neg60.5%
div-sub60.5%
distribute-rgt1-in81.9%
associate-/l/81.9%
Simplified96.2%
+-commutative96.2%
div-inv96.3%
fma-def96.6%
Applied egg-rr96.6%
if -7.49999999999999933e-9 < wj < 7.99999999999999952e-11Initial program 82.3%
sub-neg82.3%
div-sub82.3%
sub-neg82.3%
+-commutative82.3%
distribute-neg-in82.3%
remove-double-neg82.3%
sub-neg82.3%
div-sub82.3%
distribute-rgt1-in82.3%
associate-/l/82.3%
Simplified82.3%
Taylor expanded in wj around 0 99.9%
Taylor expanded in x around 0 99.9%
unpow299.9%
Simplified99.9%
Final simplification99.7%
(FPCore (wj x)
:precision binary64
(let* ((t_0 (+ (* x -4.0) (* x 1.5))) (t_1 (* wj (exp wj))))
(if (<= (+ wj (/ (- x t_1) (+ (exp wj) t_1))) 4e-16)
(+
(*
(pow wj 3.0)
(- (- (- -1.0 (* -2.0 t_0)) (* x -3.0)) (* x 0.6666666666666666)))
(+ (* (- 1.0 t_0) (pow wj 2.0)) (+ x (* -2.0 (* wj x)))))
(fma (- (/ x (exp wj)) wj) (/ 1.0 (+ wj 1.0)) wj))))
double code(double wj, double x) {
double t_0 = (x * -4.0) + (x * 1.5);
double t_1 = wj * exp(wj);
double tmp;
if ((wj + ((x - t_1) / (exp(wj) + t_1))) <= 4e-16) {
tmp = (pow(wj, 3.0) * (((-1.0 - (-2.0 * t_0)) - (x * -3.0)) - (x * 0.6666666666666666))) + (((1.0 - t_0) * pow(wj, 2.0)) + (x + (-2.0 * (wj * x))));
} else {
tmp = fma(((x / exp(wj)) - wj), (1.0 / (wj + 1.0)), wj);
}
return tmp;
}
function code(wj, x) t_0 = Float64(Float64(x * -4.0) + Float64(x * 1.5)) t_1 = Float64(wj * exp(wj)) tmp = 0.0 if (Float64(wj + Float64(Float64(x - t_1) / Float64(exp(wj) + t_1))) <= 4e-16) tmp = Float64(Float64((wj ^ 3.0) * Float64(Float64(Float64(-1.0 - Float64(-2.0 * t_0)) - Float64(x * -3.0)) - Float64(x * 0.6666666666666666))) + Float64(Float64(Float64(1.0 - t_0) * (wj ^ 2.0)) + Float64(x + Float64(-2.0 * Float64(wj * x))))); else tmp = fma(Float64(Float64(x / exp(wj)) - wj), Float64(1.0 / Float64(wj + 1.0)), wj); end return tmp end
code[wj_, x_] := Block[{t$95$0 = N[(N[(x * -4.0), $MachinePrecision] + N[(x * 1.5), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(wj * N[Exp[wj], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(wj + N[(N[(x - t$95$1), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 4e-16], N[(N[(N[Power[wj, 3.0], $MachinePrecision] * N[(N[(N[(-1.0 - N[(-2.0 * t$95$0), $MachinePrecision]), $MachinePrecision] - N[(x * -3.0), $MachinePrecision]), $MachinePrecision] - N[(x * 0.6666666666666666), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(N[(N[(1.0 - t$95$0), $MachinePrecision] * N[Power[wj, 2.0], $MachinePrecision]), $MachinePrecision] + N[(x + N[(-2.0 * N[(wj * x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(x / N[Exp[wj], $MachinePrecision]), $MachinePrecision] - wj), $MachinePrecision] * N[(1.0 / N[(wj + 1.0), $MachinePrecision]), $MachinePrecision] + wj), $MachinePrecision]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := x \cdot -4 + x \cdot 1.5\\
t_1 := wj \cdot e^{wj}\\
\mathbf{if}\;wj + \frac{x - t_1}{e^{wj} + t_1} \leq 4 \cdot 10^{-16}:\\
\;\;\;\;{wj}^{3} \cdot \left(\left(\left(-1 - -2 \cdot t_0\right) - x \cdot -3\right) - x \cdot 0.6666666666666666\right) + \left(\left(1 - t_0\right) \cdot {wj}^{2} + \left(x + -2 \cdot \left(wj \cdot x\right)\right)\right)\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(\frac{x}{e^{wj}} - wj, \frac{1}{wj + 1}, wj\right)\\
\end{array}
\end{array}
if (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < 3.9999999999999999e-16Initial program 74.2%
sub-neg74.2%
div-sub74.2%
sub-neg74.2%
+-commutative74.2%
distribute-neg-in74.2%
remove-double-neg74.2%
sub-neg74.2%
div-sub74.2%
distribute-rgt1-in75.8%
associate-/l/75.8%
Simplified75.8%
Taylor expanded in wj around 0 97.8%
if 3.9999999999999999e-16 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) Initial program 97.2%
sub-neg97.2%
div-sub97.2%
sub-neg97.2%
+-commutative97.2%
distribute-neg-in97.2%
remove-double-neg97.2%
sub-neg97.2%
div-sub97.2%
distribute-rgt1-in97.2%
associate-/l/97.2%
Simplified99.8%
+-commutative99.8%
div-inv99.8%
fma-def99.8%
Applied egg-rr99.8%
Final simplification98.4%
(FPCore (wj x) :precision binary64 (if (or (<= wj -6.8e-9) (not (<= wj 8e-11))) (- wj (/ (- wj (/ x (exp wj))) (+ wj 1.0))) (+ (+ x (* -2.0 (* wj x))) (* wj wj))))
double code(double wj, double x) {
double tmp;
if ((wj <= -6.8e-9) || !(wj <= 8e-11)) {
tmp = wj - ((wj - (x / exp(wj))) / (wj + 1.0));
} else {
tmp = (x + (-2.0 * (wj * x))) + (wj * wj);
}
return tmp;
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
real(8) :: tmp
if ((wj <= (-6.8d-9)) .or. (.not. (wj <= 8d-11))) then
tmp = wj - ((wj - (x / exp(wj))) / (wj + 1.0d0))
else
tmp = (x + ((-2.0d0) * (wj * x))) + (wj * wj)
end if
code = tmp
end function
public static double code(double wj, double x) {
double tmp;
if ((wj <= -6.8e-9) || !(wj <= 8e-11)) {
tmp = wj - ((wj - (x / Math.exp(wj))) / (wj + 1.0));
} else {
tmp = (x + (-2.0 * (wj * x))) + (wj * wj);
}
return tmp;
}
def code(wj, x): tmp = 0 if (wj <= -6.8e-9) or not (wj <= 8e-11): tmp = wj - ((wj - (x / math.exp(wj))) / (wj + 1.0)) else: tmp = (x + (-2.0 * (wj * x))) + (wj * wj) return tmp
function code(wj, x) tmp = 0.0 if ((wj <= -6.8e-9) || !(wj <= 8e-11)) tmp = Float64(wj - Float64(Float64(wj - Float64(x / exp(wj))) / Float64(wj + 1.0))); else tmp = Float64(Float64(x + Float64(-2.0 * Float64(wj * x))) + Float64(wj * wj)); end return tmp end
function tmp_2 = code(wj, x) tmp = 0.0; if ((wj <= -6.8e-9) || ~((wj <= 8e-11))) tmp = wj - ((wj - (x / exp(wj))) / (wj + 1.0)); else tmp = (x + (-2.0 * (wj * x))) + (wj * wj); end tmp_2 = tmp; end
code[wj_, x_] := If[Or[LessEqual[wj, -6.8e-9], N[Not[LessEqual[wj, 8e-11]], $MachinePrecision]], N[(wj - N[(N[(wj - N[(x / N[Exp[wj], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(wj + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(x + N[(-2.0 * N[(wj * x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(wj * wj), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;wj \leq -6.8 \cdot 10^{-9} \lor \neg \left(wj \leq 8 \cdot 10^{-11}\right):\\
\;\;\;\;wj - \frac{wj - \frac{x}{e^{wj}}}{wj + 1}\\
\mathbf{else}:\\
\;\;\;\;\left(x + -2 \cdot \left(wj \cdot x\right)\right) + wj \cdot wj\\
\end{array}
\end{array}
if wj < -6.7999999999999997e-9 or 7.99999999999999952e-11 < wj Initial program 60.5%
sub-neg60.5%
div-sub60.5%
sub-neg60.5%
+-commutative60.5%
distribute-neg-in60.5%
remove-double-neg60.5%
sub-neg60.5%
div-sub60.5%
distribute-rgt1-in81.9%
associate-/l/81.9%
Simplified96.2%
if -6.7999999999999997e-9 < wj < 7.99999999999999952e-11Initial program 82.3%
sub-neg82.3%
div-sub82.3%
sub-neg82.3%
+-commutative82.3%
distribute-neg-in82.3%
remove-double-neg82.3%
sub-neg82.3%
div-sub82.3%
distribute-rgt1-in82.3%
associate-/l/82.3%
Simplified82.3%
Taylor expanded in wj around 0 99.9%
Taylor expanded in x around 0 99.9%
unpow299.9%
Simplified99.9%
Final simplification99.7%
(FPCore (wj x) :precision binary64 (+ (+ x (* -2.0 (* wj x))) (* wj wj)))
double code(double wj, double x) {
return (x + (-2.0 * (wj * x))) + (wj * wj);
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = (x + ((-2.0d0) * (wj * x))) + (wj * wj)
end function
public static double code(double wj, double x) {
return (x + (-2.0 * (wj * x))) + (wj * wj);
}
def code(wj, x): return (x + (-2.0 * (wj * x))) + (wj * wj)
function code(wj, x) return Float64(Float64(x + Float64(-2.0 * Float64(wj * x))) + Float64(wj * wj)) end
function tmp = code(wj, x) tmp = (x + (-2.0 * (wj * x))) + (wj * wj); end
code[wj_, x_] := N[(N[(x + N[(-2.0 * N[(wj * x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(wj * wj), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\left(x + -2 \cdot \left(wj \cdot x\right)\right) + wj \cdot wj
\end{array}
Initial program 81.1%
sub-neg81.1%
div-sub81.1%
sub-neg81.1%
+-commutative81.1%
distribute-neg-in81.1%
remove-double-neg81.1%
sub-neg81.1%
div-sub81.1%
distribute-rgt1-in82.3%
associate-/l/82.2%
Simplified83.0%
Taylor expanded in wj around 0 95.9%
Taylor expanded in x around 0 95.8%
unpow295.8%
Simplified95.8%
Final simplification95.8%
(FPCore (wj x) :precision binary64 (+ x (* wj wj)))
double code(double wj, double x) {
return x + (wj * wj);
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = x + (wj * wj)
end function
public static double code(double wj, double x) {
return x + (wj * wj);
}
def code(wj, x): return x + (wj * wj)
function code(wj, x) return Float64(x + Float64(wj * wj)) end
function tmp = code(wj, x) tmp = x + (wj * wj); end
code[wj_, x_] := N[(x + N[(wj * wj), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
x + wj \cdot wj
\end{array}
Initial program 81.1%
sub-neg81.1%
div-sub81.1%
sub-neg81.1%
+-commutative81.1%
distribute-neg-in81.1%
remove-double-neg81.1%
sub-neg81.1%
div-sub81.1%
distribute-rgt1-in82.3%
associate-/l/82.2%
Simplified83.0%
Taylor expanded in wj around 0 95.9%
Taylor expanded in x around 0 95.8%
unpow295.8%
Simplified95.8%
Taylor expanded in wj around 0 95.1%
Final simplification95.1%
(FPCore (wj x) :precision binary64 wj)
double code(double wj, double x) {
return wj;
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = wj
end function
public static double code(double wj, double x) {
return wj;
}
def code(wj, x): return wj
function code(wj, x) return wj end
function tmp = code(wj, x) tmp = wj; end
code[wj_, x_] := wj
\begin{array}{l}
\\
wj
\end{array}
Initial program 81.1%
sub-neg81.1%
div-sub81.1%
sub-neg81.1%
+-commutative81.1%
distribute-neg-in81.1%
remove-double-neg81.1%
sub-neg81.1%
div-sub81.1%
distribute-rgt1-in82.3%
associate-/l/82.2%
Simplified83.0%
Taylor expanded in wj around inf 4.2%
Final simplification4.2%
(FPCore (wj x) :precision binary64 x)
double code(double wj, double x) {
return x;
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = x
end function
public static double code(double wj, double x) {
return x;
}
def code(wj, x): return x
function code(wj, x) return x end
function tmp = code(wj, x) tmp = x; end
code[wj_, x_] := x
\begin{array}{l}
\\
x
\end{array}
Initial program 81.1%
sub-neg81.1%
div-sub81.1%
sub-neg81.1%
+-commutative81.1%
distribute-neg-in81.1%
remove-double-neg81.1%
sub-neg81.1%
div-sub81.1%
distribute-rgt1-in82.3%
associate-/l/82.2%
Simplified83.0%
Taylor expanded in wj around 0 84.7%
Final simplification84.7%
(FPCore (wj x) :precision binary64 (- wj (- (/ wj (+ wj 1.0)) (/ x (+ (exp wj) (* wj (exp wj)))))))
double code(double wj, double x) {
return wj - ((wj / (wj + 1.0)) - (x / (exp(wj) + (wj * exp(wj)))));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = wj - ((wj / (wj + 1.0d0)) - (x / (exp(wj) + (wj * exp(wj)))))
end function
public static double code(double wj, double x) {
return wj - ((wj / (wj + 1.0)) - (x / (Math.exp(wj) + (wj * Math.exp(wj)))));
}
def code(wj, x): return wj - ((wj / (wj + 1.0)) - (x / (math.exp(wj) + (wj * math.exp(wj)))))
function code(wj, x) return Float64(wj - Float64(Float64(wj / Float64(wj + 1.0)) - Float64(x / Float64(exp(wj) + Float64(wj * exp(wj)))))) end
function tmp = code(wj, x) tmp = wj - ((wj / (wj + 1.0)) - (x / (exp(wj) + (wj * exp(wj))))); end
code[wj_, x_] := N[(wj - N[(N[(wj / N[(wj + 1.0), $MachinePrecision]), $MachinePrecision] - N[(x / N[(N[Exp[wj], $MachinePrecision] + N[(wj * N[Exp[wj], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
wj - \left(\frac{wj}{wj + 1} - \frac{x}{e^{wj} + wj \cdot e^{wj}}\right)
\end{array}
herbie shell --seed 2023188
(FPCore (wj x)
:name "Jmat.Real.lambertw, newton loop step"
:precision binary64
:herbie-target
(- wj (- (/ wj (+ wj 1.0)) (/ x (+ (exp wj) (* wj (exp wj))))))
(- wj (/ (- (* wj (exp wj)) x) (+ (exp wj) (* wj (exp wj))))))