Jmat.Real.lambertw, newton loop step

Average Error: 13.6 → 0.3

Time: 1.4m

Precision: 64

• could not determine a ground truth for program body

  1. wj = -7.594672547429651e+73
  2. x = -2.4734580902544025e+122

Math

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

↓

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

C

double f(double wj, double x) {
        double r25031110 = wj;
        double r25031111 = exp(r25031110);
        double r25031112 = r25031110 * r25031111;
        double r25031113 = x;
        double r25031114 = r25031112 - r25031113;
        double r25031115 = r25031111 + r25031112;
        double r25031116 = r25031114 / r25031115;
        double r25031117 = r25031110 - r25031116;
        return r25031117;
}

↓

double f(double wj, double x) {
        double r25031118 = wj;
        double r25031119 = 7.273522667070596e-05;
        bool r25031120 = r25031118 <= r25031119;
        double r25031121 = r25031118 * r25031118;
        double r25031122 = r25031121 - r25031118;
        double r25031123 = fma(r25031121, r25031122, r25031121);
        double r25031124 = sqrt(r25031123);
        double r25031125 = r25031124 * r25031124;
        double r25031126 = x;
        double r25031127 = exp(r25031118);
        double r25031128 = r25031127 * r25031118;
        double r25031129 = r25031128 + r25031127;
        double r25031130 = r25031126 / r25031129;
        double r25031131 = r25031125 + r25031130;
        double r25031132 = r25031126 / r25031127;
        double r25031133 = r25031118 - r25031132;
        double r25031134 = 1.0;
        double r25031135 = r25031134 + r25031118;
        double r25031136 = r25031133 / r25031135;
        double r25031137 = r25031118 - r25031136;
        double r25031138 = r25031120 ? r25031131 : r25031137;
        return r25031138;
}

Error

Bits error versus wj

Bits error versus x

Target

Original	13.6
Target	13.1
Herbie	0.3

\[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 < 7.273522667070596e-05

   1. Initial program 13.4

      \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}}\]
   2. Using strategy rm

   3. Applied div-sub 13.4

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

   4. Applied associate--r- 6.8

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

   5. Taylor expanded around 0 0.3

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

   6. Simplified 0.3

      \[\leadsto \color{blue}{\mathsf{fma}\left(\left(wj \cdot wj\right), \left(wj \cdot wj - wj\right), \left(wj \cdot wj\right)\right)} + \frac{x}{e^{wj} + wj \cdot e^{wj}}\]
   7. Using strategy rm

   8. Applied add-sqr-sqrt 0.3

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

   if 7.273522667070596e-05 < wj

   1. Initial program 23.5

      \[wj - \frac{wj \cdot e^{wj} - x}{e^{wj} + wj \cdot e^{wj}}\]
   2. Using strategy rm

   3. Applied distribute-rgt1-in 23.5

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

   4. Applied *-un-lft-identity 23.5

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

   5. Applied times-frac 23.5

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

   6. Applied add-cube-cbrt 24.2

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

   7. Applied prod-diff 24.2

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

   8. Simplified 23.5

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

   9. Simplified 1.2

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

4. Final simplification 0.3

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

Reproduce

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

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

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