
(FPCore (a b c d) :precision binary64 (/ (+ (* a c) (* b d)) (+ (* c c) (* d d))))
double code(double a, double b, double c, double d) {
return ((a * c) + (b * d)) / ((c * c) + (d * d));
}
real(8) function code(a, b, c, d)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8), intent (in) :: d
code = ((a * c) + (b * d)) / ((c * c) + (d * d))
end function
public static double code(double a, double b, double c, double d) {
return ((a * c) + (b * d)) / ((c * c) + (d * d));
}
def code(a, b, c, d): return ((a * c) + (b * d)) / ((c * c) + (d * d))
function code(a, b, c, d) return Float64(Float64(Float64(a * c) + Float64(b * d)) / Float64(Float64(c * c) + Float64(d * d))) end
function tmp = code(a, b, c, d) tmp = ((a * c) + (b * d)) / ((c * c) + (d * d)); end
code[a_, b_, c_, d_] := N[(N[(N[(a * c), $MachinePrecision] + N[(b * d), $MachinePrecision]), $MachinePrecision] / N[(N[(c * c), $MachinePrecision] + N[(d * d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{a \cdot c + b \cdot d}{c \cdot c + d \cdot d}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 6 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b c d) :precision binary64 (/ (+ (* a c) (* b d)) (+ (* c c) (* d d))))
double code(double a, double b, double c, double d) {
return ((a * c) + (b * d)) / ((c * c) + (d * d));
}
real(8) function code(a, b, c, d)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8), intent (in) :: d
code = ((a * c) + (b * d)) / ((c * c) + (d * d))
end function
public static double code(double a, double b, double c, double d) {
return ((a * c) + (b * d)) / ((c * c) + (d * d));
}
def code(a, b, c, d): return ((a * c) + (b * d)) / ((c * c) + (d * d))
function code(a, b, c, d) return Float64(Float64(Float64(a * c) + Float64(b * d)) / Float64(Float64(c * c) + Float64(d * d))) end
function tmp = code(a, b, c, d) tmp = ((a * c) + (b * d)) / ((c * c) + (d * d)); end
code[a_, b_, c_, d_] := N[(N[(N[(a * c), $MachinePrecision] + N[(b * d), $MachinePrecision]), $MachinePrecision] / N[(N[(c * c), $MachinePrecision] + N[(d * d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{a \cdot c + b \cdot d}{c \cdot c + d \cdot d}
\end{array}
(FPCore (a b c d)
:precision binary64
(let* ((t_0 (/ (+ (* a c) (* b d)) (fma c c (* d d))))
(t_1 (/ (fma (/ b c) d a) c)))
(if (<= c -2.6e+110)
t_1
(if (<= c -3.5e-99)
t_0
(if (<= c 5.4e-68)
(/ (fma (/ c d) a b) d)
(if (<= c 2.6e+108) t_0 t_1))))))
double code(double a, double b, double c, double d) {
double t_0 = ((a * c) + (b * d)) / fma(c, c, (d * d));
double t_1 = fma((b / c), d, a) / c;
double tmp;
if (c <= -2.6e+110) {
tmp = t_1;
} else if (c <= -3.5e-99) {
tmp = t_0;
} else if (c <= 5.4e-68) {
tmp = fma((c / d), a, b) / d;
} else if (c <= 2.6e+108) {
tmp = t_0;
} else {
tmp = t_1;
}
return tmp;
}
function code(a, b, c, d) t_0 = Float64(Float64(Float64(a * c) + Float64(b * d)) / fma(c, c, Float64(d * d))) t_1 = Float64(fma(Float64(b / c), d, a) / c) tmp = 0.0 if (c <= -2.6e+110) tmp = t_1; elseif (c <= -3.5e-99) tmp = t_0; elseif (c <= 5.4e-68) tmp = Float64(fma(Float64(c / d), a, b) / d); elseif (c <= 2.6e+108) tmp = t_0; else tmp = t_1; end return tmp end
code[a_, b_, c_, d_] := Block[{t$95$0 = N[(N[(N[(a * c), $MachinePrecision] + N[(b * d), $MachinePrecision]), $MachinePrecision] / N[(c * c + N[(d * d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(N[(b / c), $MachinePrecision] * d + a), $MachinePrecision] / c), $MachinePrecision]}, If[LessEqual[c, -2.6e+110], t$95$1, If[LessEqual[c, -3.5e-99], t$95$0, If[LessEqual[c, 5.4e-68], N[(N[(N[(c / d), $MachinePrecision] * a + b), $MachinePrecision] / d), $MachinePrecision], If[LessEqual[c, 2.6e+108], t$95$0, t$95$1]]]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \frac{a \cdot c + b \cdot d}{\mathsf{fma}\left(c, c, d \cdot d\right)}\\
t_1 := \frac{\mathsf{fma}\left(\frac{b}{c}, d, a\right)}{c}\\
\mathbf{if}\;c \leq -2.6 \cdot 10^{+110}:\\
\;\;\;\;t\_1\\
\mathbf{elif}\;c \leq -3.5 \cdot 10^{-99}:\\
\;\;\;\;t\_0\\
\mathbf{elif}\;c \leq 5.4 \cdot 10^{-68}:\\
\;\;\;\;\frac{\mathsf{fma}\left(\frac{c}{d}, a, b\right)}{d}\\
\mathbf{elif}\;c \leq 2.6 \cdot 10^{+108}:\\
\;\;\;\;t\_0\\
\mathbf{else}:\\
\;\;\;\;t\_1\\
\end{array}
\end{array}
if c < -2.6e110 or 2.6000000000000002e108 < c Initial program 35.9%
Taylor expanded in c around inf
lower-/.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
*-commutativeN/A
lower-fma.f64N/A
lower-/.f6484.2
Applied rewrites84.2%
if -2.6e110 < c < -3.4999999999999999e-99 or 5.4000000000000003e-68 < c < 2.6000000000000002e108Initial program 85.6%
lift-+.f64N/A
lift-*.f64N/A
lower-fma.f6485.7
Applied rewrites85.7%
if -3.4999999999999999e-99 < c < 5.4000000000000003e-68Initial program 57.8%
lift-+.f64N/A
lift-*.f64N/A
lower-fma.f6457.8
Applied rewrites57.8%
Taylor expanded in c around 0
unpow2N/A
associate-/r*N/A
div-addN/A
lower-/.f64N/A
+-commutativeN/A
associate-/l*N/A
*-commutativeN/A
lower-fma.f64N/A
lower-/.f6489.2
Applied rewrites89.2%
(FPCore (a b c d) :precision binary64 (if (or (<= c -9e+78) (not (<= c 40000.0))) (/ (fma (/ b c) d a) c) (/ (fma (/ c d) a b) d)))
double code(double a, double b, double c, double d) {
double tmp;
if ((c <= -9e+78) || !(c <= 40000.0)) {
tmp = fma((b / c), d, a) / c;
} else {
tmp = fma((c / d), a, b) / d;
}
return tmp;
}
function code(a, b, c, d) tmp = 0.0 if ((c <= -9e+78) || !(c <= 40000.0)) tmp = Float64(fma(Float64(b / c), d, a) / c); else tmp = Float64(fma(Float64(c / d), a, b) / d); end return tmp end
code[a_, b_, c_, d_] := If[Or[LessEqual[c, -9e+78], N[Not[LessEqual[c, 40000.0]], $MachinePrecision]], N[(N[(N[(b / c), $MachinePrecision] * d + a), $MachinePrecision] / c), $MachinePrecision], N[(N[(N[(c / d), $MachinePrecision] * a + b), $MachinePrecision] / d), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;c \leq -9 \cdot 10^{+78} \lor \neg \left(c \leq 40000\right):\\
\;\;\;\;\frac{\mathsf{fma}\left(\frac{b}{c}, d, a\right)}{c}\\
\mathbf{else}:\\
\;\;\;\;\frac{\mathsf{fma}\left(\frac{c}{d}, a, b\right)}{d}\\
\end{array}
\end{array}
if c < -8.9999999999999999e78 or 4e4 < c Initial program 49.7%
Taylor expanded in c around inf
lower-/.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
*-commutativeN/A
lower-fma.f64N/A
lower-/.f6483.9
Applied rewrites83.9%
if -8.9999999999999999e78 < c < 4e4Initial program 65.9%
lift-+.f64N/A
lift-*.f64N/A
lower-fma.f6465.9
Applied rewrites65.9%
Taylor expanded in c around 0
unpow2N/A
associate-/r*N/A
div-addN/A
lower-/.f64N/A
+-commutativeN/A
associate-/l*N/A
*-commutativeN/A
lower-fma.f64N/A
lower-/.f6481.4
Applied rewrites81.4%
Final simplification82.4%
(FPCore (a b c d) :precision binary64 (if (or (<= c -9e+78) (not (<= c 1.5e-54))) (/ (fma (/ b c) d a) c) (/ (fma (/ a d) c b) d)))
double code(double a, double b, double c, double d) {
double tmp;
if ((c <= -9e+78) || !(c <= 1.5e-54)) {
tmp = fma((b / c), d, a) / c;
} else {
tmp = fma((a / d), c, b) / d;
}
return tmp;
}
function code(a, b, c, d) tmp = 0.0 if ((c <= -9e+78) || !(c <= 1.5e-54)) tmp = Float64(fma(Float64(b / c), d, a) / c); else tmp = Float64(fma(Float64(a / d), c, b) / d); end return tmp end
code[a_, b_, c_, d_] := If[Or[LessEqual[c, -9e+78], N[Not[LessEqual[c, 1.5e-54]], $MachinePrecision]], N[(N[(N[(b / c), $MachinePrecision] * d + a), $MachinePrecision] / c), $MachinePrecision], N[(N[(N[(a / d), $MachinePrecision] * c + b), $MachinePrecision] / d), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;c \leq -9 \cdot 10^{+78} \lor \neg \left(c \leq 1.5 \cdot 10^{-54}\right):\\
\;\;\;\;\frac{\mathsf{fma}\left(\frac{b}{c}, d, a\right)}{c}\\
\mathbf{else}:\\
\;\;\;\;\frac{\mathsf{fma}\left(\frac{a}{d}, c, b\right)}{d}\\
\end{array}
\end{array}
if c < -8.9999999999999999e78 or 1.50000000000000005e-54 < c Initial program 53.9%
Taylor expanded in c around inf
lower-/.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
*-commutativeN/A
lower-fma.f64N/A
lower-/.f6480.1
Applied rewrites80.1%
if -8.9999999999999999e78 < c < 1.50000000000000005e-54Initial program 63.6%
Taylor expanded in c around 0
unpow2N/A
associate-/r*N/A
div-addN/A
lower-/.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
*-commutativeN/A
lower-fma.f64N/A
lower-/.f6483.2
Applied rewrites83.2%
Final simplification81.7%
(FPCore (a b c d) :precision binary64 (if (or (<= c -8.4e+82) (not (<= c 1.5e-54))) (/ a c) (/ (fma (/ a d) c b) d)))
double code(double a, double b, double c, double d) {
double tmp;
if ((c <= -8.4e+82) || !(c <= 1.5e-54)) {
tmp = a / c;
} else {
tmp = fma((a / d), c, b) / d;
}
return tmp;
}
function code(a, b, c, d) tmp = 0.0 if ((c <= -8.4e+82) || !(c <= 1.5e-54)) tmp = Float64(a / c); else tmp = Float64(fma(Float64(a / d), c, b) / d); end return tmp end
code[a_, b_, c_, d_] := If[Or[LessEqual[c, -8.4e+82], N[Not[LessEqual[c, 1.5e-54]], $MachinePrecision]], N[(a / c), $MachinePrecision], N[(N[(N[(a / d), $MachinePrecision] * c + b), $MachinePrecision] / d), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;c \leq -8.4 \cdot 10^{+82} \lor \neg \left(c \leq 1.5 \cdot 10^{-54}\right):\\
\;\;\;\;\frac{a}{c}\\
\mathbf{else}:\\
\;\;\;\;\frac{\mathsf{fma}\left(\frac{a}{d}, c, b\right)}{d}\\
\end{array}
\end{array}
if c < -8.4000000000000001e82 or 1.50000000000000005e-54 < c Initial program 53.9%
Taylor expanded in c around inf
lower-/.f6468.0
Applied rewrites68.0%
if -8.4000000000000001e82 < c < 1.50000000000000005e-54Initial program 63.6%
Taylor expanded in c around 0
unpow2N/A
associate-/r*N/A
div-addN/A
lower-/.f64N/A
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
*-commutativeN/A
lower-fma.f64N/A
lower-/.f6483.2
Applied rewrites83.2%
Final simplification76.1%
(FPCore (a b c d) :precision binary64 (if (or (<= d -9.5e-50) (not (<= d 7.5e+23))) (/ b d) (/ a c)))
double code(double a, double b, double c, double d) {
double tmp;
if ((d <= -9.5e-50) || !(d <= 7.5e+23)) {
tmp = b / d;
} else {
tmp = a / c;
}
return tmp;
}
real(8) function code(a, b, c, d)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8), intent (in) :: d
real(8) :: tmp
if ((d <= (-9.5d-50)) .or. (.not. (d <= 7.5d+23))) then
tmp = b / d
else
tmp = a / c
end if
code = tmp
end function
public static double code(double a, double b, double c, double d) {
double tmp;
if ((d <= -9.5e-50) || !(d <= 7.5e+23)) {
tmp = b / d;
} else {
tmp = a / c;
}
return tmp;
}
def code(a, b, c, d): tmp = 0 if (d <= -9.5e-50) or not (d <= 7.5e+23): tmp = b / d else: tmp = a / c return tmp
function code(a, b, c, d) tmp = 0.0 if ((d <= -9.5e-50) || !(d <= 7.5e+23)) tmp = Float64(b / d); else tmp = Float64(a / c); end return tmp end
function tmp_2 = code(a, b, c, d) tmp = 0.0; if ((d <= -9.5e-50) || ~((d <= 7.5e+23))) tmp = b / d; else tmp = a / c; end tmp_2 = tmp; end
code[a_, b_, c_, d_] := If[Or[LessEqual[d, -9.5e-50], N[Not[LessEqual[d, 7.5e+23]], $MachinePrecision]], N[(b / d), $MachinePrecision], N[(a / c), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;d \leq -9.5 \cdot 10^{-50} \lor \neg \left(d \leq 7.5 \cdot 10^{+23}\right):\\
\;\;\;\;\frac{b}{d}\\
\mathbf{else}:\\
\;\;\;\;\frac{a}{c}\\
\end{array}
\end{array}
if d < -9.4999999999999993e-50 or 7.49999999999999987e23 < d Initial program 54.7%
Taylor expanded in c around 0
lower-/.f6465.2
Applied rewrites65.2%
if -9.4999999999999993e-50 < d < 7.49999999999999987e23Initial program 64.4%
Taylor expanded in c around inf
lower-/.f6470.7
Applied rewrites70.7%
Final simplification67.7%
(FPCore (a b c d) :precision binary64 (/ a c))
double code(double a, double b, double c, double d) {
return a / c;
}
real(8) function code(a, b, c, d)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8), intent (in) :: d
code = a / c
end function
public static double code(double a, double b, double c, double d) {
return a / c;
}
def code(a, b, c, d): return a / c
function code(a, b, c, d) return Float64(a / c) end
function tmp = code(a, b, c, d) tmp = a / c; end
code[a_, b_, c_, d_] := N[(a / c), $MachinePrecision]
\begin{array}{l}
\\
\frac{a}{c}
\end{array}
Initial program 59.1%
Taylor expanded in c around inf
lower-/.f6441.5
Applied rewrites41.5%
(FPCore (a b c d) :precision binary64 (if (< (fabs d) (fabs c)) (/ (+ a (* b (/ d c))) (+ c (* d (/ d c)))) (/ (+ b (* a (/ c d))) (+ d (* c (/ c d))))))
double code(double a, double b, double c, double d) {
double tmp;
if (fabs(d) < fabs(c)) {
tmp = (a + (b * (d / c))) / (c + (d * (d / c)));
} else {
tmp = (b + (a * (c / d))) / (d + (c * (c / d)));
}
return tmp;
}
real(8) function code(a, b, c, d)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8), intent (in) :: d
real(8) :: tmp
if (abs(d) < abs(c)) then
tmp = (a + (b * (d / c))) / (c + (d * (d / c)))
else
tmp = (b + (a * (c / d))) / (d + (c * (c / d)))
end if
code = tmp
end function
public static double code(double a, double b, double c, double d) {
double tmp;
if (Math.abs(d) < Math.abs(c)) {
tmp = (a + (b * (d / c))) / (c + (d * (d / c)));
} else {
tmp = (b + (a * (c / d))) / (d + (c * (c / d)));
}
return tmp;
}
def code(a, b, c, d): tmp = 0 if math.fabs(d) < math.fabs(c): tmp = (a + (b * (d / c))) / (c + (d * (d / c))) else: tmp = (b + (a * (c / d))) / (d + (c * (c / d))) return tmp
function code(a, b, c, d) tmp = 0.0 if (abs(d) < abs(c)) tmp = Float64(Float64(a + Float64(b * Float64(d / c))) / Float64(c + Float64(d * Float64(d / c)))); else tmp = Float64(Float64(b + Float64(a * Float64(c / d))) / Float64(d + Float64(c * Float64(c / d)))); end return tmp end
function tmp_2 = code(a, b, c, d) tmp = 0.0; if (abs(d) < abs(c)) tmp = (a + (b * (d / c))) / (c + (d * (d / c))); else tmp = (b + (a * (c / d))) / (d + (c * (c / d))); end tmp_2 = tmp; end
code[a_, b_, c_, d_] := If[Less[N[Abs[d], $MachinePrecision], N[Abs[c], $MachinePrecision]], N[(N[(a + N[(b * N[(d / c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(c + N[(d * N[(d / c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(b + N[(a * N[(c / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(d + N[(c * N[(c / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;\left|d\right| < \left|c\right|:\\
\;\;\;\;\frac{a + b \cdot \frac{d}{c}}{c + d \cdot \frac{d}{c}}\\
\mathbf{else}:\\
\;\;\;\;\frac{b + a \cdot \frac{c}{d}}{d + c \cdot \frac{c}{d}}\\
\end{array}
\end{array}
herbie shell --seed 2024326
(FPCore (a b c d)
:name "Complex division, real part"
:precision binary64
:alt
(! :herbie-platform default (if (< (fabs d) (fabs c)) (/ (+ a (* b (/ d c))) (+ c (* d (/ d c)))) (/ (+ b (* a (/ c d))) (+ d (* c (/ c d))))))
(/ (+ (* a c) (* b d)) (+ (* c c) (* d d))))