Jmat.Real.lambertw, newton loop step

Percentage Accurate: 78.2% → 96.8%
Time: 8.3s
Alternatives: 8
Speedup: 0.9×

Specification

?
\[\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)))) (- 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:

Local Percentage Accuracy vs ?

The average percentage accuracy by input value. Horizontal axis shows value of an input variable; the variable is choosen in the title. Vertical axis is accuracy; higher is better. Red represent the original program, while blue represents Herbie's suggestion. These can be toggled with buttons below the plot. The line is an average while dots represent individual samples.

Accuracy vs Speed?

Herbie found 8 alternatives:

AlternativeAccuracySpeedup
The accuracy (vertical axis) and speed (horizontal axis) of each alternatives. Up and to the right is better. The red square shows the initial program, and each blue circle shows an alternative.The line shows the best available speed-accuracy tradeoffs.

Initial Program: 78.2% accurate, 1.0× speedup?

\[\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)))) (- 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}

Alternative 1: 96.8% accurate, 0.9× speedup?

\[\begin{array}{l} \\ \begin{array}{l} t_0 := e^{wj} \cdot wj\\ \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(-2.6666666666666665, wj, 2.5\right), wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot x\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{1 + wj}\\ \end{array} \end{array} \]
(FPCore (wj x)
 :precision binary64
 (let* ((t_0 (* (exp wj) wj)))
   (if (<= (- wj (/ (- t_0 x) (+ (exp wj) t_0))) 5e+302)
     (*
      (fma
       (fma (fma -2.6666666666666665 wj 2.5) wj -2.0)
       wj
       (fma (/ (- 1.0 wj) x) (* wj wj) 1.0))
      x)
     (- wj (/ wj (+ 1.0 wj))))))
double code(double wj, double x) {
	double t_0 = exp(wj) * wj;
	double tmp;
	if ((wj - ((t_0 - x) / (exp(wj) + t_0))) <= 5e+302) {
		tmp = fma(fma(fma(-2.6666666666666665, wj, 2.5), wj, -2.0), wj, fma(((1.0 - wj) / x), (wj * wj), 1.0)) * x;
	} else {
		tmp = wj - (wj / (1.0 + wj));
	}
	return tmp;
}
function code(wj, x)
	t_0 = Float64(exp(wj) * wj)
	tmp = 0.0
	if (Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0))) <= 5e+302)
		tmp = Float64(fma(fma(fma(-2.6666666666666665, wj, 2.5), wj, -2.0), wj, fma(Float64(Float64(1.0 - wj) / x), Float64(wj * wj), 1.0)) * x);
	else
		tmp = Float64(wj - Float64(wj / Float64(1.0 + wj)));
	end
	return tmp
end
code[wj_, x_] := Block[{t$95$0 = N[(N[Exp[wj], $MachinePrecision] * wj), $MachinePrecision]}, If[LessEqual[N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 5e+302], N[(N[(N[(N[(-2.6666666666666665 * wj + 2.5), $MachinePrecision] * wj + -2.0), $MachinePrecision] * wj + N[(N[(N[(1.0 - wj), $MachinePrecision] / x), $MachinePrecision] * N[(wj * wj), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] * x), $MachinePrecision], N[(wj - N[(wj / N[(1.0 + wj), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}

\\
\begin{array}{l}
t_0 := e^{wj} \cdot wj\\
\mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\
\;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(-2.6666666666666665, wj, 2.5\right), wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot x\\

\mathbf{else}:\\
\;\;\;\;wj - \frac{wj}{1 + wj}\\


\end{array}
\end{array}
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))))) < 5e302

    1. Initial program 74.6%

      \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
    2. Add Preprocessing
    3. Taylor expanded in wj around 0

      \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
    4. Applied rewrites99.3%

      \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
    5. Taylor expanded in x around inf

      \[\leadsto x \cdot \color{blue}{\left(1 + \left(wj \cdot \left(wj \cdot \left(\frac{5}{2} - \frac{8}{3} \cdot wj\right) - 2\right) + \frac{{wj}^{2} \cdot \left(1 - wj\right)}{x}\right)\right)} \]
    6. Step-by-step derivation
      1. Applied rewrites99.7%

        \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(-2.6666666666666665, wj, 2.5\right), wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot \color{blue}{x} \]

      if 5e302 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

      1. Initial program 0.0%

        \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
      2. Add Preprocessing
      3. Taylor expanded in x around 0

        \[\leadsto wj - \color{blue}{\frac{wj \cdot e^{wj}}{e^{wj} + wj \cdot e^{wj}}} \]
      4. Step-by-step derivation
        1. distribute-rgt1-inN/A

          \[\leadsto wj - \frac{wj \cdot e^{wj}}{\color{blue}{\left(wj + 1\right) \cdot e^{wj}}} \]
        2. +-commutativeN/A

          \[\leadsto wj - \frac{wj \cdot e^{wj}}{\color{blue}{\left(1 + wj\right)} \cdot e^{wj}} \]
        3. times-fracN/A

          \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj} \cdot \frac{e^{wj}}{e^{wj}}} \]
        4. *-inversesN/A

          \[\leadsto wj - \frac{wj}{1 + wj} \cdot \color{blue}{1} \]
        5. associate-*l/N/A

          \[\leadsto wj - \color{blue}{\frac{wj \cdot 1}{1 + wj}} \]
        6. *-rgt-identityN/A

          \[\leadsto wj - \frac{\color{blue}{wj}}{1 + wj} \]
        7. lower-/.f64N/A

          \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj}} \]
        8. lower-+.f6483.9

          \[\leadsto wj - \frac{wj}{\color{blue}{1 + wj}} \]
      5. Applied rewrites83.9%

        \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj}} \]
    7. Recombined 2 regimes into one program.
    8. Final simplification98.9%

      \[\leadsto \begin{array}{l} \mathbf{if}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(-2.6666666666666665, wj, 2.5\right), wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot x\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{1 + wj}\\ \end{array} \]
    9. Add Preprocessing

    Alternative 2: 85.0% accurate, 0.3× speedup?

    \[\begin{array}{l} \\ \begin{array}{l} t_0 := e^{wj} \cdot wj\\ t_1 := wj - \frac{t\_0 - x}{e^{wj} + t\_0}\\ \mathbf{if}\;t\_1 \leq -1 \cdot 10^{-308}:\\ \;\;\;\;1 \cdot x\\ \mathbf{elif}\;t\_1 \leq 0:\\ \;\;\;\;wj \cdot wj\\ \mathbf{elif}\;t\_1 \leq 5 \cdot 10^{+302}:\\ \;\;\;\;\mathsf{fma}\left(-2, wj, 1\right) \cdot x\\ \mathbf{else}:\\ \;\;\;\;wj - 1\\ \end{array} \end{array} \]
    (FPCore (wj x)
     :precision binary64
     (let* ((t_0 (* (exp wj) wj)) (t_1 (- wj (/ (- t_0 x) (+ (exp wj) t_0)))))
       (if (<= t_1 -1e-308)
         (* 1.0 x)
         (if (<= t_1 0.0)
           (* wj wj)
           (if (<= t_1 5e+302) (* (fma -2.0 wj 1.0) x) (- wj 1.0))))))
    double code(double wj, double x) {
    	double t_0 = exp(wj) * wj;
    	double t_1 = wj - ((t_0 - x) / (exp(wj) + t_0));
    	double tmp;
    	if (t_1 <= -1e-308) {
    		tmp = 1.0 * x;
    	} else if (t_1 <= 0.0) {
    		tmp = wj * wj;
    	} else if (t_1 <= 5e+302) {
    		tmp = fma(-2.0, wj, 1.0) * x;
    	} else {
    		tmp = wj - 1.0;
    	}
    	return tmp;
    }
    
    function code(wj, x)
    	t_0 = Float64(exp(wj) * wj)
    	t_1 = Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0)))
    	tmp = 0.0
    	if (t_1 <= -1e-308)
    		tmp = Float64(1.0 * x);
    	elseif (t_1 <= 0.0)
    		tmp = Float64(wj * wj);
    	elseif (t_1 <= 5e+302)
    		tmp = Float64(fma(-2.0, wj, 1.0) * x);
    	else
    		tmp = Float64(wj - 1.0);
    	end
    	return tmp
    end
    
    code[wj_, x_] := Block[{t$95$0 = N[(N[Exp[wj], $MachinePrecision] * wj), $MachinePrecision]}, Block[{t$95$1 = N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$1, -1e-308], N[(1.0 * x), $MachinePrecision], If[LessEqual[t$95$1, 0.0], N[(wj * wj), $MachinePrecision], If[LessEqual[t$95$1, 5e+302], N[(N[(-2.0 * wj + 1.0), $MachinePrecision] * x), $MachinePrecision], N[(wj - 1.0), $MachinePrecision]]]]]]
    
    \begin{array}{l}
    
    \\
    \begin{array}{l}
    t_0 := e^{wj} \cdot wj\\
    t_1 := wj - \frac{t\_0 - x}{e^{wj} + t\_0}\\
    \mathbf{if}\;t\_1 \leq -1 \cdot 10^{-308}:\\
    \;\;\;\;1 \cdot x\\
    
    \mathbf{elif}\;t\_1 \leq 0:\\
    \;\;\;\;wj \cdot wj\\
    
    \mathbf{elif}\;t\_1 \leq 5 \cdot 10^{+302}:\\
    \;\;\;\;\mathsf{fma}\left(-2, wj, 1\right) \cdot x\\
    
    \mathbf{else}:\\
    \;\;\;\;wj - 1\\
    
    
    \end{array}
    \end{array}
    
    Derivation
    1. Split input into 4 regimes
    2. if (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < -9.9999999999999991e-309

      1. Initial program 98.7%

        \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
      2. Add Preprocessing
      3. Taylor expanded in wj around 0

        \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
      4. Applied rewrites99.0%

        \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
      5. Taylor expanded in x around inf

        \[\leadsto x \cdot \color{blue}{\left(1 + \left(wj \cdot \left(wj \cdot \left(\frac{5}{2} - \frac{8}{3} \cdot wj\right) - 2\right) + \frac{{wj}^{2} \cdot \left(1 - wj\right)}{x}\right)\right)} \]
      6. Step-by-step derivation
        1. Applied rewrites100.0%

          \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(-2.6666666666666665, wj, 2.5\right), wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot \color{blue}{x} \]
        2. Taylor expanded in wj around 0

          \[\leadsto 1 \cdot x \]
        3. Step-by-step derivation
          1. Applied rewrites99.0%

            \[\leadsto 1 \cdot x \]

          if -9.9999999999999991e-309 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < 0.0

          1. Initial program 5.2%

            \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
          2. Add Preprocessing
          3. Taylor expanded in wj around 0

            \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
          4. Applied rewrites100.0%

            \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
          5. Taylor expanded in x around 0

            \[\leadsto {wj}^{2} \cdot \color{blue}{\left(1 - wj\right)} \]
          6. Step-by-step derivation
            1. Applied rewrites63.8%

              \[\leadsto \left(\left(1 - wj\right) \cdot wj\right) \cdot \color{blue}{wj} \]
            2. Taylor expanded in wj around 0

              \[\leadsto {wj}^{2} \]
            3. Step-by-step derivation
              1. Applied rewrites63.8%

                \[\leadsto wj \cdot wj \]

              if 0.0 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < 5e302

              1. Initial program 95.7%

                \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
              2. Add Preprocessing
              3. Taylor expanded in wj around 0

                \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
              4. Applied rewrites99.3%

                \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
              5. Taylor expanded in wj around 0

                \[\leadsto \color{blue}{x + -2 \cdot \left(wj \cdot x\right)} \]
              6. Step-by-step derivation
                1. associate-*r*N/A

                  \[\leadsto x + \color{blue}{\left(-2 \cdot wj\right) \cdot x} \]
                2. distribute-rgt1-inN/A

                  \[\leadsto \color{blue}{\left(-2 \cdot wj + 1\right) \cdot x} \]
                3. +-commutativeN/A

                  \[\leadsto \color{blue}{\left(1 + -2 \cdot wj\right)} \cdot x \]
                4. lower-*.f64N/A

                  \[\leadsto \color{blue}{\left(1 + -2 \cdot wj\right) \cdot x} \]
                5. +-commutativeN/A

                  \[\leadsto \color{blue}{\left(-2 \cdot wj + 1\right)} \cdot x \]
                6. lower-fma.f6491.9

                  \[\leadsto \color{blue}{\mathsf{fma}\left(-2, wj, 1\right)} \cdot x \]
              7. Applied rewrites91.9%

                \[\leadsto \color{blue}{\mathsf{fma}\left(-2, wj, 1\right) \cdot x} \]

              if 5e302 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

              1. Initial program 0.0%

                \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
              2. Add Preprocessing
              3. Taylor expanded in wj around inf

                \[\leadsto wj - \color{blue}{1} \]
              4. Step-by-step derivation
                1. Applied rewrites64.5%

                  \[\leadsto wj - \color{blue}{1} \]
              5. Recombined 4 regimes into one program.
              6. Final simplification86.7%

                \[\leadsto \begin{array}{l} \mathbf{if}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq -1 \cdot 10^{-308}:\\ \;\;\;\;1 \cdot x\\ \mathbf{elif}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq 0:\\ \;\;\;\;wj \cdot wj\\ \mathbf{elif}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;\mathsf{fma}\left(-2, wj, 1\right) \cdot x\\ \mathbf{else}:\\ \;\;\;\;wj - 1\\ \end{array} \]
              7. Add Preprocessing

              Alternative 3: 84.7% accurate, 0.3× speedup?

              \[\begin{array}{l} \\ \begin{array}{l} t_0 := e^{wj} \cdot wj\\ t_1 := wj - \frac{t\_0 - x}{e^{wj} + t\_0}\\ \mathbf{if}\;t\_1 \leq -1 \cdot 10^{-308}:\\ \;\;\;\;1 \cdot x\\ \mathbf{elif}\;t\_1 \leq 0:\\ \;\;\;\;wj \cdot wj\\ \mathbf{elif}\;t\_1 \leq 5 \cdot 10^{+302}:\\ \;\;\;\;1 \cdot x\\ \mathbf{else}:\\ \;\;\;\;wj - 1\\ \end{array} \end{array} \]
              (FPCore (wj x)
               :precision binary64
               (let* ((t_0 (* (exp wj) wj)) (t_1 (- wj (/ (- t_0 x) (+ (exp wj) t_0)))))
                 (if (<= t_1 -1e-308)
                   (* 1.0 x)
                   (if (<= t_1 0.0) (* wj wj) (if (<= t_1 5e+302) (* 1.0 x) (- wj 1.0))))))
              double code(double wj, double x) {
              	double t_0 = exp(wj) * wj;
              	double t_1 = wj - ((t_0 - x) / (exp(wj) + t_0));
              	double tmp;
              	if (t_1 <= -1e-308) {
              		tmp = 1.0 * x;
              	} else if (t_1 <= 0.0) {
              		tmp = wj * wj;
              	} else if (t_1 <= 5e+302) {
              		tmp = 1.0 * x;
              	} else {
              		tmp = wj - 1.0;
              	}
              	return tmp;
              }
              
              real(8) function code(wj, x)
                  real(8), intent (in) :: wj
                  real(8), intent (in) :: x
                  real(8) :: t_0
                  real(8) :: t_1
                  real(8) :: tmp
                  t_0 = exp(wj) * wj
                  t_1 = wj - ((t_0 - x) / (exp(wj) + t_0))
                  if (t_1 <= (-1d-308)) then
                      tmp = 1.0d0 * x
                  else if (t_1 <= 0.0d0) then
                      tmp = wj * wj
                  else if (t_1 <= 5d+302) then
                      tmp = 1.0d0 * x
                  else
                      tmp = wj - 1.0d0
                  end if
                  code = tmp
              end function
              
              public static double code(double wj, double x) {
              	double t_0 = Math.exp(wj) * wj;
              	double t_1 = wj - ((t_0 - x) / (Math.exp(wj) + t_0));
              	double tmp;
              	if (t_1 <= -1e-308) {
              		tmp = 1.0 * x;
              	} else if (t_1 <= 0.0) {
              		tmp = wj * wj;
              	} else if (t_1 <= 5e+302) {
              		tmp = 1.0 * x;
              	} else {
              		tmp = wj - 1.0;
              	}
              	return tmp;
              }
              
              def code(wj, x):
              	t_0 = math.exp(wj) * wj
              	t_1 = wj - ((t_0 - x) / (math.exp(wj) + t_0))
              	tmp = 0
              	if t_1 <= -1e-308:
              		tmp = 1.0 * x
              	elif t_1 <= 0.0:
              		tmp = wj * wj
              	elif t_1 <= 5e+302:
              		tmp = 1.0 * x
              	else:
              		tmp = wj - 1.0
              	return tmp
              
              function code(wj, x)
              	t_0 = Float64(exp(wj) * wj)
              	t_1 = Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0)))
              	tmp = 0.0
              	if (t_1 <= -1e-308)
              		tmp = Float64(1.0 * x);
              	elseif (t_1 <= 0.0)
              		tmp = Float64(wj * wj);
              	elseif (t_1 <= 5e+302)
              		tmp = Float64(1.0 * x);
              	else
              		tmp = Float64(wj - 1.0);
              	end
              	return tmp
              end
              
              function tmp_2 = code(wj, x)
              	t_0 = exp(wj) * wj;
              	t_1 = wj - ((t_0 - x) / (exp(wj) + t_0));
              	tmp = 0.0;
              	if (t_1 <= -1e-308)
              		tmp = 1.0 * x;
              	elseif (t_1 <= 0.0)
              		tmp = wj * wj;
              	elseif (t_1 <= 5e+302)
              		tmp = 1.0 * x;
              	else
              		tmp = wj - 1.0;
              	end
              	tmp_2 = tmp;
              end
              
              code[wj_, x_] := Block[{t$95$0 = N[(N[Exp[wj], $MachinePrecision] * wj), $MachinePrecision]}, Block[{t$95$1 = N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$1, -1e-308], N[(1.0 * x), $MachinePrecision], If[LessEqual[t$95$1, 0.0], N[(wj * wj), $MachinePrecision], If[LessEqual[t$95$1, 5e+302], N[(1.0 * x), $MachinePrecision], N[(wj - 1.0), $MachinePrecision]]]]]]
              
              \begin{array}{l}
              
              \\
              \begin{array}{l}
              t_0 := e^{wj} \cdot wj\\
              t_1 := wj - \frac{t\_0 - x}{e^{wj} + t\_0}\\
              \mathbf{if}\;t\_1 \leq -1 \cdot 10^{-308}:\\
              \;\;\;\;1 \cdot x\\
              
              \mathbf{elif}\;t\_1 \leq 0:\\
              \;\;\;\;wj \cdot wj\\
              
              \mathbf{elif}\;t\_1 \leq 5 \cdot 10^{+302}:\\
              \;\;\;\;1 \cdot x\\
              
              \mathbf{else}:\\
              \;\;\;\;wj - 1\\
              
              
              \end{array}
              \end{array}
              
              Derivation
              1. Split input into 3 regimes
              2. if (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < -9.9999999999999991e-309 or 0.0 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < 5e302

                1. Initial program 97.3%

                  \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                2. Add Preprocessing
                3. Taylor expanded in wj around 0

                  \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
                4. Applied rewrites99.1%

                  \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
                5. Taylor expanded in x around inf

                  \[\leadsto x \cdot \color{blue}{\left(1 + \left(wj \cdot \left(wj \cdot \left(\frac{5}{2} - \frac{8}{3} \cdot wj\right) - 2\right) + \frac{{wj}^{2} \cdot \left(1 - wj\right)}{x}\right)\right)} \]
                6. Step-by-step derivation
                  1. Applied rewrites99.7%

                    \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(-2.6666666666666665, wj, 2.5\right), wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot \color{blue}{x} \]
                  2. Taylor expanded in wj around 0

                    \[\leadsto 1 \cdot x \]
                  3. Step-by-step derivation
                    1. Applied rewrites95.5%

                      \[\leadsto 1 \cdot x \]

                    if -9.9999999999999991e-309 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj))))) < 0.0

                    1. Initial program 5.2%

                      \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                    2. Add Preprocessing
                    3. Taylor expanded in wj around 0

                      \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
                    4. Applied rewrites100.0%

                      \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
                    5. Taylor expanded in x around 0

                      \[\leadsto {wj}^{2} \cdot \color{blue}{\left(1 - wj\right)} \]
                    6. Step-by-step derivation
                      1. Applied rewrites63.8%

                        \[\leadsto \left(\left(1 - wj\right) \cdot wj\right) \cdot \color{blue}{wj} \]
                      2. Taylor expanded in wj around 0

                        \[\leadsto {wj}^{2} \]
                      3. Step-by-step derivation
                        1. Applied rewrites63.8%

                          \[\leadsto wj \cdot wj \]

                        if 5e302 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

                        1. Initial program 0.0%

                          \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                        2. Add Preprocessing
                        3. Taylor expanded in wj around inf

                          \[\leadsto wj - \color{blue}{1} \]
                        4. Step-by-step derivation
                          1. Applied rewrites64.5%

                            \[\leadsto wj - \color{blue}{1} \]
                        5. Recombined 3 regimes into one program.
                        6. Final simplification86.6%

                          \[\leadsto \begin{array}{l} \mathbf{if}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq -1 \cdot 10^{-308}:\\ \;\;\;\;1 \cdot x\\ \mathbf{elif}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq 0:\\ \;\;\;\;wj \cdot wj\\ \mathbf{elif}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;1 \cdot x\\ \mathbf{else}:\\ \;\;\;\;wj - 1\\ \end{array} \]
                        7. Add Preprocessing

                        Alternative 4: 96.8% accurate, 0.9× speedup?

                        \[\begin{array}{l} \\ \begin{array}{l} t_0 := e^{wj} \cdot wj\\ \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;\mathsf{fma}\left(\left(\frac{1 - wj}{x} \cdot wj + \mathsf{fma}\left(2.5, wj, -2\right)\right) \cdot wj, x, x\right)\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{1 + wj}\\ \end{array} \end{array} \]
                        (FPCore (wj x)
                         :precision binary64
                         (let* ((t_0 (* (exp wj) wj)))
                           (if (<= (- wj (/ (- t_0 x) (+ (exp wj) t_0))) 5e+302)
                             (fma (* (+ (* (/ (- 1.0 wj) x) wj) (fma 2.5 wj -2.0)) wj) x x)
                             (- wj (/ wj (+ 1.0 wj))))))
                        double code(double wj, double x) {
                        	double t_0 = exp(wj) * wj;
                        	double tmp;
                        	if ((wj - ((t_0 - x) / (exp(wj) + t_0))) <= 5e+302) {
                        		tmp = fma((((((1.0 - wj) / x) * wj) + fma(2.5, wj, -2.0)) * wj), x, x);
                        	} else {
                        		tmp = wj - (wj / (1.0 + wj));
                        	}
                        	return tmp;
                        }
                        
                        function code(wj, x)
                        	t_0 = Float64(exp(wj) * wj)
                        	tmp = 0.0
                        	if (Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0))) <= 5e+302)
                        		tmp = fma(Float64(Float64(Float64(Float64(Float64(1.0 - wj) / x) * wj) + fma(2.5, wj, -2.0)) * wj), x, x);
                        	else
                        		tmp = Float64(wj - Float64(wj / Float64(1.0 + wj)));
                        	end
                        	return tmp
                        end
                        
                        code[wj_, x_] := Block[{t$95$0 = N[(N[Exp[wj], $MachinePrecision] * wj), $MachinePrecision]}, If[LessEqual[N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 5e+302], N[(N[(N[(N[(N[(N[(1.0 - wj), $MachinePrecision] / x), $MachinePrecision] * wj), $MachinePrecision] + N[(2.5 * wj + -2.0), $MachinePrecision]), $MachinePrecision] * wj), $MachinePrecision] * x + x), $MachinePrecision], N[(wj - N[(wj / N[(1.0 + wj), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]
                        
                        \begin{array}{l}
                        
                        \\
                        \begin{array}{l}
                        t_0 := e^{wj} \cdot wj\\
                        \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\
                        \;\;\;\;\mathsf{fma}\left(\left(\frac{1 - wj}{x} \cdot wj + \mathsf{fma}\left(2.5, wj, -2\right)\right) \cdot wj, x, x\right)\\
                        
                        \mathbf{else}:\\
                        \;\;\;\;wj - \frac{wj}{1 + wj}\\
                        
                        
                        \end{array}
                        \end{array}
                        
                        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))))) < 5e302

                          1. Initial program 74.6%

                            \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                          2. Add Preprocessing
                          3. Taylor expanded in wj around 0

                            \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
                          4. Applied rewrites99.3%

                            \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
                          5. Taylor expanded in x around inf

                            \[\leadsto x \cdot \color{blue}{\left(1 + \left(wj \cdot \left(wj \cdot \left(\frac{5}{2} - \frac{8}{3} \cdot wj\right) - 2\right) + \frac{{wj}^{2} \cdot \left(1 - wj\right)}{x}\right)\right)} \]
                          6. Step-by-step derivation
                            1. Applied rewrites99.7%

                              \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(-2.6666666666666665, wj, 2.5\right), wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot \color{blue}{x} \]
                            2. Taylor expanded in wj around 0

                              \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\frac{5}{2}, wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot x \]
                            3. Step-by-step derivation
                              1. Applied rewrites99.6%

                                \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(2.5, wj, -2\right), wj, \mathsf{fma}\left(\frac{1 - wj}{x}, wj \cdot wj, 1\right)\right) \cdot x \]
                              2. Step-by-step derivation
                                1. Applied rewrites99.6%

                                  \[\leadsto \mathsf{fma}\left(wj \cdot \left(\mathsf{fma}\left(2.5, wj, -2\right) + \frac{1 - wj}{x} \cdot wj\right), x, x\right) \]

                                if 5e302 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

                                1. Initial program 0.0%

                                  \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                2. Add Preprocessing
                                3. Taylor expanded in x around 0

                                  \[\leadsto wj - \color{blue}{\frac{wj \cdot e^{wj}}{e^{wj} + wj \cdot e^{wj}}} \]
                                4. Step-by-step derivation
                                  1. distribute-rgt1-inN/A

                                    \[\leadsto wj - \frac{wj \cdot e^{wj}}{\color{blue}{\left(wj + 1\right) \cdot e^{wj}}} \]
                                  2. +-commutativeN/A

                                    \[\leadsto wj - \frac{wj \cdot e^{wj}}{\color{blue}{\left(1 + wj\right)} \cdot e^{wj}} \]
                                  3. times-fracN/A

                                    \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj} \cdot \frac{e^{wj}}{e^{wj}}} \]
                                  4. *-inversesN/A

                                    \[\leadsto wj - \frac{wj}{1 + wj} \cdot \color{blue}{1} \]
                                  5. associate-*l/N/A

                                    \[\leadsto wj - \color{blue}{\frac{wj \cdot 1}{1 + wj}} \]
                                  6. *-rgt-identityN/A

                                    \[\leadsto wj - \frac{\color{blue}{wj}}{1 + wj} \]
                                  7. lower-/.f64N/A

                                    \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj}} \]
                                  8. lower-+.f6483.9

                                    \[\leadsto wj - \frac{wj}{\color{blue}{1 + wj}} \]
                                5. Applied rewrites83.9%

                                  \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj}} \]
                              3. Recombined 2 regimes into one program.
                              4. Final simplification98.9%

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

                              Alternative 5: 96.0% accurate, 0.9× speedup?

                              \[\begin{array}{l} \\ \begin{array}{l} t_0 := e^{wj} \cdot wj\\ \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;\mathsf{fma}\left(\left(1 - wj\right) \cdot wj, wj, x\right)\\ \mathbf{else}:\\ \;\;\;\;wj - \frac{wj}{1 + wj}\\ \end{array} \end{array} \]
                              (FPCore (wj x)
                               :precision binary64
                               (let* ((t_0 (* (exp wj) wj)))
                                 (if (<= (- wj (/ (- t_0 x) (+ (exp wj) t_0))) 5e+302)
                                   (fma (* (- 1.0 wj) wj) wj x)
                                   (- wj (/ wj (+ 1.0 wj))))))
                              double code(double wj, double x) {
                              	double t_0 = exp(wj) * wj;
                              	double tmp;
                              	if ((wj - ((t_0 - x) / (exp(wj) + t_0))) <= 5e+302) {
                              		tmp = fma(((1.0 - wj) * wj), wj, x);
                              	} else {
                              		tmp = wj - (wj / (1.0 + wj));
                              	}
                              	return tmp;
                              }
                              
                              function code(wj, x)
                              	t_0 = Float64(exp(wj) * wj)
                              	tmp = 0.0
                              	if (Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0))) <= 5e+302)
                              		tmp = fma(Float64(Float64(1.0 - wj) * wj), wj, x);
                              	else
                              		tmp = Float64(wj - Float64(wj / Float64(1.0 + wj)));
                              	end
                              	return tmp
                              end
                              
                              code[wj_, x_] := Block[{t$95$0 = N[(N[Exp[wj], $MachinePrecision] * wj), $MachinePrecision]}, If[LessEqual[N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 5e+302], N[(N[(N[(1.0 - wj), $MachinePrecision] * wj), $MachinePrecision] * wj + x), $MachinePrecision], N[(wj - N[(wj / N[(1.0 + wj), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]
                              
                              \begin{array}{l}
                              
                              \\
                              \begin{array}{l}
                              t_0 := e^{wj} \cdot wj\\
                              \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\
                              \;\;\;\;\mathsf{fma}\left(\left(1 - wj\right) \cdot wj, wj, x\right)\\
                              
                              \mathbf{else}:\\
                              \;\;\;\;wj - \frac{wj}{1 + wj}\\
                              
                              
                              \end{array}
                              \end{array}
                              
                              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))))) < 5e302

                                1. Initial program 74.6%

                                  \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                2. Add Preprocessing
                                3. Taylor expanded in wj around 0

                                  \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
                                4. Applied rewrites99.3%

                                  \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
                                5. Taylor expanded in x around 0

                                  \[\leadsto \mathsf{fma}\left(wj \cdot \left(1 - wj\right), wj, x\right) \]
                                6. Step-by-step derivation
                                  1. Applied rewrites99.1%

                                    \[\leadsto \mathsf{fma}\left(\left(1 - wj\right) \cdot wj, wj, x\right) \]

                                  if 5e302 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

                                  1. Initial program 0.0%

                                    \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                  2. Add Preprocessing
                                  3. Taylor expanded in x around 0

                                    \[\leadsto wj - \color{blue}{\frac{wj \cdot e^{wj}}{e^{wj} + wj \cdot e^{wj}}} \]
                                  4. Step-by-step derivation
                                    1. distribute-rgt1-inN/A

                                      \[\leadsto wj - \frac{wj \cdot e^{wj}}{\color{blue}{\left(wj + 1\right) \cdot e^{wj}}} \]
                                    2. +-commutativeN/A

                                      \[\leadsto wj - \frac{wj \cdot e^{wj}}{\color{blue}{\left(1 + wj\right)} \cdot e^{wj}} \]
                                    3. times-fracN/A

                                      \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj} \cdot \frac{e^{wj}}{e^{wj}}} \]
                                    4. *-inversesN/A

                                      \[\leadsto wj - \frac{wj}{1 + wj} \cdot \color{blue}{1} \]
                                    5. associate-*l/N/A

                                      \[\leadsto wj - \color{blue}{\frac{wj \cdot 1}{1 + wj}} \]
                                    6. *-rgt-identityN/A

                                      \[\leadsto wj - \frac{\color{blue}{wj}}{1 + wj} \]
                                    7. lower-/.f64N/A

                                      \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj}} \]
                                    8. lower-+.f6483.9

                                      \[\leadsto wj - \frac{wj}{\color{blue}{1 + wj}} \]
                                  5. Applied rewrites83.9%

                                    \[\leadsto wj - \color{blue}{\frac{wj}{1 + wj}} \]
                                7. Recombined 2 regimes into one program.
                                8. Final simplification98.4%

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

                                Alternative 6: 95.7% accurate, 0.9× speedup?

                                \[\begin{array}{l} \\ \begin{array}{l} t_0 := e^{wj} \cdot wj\\ \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;\mathsf{fma}\left(\left(1 - wj\right) \cdot wj, wj, x\right)\\ \mathbf{else}:\\ \;\;\;\;wj - 1\\ \end{array} \end{array} \]
                                (FPCore (wj x)
                                 :precision binary64
                                 (let* ((t_0 (* (exp wj) wj)))
                                   (if (<= (- wj (/ (- t_0 x) (+ (exp wj) t_0))) 5e+302)
                                     (fma (* (- 1.0 wj) wj) wj x)
                                     (- wj 1.0))))
                                double code(double wj, double x) {
                                	double t_0 = exp(wj) * wj;
                                	double tmp;
                                	if ((wj - ((t_0 - x) / (exp(wj) + t_0))) <= 5e+302) {
                                		tmp = fma(((1.0 - wj) * wj), wj, x);
                                	} else {
                                		tmp = wj - 1.0;
                                	}
                                	return tmp;
                                }
                                
                                function code(wj, x)
                                	t_0 = Float64(exp(wj) * wj)
                                	tmp = 0.0
                                	if (Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0))) <= 5e+302)
                                		tmp = fma(Float64(Float64(1.0 - wj) * wj), wj, x);
                                	else
                                		tmp = Float64(wj - 1.0);
                                	end
                                	return tmp
                                end
                                
                                code[wj_, x_] := Block[{t$95$0 = N[(N[Exp[wj], $MachinePrecision] * wj), $MachinePrecision]}, If[LessEqual[N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 5e+302], N[(N[(N[(1.0 - wj), $MachinePrecision] * wj), $MachinePrecision] * wj + x), $MachinePrecision], N[(wj - 1.0), $MachinePrecision]]]
                                
                                \begin{array}{l}
                                
                                \\
                                \begin{array}{l}
                                t_0 := e^{wj} \cdot wj\\
                                \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\
                                \;\;\;\;\mathsf{fma}\left(\left(1 - wj\right) \cdot wj, wj, x\right)\\
                                
                                \mathbf{else}:\\
                                \;\;\;\;wj - 1\\
                                
                                
                                \end{array}
                                \end{array}
                                
                                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))))) < 5e302

                                  1. Initial program 74.6%

                                    \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                  2. Add Preprocessing
                                  3. Taylor expanded in wj around 0

                                    \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
                                  4. Applied rewrites99.3%

                                    \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
                                  5. Taylor expanded in x around 0

                                    \[\leadsto \mathsf{fma}\left(wj \cdot \left(1 - wj\right), wj, x\right) \]
                                  6. Step-by-step derivation
                                    1. Applied rewrites99.1%

                                      \[\leadsto \mathsf{fma}\left(\left(1 - wj\right) \cdot wj, wj, x\right) \]

                                    if 5e302 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

                                    1. Initial program 0.0%

                                      \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                    2. Add Preprocessing
                                    3. Taylor expanded in wj around inf

                                      \[\leadsto wj - \color{blue}{1} \]
                                    4. Step-by-step derivation
                                      1. Applied rewrites64.5%

                                        \[\leadsto wj - \color{blue}{1} \]
                                    5. Recombined 2 regimes into one program.
                                    6. Final simplification97.5%

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

                                    Alternative 7: 14.5% accurate, 1.0× speedup?

                                    \[\begin{array}{l} \\ \begin{array}{l} t_0 := e^{wj} \cdot wj\\ \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;wj \cdot wj\\ \mathbf{else}:\\ \;\;\;\;wj - 1\\ \end{array} \end{array} \]
                                    (FPCore (wj x)
                                     :precision binary64
                                     (let* ((t_0 (* (exp wj) wj)))
                                       (if (<= (- wj (/ (- t_0 x) (+ (exp wj) t_0))) 5e+302)
                                         (* wj wj)
                                         (- wj 1.0))))
                                    double code(double wj, double x) {
                                    	double t_0 = exp(wj) * wj;
                                    	double tmp;
                                    	if ((wj - ((t_0 - x) / (exp(wj) + t_0))) <= 5e+302) {
                                    		tmp = wj * wj;
                                    	} else {
                                    		tmp = wj - 1.0;
                                    	}
                                    	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 = exp(wj) * wj
                                        if ((wj - ((t_0 - x) / (exp(wj) + t_0))) <= 5d+302) then
                                            tmp = wj * wj
                                        else
                                            tmp = wj - 1.0d0
                                        end if
                                        code = tmp
                                    end function
                                    
                                    public static double code(double wj, double x) {
                                    	double t_0 = Math.exp(wj) * wj;
                                    	double tmp;
                                    	if ((wj - ((t_0 - x) / (Math.exp(wj) + t_0))) <= 5e+302) {
                                    		tmp = wj * wj;
                                    	} else {
                                    		tmp = wj - 1.0;
                                    	}
                                    	return tmp;
                                    }
                                    
                                    def code(wj, x):
                                    	t_0 = math.exp(wj) * wj
                                    	tmp = 0
                                    	if (wj - ((t_0 - x) / (math.exp(wj) + t_0))) <= 5e+302:
                                    		tmp = wj * wj
                                    	else:
                                    		tmp = wj - 1.0
                                    	return tmp
                                    
                                    function code(wj, x)
                                    	t_0 = Float64(exp(wj) * wj)
                                    	tmp = 0.0
                                    	if (Float64(wj - Float64(Float64(t_0 - x) / Float64(exp(wj) + t_0))) <= 5e+302)
                                    		tmp = Float64(wj * wj);
                                    	else
                                    		tmp = Float64(wj - 1.0);
                                    	end
                                    	return tmp
                                    end
                                    
                                    function tmp_2 = code(wj, x)
                                    	t_0 = exp(wj) * wj;
                                    	tmp = 0.0;
                                    	if ((wj - ((t_0 - x) / (exp(wj) + t_0))) <= 5e+302)
                                    		tmp = wj * wj;
                                    	else
                                    		tmp = wj - 1.0;
                                    	end
                                    	tmp_2 = tmp;
                                    end
                                    
                                    code[wj_, x_] := Block[{t$95$0 = N[(N[Exp[wj], $MachinePrecision] * wj), $MachinePrecision]}, If[LessEqual[N[(wj - N[(N[(t$95$0 - x), $MachinePrecision] / N[(N[Exp[wj], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 5e+302], N[(wj * wj), $MachinePrecision], N[(wj - 1.0), $MachinePrecision]]]
                                    
                                    \begin{array}{l}
                                    
                                    \\
                                    \begin{array}{l}
                                    t_0 := e^{wj} \cdot wj\\
                                    \mathbf{if}\;wj - \frac{t\_0 - x}{e^{wj} + t\_0} \leq 5 \cdot 10^{+302}:\\
                                    \;\;\;\;wj \cdot wj\\
                                    
                                    \mathbf{else}:\\
                                    \;\;\;\;wj - 1\\
                                    
                                    
                                    \end{array}
                                    \end{array}
                                    
                                    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))))) < 5e302

                                      1. Initial program 74.6%

                                        \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                      2. Add Preprocessing
                                      3. Taylor expanded in wj around 0

                                        \[\leadsto \color{blue}{x + wj \cdot \left(wj \cdot \left(\left(1 + -1 \cdot \left(wj \cdot \left(1 + \left(-3 \cdot x + \left(-2 \cdot \left(-4 \cdot x + \frac{3}{2} \cdot x\right) + \frac{2}{3} \cdot x\right)\right)\right)\right)\right) - \left(-4 \cdot x + \frac{3}{2} \cdot x\right)\right) - 2 \cdot x\right)} \]
                                      4. Applied rewrites99.3%

                                        \[\leadsto \color{blue}{\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2.5, x, 1 - \mathsf{fma}\left(0.6666666666666666, x, \mathsf{fma}\left(2, x, 1\right)\right) \cdot wj\right), wj, -2 \cdot x\right), wj, x\right)} \]
                                      5. Taylor expanded in x around 0

                                        \[\leadsto {wj}^{2} \cdot \color{blue}{\left(1 - wj\right)} \]
                                      6. Step-by-step derivation
                                        1. Applied rewrites20.7%

                                          \[\leadsto \left(\left(1 - wj\right) \cdot wj\right) \cdot \color{blue}{wj} \]
                                        2. Taylor expanded in wj around 0

                                          \[\leadsto {wj}^{2} \]
                                        3. Step-by-step derivation
                                          1. Applied rewrites20.0%

                                            \[\leadsto wj \cdot wj \]

                                          if 5e302 < (-.f64 wj (/.f64 (-.f64 (*.f64 wj (exp.f64 wj)) x) (+.f64 (exp.f64 wj) (*.f64 wj (exp.f64 wj)))))

                                          1. Initial program 0.0%

                                            \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                          2. Add Preprocessing
                                          3. Taylor expanded in wj around inf

                                            \[\leadsto wj - \color{blue}{1} \]
                                          4. Step-by-step derivation
                                            1. Applied rewrites64.5%

                                              \[\leadsto wj - \color{blue}{1} \]
                                          5. Recombined 2 regimes into one program.
                                          6. Final simplification22.1%

                                            \[\leadsto \begin{array}{l} \mathbf{if}\;wj - \frac{e^{wj} \cdot wj - x}{e^{wj} + e^{wj} \cdot wj} \leq 5 \cdot 10^{+302}:\\ \;\;\;\;wj \cdot wj\\ \mathbf{else}:\\ \;\;\;\;wj - 1\\ \end{array} \]
                                          7. Add Preprocessing

                                          Alternative 8: 4.1% accurate, 82.8× speedup?

                                          \[\begin{array}{l} \\ wj - 1 \end{array} \]
                                          (FPCore (wj x) :precision binary64 (- wj 1.0))
                                          double code(double wj, double x) {
                                          	return wj - 1.0;
                                          }
                                          
                                          real(8) function code(wj, x)
                                              real(8), intent (in) :: wj
                                              real(8), intent (in) :: x
                                              code = wj - 1.0d0
                                          end function
                                          
                                          public static double code(double wj, double x) {
                                          	return wj - 1.0;
                                          }
                                          
                                          def code(wj, x):
                                          	return wj - 1.0
                                          
                                          function code(wj, x)
                                          	return Float64(wj - 1.0)
                                          end
                                          
                                          function tmp = code(wj, x)
                                          	tmp = wj - 1.0;
                                          end
                                          
                                          code[wj_, x_] := N[(wj - 1.0), $MachinePrecision]
                                          
                                          \begin{array}{l}
                                          
                                          \\
                                          wj - 1
                                          \end{array}
                                          
                                          Derivation
                                          1. Initial program 71.1%

                                            \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}} \]
                                          2. Add Preprocessing
                                          3. Taylor expanded in wj around inf

                                            \[\leadsto wj - \color{blue}{1} \]
                                          4. Step-by-step derivation
                                            1. Applied rewrites6.0%

                                              \[\leadsto wj - \color{blue}{1} \]
                                            2. Add Preprocessing

                                            Developer Target 1: 79.2% accurate, 1.4× speedup?

                                            \[\begin{array}{l} \\ wj - \left(\frac{wj}{wj + 1} - \frac{x}{e^{wj} + wj \cdot e^{wj}}\right) \end{array} \]
                                            (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}
                                            

                                            Reproduce

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