
(FPCore (a b angle x-scale y-scale)
:precision binary64
(let* ((t_0 (* (/ angle 180.0) (PI)))
(t_1 (sin t_0))
(t_2 (cos t_0))
(t_3
(/
(/ (* (* (* 2.0 (- (pow b 2.0) (pow a 2.0))) t_1) t_2) x-scale)
y-scale)))
(-
(* t_3 t_3)
(*
(*
4.0
(/ (/ (+ (pow (* a t_1) 2.0) (pow (* b t_2) 2.0)) x-scale) x-scale))
(/ (/ (+ (pow (* a t_2) 2.0) (pow (* b t_1) 2.0)) y-scale) y-scale)))))\begin{array}{l}
\\
\begin{array}{l}
t_0 := \frac{angle}{180} \cdot \mathsf{PI}\left(\right)\\
t_1 := \sin t\_0\\
t_2 := \cos t\_0\\
t_3 := \frac{\frac{\left(\left(2 \cdot \left({b}^{2} - {a}^{2}\right)\right) \cdot t\_1\right) \cdot t\_2}{x-scale}}{y-scale}\\
t\_3 \cdot t\_3 - \left(4 \cdot \frac{\frac{{\left(a \cdot t\_1\right)}^{2} + {\left(b \cdot t\_2\right)}^{2}}{x-scale}}{x-scale}\right) \cdot \frac{\frac{{\left(a \cdot t\_2\right)}^{2} + {\left(b \cdot t\_1\right)}^{2}}{y-scale}}{y-scale}
\end{array}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 6 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b angle x-scale y-scale)
:precision binary64
(let* ((t_0 (* (/ angle 180.0) (PI)))
(t_1 (sin t_0))
(t_2 (cos t_0))
(t_3
(/
(/ (* (* (* 2.0 (- (pow b 2.0) (pow a 2.0))) t_1) t_2) x-scale)
y-scale)))
(-
(* t_3 t_3)
(*
(*
4.0
(/ (/ (+ (pow (* a t_1) 2.0) (pow (* b t_2) 2.0)) x-scale) x-scale))
(/ (/ (+ (pow (* a t_2) 2.0) (pow (* b t_1) 2.0)) y-scale) y-scale)))))\begin{array}{l}
\\
\begin{array}{l}
t_0 := \frac{angle}{180} \cdot \mathsf{PI}\left(\right)\\
t_1 := \sin t\_0\\
t_2 := \cos t\_0\\
t_3 := \frac{\frac{\left(\left(2 \cdot \left({b}^{2} - {a}^{2}\right)\right) \cdot t\_1\right) \cdot t\_2}{x-scale}}{y-scale}\\
t\_3 \cdot t\_3 - \left(4 \cdot \frac{\frac{{\left(a \cdot t\_1\right)}^{2} + {\left(b \cdot t\_2\right)}^{2}}{x-scale}}{x-scale}\right) \cdot \frac{\frac{{\left(a \cdot t\_2\right)}^{2} + {\left(b \cdot t\_1\right)}^{2}}{y-scale}}{y-scale}
\end{array}
\end{array}
(FPCore (a b angle x-scale y-scale) :precision binary64 (let* ((t_0 (/ (* a b) (* y-scale x-scale)))) (* -4.0 (* t_0 t_0))))
double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double t_0 = (a * b) / (y_45_scale * x_45_scale);
return -4.0 * (t_0 * t_0);
}
real(8) function code(a, b, angle, x_45scale, y_45scale)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: angle
real(8), intent (in) :: x_45scale
real(8), intent (in) :: y_45scale
real(8) :: t_0
t_0 = (a * b) / (y_45scale * x_45scale)
code = (-4.0d0) * (t_0 * t_0)
end function
public static double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double t_0 = (a * b) / (y_45_scale * x_45_scale);
return -4.0 * (t_0 * t_0);
}
def code(a, b, angle, x_45_scale, y_45_scale): t_0 = (a * b) / (y_45_scale * x_45_scale) return -4.0 * (t_0 * t_0)
function code(a, b, angle, x_45_scale, y_45_scale) t_0 = Float64(Float64(a * b) / Float64(y_45_scale * x_45_scale)) return Float64(-4.0 * Float64(t_0 * t_0)) end
function tmp = code(a, b, angle, x_45_scale, y_45_scale) t_0 = (a * b) / (y_45_scale * x_45_scale); tmp = -4.0 * (t_0 * t_0); end
code[a_, b_, angle_, x$45$scale_, y$45$scale_] := Block[{t$95$0 = N[(N[(a * b), $MachinePrecision] / N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision]}, N[(-4.0 * N[(t$95$0 * t$95$0), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \frac{a \cdot b}{y-scale \cdot x-scale}\\
-4 \cdot \left(t\_0 \cdot t\_0\right)
\end{array}
\end{array}
Initial program 28.1%
Taylor expanded in angle around 0
associate-/l*N/A
associate-*r*N/A
lower-*.f64N/A
lower-*.f64N/A
unpow2N/A
lower-*.f64N/A
unpow2N/A
*-commutativeN/A
times-fracN/A
lower-*.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f6457.6
Applied rewrites57.6%
Applied rewrites83.4%
Applied rewrites96.3%
Final simplification96.3%
(FPCore (a b angle x-scale y-scale)
:precision binary64
(let* ((t_0
(*
(* b b)
(* (/ a (* y-scale x-scale)) (/ (* -4.0 a) (* y-scale x-scale)))))
(t_1 (/ b (* y-scale x-scale))))
(if (<= a 1.9e-163)
t_0
(if (<= a 7e+153) (* (* t_1 t_1) (* (* a a) -4.0)) t_0))))
double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double t_0 = (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale)));
double t_1 = b / (y_45_scale * x_45_scale);
double tmp;
if (a <= 1.9e-163) {
tmp = t_0;
} else if (a <= 7e+153) {
tmp = (t_1 * t_1) * ((a * a) * -4.0);
} else {
tmp = t_0;
}
return tmp;
}
real(8) function code(a, b, angle, x_45scale, y_45scale)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: angle
real(8), intent (in) :: x_45scale
real(8), intent (in) :: y_45scale
real(8) :: t_0
real(8) :: t_1
real(8) :: tmp
t_0 = (b * b) * ((a / (y_45scale * x_45scale)) * (((-4.0d0) * a) / (y_45scale * x_45scale)))
t_1 = b / (y_45scale * x_45scale)
if (a <= 1.9d-163) then
tmp = t_0
else if (a <= 7d+153) then
tmp = (t_1 * t_1) * ((a * a) * (-4.0d0))
else
tmp = t_0
end if
code = tmp
end function
public static double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double t_0 = (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale)));
double t_1 = b / (y_45_scale * x_45_scale);
double tmp;
if (a <= 1.9e-163) {
tmp = t_0;
} else if (a <= 7e+153) {
tmp = (t_1 * t_1) * ((a * a) * -4.0);
} else {
tmp = t_0;
}
return tmp;
}
def code(a, b, angle, x_45_scale, y_45_scale): t_0 = (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale))) t_1 = b / (y_45_scale * x_45_scale) tmp = 0 if a <= 1.9e-163: tmp = t_0 elif a <= 7e+153: tmp = (t_1 * t_1) * ((a * a) * -4.0) else: tmp = t_0 return tmp
function code(a, b, angle, x_45_scale, y_45_scale) t_0 = Float64(Float64(b * b) * Float64(Float64(a / Float64(y_45_scale * x_45_scale)) * Float64(Float64(-4.0 * a) / Float64(y_45_scale * x_45_scale)))) t_1 = Float64(b / Float64(y_45_scale * x_45_scale)) tmp = 0.0 if (a <= 1.9e-163) tmp = t_0; elseif (a <= 7e+153) tmp = Float64(Float64(t_1 * t_1) * Float64(Float64(a * a) * -4.0)); else tmp = t_0; end return tmp end
function tmp_2 = code(a, b, angle, x_45_scale, y_45_scale) t_0 = (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale))); t_1 = b / (y_45_scale * x_45_scale); tmp = 0.0; if (a <= 1.9e-163) tmp = t_0; elseif (a <= 7e+153) tmp = (t_1 * t_1) * ((a * a) * -4.0); else tmp = t_0; end tmp_2 = tmp; end
code[a_, b_, angle_, x$45$scale_, y$45$scale_] := Block[{t$95$0 = N[(N[(b * b), $MachinePrecision] * N[(N[(a / N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision] * N[(N[(-4.0 * a), $MachinePrecision] / N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(b / N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[a, 1.9e-163], t$95$0, If[LessEqual[a, 7e+153], N[(N[(t$95$1 * t$95$1), $MachinePrecision] * N[(N[(a * a), $MachinePrecision] * -4.0), $MachinePrecision]), $MachinePrecision], t$95$0]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \left(b \cdot b\right) \cdot \left(\frac{a}{y-scale \cdot x-scale} \cdot \frac{-4 \cdot a}{y-scale \cdot x-scale}\right)\\
t_1 := \frac{b}{y-scale \cdot x-scale}\\
\mathbf{if}\;a \leq 1.9 \cdot 10^{-163}:\\
\;\;\;\;t\_0\\
\mathbf{elif}\;a \leq 7 \cdot 10^{+153}:\\
\;\;\;\;\left(t\_1 \cdot t\_1\right) \cdot \left(\left(a \cdot a\right) \cdot -4\right)\\
\mathbf{else}:\\
\;\;\;\;t\_0\\
\end{array}
\end{array}
if a < 1.9e-163 or 6.9999999999999998e153 < a Initial program 25.8%
Taylor expanded in b around 0
Applied rewrites48.6%
Taylor expanded in angle around 0
Applied rewrites56.8%
Applied rewrites74.9%
if 1.9e-163 < a < 6.9999999999999998e153Initial program 36.1%
Taylor expanded in angle around 0
associate-/l*N/A
associate-*r*N/A
lower-*.f64N/A
lower-*.f64N/A
unpow2N/A
lower-*.f64N/A
unpow2N/A
*-commutativeN/A
times-fracN/A
lower-*.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f6474.2
Applied rewrites74.2%
Applied rewrites95.6%
Final simplification79.4%
(FPCore (a b angle x-scale y-scale)
:precision binary64
(if (<= b 4.8e-125)
0.0
(*
(/ (* b b) (* (* y-scale x-scale) (* y-scale x-scale)))
(* (* a a) -4.0))))
double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double tmp;
if (b <= 4.8e-125) {
tmp = 0.0;
} else {
tmp = ((b * b) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * ((a * a) * -4.0);
}
return tmp;
}
real(8) function code(a, b, angle, x_45scale, y_45scale)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: angle
real(8), intent (in) :: x_45scale
real(8), intent (in) :: y_45scale
real(8) :: tmp
if (b <= 4.8d-125) then
tmp = 0.0d0
else
tmp = ((b * b) / ((y_45scale * x_45scale) * (y_45scale * x_45scale))) * ((a * a) * (-4.0d0))
end if
code = tmp
end function
public static double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double tmp;
if (b <= 4.8e-125) {
tmp = 0.0;
} else {
tmp = ((b * b) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * ((a * a) * -4.0);
}
return tmp;
}
def code(a, b, angle, x_45_scale, y_45_scale): tmp = 0 if b <= 4.8e-125: tmp = 0.0 else: tmp = ((b * b) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * ((a * a) * -4.0) return tmp
function code(a, b, angle, x_45_scale, y_45_scale) tmp = 0.0 if (b <= 4.8e-125) tmp = 0.0; else tmp = Float64(Float64(Float64(b * b) / Float64(Float64(y_45_scale * x_45_scale) * Float64(y_45_scale * x_45_scale))) * Float64(Float64(a * a) * -4.0)); end return tmp end
function tmp_2 = code(a, b, angle, x_45_scale, y_45_scale) tmp = 0.0; if (b <= 4.8e-125) tmp = 0.0; else tmp = ((b * b) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * ((a * a) * -4.0); end tmp_2 = tmp; end
code[a_, b_, angle_, x$45$scale_, y$45$scale_] := If[LessEqual[b, 4.8e-125], 0.0, N[(N[(N[(b * b), $MachinePrecision] / N[(N[(y$45$scale * x$45$scale), $MachinePrecision] * N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(N[(a * a), $MachinePrecision] * -4.0), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \leq 4.8 \cdot 10^{-125}:\\
\;\;\;\;0\\
\mathbf{else}:\\
\;\;\;\;\frac{b \cdot b}{\left(y-scale \cdot x-scale\right) \cdot \left(y-scale \cdot x-scale\right)} \cdot \left(\left(a \cdot a\right) \cdot -4\right)\\
\end{array}
\end{array}
if b < 4.8000000000000003e-125Initial program 29.2%
Applied rewrites26.1%
Taylor expanded in a around 0
cancel-sign-sub-invN/A
metadata-evalN/A
+-inverses43.3
Applied rewrites43.3%
if 4.8000000000000003e-125 < b Initial program 26.2%
Taylor expanded in angle around 0
associate-/l*N/A
associate-*r*N/A
lower-*.f64N/A
lower-*.f64N/A
unpow2N/A
lower-*.f64N/A
unpow2N/A
*-commutativeN/A
times-fracN/A
lower-*.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f6460.0
Applied rewrites60.0%
Taylor expanded in b around 0
Applied rewrites72.4%
Final simplification54.4%
(FPCore (a b angle x-scale y-scale)
:precision binary64
(if (<= b 4.8e-125)
0.0
(*
(/ (* (* a a) -4.0) (* (* y-scale x-scale) (* y-scale x-scale)))
(* b b))))
double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double tmp;
if (b <= 4.8e-125) {
tmp = 0.0;
} else {
tmp = (((a * a) * -4.0) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * (b * b);
}
return tmp;
}
real(8) function code(a, b, angle, x_45scale, y_45scale)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: angle
real(8), intent (in) :: x_45scale
real(8), intent (in) :: y_45scale
real(8) :: tmp
if (b <= 4.8d-125) then
tmp = 0.0d0
else
tmp = (((a * a) * (-4.0d0)) / ((y_45scale * x_45scale) * (y_45scale * x_45scale))) * (b * b)
end if
code = tmp
end function
public static double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
double tmp;
if (b <= 4.8e-125) {
tmp = 0.0;
} else {
tmp = (((a * a) * -4.0) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * (b * b);
}
return tmp;
}
def code(a, b, angle, x_45_scale, y_45_scale): tmp = 0 if b <= 4.8e-125: tmp = 0.0 else: tmp = (((a * a) * -4.0) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * (b * b) return tmp
function code(a, b, angle, x_45_scale, y_45_scale) tmp = 0.0 if (b <= 4.8e-125) tmp = 0.0; else tmp = Float64(Float64(Float64(Float64(a * a) * -4.0) / Float64(Float64(y_45_scale * x_45_scale) * Float64(y_45_scale * x_45_scale))) * Float64(b * b)); end return tmp end
function tmp_2 = code(a, b, angle, x_45_scale, y_45_scale) tmp = 0.0; if (b <= 4.8e-125) tmp = 0.0; else tmp = (((a * a) * -4.0) / ((y_45_scale * x_45_scale) * (y_45_scale * x_45_scale))) * (b * b); end tmp_2 = tmp; end
code[a_, b_, angle_, x$45$scale_, y$45$scale_] := If[LessEqual[b, 4.8e-125], 0.0, N[(N[(N[(N[(a * a), $MachinePrecision] * -4.0), $MachinePrecision] / N[(N[(y$45$scale * x$45$scale), $MachinePrecision] * N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(b * b), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \leq 4.8 \cdot 10^{-125}:\\
\;\;\;\;0\\
\mathbf{else}:\\
\;\;\;\;\frac{\left(a \cdot a\right) \cdot -4}{\left(y-scale \cdot x-scale\right) \cdot \left(y-scale \cdot x-scale\right)} \cdot \left(b \cdot b\right)\\
\end{array}
\end{array}
if b < 4.8000000000000003e-125Initial program 29.2%
Applied rewrites26.1%
Taylor expanded in a around 0
cancel-sign-sub-invN/A
metadata-evalN/A
+-inverses43.3
Applied rewrites43.3%
if 4.8000000000000003e-125 < b Initial program 26.2%
Taylor expanded in b around 0
Applied rewrites60.9%
Taylor expanded in angle around 0
Applied rewrites71.6%
(FPCore (a b angle x-scale y-scale) :precision binary64 (* (* b b) (* (/ a (* y-scale x-scale)) (/ (* -4.0 a) (* y-scale x-scale)))))
double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
return (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale)));
}
real(8) function code(a, b, angle, x_45scale, y_45scale)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: angle
real(8), intent (in) :: x_45scale
real(8), intent (in) :: y_45scale
code = (b * b) * ((a / (y_45scale * x_45scale)) * (((-4.0d0) * a) / (y_45scale * x_45scale)))
end function
public static double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
return (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale)));
}
def code(a, b, angle, x_45_scale, y_45_scale): return (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale)))
function code(a, b, angle, x_45_scale, y_45_scale) return Float64(Float64(b * b) * Float64(Float64(a / Float64(y_45_scale * x_45_scale)) * Float64(Float64(-4.0 * a) / Float64(y_45_scale * x_45_scale)))) end
function tmp = code(a, b, angle, x_45_scale, y_45_scale) tmp = (b * b) * ((a / (y_45_scale * x_45_scale)) * ((-4.0 * a) / (y_45_scale * x_45_scale))); end
code[a_, b_, angle_, x$45$scale_, y$45$scale_] := N[(N[(b * b), $MachinePrecision] * N[(N[(a / N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision] * N[(N[(-4.0 * a), $MachinePrecision] / N[(y$45$scale * x$45$scale), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\left(b \cdot b\right) \cdot \left(\frac{a}{y-scale \cdot x-scale} \cdot \frac{-4 \cdot a}{y-scale \cdot x-scale}\right)
\end{array}
Initial program 28.1%
Taylor expanded in b around 0
Applied rewrites54.1%
Taylor expanded in angle around 0
Applied rewrites62.1%
Applied rewrites76.6%
Final simplification76.6%
(FPCore (a b angle x-scale y-scale) :precision binary64 0.0)
double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
return 0.0;
}
real(8) function code(a, b, angle, x_45scale, y_45scale)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: angle
real(8), intent (in) :: x_45scale
real(8), intent (in) :: y_45scale
code = 0.0d0
end function
public static double code(double a, double b, double angle, double x_45_scale, double y_45_scale) {
return 0.0;
}
def code(a, b, angle, x_45_scale, y_45_scale): return 0.0
function code(a, b, angle, x_45_scale, y_45_scale) return 0.0 end
function tmp = code(a, b, angle, x_45_scale, y_45_scale) tmp = 0.0; end
code[a_, b_, angle_, x$45$scale_, y$45$scale_] := 0.0
\begin{array}{l}
\\
0
\end{array}
Initial program 28.1%
Applied rewrites25.0%
Taylor expanded in a around 0
cancel-sign-sub-invN/A
metadata-evalN/A
+-inverses36.8
Applied rewrites36.8%
herbie shell --seed 2024312
(FPCore (a b angle x-scale y-scale)
:name "Simplification of discriminant from scale-rotated-ellipse"
:precision binary64
(- (* (/ (/ (* (* (* 2.0 (- (pow b 2.0) (pow a 2.0))) (sin (* (/ angle 180.0) (PI)))) (cos (* (/ angle 180.0) (PI)))) x-scale) y-scale) (/ (/ (* (* (* 2.0 (- (pow b 2.0) (pow a 2.0))) (sin (* (/ angle 180.0) (PI)))) (cos (* (/ angle 180.0) (PI)))) x-scale) y-scale)) (* (* 4.0 (/ (/ (+ (pow (* a (sin (* (/ angle 180.0) (PI)))) 2.0) (pow (* b (cos (* (/ angle 180.0) (PI)))) 2.0)) x-scale) x-scale)) (/ (/ (+ (pow (* a (cos (* (/ angle 180.0) (PI)))) 2.0) (pow (* b (sin (* (/ angle 180.0) (PI)))) 2.0)) y-scale) y-scale))))