
(FPCore (a1 a2 th) :precision binary64 (let* ((t_1 (/ (cos th) (sqrt 2.0)))) (+ (* t_1 (* a1 a1)) (* t_1 (* a2 a2)))))
double code(double a1, double a2, double th) {
double t_1 = cos(th) / sqrt(2.0);
return (t_1 * (a1 * a1)) + (t_1 * (a2 * a2));
}
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(a1, a2, th)
use fmin_fmax_functions
real(8), intent (in) :: a1
real(8), intent (in) :: a2
real(8), intent (in) :: th
real(8) :: t_1
t_1 = cos(th) / sqrt(2.0d0)
code = (t_1 * (a1 * a1)) + (t_1 * (a2 * a2))
end function
public static double code(double a1, double a2, double th) {
double t_1 = Math.cos(th) / Math.sqrt(2.0);
return (t_1 * (a1 * a1)) + (t_1 * (a2 * a2));
}
def code(a1, a2, th): t_1 = math.cos(th) / math.sqrt(2.0) return (t_1 * (a1 * a1)) + (t_1 * (a2 * a2))
function code(a1, a2, th) t_1 = Float64(cos(th) / sqrt(2.0)) return Float64(Float64(t_1 * Float64(a1 * a1)) + Float64(t_1 * Float64(a2 * a2))) end
function tmp = code(a1, a2, th) t_1 = cos(th) / sqrt(2.0); tmp = (t_1 * (a1 * a1)) + (t_1 * (a2 * a2)); end
code[a1_, a2_, th_] := Block[{t$95$1 = N[(N[Cos[th], $MachinePrecision] / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]}, N[(N[(t$95$1 * N[(a1 * a1), $MachinePrecision]), $MachinePrecision] + N[(t$95$1 * N[(a2 * a2), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
t_1 := \frac{\cos th}{\sqrt{2}}\\
t\_1 \cdot \left(a1 \cdot a1\right) + t\_1 \cdot \left(a2 \cdot a2\right)
\end{array}
Herbie found 6 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a1 a2 th) :precision binary64 (let* ((t_1 (/ (cos th) (sqrt 2.0)))) (+ (* t_1 (* a1 a1)) (* t_1 (* a2 a2)))))
double code(double a1, double a2, double th) {
double t_1 = cos(th) / sqrt(2.0);
return (t_1 * (a1 * a1)) + (t_1 * (a2 * a2));
}
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(a1, a2, th)
use fmin_fmax_functions
real(8), intent (in) :: a1
real(8), intent (in) :: a2
real(8), intent (in) :: th
real(8) :: t_1
t_1 = cos(th) / sqrt(2.0d0)
code = (t_1 * (a1 * a1)) + (t_1 * (a2 * a2))
end function
public static double code(double a1, double a2, double th) {
double t_1 = Math.cos(th) / Math.sqrt(2.0);
return (t_1 * (a1 * a1)) + (t_1 * (a2 * a2));
}
def code(a1, a2, th): t_1 = math.cos(th) / math.sqrt(2.0) return (t_1 * (a1 * a1)) + (t_1 * (a2 * a2))
function code(a1, a2, th) t_1 = Float64(cos(th) / sqrt(2.0)) return Float64(Float64(t_1 * Float64(a1 * a1)) + Float64(t_1 * Float64(a2 * a2))) end
function tmp = code(a1, a2, th) t_1 = cos(th) / sqrt(2.0); tmp = (t_1 * (a1 * a1)) + (t_1 * (a2 * a2)); end
code[a1_, a2_, th_] := Block[{t$95$1 = N[(N[Cos[th], $MachinePrecision] / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]}, N[(N[(t$95$1 * N[(a1 * a1), $MachinePrecision]), $MachinePrecision] + N[(t$95$1 * N[(a2 * a2), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
t_1 := \frac{\cos th}{\sqrt{2}}\\
t\_1 \cdot \left(a1 \cdot a1\right) + t\_1 \cdot \left(a2 \cdot a2\right)
\end{array}
(FPCore (a1 a2 th) :precision binary64 (* (cos th) (/ (fma a2 a2 (* a1 a1)) 1.4142135623730951)))
double code(double a1, double a2, double th) {
return cos(th) * (fma(a2, a2, (a1 * a1)) / 1.4142135623730951);
}
function code(a1, a2, th) return Float64(cos(th) * Float64(fma(a2, a2, Float64(a1 * a1)) / 1.4142135623730951)) end
code[a1_, a2_, th_] := N[(N[Cos[th], $MachinePrecision] * N[(N[(a2 * a2 + N[(a1 * a1), $MachinePrecision]), $MachinePrecision] / 1.4142135623730951), $MachinePrecision]), $MachinePrecision]
\cos th \cdot \frac{\mathsf{fma}\left(a2, a2, a1 \cdot a1\right)}{1.4142135623730951}
Initial program 99.6%
lift-+.f64N/A
lift-*.f64N/A
lift-*.f64N/A
distribute-lft-outN/A
lift-/.f64N/A
associate-*l/N/A
associate-/l*N/A
lower-*.f64N/A
lower-/.f64N/A
+-commutativeN/A
lift-*.f64N/A
lower-fma.f6499.6%
Applied rewrites99.6%
Evaluated real constant99.6%
(FPCore (a1 a2 th) :precision binary64 (* (* (cos th) (fma a2 a2 (* a1 a1))) 0.7071067811865475))
double code(double a1, double a2, double th) {
return (cos(th) * fma(a2, a2, (a1 * a1))) * 0.7071067811865475;
}
function code(a1, a2, th) return Float64(Float64(cos(th) * fma(a2, a2, Float64(a1 * a1))) * 0.7071067811865475) end
code[a1_, a2_, th_] := N[(N[(N[Cos[th], $MachinePrecision] * N[(a2 * a2 + N[(a1 * a1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * 0.7071067811865475), $MachinePrecision]
\left(\cos th \cdot \mathsf{fma}\left(a2, a2, a1 \cdot a1\right)\right) \cdot 0.7071067811865475
Initial program 99.6%
lift-+.f64N/A
lift-*.f64N/A
lift-*.f64N/A
distribute-lft-outN/A
lift-/.f64N/A
associate-*l/N/A
associate-/l*N/A
lower-*.f64N/A
lower-/.f64N/A
+-commutativeN/A
lift-*.f64N/A
lower-fma.f6499.6%
Applied rewrites99.6%
Evaluated real constant99.6%
lift-cos.f64N/A
lift-*.f64N/A
lift-/.f64N/A
associate-*r/N/A
mult-flipN/A
lower-*.f64N/A
lower-*.f64N/A
lift-cos.f64N/A
metadata-eval99.5%
Applied rewrites99.5%
(FPCore (a1 a2 th) :precision binary64 (* (* (cos th) 0.7071067811865475) (fma a2 a2 (* a1 a1))))
double code(double a1, double a2, double th) {
return (cos(th) * 0.7071067811865475) * fma(a2, a2, (a1 * a1));
}
function code(a1, a2, th) return Float64(Float64(cos(th) * 0.7071067811865475) * fma(a2, a2, Float64(a1 * a1))) end
code[a1_, a2_, th_] := N[(N[(N[Cos[th], $MachinePrecision] * 0.7071067811865475), $MachinePrecision] * N[(a2 * a2 + N[(a1 * a1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\left(\cos th \cdot 0.7071067811865475\right) \cdot \mathsf{fma}\left(a2, a2, a1 \cdot a1\right)
Initial program 99.6%
lift-+.f64N/A
lift-*.f64N/A
lift-*.f64N/A
distribute-lft-outN/A
lift-/.f64N/A
associate-*l/N/A
associate-/l*N/A
lower-*.f64N/A
lower-/.f64N/A
+-commutativeN/A
lift-*.f64N/A
lower-fma.f6499.6%
Applied rewrites99.6%
Evaluated real constant99.6%
lift-/.f64N/A
lift-fma.f64N/A
div-addN/A
associate-/l*N/A
lower-fma.f64N/A
lower-/.f64N/A
mult-flipN/A
lower-*.f64N/A
metadata-eval99.6%
Applied rewrites99.6%
lift-fma.f64N/A
+-commutativeN/A
lift-*.f64N/A
lift-/.f64N/A
associate-*r/N/A
mult-flip-revN/A
metadata-evalN/A
distribute-rgt-inN/A
+-commutativeN/A
lift-fma.f64N/A
lift-*.f6499.5%
lower-*.f64N/A
lift-*.f64N/A
associate-*r*N/A
lower-*.f64N/A
lower-*.f6499.5%
Applied rewrites99.5%
(FPCore (a1 a2 th)
:precision binary64
(let* ((t_1 (fma a2 a2 (* a1 a1))) (t_2 (/ (cos th) (sqrt 2.0))))
(if (<= (+ (* t_2 (* a1 a1)) (* t_2 (* a2 a2))) -5e-160)
(* (* 0.7071067811865475 t_1) (+ 1.0 (* -0.5 (pow th 2.0))))
(/ t_1 (sqrt 2.0)))))double code(double a1, double a2, double th) {
double t_1 = fma(a2, a2, (a1 * a1));
double t_2 = cos(th) / sqrt(2.0);
double tmp;
if (((t_2 * (a1 * a1)) + (t_2 * (a2 * a2))) <= -5e-160) {
tmp = (0.7071067811865475 * t_1) * (1.0 + (-0.5 * pow(th, 2.0)));
} else {
tmp = t_1 / sqrt(2.0);
}
return tmp;
}
function code(a1, a2, th) t_1 = fma(a2, a2, Float64(a1 * a1)) t_2 = Float64(cos(th) / sqrt(2.0)) tmp = 0.0 if (Float64(Float64(t_2 * Float64(a1 * a1)) + Float64(t_2 * Float64(a2 * a2))) <= -5e-160) tmp = Float64(Float64(0.7071067811865475 * t_1) * Float64(1.0 + Float64(-0.5 * (th ^ 2.0)))); else tmp = Float64(t_1 / sqrt(2.0)); end return tmp end
code[a1_, a2_, th_] := Block[{t$95$1 = N[(a2 * a2 + N[(a1 * a1), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(N[Cos[th], $MachinePrecision] / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(N[(t$95$2 * N[(a1 * a1), $MachinePrecision]), $MachinePrecision] + N[(t$95$2 * N[(a2 * a2), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], -5e-160], N[(N[(0.7071067811865475 * t$95$1), $MachinePrecision] * N[(1.0 + N[(-0.5 * N[Power[th, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(t$95$1 / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]]]]
\begin{array}{l}
t_1 := \mathsf{fma}\left(a2, a2, a1 \cdot a1\right)\\
t_2 := \frac{\cos th}{\sqrt{2}}\\
\mathbf{if}\;t\_2 \cdot \left(a1 \cdot a1\right) + t\_2 \cdot \left(a2 \cdot a2\right) \leq -5 \cdot 10^{-160}:\\
\;\;\;\;\left(0.7071067811865475 \cdot t\_1\right) \cdot \left(1 + -0.5 \cdot {th}^{2}\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{t\_1}{\sqrt{2}}\\
\end{array}
if (+.f64 (*.f64 (/.f64 (cos.f64 th) (sqrt.f64 #s(literal 2 binary64))) (*.f64 a1 a1)) (*.f64 (/.f64 (cos.f64 th) (sqrt.f64 #s(literal 2 binary64))) (*.f64 a2 a2))) < -4.9999999999999999e-160Initial program 99.6%
lift-+.f64N/A
lift-*.f64N/A
lift-*.f64N/A
distribute-lft-outN/A
lift-/.f64N/A
associate-*l/N/A
associate-/l*N/A
lower-*.f64N/A
lower-/.f64N/A
+-commutativeN/A
lift-*.f64N/A
lower-fma.f6499.6%
Applied rewrites99.6%
Evaluated real constant99.6%
lift-cos.f64N/A
lift-*.f64N/A
*-commutativeN/A
lower-*.f64N/A
lift-/.f64N/A
mult-flipN/A
*-commutativeN/A
lower-*.f64N/A
metadata-evalN/A
lift-cos.f6499.5%
Applied rewrites99.5%
Taylor expanded in th around 0
lower-+.f64N/A
lower-*.f64N/A
lower-pow.f6463.4%
Applied rewrites63.4%
if -4.9999999999999999e-160 < (+.f64 (*.f64 (/.f64 (cos.f64 th) (sqrt.f64 #s(literal 2 binary64))) (*.f64 a1 a1)) (*.f64 (/.f64 (cos.f64 th) (sqrt.f64 #s(literal 2 binary64))) (*.f64 a2 a2))) Initial program 99.6%
Taylor expanded in th around 0
lower-+.f64N/A
lower-/.f64N/A
lower-pow.f64N/A
lower-sqrt.f64N/A
lower-/.f64N/A
lower-pow.f64N/A
lower-sqrt.f6466.2%
Applied rewrites66.2%
lift-+.f64N/A
+-commutativeN/A
lift-/.f64N/A
lift-pow.f64N/A
pow2N/A
lift-*.f64N/A
lift-/.f64N/A
lift-pow.f64N/A
pow2N/A
lift-*.f64N/A
div-addN/A
lift-*.f64N/A
lift-fma.f64N/A
lift-/.f6466.2%
Applied rewrites66.2%
(FPCore (a1 a2 th) :precision binary64 (/ (fma a2 a2 (* a1 a1)) (sqrt 2.0)))
double code(double a1, double a2, double th) {
return fma(a2, a2, (a1 * a1)) / sqrt(2.0);
}
function code(a1, a2, th) return Float64(fma(a2, a2, Float64(a1 * a1)) / sqrt(2.0)) end
code[a1_, a2_, th_] := N[(N[(a2 * a2 + N[(a1 * a1), $MachinePrecision]), $MachinePrecision] / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]
\frac{\mathsf{fma}\left(a2, a2, a1 \cdot a1\right)}{\sqrt{2}}
Initial program 99.6%
Taylor expanded in th around 0
lower-+.f64N/A
lower-/.f64N/A
lower-pow.f64N/A
lower-sqrt.f64N/A
lower-/.f64N/A
lower-pow.f64N/A
lower-sqrt.f6466.2%
Applied rewrites66.2%
lift-+.f64N/A
+-commutativeN/A
lift-/.f64N/A
lift-pow.f64N/A
pow2N/A
lift-*.f64N/A
lift-/.f64N/A
lift-pow.f64N/A
pow2N/A
lift-*.f64N/A
div-addN/A
lift-*.f64N/A
lift-fma.f64N/A
lift-/.f6466.2%
Applied rewrites66.2%
(FPCore (a1 a2 th) :precision binary64 (/ (* 1.0 (* a1 a1)) (sqrt 2.0)))
double code(double a1, double a2, double th) {
return (1.0 * (a1 * a1)) / sqrt(2.0);
}
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(a1, a2, th)
use fmin_fmax_functions
real(8), intent (in) :: a1
real(8), intent (in) :: a2
real(8), intent (in) :: th
code = (1.0d0 * (a1 * a1)) / sqrt(2.0d0)
end function
public static double code(double a1, double a2, double th) {
return (1.0 * (a1 * a1)) / Math.sqrt(2.0);
}
def code(a1, a2, th): return (1.0 * (a1 * a1)) / math.sqrt(2.0)
function code(a1, a2, th) return Float64(Float64(1.0 * Float64(a1 * a1)) / sqrt(2.0)) end
function tmp = code(a1, a2, th) tmp = (1.0 * (a1 * a1)) / sqrt(2.0); end
code[a1_, a2_, th_] := N[(N[(1.0 * N[(a1 * a1), $MachinePrecision]), $MachinePrecision] / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]
\frac{1 \cdot \left(a1 \cdot a1\right)}{\sqrt{2}}
Initial program 99.6%
Taylor expanded in th around 0
lower-+.f64N/A
lower-/.f64N/A
lower-pow.f64N/A
lower-sqrt.f64N/A
lower-/.f64N/A
lower-pow.f64N/A
lower-sqrt.f6466.2%
Applied rewrites66.2%
lift-+.f64N/A
+-commutativeN/A
lift-/.f64N/A
lift-pow.f64N/A
pow2N/A
lift-*.f64N/A
lift-/.f64N/A
lift-pow.f64N/A
pow2N/A
lift-*.f64N/A
div-addN/A
lift-*.f64N/A
lift-fma.f64N/A
lift-/.f6466.2%
Applied rewrites66.2%
lift-fma.f64N/A
+-commutativeN/A
sum-to-multN/A
lower-unsound-*.f64N/A
lower-unsound-+.f64N/A
lower-unsound-/.f64N/A
lower-*.f6443.3%
Applied rewrites43.3%
Taylor expanded in a1 around inf
Applied rewrites39.4%
herbie shell --seed 2025193
(FPCore (a1 a2 th)
:name "Migdal et al, Equation (64)"
:precision binary64
(+ (* (/ (cos th) (sqrt 2.0)) (* a1 a1)) (* (/ (cos th) (sqrt 2.0)) (* a2 a2))))