Jmat.Real.lambertw, newton loop step

Average Error: 14.2 → 0.8

Time: 14.1s

Precision: binary64

Cost: 33476

Math

\[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]

↓

\[\begin{array}{l} t_0 := wj \cdot e^{wj}\\ \mathbf{if}\;wj + \frac{x - t_0}{e^{wj} + t_0} \leq 10^{-30}:\\ \;\;\;\;\mathsf{fma}\left(wj, wj + x \cdot \mathsf{fma}\left(wj, 2, -2\right), x\right)\\ \mathbf{else}:\\ \;\;\;\;wj + \frac{\frac{x}{e^{wj}} - wj}{wj + 1}\\ \end{array} \]

FPCore

(FPCore (wj x)
 :precision binary64
 (- wj (/ (- (* wj (exp wj)) x) (+ (exp wj) (* wj (exp wj))))))

↓

(FPCore (wj x)
 :precision binary64
 (let* ((t_0 (* wj (exp wj))))
   (if (<= (+ wj (/ (- x t_0) (+ (exp wj) t_0))) 1e-30)
     (fma wj (+ wj (* x (fma wj 2.0 -2.0))) x)
     (+ wj (/ (- (/ x (exp wj)) wj) (+ wj 1.0))))))

C

double code(double wj, double x) {
	return wj - (((wj * exp(wj)) - x) / (exp(wj) + (wj * exp(wj))));
}

↓

double code(double wj, double x) {
	double t_0 = wj * exp(wj);
	double tmp;
	if ((wj + ((x - t_0) / (exp(wj) + t_0))) <= 1e-30) {
		tmp = fma(wj, (wj + (x * fma(wj, 2.0, -2.0))), x);
	} else {
		tmp = wj + (((x / exp(wj)) - wj) / (wj + 1.0));
	}
	return tmp;
}

Julia

function code(wj, x)
	return Float64(wj - Float64(Float64(Float64(wj * exp(wj)) - x) / Float64(exp(wj) + Float64(wj * exp(wj)))))
end

↓

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))) <= 1e-30)
		tmp = fma(wj, Float64(wj + Float64(x * fma(wj, 2.0, -2.0))), x);
	else
		tmp = Float64(wj + Float64(Float64(Float64(x / exp(wj)) - wj) / Float64(wj + 1.0)));
	end
	return tmp
end

Wolfram

code[wj_, x_] := N[(wj - N[(N[(N[(wj * N[Exp[wj], $MachinePrecision]), $MachinePrecision] - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + N[(wj * N[Exp[wj], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

↓

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], 1e-30], N[(wj * N[(wj + N[(x * N[(wj * 2.0 + -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + x), $MachinePrecision], N[(wj + N[(N[(N[(x / N[Exp[wj], $MachinePrecision]), $MachinePrecision] - wj), $MachinePrecision] / N[(wj + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Target

Original	14.2
Target	13.5
Herbie	0.8

\[wj - \left(\frac{wj}{wj + 1} - \frac{x}{e^{wj} + wj \cdot e^{wj}}\right) \]

Derivation

1. Split input into 2 regimes

2. if (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < 1e-30

   1. Initial program 18.7

      \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]

   2. Simplified 18.7

      \[\leadsto \color{blue}{wj - \frac{wj - \frac{x}{e^{wj}}}{wj + 1}} \]
      Proof
      (-.f64 wj (/.f64 (-.f64 wj (/.f64 x (exp.f64 wj))) (+.f64 wj 1))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (Rewrite=> div-sub_binary64 (-.f64 (/.f64 wj (+.f64 wj 1)) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1))))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (-.f64 (Rewrite<= *-rgt-identity_binary64 (*.f64 (/.f64 wj (+.f64 wj 1)) 1)) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (-.f64 (*.f64 (/.f64 wj (+.f64 wj 1)) (Rewrite<= *-inverses_binary64 (/.f64 (exp.f64 wj) (exp.f64 wj)))) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 3 points increase in error, 0 points decrease in error
      (-.f64 wj (-.f64 (Rewrite<= times-frac_binary64 (/.f64 (*.f64 wj (exp.f64 wj)) (*.f64 (+.f64 wj 1) (exp.f64 wj)))) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 3 points increase in error, 1 points decrease in error
      (-.f64 wj (-.f64 (/.f64 (*.f64 wj (exp.f64 wj)) (Rewrite<= distribute-rgt1-in_binary64 (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 0 points increase in error, 1 points decrease in error
      (-.f64 wj (-.f64 (/.f64 (*.f64 wj (exp.f64 wj)) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))) (Rewrite=> associate-/l/_binary64 (/.f64 x (*.f64 (+.f64 wj 1) (exp.f64 wj)))))): 0 points increase in error, 1 points decrease in error
      (-.f64 wj (-.f64 (/.f64 (*.f64 wj (exp.f64 wj)) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))) (/.f64 x (Rewrite<= distribute-rgt1-in_binary64 (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (Rewrite<= div-sub_binary64 (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))): 0 points increase in error, 0 points decrease in error

   3. Taylor expanded in wj around 0 19.0

      \[\leadsto wj - \frac{wj - \color{blue}{\left(-1 \cdot \left(wj \cdot x\right) + x\right)}}{wj + 1} \]

   4. Simplified 19.0

      \[\leadsto wj - \frac{wj - \color{blue}{\left(x - x \cdot wj\right)}}{wj + 1} \]
      Proof
      (-.f64 x (*.f64 x wj)): 0 points increase in error, 0 points decrease in error
      (-.f64 x (Rewrite<= *-commutative_binary64 (*.f64 wj x))): 0 points increase in error, 0 points decrease in error
      (Rewrite<= unsub-neg_binary64 (+.f64 x (neg.f64 (*.f64 wj x)))): 0 points increase in error, 0 points decrease in error
      (+.f64 x (Rewrite<= mul-1-neg_binary64 (*.f64 -1 (*.f64 wj x)))): 0 points increase in error, 0 points decrease in error
      (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 -1 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error

   5. Taylor expanded in wj around 0 0.5

      \[\leadsto \color{blue}{{wj}^{2} \cdot \left(\left(1 + x\right) - -1 \cdot x\right) + \left(-2 \cdot \left(wj \cdot x\right) + x\right)} \]

   6. Taylor expanded in x around 0 0.5

      \[\leadsto \color{blue}{{wj}^{2} + \left(2 \cdot {wj}^{2} + \left(1 + -2 \cdot wj\right)\right) \cdot x} \]

   7. Simplified 0.5

      \[\leadsto \color{blue}{\mathsf{fma}\left(wj, wj + x \cdot \mathsf{fma}\left(wj, 2, -2\right), x\right)} \]
      Proof
      (fma.f64 wj (+.f64 wj (*.f64 x (fma.f64 wj 2 -2))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 wj (*.f64 x (Rewrite<= fma-def_binary64 (+.f64 (*.f64 wj 2) -2)))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 wj (*.f64 x (+.f64 (Rewrite<= *-commutative_binary64 (*.f64 2 wj)) -2))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 wj (*.f64 x (Rewrite<= +-commutative_binary64 (+.f64 -2 (*.f64 2 wj))))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 wj (Rewrite<= distribute-rgt-out_binary64 (+.f64 (*.f64 -2 x) (*.f64 (*.f64 2 wj) x)))) x): 1 points increase in error, 0 points decrease in error
      (fma.f64 wj (Rewrite=> +-commutative_binary64 (+.f64 (+.f64 (*.f64 -2 x) (*.f64 (*.f64 2 wj) x)) wj)) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (Rewrite<= associate-+r+_binary64 (+.f64 (*.f64 -2 x) (+.f64 (*.f64 (*.f64 2 wj) x) wj))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (+.f64 (*.f64 (Rewrite=> *-commutative_binary64 (*.f64 wj 2)) x) wj)) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (+.f64 (Rewrite=> associate-*l*_binary64 (*.f64 wj (*.f64 2 x))) wj)) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (+.f64 (*.f64 wj (*.f64 2 x)) (Rewrite<= *-rgt-identity_binary64 (*.f64 wj 1)))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (Rewrite<= distribute-lft-in_binary64 (*.f64 wj (+.f64 (*.f64 2 x) 1)))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (*.f64 wj (+.f64 (Rewrite<= count-2_binary64 (+.f64 x x)) 1))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (*.f64 wj (Rewrite=> associate-+l+_binary64 (+.f64 x (+.f64 x 1))))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (*.f64 wj (+.f64 x (Rewrite<= +-commutative_binary64 (+.f64 1 x))))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (*.f64 wj (Rewrite<= +-commutative_binary64 (+.f64 (+.f64 1 x) x)))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (*.f64 wj (+.f64 (+.f64 1 x) (Rewrite<= remove-double-neg_binary64 (neg.f64 (neg.f64 x)))))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (*.f64 wj (+.f64 (+.f64 1 x) (neg.f64 (Rewrite<= mul-1-neg_binary64 (*.f64 -1 x)))))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (+.f64 (*.f64 -2 x) (*.f64 wj (Rewrite<= sub-neg_binary64 (-.f64 (+.f64 1 x) (*.f64 -1 x))))) x): 0 points increase in error, 0 points decrease in error
      (fma.f64 wj (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 wj (-.f64 (+.f64 1 x) (*.f64 -1 x))) (*.f64 -2 x))) x): 0 points increase in error, 0 points decrease in error
      (Rewrite<= fma-def_binary64 (+.f64 (*.f64 wj (+.f64 (*.f64 wj (-.f64 (+.f64 1 x) (*.f64 -1 x))) (*.f64 -2 x))) x)): 1 points increase in error, 2 points decrease in error
      (+.f64 (Rewrite<= distribute-lft-out_binary64 (+.f64 (*.f64 wj (*.f64 wj (-.f64 (+.f64 1 x) (*.f64 -1 x)))) (*.f64 wj (*.f64 -2 x)))) x): 2 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 wj wj) (-.f64 (+.f64 1 x) (*.f64 -1 x)))) (*.f64 wj (*.f64 -2 x))) x): 0 points increase in error, 1 points decrease in error
      (+.f64 (+.f64 (*.f64 (Rewrite<= unpow2_binary64 (pow.f64 wj 2)) (-.f64 (+.f64 1 x) (*.f64 -1 x))) (*.f64 wj (*.f64 -2 x))) x): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 (pow.f64 wj 2) (-.f64 (+.f64 1 x) (*.f64 -1 x))) (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 wj -2) x))) x): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 (pow.f64 wj 2) (-.f64 (+.f64 1 x) (*.f64 -1 x))) (*.f64 (Rewrite<= *-commutative_binary64 (*.f64 -2 wj)) x)) x): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 (pow.f64 wj 2) (-.f64 (+.f64 1 x) (*.f64 -1 x))) (Rewrite<= associate-*r*_binary64 (*.f64 -2 (*.f64 wj x)))) x): 0 points increase in error, 0 points decrease in error
      (Rewrite=> associate-+l+_binary64 (+.f64 (*.f64 (pow.f64 wj 2) (-.f64 (+.f64 1 x) (*.f64 -1 x))) (+.f64 (*.f64 -2 (*.f64 wj x)) x))): 3 points increase in error, 2 points decrease in error
      (+.f64 (*.f64 (pow.f64 wj 2) (Rewrite=> associate--l+_binary64 (+.f64 1 (-.f64 x (*.f64 -1 x))))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (Rewrite=> distribute-rgt-in_binary64 (+.f64 (*.f64 1 (pow.f64 wj 2)) (*.f64 (-.f64 x (*.f64 -1 x)) (pow.f64 wj 2)))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 1 (pow.f64 wj 2)) (*.f64 (-.f64 (Rewrite<= *-lft-identity_binary64 (*.f64 1 x)) (*.f64 -1 x)) (pow.f64 wj 2))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 1 (pow.f64 wj 2)) (*.f64 (Rewrite=> distribute-rgt-out--_binary64 (*.f64 x (-.f64 1 -1))) (pow.f64 wj 2))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 1 (pow.f64 wj 2)) (*.f64 (*.f64 x (Rewrite=> metadata-eval 2)) (pow.f64 wj 2))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 1 (pow.f64 wj 2)) (*.f64 (Rewrite<= *-commutative_binary64 (*.f64 2 x)) (pow.f64 wj 2))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 1 (pow.f64 wj 2)) (Rewrite<= associate-*r*_binary64 (*.f64 2 (*.f64 x (pow.f64 wj 2))))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (*.f64 1 (pow.f64 wj 2)) (*.f64 2 (Rewrite<= *-commutative_binary64 (*.f64 (pow.f64 wj 2) x)))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (+.f64 (+.f64 (Rewrite=> *-lft-identity_binary64 (pow.f64 wj 2)) (*.f64 2 (*.f64 (pow.f64 wj 2) x))) (+.f64 (*.f64 -2 (*.f64 wj x)) x)): 0 points increase in error, 0 points decrease in error
      (Rewrite<= associate-+r+_binary64 (+.f64 (pow.f64 wj 2) (+.f64 (*.f64 2 (*.f64 (pow.f64 wj 2) x)) (+.f64 (*.f64 -2 (*.f64 wj x)) x)))): 0 points increase in error, 0 points decrease in error
      (+.f64 (pow.f64 wj 2) (+.f64 (Rewrite=> associate-*r*_binary64 (*.f64 (*.f64 2 (pow.f64 wj 2)) x)) (+.f64 (*.f64 -2 (*.f64 wj x)) x))): 0 points increase in error, 0 points decrease in error
      (+.f64 (pow.f64 wj 2) (+.f64 (*.f64 (*.f64 2 (pow.f64 wj 2)) x) (+.f64 (Rewrite=> associate-*r*_binary64 (*.f64 (*.f64 -2 wj) x)) x))): 0 points increase in error, 0 points decrease in error
      (+.f64 (pow.f64 wj 2) (+.f64 (*.f64 (*.f64 2 (pow.f64 wj 2)) x) (Rewrite=> distribute-lft1-in_binary64 (*.f64 (+.f64 (*.f64 -2 wj) 1) x)))): 1 points increase in error, 1 points decrease in error
      (+.f64 (pow.f64 wj 2) (Rewrite<= distribute-rgt-in_binary64 (*.f64 x (+.f64 (*.f64 2 (pow.f64 wj 2)) (+.f64 (*.f64 -2 wj) 1))))): 0 points increase in error, 2 points decrease in error
      (+.f64 (pow.f64 wj 2) (*.f64 x (+.f64 (*.f64 2 (pow.f64 wj 2)) (Rewrite<= +-commutative_binary64 (+.f64 1 (*.f64 -2 wj)))))): 0 points increase in error, 0 points decrease in error
      (+.f64 (pow.f64 wj 2) (Rewrite<= *-commutative_binary64 (*.f64 (+.f64 (*.f64 2 (pow.f64 wj 2)) (+.f64 1 (*.f64 -2 wj))) x))): 0 points increase in error, 0 points decrease in error

   if 1e-30 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

   1. Initial program 3.6

      \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]

   2. Simplified 1.5

      \[\leadsto \color{blue}{wj - \frac{wj - \frac{x}{e^{wj}}}{wj + 1}} \]
      Proof
      (-.f64 wj (/.f64 (-.f64 wj (/.f64 x (exp.f64 wj))) (+.f64 wj 1))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (Rewrite=> div-sub_binary64 (-.f64 (/.f64 wj (+.f64 wj 1)) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1))))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (-.f64 (Rewrite<= *-rgt-identity_binary64 (*.f64 (/.f64 wj (+.f64 wj 1)) 1)) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (-.f64 (*.f64 (/.f64 wj (+.f64 wj 1)) (Rewrite<= *-inverses_binary64 (/.f64 (exp.f64 wj) (exp.f64 wj)))) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 3 points increase in error, 0 points decrease in error
      (-.f64 wj (-.f64 (Rewrite<= times-frac_binary64 (/.f64 (*.f64 wj (exp.f64 wj)) (*.f64 (+.f64 wj 1) (exp.f64 wj)))) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 3 points increase in error, 1 points decrease in error
      (-.f64 wj (-.f64 (/.f64 (*.f64 wj (exp.f64 wj)) (Rewrite<= distribute-rgt1-in_binary64 (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) (/.f64 (/.f64 x (exp.f64 wj)) (+.f64 wj 1)))): 0 points increase in error, 1 points decrease in error
      (-.f64 wj (-.f64 (/.f64 (*.f64 wj (exp.f64 wj)) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))) (Rewrite=> associate-/l/_binary64 (/.f64 x (*.f64 (+.f64 wj 1) (exp.f64 wj)))))): 0 points increase in error, 1 points decrease in error
      (-.f64 wj (-.f64 (/.f64 (*.f64 wj (exp.f64 wj)) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))) (/.f64 x (Rewrite<= distribute-rgt1-in_binary64 (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))))): 0 points increase in error, 0 points decrease in error
      (-.f64 wj (Rewrite<= div-sub_binary64 (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))): 0 points increase in error, 0 points decrease in error
3. Recombined 2 regimes into one program.

4. Final simplification 0.8

   \[\leadsto \begin{array}{l} \mathbf{if}\;wj + \frac{x - wj \cdot e^{wj}}{e^{wj} + wj \cdot e^{wj}} \leq 10^{-30}:\\ \;\;\;\;\mathsf{fma}\left(wj, wj + x \cdot \mathsf{fma}\left(wj, 2, -2\right), x\right)\\ \mathbf{else}:\\ \;\;\;\;wj + \frac{\frac{x}{e^{wj}} - wj}{wj + 1}\\ \end{array} \]

Alternatives

Alternative 1
Error	0.8
Cost	7556

\[\begin{array}{l} \mathbf{if}\;wj \leq 2.7697794317313065 \cdot 10^{-9}:\\ \;\;\;\;{wj}^{2} \cdot \left(x + \left(x + 1\right)\right) + \left(x + -2 \cdot \left(wj \cdot x\right)\right)\\ \mathbf{else}:\\ \;\;\;\;wj + \frac{\frac{x}{e^{wj}} - wj}{wj + 1}\\ \end{array} \]

Alternative 2
Error	0.9
Cost	7236

\[\begin{array}{l} \mathbf{if}\;wj \leq 2.7697794317313065 \cdot 10^{-9}:\\ \;\;\;\;\mathsf{fma}\left(wj, wj + x \cdot -2, x\right)\\ \mathbf{else}:\\ \;\;\;\;wj + \frac{\frac{x}{e^{wj}} - wj}{wj + 1}\\ \end{array} \]

Alternative 3
Error	1.3
Cost	6980

\[\begin{array}{l} \mathbf{if}\;wj \leq 2.7697794317313065 \cdot 10^{-9}:\\ \;\;\;\;\mathsf{fma}\left(wj, wj + x \cdot -2, x\right)\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{wj + 1}\\ \end{array} \]

Alternative 4
Error	1.3
Cost	836

\[\begin{array}{l} \mathbf{if}\;wj \leq 2.7697794317313065 \cdot 10^{-9}:\\ \;\;\;\;\left(x + -2 \cdot \left(wj \cdot x\right)\right) + wj \cdot wj\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{wj + 1}\\ \end{array} \]

Alternative 5
Error	8.6
Cost	708

\[\begin{array}{l} \mathbf{if}\;wj \leq 2.7697794317313065 \cdot 10^{-9}:\\ \;\;\;\;x \cdot \frac{1 - wj}{wj + 1}\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{wj + 1}\\ \end{array} \]

Alternative 6
Error	8.6
Cost	580

\[\begin{array}{l} \mathbf{if}\;wj \leq 2.7697794317313065 \cdot 10^{-9}:\\ \;\;\;\;x + -2 \cdot \left(wj \cdot x\right)\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{wj + 1}\\ \end{array} \]

Alternative 7
Error	9.4
Cost	448

\[x + -2 \cdot \left(wj \cdot x\right) \]

Alternative 8
Error	9.7
Cost	64

\[x \]

Reproduce

herbie shell --seed 2022316 
(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))))))