(FPCore (a b c d)
:precision binary64
(/ (- (* b c) (* a d)) (+ (* c c) (* d d))))
↓
(FPCore (a b c d)
:precision binary64
(let* ((t_0 (/ (* d a) (hypot d c)))
(t_1 (- (/ (* (/ (- b) d) c) (- d)) (/ a d)))
(t_2 (/ 1.0 (hypot d c))))
(if (<= d -2.5355883399872014e+249)
t_1
(if (<= d -9.766740365443902e+235)
(/ (- b (* d (/ a c))) c)
(if (<= d -4.939221699676908e+146)
(- (* c (/ (/ b d) d)) (/ a d))
(if (<= d 2.926004778278525e+141)
(+
(fma (* b (/ c (hypot d c))) t_2 (* t_0 (/ -1.0 (hypot d c))))
(fma (/ (* d (- a)) (hypot d c)) t_2 (* t_2 t_0)))
t_1))))))
double code(double a, double b, double c, double d) {
return ((b * c) - (a * d)) / ((c * c) + (d * d));
}
\[\begin{array}{l}
\mathbf{if}\;\left|d\right| < \left|c\right|:\\
\;\;\;\;\frac{b - a \cdot \frac{d}{c}}{c + d \cdot \frac{d}{c}}\\
\mathbf{else}:\\
\;\;\;\;\frac{\left(-a\right) + b \cdot \frac{c}{d}}{d + c \cdot \frac{c}{d}}\\
\end{array}
\]
Derivation
Split input into 4 regimes
if d < -2.5355883399872014e249 or 2.9260047782785248e141 < d
Initial program 42.5
\[\frac{b \cdot c - a \cdot d}{c \cdot c + d \cdot d}
\]
Simplified42.5
\[\leadsto \color{blue}{\frac{\mathsf{fma}\left(d, -a, b \cdot c\right)}{\mathsf{fma}\left(c, c, d \cdot d\right)}}
\]
Proof
(/.f64 (fma.f64 d (neg.f64 a) (*.f64 b c)) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (Rewrite<= fma-def_binary64 (+.f64 (*.f64 d (neg.f64 a)) (*.f64 b c))) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (+.f64 (Rewrite<= *-commutative_binary64 (*.f64 (neg.f64 a) d)) (*.f64 b c)) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 b c) (*.f64 (neg.f64 a) d))) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (Rewrite<= cancel-sign-sub-inv_binary64 (-.f64 (*.f64 b c) (*.f64 a d))) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (-.f64 (*.f64 b c) (*.f64 a d)) (Rewrite<= fma-def_binary64 (+.f64 (*.f64 c c) (*.f64 d d)))): 1 points increase in error, 0 points decrease in error
(-.f64 (*.f64 c (/.f64 (/.f64 b d) d)) (/.f64 a d)): 0 points increase in error, 0 points decrease in error
(-.f64 (*.f64 c (Rewrite<= associate-/r*_binary64 (/.f64 b (*.f64 d d)))) (/.f64 a d)): 32 points increase in error, 6 points decrease in error
(-.f64 (*.f64 c (/.f64 b (Rewrite<= unpow2_binary64 (pow.f64 d 2)))) (/.f64 a d)): 0 points increase in error, 0 points decrease in error
(-.f64 (Rewrite=> associate-*r/_binary64 (/.f64 (*.f64 c b) (pow.f64 d 2))) (/.f64 a d)): 24 points increase in error, 8 points decrease in error
(Rewrite<= unsub-neg_binary64 (+.f64 (/.f64 (*.f64 c b) (pow.f64 d 2)) (neg.f64 (/.f64 a d)))): 0 points increase in error, 0 points decrease in error
(+.f64 (/.f64 (*.f64 c b) (pow.f64 d 2)) (Rewrite<= mul-1-neg_binary64 (*.f64 -1 (/.f64 a d)))): 0 points increase in error, 0 points decrease in error
(Rewrite<= +-commutative_binary64 (+.f64 (*.f64 -1 (/.f64 a d)) (/.f64 (*.f64 c b) (pow.f64 d 2)))): 0 points increase in error, 0 points decrease in error
\[\leadsto \color{blue}{\frac{b}{c} - d \cdot \frac{\frac{a}{c}}{c}}
\]
Proof
(-.f64 (/.f64 b c) (*.f64 d (/.f64 (/.f64 a c) c))): 0 points increase in error, 0 points decrease in error
(-.f64 (/.f64 b c) (*.f64 d (Rewrite<= associate-/r*_binary64 (/.f64 a (*.f64 c c))))): 25 points increase in error, 9 points decrease in error
(-.f64 (/.f64 b c) (*.f64 d (/.f64 a (Rewrite<= unpow2_binary64 (pow.f64 c 2))))): 0 points increase in error, 0 points decrease in error
(-.f64 (/.f64 b c) (Rewrite<= *-commutative_binary64 (*.f64 (/.f64 a (pow.f64 c 2)) d))): 0 points increase in error, 0 points decrease in error
(-.f64 (/.f64 b c) (Rewrite=> associate-*l/_binary64 (/.f64 (*.f64 a d) (pow.f64 c 2)))): 34 points increase in error, 6 points decrease in error
(Rewrite<= unsub-neg_binary64 (+.f64 (/.f64 b c) (neg.f64 (/.f64 (*.f64 a d) (pow.f64 c 2))))): 0 points increase in error, 0 points decrease in error
(+.f64 (/.f64 b c) (Rewrite<= mul-1-neg_binary64 (*.f64 -1 (/.f64 (*.f64 a d) (pow.f64 c 2))))): 0 points increase in error, 0 points decrease in error
(Rewrite<= +-commutative_binary64 (+.f64 (*.f64 -1 (/.f64 (*.f64 a d) (pow.f64 c 2))) (/.f64 b c))): 0 points increase in error, 0 points decrease in error
Applied egg-rr55.1
\[\leadsto \color{blue}{\frac{b - d \cdot \frac{a}{c}}{c}}
\]
if -9.76674036544390209e235 < d < -4.9392216996769081e146
Initial program 47.7
\[\frac{b \cdot c - a \cdot d}{c \cdot c + d \cdot d}
\]
Simplified47.7
\[\leadsto \color{blue}{\frac{\mathsf{fma}\left(d, -a, b \cdot c\right)}{\mathsf{fma}\left(c, c, d \cdot d\right)}}
\]
Proof
(/.f64 (fma.f64 d (neg.f64 a) (*.f64 b c)) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (Rewrite<= fma-def_binary64 (+.f64 (*.f64 d (neg.f64 a)) (*.f64 b c))) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (+.f64 (Rewrite<= *-commutative_binary64 (*.f64 (neg.f64 a) d)) (*.f64 b c)) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 b c) (*.f64 (neg.f64 a) d))) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (Rewrite<= cancel-sign-sub-inv_binary64 (-.f64 (*.f64 b c) (*.f64 a d))) (fma.f64 c c (*.f64 d d))): 0 points increase in error, 0 points decrease in error
(/.f64 (-.f64 (*.f64 b c) (*.f64 a d)) (Rewrite<= fma-def_binary64 (+.f64 (*.f64 c c) (*.f64 d d)))): 1 points increase in error, 0 points decrease in error
(-.f64 (*.f64 c (/.f64 (/.f64 b d) d)) (/.f64 a d)): 0 points increase in error, 0 points decrease in error
(-.f64 (*.f64 c (Rewrite<= associate-/r*_binary64 (/.f64 b (*.f64 d d)))) (/.f64 a d)): 32 points increase in error, 6 points decrease in error
(-.f64 (*.f64 c (/.f64 b (Rewrite<= unpow2_binary64 (pow.f64 d 2)))) (/.f64 a d)): 0 points increase in error, 0 points decrease in error
(-.f64 (Rewrite=> associate-*r/_binary64 (/.f64 (*.f64 c b) (pow.f64 d 2))) (/.f64 a d)): 24 points increase in error, 8 points decrease in error
(Rewrite<= unsub-neg_binary64 (+.f64 (/.f64 (*.f64 c b) (pow.f64 d 2)) (neg.f64 (/.f64 a d)))): 0 points increase in error, 0 points decrease in error
(+.f64 (/.f64 (*.f64 c b) (pow.f64 d 2)) (Rewrite<= mul-1-neg_binary64 (*.f64 -1 (/.f64 a d)))): 0 points increase in error, 0 points decrease in error
(Rewrite<= +-commutative_binary64 (+.f64 (*.f64 -1 (/.f64 a d)) (/.f64 (*.f64 c b) (pow.f64 d 2)))): 0 points increase in error, 0 points decrease in error
if -4.9392216996769081e146 < d < 2.9260047782785248e141
Initial program 19.6
\[\frac{b \cdot c - a \cdot d}{c \cdot c + d \cdot d}
\]
herbie shell --seed 2022308
(FPCore (a b c d)
:name "Complex division, imag part"
:precision binary64
:herbie-target
(if (< (fabs d) (fabs c)) (/ (- b (* a (/ d c))) (+ c (* d (/ d c)))) (/ (+ (- a) (* b (/ c d))) (+ d (* c (/ c d)))))
(/ (- (* b c) (* a d)) (+ (* c c) (* d d))))