Jmat.Real.lambertw, newton loop step

Average Error: 13.4 → 0.9

Time: 4.8s

Precision: 64

• could not determine a ground truth for program body

  1. wj = -5.293892094996653e+190
  2. x = 5.81386847535435e+28

Math

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

↓

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

C

double f(double wj, double x) {
        double r185995 = wj;
        double r185996 = exp(r185995);
        double r185997 = r185995 * r185996;
        double r185998 = x;
        double r185999 = r185997 - r185998;
        double r186000 = r185996 + r185997;
        double r186001 = r185999 / r186000;
        double r186002 = r185995 - r186001;
        return r186002;
}

↓

double f(double wj, double x) {
        double r186003 = wj;
        double r186004 = 9.60333635941626e-09;
        bool r186005 = r186003 <= r186004;
        double r186006 = 1.0;
        double r186007 = x;
        double r186008 = fma(r186003, r186003, r186007);
        double r186009 = r186006 * r186008;
        double r186010 = 2.0;
        double r186011 = r186003 * r186007;
        double r186012 = r186010 * r186011;
        double r186013 = r186009 - r186012;
        double r186014 = r186003 + r186006;
        double r186015 = r186007 / r186014;
        double r186016 = exp(r186003);
        double r186017 = r186015 / r186016;
        double r186018 = r186003 / r186014;
        double r186019 = r186003 - r186018;
        double r186020 = r186017 + r186019;
        double r186021 = r186005 ? r186013 : r186020;
        return r186021;
}

Error

Bits error versus wj

Bits error versus x

Target

Original	13.4
Target	12.7
Herbie	0.9

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

Derivation

1. Split input into 2 regimes

2. if wj < 9.60333635941626e-09

   1. Initial program 13.0

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

   2. Simplified 13.0

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

   3. Taylor expanded around 0 0.8

      \[\leadsto \color{blue}{\left(x + {wj}^{2}\right) - 2 \cdot \left(wj \cdot x\right)}\]
   4. Using strategy rm

   5. Applied *-un-lft-identity 0.8

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

   6. Applied *-un-lft-identity 0.8

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

   7. Applied distribute-lft-out 0.8

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

   8. Simplified 0.8

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

   if 9.60333635941626e-09 < wj

   1. Initial program 26.2

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

   2. Simplified 3.2

      \[\leadsto \color{blue}{\left(\frac{\frac{x}{wj + 1}}{e^{wj}} + wj\right) - \frac{wj}{wj + 1}}\]
   3. Using strategy rm

   4. Applied associate--l+ 3.2

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

4. Final simplification 0.9

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

Reproduce

herbie shell --seed 2020027 +o rules:numerics
(FPCore (wj x)
  :name "Jmat.Real.lambertw, newton loop step"
  :precision binary64

  :herbie-target
  (- wj (- (/ wj (+ wj 1)) (/ x (+ (exp wj) (* wj (exp wj))))))

  (- wj (/ (- (* wj (exp wj)) x) (+ (exp wj) (* wj (exp wj))))))