Complex division, imag part

Average Error: 26.0 → 5.9

Time: 10.8s

Precision: binary64

Cost: 33480

Math

\[\frac{b \cdot c - a \cdot d}{c \cdot c + d \cdot d} \]

↓

\[\begin{array}{l} t_0 := \frac{b}{c} - \frac{a}{c} \cdot \frac{d}{c}\\ \mathbf{if}\;c \leq -1.25 \cdot 10^{+113}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;c \leq 1.7 \cdot 10^{+115}:\\ \;\;\;\;\frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \left(\frac{{\left(\sqrt[3]{c \cdot b}\right)}^{3}}{\mathsf{hypot}\left(c, d\right)} - d \cdot \frac{a}{\mathsf{hypot}\left(c, d\right)}\right)\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

FPCore

(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 (- (/ b c) (* (/ a c) (/ d c)))))
   (if (<= c -1.25e+113)
     t_0
     (if (<= c 1.7e+115)
       (*
        (/ 1.0 (hypot c d))
        (- (/ (pow (cbrt (* c b)) 3.0) (hypot c d)) (* d (/ a (hypot c d)))))
       t_0))))

C

double code(double a, double b, double c, double d) {
	return ((b * c) - (a * d)) / ((c * c) + (d * d));
}

↓

double code(double a, double b, double c, double d) {
	double t_0 = (b / c) - ((a / c) * (d / c));
	double tmp;
	if (c <= -1.25e+113) {
		tmp = t_0;
	} else if (c <= 1.7e+115) {
		tmp = (1.0 / hypot(c, d)) * ((pow(cbrt((c * b)), 3.0) / hypot(c, d)) - (d * (a / hypot(c, d))));
	} else {
		tmp = t_0;
	}
	return tmp;
}

Java

public static double code(double a, double b, double c, double d) {
	return ((b * c) - (a * d)) / ((c * c) + (d * d));
}

↓

public static double code(double a, double b, double c, double d) {
	double t_0 = (b / c) - ((a / c) * (d / c));
	double tmp;
	if (c <= -1.25e+113) {
		tmp = t_0;
	} else if (c <= 1.7e+115) {
		tmp = (1.0 / Math.hypot(c, d)) * ((Math.pow(Math.cbrt((c * b)), 3.0) / Math.hypot(c, d)) - (d * (a / Math.hypot(c, d))));
	} else {
		tmp = t_0;
	}
	return tmp;
}

Julia

function code(a, b, c, d)
	return Float64(Float64(Float64(b * c) - Float64(a * d)) / Float64(Float64(c * c) + Float64(d * d)))
end

↓

function code(a, b, c, d)
	t_0 = Float64(Float64(b / c) - Float64(Float64(a / c) * Float64(d / c)))
	tmp = 0.0
	if (c <= -1.25e+113)
		tmp = t_0;
	elseif (c <= 1.7e+115)
		tmp = Float64(Float64(1.0 / hypot(c, d)) * Float64(Float64((cbrt(Float64(c * b)) ^ 3.0) / hypot(c, d)) - Float64(d * Float64(a / hypot(c, d)))));
	else
		tmp = t_0;
	end
	return tmp
end

Wolfram

code[a_, b_, c_, d_] := N[(N[(N[(b * c), $MachinePrecision] - N[(a * d), $MachinePrecision]), $MachinePrecision] / N[(N[(c * c), $MachinePrecision] + N[(d * d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

↓

code[a_, b_, c_, d_] := Block[{t$95$0 = N[(N[(b / c), $MachinePrecision] - N[(N[(a / c), $MachinePrecision] * N[(d / c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[c, -1.25e+113], t$95$0, If[LessEqual[c, 1.7e+115], N[(N[(1.0 / N[Sqrt[c ^ 2 + d ^ 2], $MachinePrecision]), $MachinePrecision] * N[(N[(N[Power[N[Power[N[(c * b), $MachinePrecision], 1/3], $MachinePrecision], 3.0], $MachinePrecision] / N[Sqrt[c ^ 2 + d ^ 2], $MachinePrecision]), $MachinePrecision] - N[(d * N[(a / N[Sqrt[c ^ 2 + d ^ 2], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], t$95$0]]]

Target

Original	26.0
Target	0.4
Herbie	5.9

\[\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

1. Split input into 2 regimes

2. if c < -1.25e113 or 1.7e115 < c

   1. Initial program 40.9

      \[\frac{b \cdot c - a \cdot d}{c \cdot c + d \cdot d} \]

   2. Applied egg-rr 26.7

      \[\leadsto \color{blue}{\frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \frac{b \cdot c - a \cdot d}{\mathsf{hypot}\left(c, d\right)}} \]

   3. Taylor expanded in c around -inf 31.0

      \[\leadsto \frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \color{blue}{\left(\frac{a \cdot d}{c} + -1 \cdot b\right)} \]

   4. Simplified 28.9

      \[\leadsto \frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \color{blue}{\left(\frac{a}{c} \cdot d - b\right)} \]
      Proof
      (-.f64 (*.f64 (/.f64 a c) d) b): 0 points increase in error, 0 points decrease in error
      (-.f64 (Rewrite<= associate-/r/_binary64 (/.f64 a (/.f64 c d))) b): 19 points increase in error, 23 points decrease in error
      (-.f64 (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 a d) c)) b): 26 points increase in error, 17 points decrease in error
      (Rewrite<= unsub-neg_binary64 (+.f64 (/.f64 (*.f64 a d) c) (neg.f64 b))): 0 points increase in error, 0 points decrease in error
      (+.f64 (/.f64 (*.f64 a d) c) (Rewrite<= mul-1-neg_binary64 (*.f64 -1 b))): 0 points increase in error, 0 points decrease in error

   5. Taylor expanded in c around -inf 16.9

      \[\leadsto \color{blue}{-1 \cdot \frac{a \cdot d}{{c}^{2}} + \frac{b}{c}} \]

   6. Simplified 10.0

      \[\leadsto \color{blue}{\frac{b}{c} - \frac{a}{c} \cdot \frac{d}{c}} \]
      Proof
      (-.f64 (/.f64 b c) (*.f64 (/.f64 a c) (/.f64 d c))): 0 points increase in error, 0 points decrease in error
      (-.f64 (/.f64 b c) (Rewrite<= times-frac_binary64 (/.f64 (*.f64 a d) (*.f64 c c)))): 31 points increase in error, 16 points decrease in error
      (-.f64 (/.f64 b c) (/.f64 (*.f64 a d) (Rewrite<= unpow2_binary64 (pow.f64 c 2)))): 0 points increase in error, 0 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

   if -1.25e113 < c < 1.7e115

   1. Initial program 18.6

      \[\frac{b \cdot c - a \cdot d}{c \cdot c + d \cdot d} \]

   2. Applied egg-rr 11.6

      \[\leadsto \color{blue}{\frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \frac{b \cdot c - a \cdot d}{\mathsf{hypot}\left(c, d\right)}} \]

   3. Applied egg-rr 11.9

      \[\leadsto \frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \color{blue}{\mathsf{fma}\left(\frac{{\left(\sqrt[3]{b \cdot c}\right)}^{2}}{1}, \frac{\sqrt[3]{b \cdot c}}{\mathsf{hypot}\left(c, d\right)}, -\frac{a \cdot d}{\mathsf{hypot}\left(c, d\right)}\right)} \]

   4. Simplified 3.9

      \[\leadsto \frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \color{blue}{\left(\frac{{\left(\sqrt[3]{c \cdot b}\right)}^{3}}{\mathsf{hypot}\left(c, d\right)} - \frac{a}{\mathsf{hypot}\left(c, d\right)} \cdot d\right)} \]
      Proof
      (-.f64 (/.f64 (pow.f64 (cbrt.f64 (*.f64 c b)) 3) (hypot.f64 c d)) (*.f64 (/.f64 a (hypot.f64 c d)) d)): 0 points increase in error, 0 points decrease in error
      (-.f64 (/.f64 (pow.f64 (cbrt.f64 (Rewrite<= *-commutative_binary64 (*.f64 b c))) 3) (hypot.f64 c d)) (*.f64 (/.f64 a (hypot.f64 c d)) d)): 0 points increase in error, 0 points decrease in error
      (-.f64 (/.f64 (pow.f64 (cbrt.f64 (*.f64 b c)) (Rewrite<= metadata-eval (+.f64 2 1))) (hypot.f64 c d)) (*.f64 (/.f64 a (hypot.f64 c d)) d)): 0 points increase in error, 0 points decrease in error
      (-.f64 (/.f64 (Rewrite<= pow-plus_binary64 (*.f64 (pow.f64 (cbrt.f64 (*.f64 b c)) 2) (cbrt.f64 (*.f64 b c)))) (hypot.f64 c d)) (*.f64 (/.f64 a (hypot.f64 c d)) d)): 11 points increase in error, 11 points decrease in error
      (-.f64 (/.f64 (*.f64 (Rewrite<= /-rgt-identity_binary64 (/.f64 (pow.f64 (cbrt.f64 (*.f64 b c)) 2) 1)) (cbrt.f64 (*.f64 b c))) (hypot.f64 c d)) (*.f64 (/.f64 a (hypot.f64 c d)) d)): 0 points increase in error, 0 points decrease in error
      (-.f64 (Rewrite<= associate-*r/_binary64 (*.f64 (/.f64 (pow.f64 (cbrt.f64 (*.f64 b c)) 2) 1) (/.f64 (cbrt.f64 (*.f64 b c)) (hypot.f64 c d)))) (*.f64 (/.f64 a (hypot.f64 c d)) d)): 18 points increase in error, 13 points decrease in error
      (-.f64 (*.f64 (/.f64 (pow.f64 (cbrt.f64 (*.f64 b c)) 2) 1) (/.f64 (cbrt.f64 (*.f64 b c)) (hypot.f64 c d))) (Rewrite<= associate-/r/_binary64 (/.f64 a (/.f64 (hypot.f64 c d) d)))): 24 points increase in error, 19 points decrease in error
      (-.f64 (*.f64 (/.f64 (pow.f64 (cbrt.f64 (*.f64 b c)) 2) 1) (/.f64 (cbrt.f64 (*.f64 b c)) (hypot.f64 c d))) (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 a d) (hypot.f64 c d)))): 40 points increase in error, 8 points decrease in error
      (Rewrite=> fma-neg_binary64 (fma.f64 (/.f64 (pow.f64 (cbrt.f64 (*.f64 b c)) 2) 1) (/.f64 (cbrt.f64 (*.f64 b c)) (hypot.f64 c d)) (neg.f64 (/.f64 (*.f64 a d) (hypot.f64 c d))))): 1 points increase in error, 1 points decrease in error
3. Recombined 2 regimes into one program.

4. Final simplification 5.9

   \[\leadsto \begin{array}{l} \mathbf{if}\;c \leq -1.25 \cdot 10^{+113}:\\ \;\;\;\;\frac{b}{c} - \frac{a}{c} \cdot \frac{d}{c}\\ \mathbf{elif}\;c \leq 1.7 \cdot 10^{+115}:\\ \;\;\;\;\frac{1}{\mathsf{hypot}\left(c, d\right)} \cdot \left(\frac{{\left(\sqrt[3]{c \cdot b}\right)}^{3}}{\mathsf{hypot}\left(c, d\right)} - d \cdot \frac{a}{\mathsf{hypot}\left(c, d\right)}\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{b}{c} - \frac{a}{c} \cdot \frac{d}{c}\\ \end{array} \]

Alternatives

Alternative 1
Error	10.1
Cost	14160

\[\begin{array}{l} t_0 := \frac{\frac{c \cdot b - a \cdot d}{\mathsf{hypot}\left(c, d\right)}}{\mathsf{hypot}\left(c, d\right)}\\ t_1 := \frac{c \cdot \frac{b}{d} - a}{d}\\ \mathbf{if}\;d \leq -1.15 \cdot 10^{+78}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;d \leq -1.42 \cdot 10^{-67}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;d \leq 2.05 \cdot 10^{-114}:\\ \;\;\;\;\frac{\frac{a}{c} \cdot d - b}{-c}\\ \mathbf{elif}\;d \leq 1.1 \cdot 10^{+132}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;t_1\\ \end{array} \]

Alternative 2
Error	11.9
Cost	1488

\[\begin{array}{l} t_0 := \frac{c \cdot b - a \cdot d}{c \cdot c + d \cdot d}\\ t_1 := \frac{c \cdot \frac{b}{d} - a}{d}\\ \mathbf{if}\;d \leq -1.6 \cdot 10^{+78}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-66}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;d \leq 10^{-114}:\\ \;\;\;\;\frac{\frac{a}{c} \cdot d - b}{-c}\\ \mathbf{elif}\;d \leq 8 \cdot 10^{+131}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;t_1\\ \end{array} \]

Alternative 3
Error	20.8
Cost	1300

\[\begin{array}{l} t_0 := \frac{c \cdot \frac{b}{d} - a}{d}\\ t_1 := \frac{-a}{d + c \cdot \frac{c}{d}}\\ \mathbf{if}\;d \leq -2.1 \cdot 10^{-125}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;d \leq 8.2 \cdot 10^{-41}:\\ \;\;\;\;\frac{b}{c}\\ \mathbf{elif}\;d \leq 5.8 \cdot 10^{+63}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;d \leq 4.5 \cdot 10^{+124}:\\ \;\;\;\;\frac{b}{c}\\ \mathbf{elif}\;d \leq 3.6 \cdot 10^{+131}:\\ \;\;\;\;t_1\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

Alternative 4
Error	16.8
Cost	1168

\[\begin{array}{l} t_0 := \frac{\frac{a}{c} \cdot d - b}{-c}\\ t_1 := \frac{c \cdot \frac{b}{d} - a}{d}\\ \mathbf{if}\;d \leq -2.85 \cdot 10^{+78}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;d \leq -1.7 \cdot 10^{+43}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;d \leq -1.86 \cdot 10^{-28}:\\ \;\;\;\;\frac{-a}{d + c \cdot \frac{c}{d}}\\ \mathbf{elif}\;d \leq 4.5 \cdot 10^{+124}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;t_1\\ \end{array} \]

Alternative 5
Error	19.5
Cost	1104

\[\begin{array}{l} t_0 := \frac{c \cdot \frac{b}{d} - a}{d}\\ \mathbf{if}\;d \leq -1.12 \cdot 10^{+98}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;d \leq -3.3 \cdot 10^{+38}:\\ \;\;\;\;\frac{b}{c}\\ \mathbf{elif}\;d \leq -9.5 \cdot 10^{-25}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;d \leq 1.8 \cdot 10^{-41}:\\ \;\;\;\;\frac{b}{c}\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

Alternative 6
Error	24.2
Cost	520

\[\begin{array}{l} t_0 := \frac{-a}{d}\\ \mathbf{if}\;d \leq -1.45 \cdot 10^{+53}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;d \leq 4.5 \cdot 10^{+124}:\\ \;\;\;\;\frac{b}{c}\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

Alternative 7
Error	58.7
Cost	192

\[\frac{a}{c} \]

Alternative 8
Error	37.7
Cost	192

\[\frac{b}{c} \]

Reproduce

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