
(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 6 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
(let* ((t_0 (* wj (exp wj))))
(if (<= (+ wj (/ (- x t_0) (+ (exp wj) t_0))) 2e-18)
(+ x (* wj (+ wj (* x (+ (* wj 2.5) -2.0)))))
(fma (/ -1.0 (+ wj 1.0)) (- wj (/ x (exp wj))) wj))))
double code(double wj, double x) {
double t_0 = wj * exp(wj);
double tmp;
if ((wj + ((x - t_0) / (exp(wj) + t_0))) <= 2e-18) {
tmp = x + (wj * (wj + (x * ((wj * 2.5) + -2.0))));
} else {
tmp = fma((-1.0 / (wj + 1.0)), (wj - (x / exp(wj))), wj);
}
return tmp;
}
function code(wj, x) t_0 = Float64(wj * exp(wj)) tmp = 0.0 if (Float64(wj + Float64(Float64(x - t_0) / Float64(exp(wj) + t_0))) <= 2e-18) tmp = Float64(x + Float64(wj * Float64(wj + Float64(x * Float64(Float64(wj * 2.5) + -2.0))))); else tmp = fma(Float64(-1.0 / Float64(wj + 1.0)), Float64(wj - Float64(x / exp(wj))), wj); end return tmp end
code[wj_, x_] := Block[{t$95$0 = N[(wj * N[Exp[wj], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(wj + N[(N[(x - t$95$0), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2e-18], N[(x + N[(wj * N[(wj + N[(x * N[(N[(wj * 2.5), $MachinePrecision] + -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(-1.0 / N[(wj + 1.0), $MachinePrecision]), $MachinePrecision] * N[(wj - N[(x / N[Exp[wj], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + wj), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := wj \cdot e^{wj}\\
\mathbf{if}\;wj + \frac{x - t\_0}{e^{wj} + t\_0} \leq 2 \cdot 10^{-18}:\\
\;\;\;\;x + wj \cdot \left(wj + x \cdot \left(wj \cdot 2.5 + -2\right)\right)\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(\frac{-1}{wj + 1}, wj - \frac{x}{e^{wj}}, 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))))) < 2.0000000000000001e-18Initial program 78.0%
sub-negN/A
+-lowering-+.f64N/A
distribute-rgt1-inN/A
associate-/l/N/A
distribute-neg-frac2N/A
/-lowering-/.f64N/A
div-subN/A
associate-/l*N/A
*-inversesN/A
*-rgt-identityN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f64N/A
+-commutativeN/A
*-inversesN/A
distribute-neg-inN/A
Simplified78.0%
Taylor expanded in wj around 0
Simplified99.1%
Taylor expanded in x around 0
Simplified99.1%
if 2.0000000000000001e-18 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) Initial program 96.8%
sub-negN/A
+-lowering-+.f64N/A
distribute-rgt1-inN/A
associate-/l/N/A
distribute-neg-frac2N/A
/-lowering-/.f64N/A
div-subN/A
associate-/l*N/A
*-inversesN/A
*-rgt-identityN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f64N/A
+-commutativeN/A
*-inversesN/A
distribute-neg-inN/A
Simplified99.4%
+-commutativeN/A
div-invN/A
*-commutativeN/A
fma-defineN/A
fma-lowering-fma.f64N/A
metadata-evalN/A
remove-double-negN/A
frac-2negN/A
/-lowering-/.f64N/A
sub-negN/A
+-commutativeN/A
distribute-neg-inN/A
remove-double-negN/A
metadata-evalN/A
metadata-evalN/A
+-lowering-+.f64N/A
metadata-evalN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f6499.5%
Applied egg-rr99.5%
Final simplification99.2%
(FPCore (wj x) :precision binary64 (+ x (* wj (+ wj (* x (+ (* wj 2.5) -2.0))))))
double code(double wj, double x) {
return x + (wj * (wj + (x * ((wj * 2.5) + -2.0))));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = x + (wj * (wj + (x * ((wj * 2.5d0) + (-2.0d0)))))
end function
public static double code(double wj, double x) {
return x + (wj * (wj + (x * ((wj * 2.5) + -2.0))));
}
def code(wj, x): return x + (wj * (wj + (x * ((wj * 2.5) + -2.0))))
function code(wj, x) return Float64(x + Float64(wj * Float64(wj + Float64(x * Float64(Float64(wj * 2.5) + -2.0))))) end
function tmp = code(wj, x) tmp = x + (wj * (wj + (x * ((wj * 2.5) + -2.0)))); end
code[wj_, x_] := N[(x + N[(wj * N[(wj + N[(x * N[(N[(wj * 2.5), $MachinePrecision] + -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
x + wj \cdot \left(wj + x \cdot \left(wj \cdot 2.5 + -2\right)\right)
\end{array}
Initial program 83.7%
sub-negN/A
+-lowering-+.f64N/A
distribute-rgt1-inN/A
associate-/l/N/A
distribute-neg-frac2N/A
/-lowering-/.f64N/A
div-subN/A
associate-/l*N/A
*-inversesN/A
*-rgt-identityN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f64N/A
+-commutativeN/A
*-inversesN/A
distribute-neg-inN/A
Simplified84.5%
Taylor expanded in wj around 0
Simplified97.1%
Taylor expanded in x around 0
Simplified97.2%
(FPCore (wj x) :precision binary64 (+ x (* wj (+ wj (* x -2.0)))))
double code(double wj, double x) {
return x + (wj * (wj + (x * -2.0)));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = x + (wj * (wj + (x * (-2.0d0))))
end function
public static double code(double wj, double x) {
return x + (wj * (wj + (x * -2.0)));
}
def code(wj, x): return x + (wj * (wj + (x * -2.0)))
function code(wj, x) return Float64(x + Float64(wj * Float64(wj + Float64(x * -2.0)))) end
function tmp = code(wj, x) tmp = x + (wj * (wj + (x * -2.0))); end
code[wj_, x_] := N[(x + N[(wj * N[(wj + N[(x * -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
x + wj \cdot \left(wj + x \cdot -2\right)
\end{array}
Initial program 83.7%
sub-negN/A
+-lowering-+.f64N/A
distribute-rgt1-inN/A
associate-/l/N/A
distribute-neg-frac2N/A
/-lowering-/.f64N/A
div-subN/A
associate-/l*N/A
*-inversesN/A
*-rgt-identityN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f64N/A
+-commutativeN/A
*-inversesN/A
distribute-neg-inN/A
Simplified84.5%
Taylor expanded in wj around 0
Simplified97.1%
Taylor expanded in x around 0
Simplified97.2%
Taylor expanded in wj around 0
*-commutativeN/A
*-lowering-*.f6496.8%
Simplified96.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 83.7%
sub-negN/A
+-lowering-+.f64N/A
distribute-rgt1-inN/A
associate-/l/N/A
distribute-neg-frac2N/A
/-lowering-/.f64N/A
div-subN/A
associate-/l*N/A
*-inversesN/A
*-rgt-identityN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f64N/A
+-commutativeN/A
*-inversesN/A
distribute-neg-inN/A
Simplified84.5%
Taylor expanded in wj around 0
Simplified97.1%
Taylor expanded in x around 0
unpow2N/A
*-lowering-*.f6496.2%
Simplified96.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 83.7%
sub-negN/A
+-lowering-+.f64N/A
distribute-rgt1-inN/A
associate-/l/N/A
distribute-neg-frac2N/A
/-lowering-/.f64N/A
div-subN/A
associate-/l*N/A
*-inversesN/A
*-rgt-identityN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f64N/A
+-commutativeN/A
*-inversesN/A
distribute-neg-inN/A
Simplified84.5%
Taylor expanded in wj around 0
Simplified89.8%
(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 83.7%
sub-negN/A
+-lowering-+.f64N/A
distribute-rgt1-inN/A
associate-/l/N/A
distribute-neg-frac2N/A
/-lowering-/.f64N/A
div-subN/A
associate-/l*N/A
*-inversesN/A
*-rgt-identityN/A
--lowering--.f64N/A
/-lowering-/.f64N/A
exp-lowering-exp.f64N/A
+-commutativeN/A
*-inversesN/A
distribute-neg-inN/A
Simplified84.5%
Taylor expanded in wj around inf
Simplified4.3%
(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 2024138
(FPCore (wj x)
:name "Jmat.Real.lambertw, newton loop step"
:precision binary64
:alt
(! :herbie-platform default (let ((ew (exp wj))) (- wj (- (/ wj (+ wj 1)) (/ x (+ ew (* wj ew)))))))
(- wj (/ (- (* wj (exp wj)) x) (+ (exp wj) (* wj (exp wj))))))