Jmat.Real.lambertw, newton loop step

Average Error: 14.4 → 0.4

Time: 6.9s

Precision: binary64

Math

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

↓

\[\begin{array}{l} t_0 := \frac{x}{e^{wj}} - wj\\ \mathbf{if}\;wj \leq -1.3031958808623259 \cdot 10^{-6}:\\ \;\;\;\;wj + t_0 \cdot \frac{1}{wj + 1}\\ \mathbf{elif}\;wj \leq 7.73535415808421 \cdot 10^{-12}:\\ \;\;\;\;\mathsf{fma}\left(wj, wj, x\right) + x \cdot \left(wj \cdot \mathsf{fma}\left(2.5, wj, -2\right)\right)\\ \mathbf{else}:\\ \;\;\;\;wj + \frac{t_0}{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 (- (/ x (exp wj)) wj)))
   (if (<= wj -1.3031958808623259e-6)
     (+ wj (* t_0 (/ 1.0 (+ wj 1.0))))
     (if (<= wj 7.73535415808421e-12)
       (+ (fma wj wj x) (* x (* wj (fma 2.5 wj -2.0))))
       (+ wj (/ t_0 (+ 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 = (x / exp(wj)) - wj;
	double tmp;
	if (wj <= -1.3031958808623259e-6) {
		tmp = wj + (t_0 * (1.0 / (wj + 1.0)));
	} else if (wj <= 7.73535415808421e-12) {
		tmp = fma(wj, wj, x) + (x * (wj * fma(2.5, wj, -2.0)));
	} else {
		tmp = wj + (t_0 / (wj + 1.0));
	}
	return tmp;
}

Error

Bits error versus wj

Bits error versus x

Target

Original	14.4
Target	13.7
Herbie	0.4

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

Derivation

1. Split input into 3 regimes

2. if wj < -1.30319588086232587e-6

   1. Initial program 29.3

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

   2. Simplified 1.9

      \[\leadsto \color{blue}{wj + \frac{\frac{x}{e^{wj}} - wj}{wj + 1}} \]

   3. Applied div-inv_binary64 2.1

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

   if -1.30319588086232587e-6 < wj < 7.7353541580842095e-12

   1. Initial program 13.6

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

   2. Simplified 13.6

      \[\leadsto \color{blue}{wj + \frac{\frac{x}{e^{wj}} - wj}{wj + 1}} \]

   3. Taylor expanded in wj around 0 0.2

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

   4. Simplified 0.2

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

   5. Taylor expanded in wj around 0 0.2

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

   6. Simplified 0.2

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

   if 7.7353541580842095e-12 < wj

   1. Initial program 27.5

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

   2. Simplified 5.2

      \[\leadsto \color{blue}{wj + \frac{\frac{x}{e^{wj}} - wj}{wj + 1}} \]
3. Recombined 3 regimes into one program.

4. Final simplification 0.4

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

Reproduce

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