Average Error: 59.3 → 34.5
Time: 33.6s
Precision: 64
\[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right)\]
\[\begin{array}{l} \mathbf{if}\;d \le 1.028014946374977112793459567345674978253 \cdot 10^{-41} \lor \neg \left(d \le 1.879132275263087348510262122271649682884 \cdot 10^{52}\right):\\ \;\;\;\;\mathsf{log1p}\left(\mathsf{expm1}\left(0\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{c0}{2 \cdot w} \cdot \left(2 \cdot \frac{{d}^{2} \cdot c0}{w \cdot \left({D}^{2} \cdot h\right)}\right)\\ \end{array}\]
\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right)
\begin{array}{l}
\mathbf{if}\;d \le 1.028014946374977112793459567345674978253 \cdot 10^{-41} \lor \neg \left(d \le 1.879132275263087348510262122271649682884 \cdot 10^{52}\right):\\
\;\;\;\;\mathsf{log1p}\left(\mathsf{expm1}\left(0\right)\right)\\

\mathbf{else}:\\
\;\;\;\;\frac{c0}{2 \cdot w} \cdot \left(2 \cdot \frac{{d}^{2} \cdot c0}{w \cdot \left({D}^{2} \cdot h\right)}\right)\\

\end{array}
double f(double c0, double w, double h, double D, double d, double M) {
        double r139439 = c0;
        double r139440 = 2.0;
        double r139441 = w;
        double r139442 = r139440 * r139441;
        double r139443 = r139439 / r139442;
        double r139444 = d;
        double r139445 = r139444 * r139444;
        double r139446 = r139439 * r139445;
        double r139447 = h;
        double r139448 = r139441 * r139447;
        double r139449 = D;
        double r139450 = r139449 * r139449;
        double r139451 = r139448 * r139450;
        double r139452 = r139446 / r139451;
        double r139453 = r139452 * r139452;
        double r139454 = M;
        double r139455 = r139454 * r139454;
        double r139456 = r139453 - r139455;
        double r139457 = sqrt(r139456);
        double r139458 = r139452 + r139457;
        double r139459 = r139443 * r139458;
        return r139459;
}

double f(double c0, double w, double h, double D, double d, double __attribute__((unused)) M) {
        double r139460 = d;
        double r139461 = 1.0280149463749771e-41;
        bool r139462 = r139460 <= r139461;
        double r139463 = 1.8791322752630873e+52;
        bool r139464 = r139460 <= r139463;
        double r139465 = !r139464;
        bool r139466 = r139462 || r139465;
        double r139467 = 0.0;
        double r139468 = expm1(r139467);
        double r139469 = log1p(r139468);
        double r139470 = c0;
        double r139471 = 2.0;
        double r139472 = w;
        double r139473 = r139471 * r139472;
        double r139474 = r139470 / r139473;
        double r139475 = 2.0;
        double r139476 = pow(r139460, r139475);
        double r139477 = r139476 * r139470;
        double r139478 = D;
        double r139479 = pow(r139478, r139475);
        double r139480 = h;
        double r139481 = r139479 * r139480;
        double r139482 = r139472 * r139481;
        double r139483 = r139477 / r139482;
        double r139484 = r139475 * r139483;
        double r139485 = r139474 * r139484;
        double r139486 = r139466 ? r139469 : r139485;
        return r139486;
}

Error

Bits error versus c0

Bits error versus w

Bits error versus h

Bits error versus D

Bits error versus d

Bits error versus M

Try it out

Your Program's Arguments

Results

Enter valid numbers for all inputs

Derivation

  1. Split input into 2 regimes
  2. if d < 1.0280149463749771e-41 or 1.8791322752630873e+52 < d

    1. Initial program 59.9

      \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right)\]
    2. Taylor expanded around inf 34.6

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{0}\]
    3. Using strategy rm
    4. Applied log1p-expm1-u34.6

      \[\leadsto \color{blue}{\mathsf{log1p}\left(\mathsf{expm1}\left(\frac{c0}{2 \cdot w} \cdot 0\right)\right)}\]
    5. Simplified32.7

      \[\leadsto \mathsf{log1p}\left(\color{blue}{\mathsf{expm1}\left(0\right)}\right)\]

    if 1.0280149463749771e-41 < d < 1.8791322752630873e+52

    1. Initial program 52.7

      \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right)\]
    2. Using strategy rm
    3. Applied associate-/l*53.8

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \left(\color{blue}{\frac{c0}{\frac{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}{d \cdot d}}} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right)\]
    4. Taylor expanded around inf 54.7

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\left(2 \cdot \frac{{d}^{2} \cdot c0}{w \cdot \left({D}^{2} \cdot h\right)}\right)}\]
  3. Recombined 2 regimes into one program.
  4. Final simplification34.5

    \[\leadsto \begin{array}{l} \mathbf{if}\;d \le 1.028014946374977112793459567345674978253 \cdot 10^{-41} \lor \neg \left(d \le 1.879132275263087348510262122271649682884 \cdot 10^{52}\right):\\ \;\;\;\;\mathsf{log1p}\left(\mathsf{expm1}\left(0\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{c0}{2 \cdot w} \cdot \left(2 \cdot \frac{{d}^{2} \cdot c0}{w \cdot \left({D}^{2} \cdot h\right)}\right)\\ \end{array}\]

Reproduce

herbie shell --seed 2019323 +o rules:numerics
(FPCore (c0 w h D d M)
  :name "Henrywood and Agarwal, Equation (13)"
  :precision binary64
  (* (/ c0 (* 2 w)) (+ (/ (* c0 (* d d)) (* (* w h) (* D D))) (sqrt (- (* (/ (* c0 (* d d)) (* (* w h) (* D D))) (/ (* c0 (* d d)) (* (* w h) (* D D)))) (* M M))))))