Migdal et al, Equation (64)

Percentage Accurate: 99.5% → 99.6%
Time: 13.4s
Alternatives: 16
Speedup: 1.9×

Specification

?
\[\begin{array}{l} \\ \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} \end{array} \]
(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));
}
real(8) function code(a1, a2, th)
    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}

\\
\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}
\end{array}

Sampling outcomes in binary64 precision:

Local Percentage Accuracy vs ?

The average percentage accuracy by input value. Horizontal axis shows value of an input variable; the variable is choosen in the title. Vertical axis is accuracy; higher is better. Red represent the original program, while blue represents Herbie's suggestion. These can be toggled with buttons below the plot. The line is an average while dots represent individual samples.

Accuracy vs Speed?

Herbie found 16 alternatives:

AlternativeAccuracySpeedup
The accuracy (vertical axis) and speed (horizontal axis) of each alternatives. Up and to the right is better. The red square shows the initial program, and each blue circle shows an alternative.The line shows the best available speed-accuracy tradeoffs.

Initial Program: 99.5% accurate, 1.0× speedup?

\[\begin{array}{l} \\ \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} \end{array} \]
(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));
}
real(8) function code(a1, a2, th)
    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}

\\
\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}
\end{array}

Alternative 1: 99.6% accurate, 1.0× speedup?

\[\begin{array}{l} a2_m = \left|a2\right| \\ [a1, a2_m, th] = \mathsf{sort}([a1, a2_m, th])\\ \\ \mathsf{fma}\left(\cos th \cdot a2\_m, \frac{a2\_m}{\sqrt{2}}, \cos th \cdot \left(a1 \cdot \frac{a1}{\sqrt{2}}\right)\right) \end{array} \]
a2_m = (fabs.f64 a2)
NOTE: a1, a2_m, and th should be sorted in increasing order before calling this function.
(FPCore (a1 a2_m th)
 :precision binary64
 (fma
  (* (cos th) a2_m)
  (/ a2_m (sqrt 2.0))
  (* (cos th) (* a1 (/ a1 (sqrt 2.0))))))
a2_m = fabs(a2);
assert(a1 < a2_m && a2_m < th);
double code(double a1, double a2_m, double th) {
	return fma((cos(th) * a2_m), (a2_m / sqrt(2.0)), (cos(th) * (a1 * (a1 / sqrt(2.0)))));
}
a2_m = abs(a2)
a1, a2_m, th = sort([a1, a2_m, th])
function code(a1, a2_m, th)
	return fma(Float64(cos(th) * a2_m), Float64(a2_m / sqrt(2.0)), Float64(cos(th) * Float64(a1 * Float64(a1 / sqrt(2.0)))))
end
a2_m = N[Abs[a2], $MachinePrecision]
NOTE: a1, a2_m, and th should be sorted in increasing order before calling this function.
code[a1_, a2$95$m_, th_] := N[(N[(N[Cos[th], $MachinePrecision] * a2$95$m), $MachinePrecision] * N[(a2$95$m / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision] + N[(N[Cos[th], $MachinePrecision] * N[(a1 * N[(a1 / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
a2_m = \left|a2\right|
\\
[a1, a2_m, th] = \mathsf{sort}([a1, a2_m, th])\\
\\
\mathsf{fma}\left(\cos th \cdot a2\_m, \frac{a2\_m}{\sqrt{2}}, \cos th \cdot \left(a1 \cdot \frac{a1}{\sqrt{2}}\right)\right)
\end{array}
Derivation
  1. Initial program 99.2%

    \[\frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) + \frac{\cos th}{\sqrt{2}} \cdot \left(a2 \cdot a2\right) \]
  2. Add Preprocessing
  3. Step-by-step derivation
    1. lift-+.f64N/A

      \[\leadsto \color{blue}{\frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) + \frac{\cos th}{\sqrt{2}} \cdot \left(a2 \cdot a2\right)} \]
    2. +-commutativeN/A

      \[\leadsto \color{blue}{\frac{\cos th}{\sqrt{2}} \cdot \left(a2 \cdot a2\right) + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right)} \]
    3. lift-*.f64N/A

      \[\leadsto \color{blue}{\frac{\cos th}{\sqrt{2}} \cdot \left(a2 \cdot a2\right)} + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    4. lift-/.f64N/A

      \[\leadsto \color{blue}{\frac{\cos th}{\sqrt{2}}} \cdot \left(a2 \cdot a2\right) + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    5. associate-*l/N/A

      \[\leadsto \color{blue}{\frac{\cos th \cdot \left(a2 \cdot a2\right)}{\sqrt{2}}} + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    6. lift-*.f64N/A

      \[\leadsto \frac{\cos th \cdot \color{blue}{\left(a2 \cdot a2\right)}}{\sqrt{2}} + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    7. associate-*r*N/A

      \[\leadsto \frac{\color{blue}{\left(\cos th \cdot a2\right) \cdot a2}}{\sqrt{2}} + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    8. associate-/l*N/A

      \[\leadsto \color{blue}{\left(\cos th \cdot a2\right) \cdot \frac{a2}{\sqrt{2}}} + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    9. *-lft-identityN/A

      \[\leadsto \left(\cos th \cdot a2\right) \cdot \frac{\color{blue}{1 \cdot a2}}{\sqrt{2}} + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    10. associate-*l/N/A

      \[\leadsto \left(\cos th \cdot a2\right) \cdot \color{blue}{\left(\frac{1}{\sqrt{2}} \cdot a2\right)} + \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) \]
    11. lower-fma.f64N/A

      \[\leadsto \color{blue}{\mathsf{fma}\left(\cos th \cdot a2, \frac{1}{\sqrt{2}} \cdot a2, \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right)\right)} \]
    12. lower-*.f64N/A

      \[\leadsto \mathsf{fma}\left(\color{blue}{\cos th \cdot a2}, \frac{1}{\sqrt{2}} \cdot a2, \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right)\right) \]
    13. associate-*l/N/A

      \[\leadsto \mathsf{fma}\left(\cos th \cdot a2, \color{blue}{\frac{1 \cdot a2}{\sqrt{2}}}, \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right)\right) \]
    14. *-lft-identityN/A

      \[\leadsto \mathsf{fma}\left(\cos th \cdot a2, \frac{\color{blue}{a2}}{\sqrt{2}}, \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right)\right) \]
    15. lower-/.f6499.6

      \[\leadsto \mathsf{fma}\left(\cos th \cdot a2, \color{blue}{\frac{a2}{\sqrt{2}}}, \frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right)\right) \]
    16. lift-*.f64N/A

      \[\leadsto \mathsf{fma}\left(\cos th \cdot a2, \frac{a2}{\sqrt{2}}, \color{blue}{\frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right)}\right) \]
    17. lift-/.f64N/A

      \[\leadsto \mathsf{fma}\left(\cos th \cdot a2, \frac{a2}{\sqrt{2}}, \color{blue}{\frac{\cos th}{\sqrt{2}}} \cdot \left(a1 \cdot a1\right)\right) \]
  4. Applied rewrites99.6%

    \[\leadsto \color{blue}{\mathsf{fma}\left(\cos th \cdot a2, \frac{a2}{\sqrt{2}}, \cos th \cdot \left(a1 \cdot \frac{a1}{\sqrt{2}}\right)\right)} \]
  5. Add Preprocessing

Alternative 2: 76.0% accurate, 0.8× speedup?

\[\begin{array}{l} a2_m = \left|a2\right| \\ [a1, a2_m, th] = \mathsf{sort}([a1, a2_m, th])\\ \\ \begin{array}{l} t_1 := \frac{\cos th}{\sqrt{2}}\\ \mathbf{if}\;t\_1 \cdot \left(a1 \cdot a1\right) + \left(a2\_m \cdot a2\_m\right) \cdot t\_1 \leq -5 \cdot 10^{-130}:\\ \;\;\;\;\frac{a2\_m \cdot \mathsf{fma}\left(th, a2\_m \cdot \left(th \cdot -0.5\right), a2\_m\right)}{\sqrt{2}}\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{1}{\frac{\sqrt{2}}{a2\_m \cdot a2\_m}}\right)\\ \end{array} \end{array} \]
a2_m = (fabs.f64 a2)
NOTE: a1, a2_m, and th should be sorted in increasing order before calling this function.
(FPCore (a1 a2_m th)
 :precision binary64
 (let* ((t_1 (/ (cos th) (sqrt 2.0))))
   (if (<= (+ (* t_1 (* a1 a1)) (* (* a2_m a2_m) t_1)) -5e-130)
     (/ (* a2_m (fma th (* a2_m (* th -0.5)) a2_m)) (sqrt 2.0))
     (fma a1 (/ a1 (sqrt 2.0)) (/ 1.0 (/ (sqrt 2.0) (* a2_m a2_m)))))))
a2_m = fabs(a2);
assert(a1 < a2_m && a2_m < th);
double code(double a1, double a2_m, double th) {
	double t_1 = cos(th) / sqrt(2.0);
	double tmp;
	if (((t_1 * (a1 * a1)) + ((a2_m * a2_m) * t_1)) <= -5e-130) {
		tmp = (a2_m * fma(th, (a2_m * (th * -0.5)), a2_m)) / sqrt(2.0);
	} else {
		tmp = fma(a1, (a1 / sqrt(2.0)), (1.0 / (sqrt(2.0) / (a2_m * a2_m))));
	}
	return tmp;
}
a2_m = abs(a2)
a1, a2_m, th = sort([a1, a2_m, th])
function code(a1, a2_m, th)
	t_1 = Float64(cos(th) / sqrt(2.0))
	tmp = 0.0
	if (Float64(Float64(t_1 * Float64(a1 * a1)) + Float64(Float64(a2_m * a2_m) * t_1)) <= -5e-130)
		tmp = Float64(Float64(a2_m * fma(th, Float64(a2_m * Float64(th * -0.5)), a2_m)) / sqrt(2.0));
	else
		tmp = fma(a1, Float64(a1 / sqrt(2.0)), Float64(1.0 / Float64(sqrt(2.0) / Float64(a2_m * a2_m))));
	end
	return tmp
end
a2_m = N[Abs[a2], $MachinePrecision]
NOTE: a1, a2_m, and th should be sorted in increasing order before calling this function.
code[a1_, a2$95$m_, th_] := Block[{t$95$1 = N[(N[Cos[th], $MachinePrecision] / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(N[(t$95$1 * N[(a1 * a1), $MachinePrecision]), $MachinePrecision] + N[(N[(a2$95$m * a2$95$m), $MachinePrecision] * t$95$1), $MachinePrecision]), $MachinePrecision], -5e-130], N[(N[(a2$95$m * N[(th * N[(a2$95$m * N[(th * -0.5), $MachinePrecision]), $MachinePrecision] + a2$95$m), $MachinePrecision]), $MachinePrecision] / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision], N[(a1 * N[(a1 / N[Sqrt[2.0], $MachinePrecision]), $MachinePrecision] + N[(1.0 / N[(N[Sqrt[2.0], $MachinePrecision] / N[(a2$95$m * a2$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
a2_m = \left|a2\right|
\\
[a1, a2_m, th] = \mathsf{sort}([a1, a2_m, th])\\
\\
\begin{array}{l}
t_1 := \frac{\cos th}{\sqrt{2}}\\
\mathbf{if}\;t\_1 \cdot \left(a1 \cdot a1\right) + \left(a2\_m \cdot a2\_m\right) \cdot t\_1 \leq -5 \cdot 10^{-130}:\\
\;\;\;\;\frac{a2\_m \cdot \mathsf{fma}\left(th, a2\_m \cdot \left(th \cdot -0.5\right), a2\_m\right)}{\sqrt{2}}\\

\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{1}{\frac{\sqrt{2}}{a2\_m \cdot a2\_m}}\right)\\


\end{array}
\end{array}
Derivation
  1. Split input into 2 regimes
  2. 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.9999999999999996e-130

    1. Initial program 99.6%

      \[\frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) + \frac{\cos th}{\sqrt{2}} \cdot \left(a2 \cdot a2\right) \]
    2. Add Preprocessing
    3. Taylor expanded in th around 0

      \[\leadsto \color{blue}{{th}^{2} \cdot \left(\frac{-1}{2} \cdot \frac{{a1}^{2}}{\sqrt{2}} + \frac{-1}{2} \cdot \frac{{a2}^{2}}{\sqrt{2}}\right) + \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right)} \]
    4. Step-by-step derivation
      1. distribute-lft-outN/A

        \[\leadsto {th}^{2} \cdot \color{blue}{\left(\frac{-1}{2} \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right)\right)} + \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      2. associate-*r*N/A

        \[\leadsto \color{blue}{\left({th}^{2} \cdot \frac{-1}{2}\right) \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right)} + \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      3. *-commutativeN/A

        \[\leadsto \color{blue}{\left(\frac{-1}{2} \cdot {th}^{2}\right)} \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) + \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      4. distribute-lft1-inN/A

        \[\leadsto \color{blue}{\left(\frac{-1}{2} \cdot {th}^{2} + 1\right) \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right)} \]
      5. +-commutativeN/A

        \[\leadsto \color{blue}{\left(1 + \frac{-1}{2} \cdot {th}^{2}\right)} \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      6. lower-*.f64N/A

        \[\leadsto \color{blue}{\left(1 + \frac{-1}{2} \cdot {th}^{2}\right) \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right)} \]
      7. +-commutativeN/A

        \[\leadsto \color{blue}{\left(\frac{-1}{2} \cdot {th}^{2} + 1\right)} \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      8. *-commutativeN/A

        \[\leadsto \left(\color{blue}{{th}^{2} \cdot \frac{-1}{2}} + 1\right) \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      9. unpow2N/A

        \[\leadsto \left(\color{blue}{\left(th \cdot th\right)} \cdot \frac{-1}{2} + 1\right) \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      10. associate-*l*N/A

        \[\leadsto \left(\color{blue}{th \cdot \left(th \cdot \frac{-1}{2}\right)} + 1\right) \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      11. lower-fma.f64N/A

        \[\leadsto \color{blue}{\mathsf{fma}\left(th, th \cdot \frac{-1}{2}, 1\right)} \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
      12. lower-*.f64N/A

        \[\leadsto \mathsf{fma}\left(th, \color{blue}{th \cdot \frac{-1}{2}}, 1\right) \cdot \left(\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}\right) \]
    5. Applied rewrites56.7%

      \[\leadsto \color{blue}{\mathsf{fma}\left(th, th \cdot -0.5, 1\right) \cdot \mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{a2 \cdot a2}{\sqrt{2}}\right)} \]
    6. Taylor expanded in a1 around 0

      \[\leadsto \frac{{a2}^{2} \cdot \left(1 + \frac{-1}{2} \cdot {th}^{2}\right)}{\color{blue}{\sqrt{2}}} \]
    7. Step-by-step derivation
      1. Applied rewrites52.4%

        \[\leadsto \frac{a2 \cdot \mathsf{fma}\left(-0.5 \cdot \left(th \cdot th\right), a2, a2\right)}{\color{blue}{\sqrt{2}}} \]
      2. Step-by-step derivation
        1. Applied rewrites50.3%

          \[\leadsto \frac{\mathsf{fma}\left(th, \left(th \cdot -0.5\right) \cdot a2, a2\right) \cdot a2}{\sqrt{2}} \]

        if -4.9999999999999996e-130 < (+.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)))

        1. Initial program 99.5%

          \[\frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) + \frac{\cos th}{\sqrt{2}} \cdot \left(a2 \cdot a2\right) \]
        2. Add Preprocessing
        3. Taylor expanded in th around 0

          \[\leadsto \color{blue}{\frac{{a1}^{2}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}}} \]
        4. Step-by-step derivation
          1. unpow2N/A

            \[\leadsto \frac{\color{blue}{a1 \cdot a1}}{\sqrt{2}} + \frac{{a2}^{2}}{\sqrt{2}} \]
          2. associate-/l*N/A

            \[\leadsto \color{blue}{a1 \cdot \frac{a1}{\sqrt{2}}} + \frac{{a2}^{2}}{\sqrt{2}} \]
          3. lower-fma.f64N/A

            \[\leadsto \color{blue}{\mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{{a2}^{2}}{\sqrt{2}}\right)} \]
          4. lower-/.f64N/A

            \[\leadsto \mathsf{fma}\left(a1, \color{blue}{\frac{a1}{\sqrt{2}}}, \frac{{a2}^{2}}{\sqrt{2}}\right) \]
          5. lower-sqrt.f64N/A

            \[\leadsto \mathsf{fma}\left(a1, \frac{a1}{\color{blue}{\sqrt{2}}}, \frac{{a2}^{2}}{\sqrt{2}}\right) \]
          6. lower-/.f64N/A

            \[\leadsto \mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \color{blue}{\frac{{a2}^{2}}{\sqrt{2}}}\right) \]
          7. unpow2N/A

            \[\leadsto \mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{\color{blue}{a2 \cdot a2}}{\sqrt{2}}\right) \]
          8. lower-*.f64N/A

            \[\leadsto \mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{\color{blue}{a2 \cdot a2}}{\sqrt{2}}\right) \]
          9. lower-sqrt.f6483.3

            \[\leadsto \mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{a2 \cdot a2}{\color{blue}{\sqrt{2}}}\right) \]
        5. Applied rewrites83.3%

          \[\leadsto \color{blue}{\mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{a2 \cdot a2}{\sqrt{2}}\right)} \]
        6. Step-by-step derivation
          1. Applied rewrites83.2%

            \[\leadsto \mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{1}{\frac{\sqrt{2}}{a2 \cdot a2}}\right) \]
        7. Recombined 2 regimes into one program.
        8. Final simplification76.0%

          \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{\cos th}{\sqrt{2}} \cdot \left(a1 \cdot a1\right) + \left(a2 \cdot a2\right) \cdot \frac{\cos th}{\sqrt{2}} \leq -5 \cdot 10^{-130}:\\ \;\;\;\;\frac{a2 \cdot \mathsf{fma}\left(th, a2 \cdot \left(th \cdot -0.5\right), a2\right)}{\sqrt{2}}\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(a1, \frac{a1}{\sqrt{2}}, \frac{1}{\frac{\sqrt{2}}{a2 \cdot a2}}\right)\\ \end{array} \]
        9. Add Preprocessing

        Reproduce

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