
(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
(let* ((t_0 (* wj (exp wj))))
(if (<= (+ wj (/ (- x t_0) (+ (exp wj) t_0))) 2e-21)
(- x (* wj (+ (* x 2.0) (* wj (+ wj -1.0)))))
(+ wj (/ (- wj (/ x (exp wj))) (- -1.0 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-21) {
tmp = x - (wj * ((x * 2.0) + (wj * (wj + -1.0))));
} else {
tmp = wj + ((wj - (x / exp(wj))) / (-1.0 - wj));
}
return tmp;
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
real(8) :: t_0
real(8) :: tmp
t_0 = wj * exp(wj)
if ((wj + ((x - t_0) / (exp(wj) + t_0))) <= 2d-21) then
tmp = x - (wj * ((x * 2.0d0) + (wj * (wj + (-1.0d0)))))
else
tmp = wj + ((wj - (x / exp(wj))) / ((-1.0d0) - wj))
end if
code = tmp
end function
public static double code(double wj, double x) {
double t_0 = wj * Math.exp(wj);
double tmp;
if ((wj + ((x - t_0) / (Math.exp(wj) + t_0))) <= 2e-21) {
tmp = x - (wj * ((x * 2.0) + (wj * (wj + -1.0))));
} else {
tmp = wj + ((wj - (x / Math.exp(wj))) / (-1.0 - wj));
}
return tmp;
}
def code(wj, x): t_0 = wj * math.exp(wj) tmp = 0 if (wj + ((x - t_0) / (math.exp(wj) + t_0))) <= 2e-21: tmp = x - (wj * ((x * 2.0) + (wj * (wj + -1.0)))) else: tmp = wj + ((wj - (x / math.exp(wj))) / (-1.0 - 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-21) tmp = Float64(x - Float64(wj * Float64(Float64(x * 2.0) + Float64(wj * Float64(wj + -1.0))))); else tmp = Float64(wj + Float64(Float64(wj - Float64(x / exp(wj))) / Float64(-1.0 - wj))); end return tmp end
function tmp_2 = code(wj, x) t_0 = wj * exp(wj); tmp = 0.0; if ((wj + ((x - t_0) / (exp(wj) + t_0))) <= 2e-21) tmp = x - (wj * ((x * 2.0) + (wj * (wj + -1.0)))); else tmp = wj + ((wj - (x / exp(wj))) / (-1.0 - wj)); end tmp_2 = 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-21], N[(x - N[(wj * N[(N[(x * 2.0), $MachinePrecision] + N[(wj * N[(wj + -1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(wj + N[(N[(wj - N[(x / N[Exp[wj], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(-1.0 - wj), $MachinePrecision]), $MachinePrecision]), $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^{-21}:\\
\;\;\;\;x - wj \cdot \left(x \cdot 2 + wj \cdot \left(wj + -1\right)\right)\\
\mathbf{else}:\\
\;\;\;\;wj + \frac{wj - \frac{x}{e^{wj}}}{-1 - wj}\\
\end{array}
\end{array}
if (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < 1.99999999999999982e-21Initial program 73.1%
distribute-rgt1-in73.7%
associate-/l/73.7%
div-sub73.1%
associate-/l*73.1%
*-inverses73.7%
*-rgt-identity73.7%
Simplified73.7%
Taylor expanded in wj around 0 98.3%
Taylor expanded in x around 0 98.3%
mul-1-neg98.3%
unsub-neg98.3%
Simplified98.3%
if 1.99999999999999982e-21 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) Initial program 96.7%
distribute-rgt1-in98.0%
associate-/l/98.0%
div-sub96.7%
associate-/l*96.7%
*-inverses99.4%
*-rgt-identity99.4%
Simplified99.4%
Final simplification98.6%
(FPCore (wj x)
:precision binary64
(let* ((t_0 (+ (* x -4.0) (* x 1.5))))
(+
x
(*
wj
(-
(*
wj
(-
(+
1.0
(*
wj
(- -1.0 (+ (* x -3.0) (+ (* -2.0 t_0) (* x 0.6666666666666666))))))
t_0))
(* x 2.0))))))
double code(double wj, double x) {
double t_0 = (x * -4.0) + (x * 1.5);
return x + (wj * ((wj * ((1.0 + (wj * (-1.0 - ((x * -3.0) + ((-2.0 * t_0) + (x * 0.6666666666666666)))))) - t_0)) - (x * 2.0)));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
real(8) :: t_0
t_0 = (x * (-4.0d0)) + (x * 1.5d0)
code = x + (wj * ((wj * ((1.0d0 + (wj * ((-1.0d0) - ((x * (-3.0d0)) + (((-2.0d0) * t_0) + (x * 0.6666666666666666d0)))))) - t_0)) - (x * 2.0d0)))
end function
public static double code(double wj, double x) {
double t_0 = (x * -4.0) + (x * 1.5);
return x + (wj * ((wj * ((1.0 + (wj * (-1.0 - ((x * -3.0) + ((-2.0 * t_0) + (x * 0.6666666666666666)))))) - t_0)) - (x * 2.0)));
}
def code(wj, x): t_0 = (x * -4.0) + (x * 1.5) return x + (wj * ((wj * ((1.0 + (wj * (-1.0 - ((x * -3.0) + ((-2.0 * t_0) + (x * 0.6666666666666666)))))) - t_0)) - (x * 2.0)))
function code(wj, x) t_0 = Float64(Float64(x * -4.0) + Float64(x * 1.5)) return Float64(x + Float64(wj * Float64(Float64(wj * Float64(Float64(1.0 + Float64(wj * Float64(-1.0 - Float64(Float64(x * -3.0) + Float64(Float64(-2.0 * t_0) + Float64(x * 0.6666666666666666)))))) - t_0)) - Float64(x * 2.0)))) end
function tmp = code(wj, x) t_0 = (x * -4.0) + (x * 1.5); tmp = x + (wj * ((wj * ((1.0 + (wj * (-1.0 - ((x * -3.0) + ((-2.0 * t_0) + (x * 0.6666666666666666)))))) - t_0)) - (x * 2.0))); end
code[wj_, x_] := Block[{t$95$0 = N[(N[(x * -4.0), $MachinePrecision] + N[(x * 1.5), $MachinePrecision]), $MachinePrecision]}, N[(x + N[(wj * N[(N[(wj * N[(N[(1.0 + N[(wj * N[(-1.0 - N[(N[(x * -3.0), $MachinePrecision] + N[(N[(-2.0 * t$95$0), $MachinePrecision] + N[(x * 0.6666666666666666), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] - t$95$0), $MachinePrecision]), $MachinePrecision] - N[(x * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := x \cdot -4 + x \cdot 1.5\\
x + wj \cdot \left(wj \cdot \left(\left(1 + wj \cdot \left(-1 - \left(x \cdot -3 + \left(-2 \cdot t\_0 + x \cdot 0.6666666666666666\right)\right)\right)\right) - t\_0\right) - x \cdot 2\right)
\end{array}
\end{array}
Initial program 80.0%
distribute-rgt1-in80.8%
associate-/l/80.8%
div-sub80.0%
associate-/l*80.0%
*-inverses81.2%
*-rgt-identity81.2%
Simplified81.2%
Taylor expanded in wj around 0 96.7%
Final simplification96.7%
(FPCore (wj x) :precision binary64 (- x (* wj (+ (* x 2.0) (* wj (+ wj -1.0))))))
double code(double wj, double x) {
return x - (wj * ((x * 2.0) + (wj * (wj + -1.0))));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = x - (wj * ((x * 2.0d0) + (wj * (wj + (-1.0d0)))))
end function
public static double code(double wj, double x) {
return x - (wj * ((x * 2.0) + (wj * (wj + -1.0))));
}
def code(wj, x): return x - (wj * ((x * 2.0) + (wj * (wj + -1.0))))
function code(wj, x) return Float64(x - Float64(wj * Float64(Float64(x * 2.0) + Float64(wj * Float64(wj + -1.0))))) end
function tmp = code(wj, x) tmp = x - (wj * ((x * 2.0) + (wj * (wj + -1.0)))); end
code[wj_, x_] := N[(x - N[(wj * N[(N[(x * 2.0), $MachinePrecision] + N[(wj * N[(wj + -1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
x - wj \cdot \left(x \cdot 2 + wj \cdot \left(wj + -1\right)\right)
\end{array}
Initial program 80.0%
distribute-rgt1-in80.8%
associate-/l/80.8%
div-sub80.0%
associate-/l*80.0%
*-inverses81.2%
*-rgt-identity81.2%
Simplified81.2%
Taylor expanded in wj around 0 96.7%
Taylor expanded in x around 0 96.4%
mul-1-neg96.4%
unsub-neg96.4%
Simplified96.4%
Final simplification96.4%
(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 80.0%
distribute-rgt1-in80.8%
associate-/l/80.8%
div-sub80.0%
associate-/l*80.0%
*-inverses81.2%
*-rgt-identity81.2%
Simplified81.2%
Taylor expanded in wj around 0 96.7%
Taylor expanded in x around 0 96.4%
mul-1-neg96.4%
unsub-neg96.4%
Simplified96.4%
Taylor expanded in wj around 0 96.2%
Final simplification96.2%
(FPCore (wj x) :precision binary64 (+ x (* -2.0 (* wj x))))
double code(double wj, double x) {
return x + (-2.0 * (wj * x));
}
real(8) function code(wj, x)
real(8), intent (in) :: wj
real(8), intent (in) :: x
code = x + ((-2.0d0) * (wj * x))
end function
public static double code(double wj, double x) {
return x + (-2.0 * (wj * x));
}
def code(wj, x): return x + (-2.0 * (wj * x))
function code(wj, x) return Float64(x + Float64(-2.0 * Float64(wj * x))) end
function tmp = code(wj, x) tmp = x + (-2.0 * (wj * x)); end
code[wj_, x_] := N[(x + N[(-2.0 * N[(wj * x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
x + -2 \cdot \left(wj \cdot x\right)
\end{array}
Initial program 80.0%
distribute-rgt1-in80.8%
associate-/l/80.8%
div-sub80.0%
associate-/l*80.0%
*-inverses81.2%
*-rgt-identity81.2%
Simplified81.2%
Taylor expanded in wj around 0 86.5%
*-commutative86.5%
Simplified86.5%
Final simplification86.5%
(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 80.0%
distribute-rgt1-in80.8%
associate-/l/80.8%
div-sub80.0%
associate-/l*80.0%
*-inverses81.2%
*-rgt-identity81.2%
Simplified81.2%
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 80.0%
distribute-rgt1-in80.8%
associate-/l/80.8%
div-sub80.0%
associate-/l*80.0%
*-inverses81.2%
*-rgt-identity81.2%
Simplified81.2%
Taylor expanded in wj around 0 86.2%
Final simplification86.2%
(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 2024076
(FPCore (wj x)
:name "Jmat.Real.lambertw, newton loop step"
:precision binary64
:alt
(- wj (- (/ wj (+ wj 1.0)) (/ x (+ (exp wj) (* wj (exp wj))))))
(- wj (/ (- (* wj (exp wj)) x) (+ (exp wj) (* wj (exp wj))))))