Henrywood and Agarwal, Equation (13)

Percentage Accurate: 24.9% → 24.9%
Time: 34.3s
Alternatives: 4
Speedup: 1.0×

Specification

?
\[\begin{array}{l} \\ \begin{array}{l} t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\ \frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right) \end{array} \end{array} \]
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (/ (* c0 (* d d)) (* (* w h) (* D D)))))
   (* (/ c0 (* 2.0 w)) (+ t_0 (sqrt (- (* t_0 t_0) (* M M)))))))
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
}
real(8) function code(c0, w, h, d, d_1, m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d
    real(8), intent (in) :: d_1
    real(8), intent (in) :: m
    real(8) :: t_0
    t_0 = (c0 * (d_1 * d_1)) / ((w * h) * (d * d))
    code = (c0 / (2.0d0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (m * m))))
end function
public static double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + Math.sqrt(((t_0 * t_0) - (M * M))));
}
def code(c0, w, h, D, d, M):
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D))
	return (c0 / (2.0 * w)) * (t_0 + math.sqrt(((t_0 * t_0) - (M * M))))
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D)))
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(t_0 + sqrt(Float64(Float64(t_0 * t_0) - Float64(M * M)))))
end
function tmp = code(c0, w, h, D, d, M)
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	tmp = (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
end
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$0 + N[Sqrt[N[(N[(t$95$0 * t$95$0), $MachinePrecision] - N[(M * M), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}

\\
\begin{array}{l}
t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\
\frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\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 4 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: 24.9% accurate, 1.0× speedup?

\[\begin{array}{l} \\ \begin{array}{l} t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\ \frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right) \end{array} \end{array} \]
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (/ (* c0 (* d d)) (* (* w h) (* D D)))))
   (* (/ c0 (* 2.0 w)) (+ t_0 (sqrt (- (* t_0 t_0) (* M M)))))))
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
}
real(8) function code(c0, w, h, d, d_1, m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d
    real(8), intent (in) :: d_1
    real(8), intent (in) :: m
    real(8) :: t_0
    t_0 = (c0 * (d_1 * d_1)) / ((w * h) * (d * d))
    code = (c0 / (2.0d0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (m * m))))
end function
public static double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + Math.sqrt(((t_0 * t_0) - (M * M))));
}
def code(c0, w, h, D, d, M):
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D))
	return (c0 / (2.0 * w)) * (t_0 + math.sqrt(((t_0 * t_0) - (M * M))))
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D)))
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(t_0 + sqrt(Float64(Float64(t_0 * t_0) - Float64(M * M)))))
end
function tmp = code(c0, w, h, D, d, M)
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	tmp = (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
end
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$0 + N[Sqrt[N[(N[(t$95$0 * t$95$0), $MachinePrecision] - N[(M * M), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}

\\
\begin{array}{l}
t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\
\frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right)
\end{array}
\end{array}

Alternative 1: 24.9% accurate, 1.0× speedup?

\[\begin{array}{l} D = |D|\\ d = |d|\\ M = |M|\\ \\ \begin{array}{l} t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\ \frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right) \end{array} \end{array} \]
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (/ (* c0 (* d d)) (* (* w h) (* D D)))))
   (* (/ c0 (* 2.0 w)) (+ t_0 (sqrt (- (* t_0 t_0) (* M M)))))))
D = abs(D);
d = abs(d);
M = abs(M);
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
}
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
real(8) function code(c0, w, h, d, d_1, m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d
    real(8), intent (in) :: d_1
    real(8), intent (in) :: m
    real(8) :: t_0
    t_0 = (c0 * (d_1 * d_1)) / ((w * h) * (d * d))
    code = (c0 / (2.0d0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (m * m))))
end function
D = Math.abs(D);
d = Math.abs(d);
M = Math.abs(M);
public static double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + Math.sqrt(((t_0 * t_0) - (M * M))));
}
D = abs(D)
d = abs(d)
M = abs(M)
def code(c0, w, h, D, d, M):
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D))
	return (c0 / (2.0 * w)) * (t_0 + math.sqrt(((t_0 * t_0) - (M * M))))
D = abs(D)
d = abs(d)
M = abs(M)
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D)))
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(t_0 + sqrt(Float64(Float64(t_0 * t_0) - Float64(M * M)))))
end
D = abs(D)
d = abs(d)
M = abs(M)
function tmp = code(c0, w, h, D, d, M)
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	tmp = (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
end
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$0 + N[Sqrt[N[(N[(t$95$0 * t$95$0), $MachinePrecision] - N[(M * M), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
D = |D|\\
d = |d|\\
M = |M|\\
\\
\begin{array}{l}
t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\
\frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right)
\end{array}
\end{array}
Derivation
  1. Initial program 24.0%

    \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right) \]

Alternative 2: 25.3% accurate, 0.6× speedup?

\[\begin{array}{l} D = |D|\\ d = |d|\\ M = |M|\\ \\ \begin{array}{l} t_0 := \frac{c0}{w \cdot h} \cdot \frac{d \cdot d}{D \cdot D}\\ \frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{\mathsf{fma}\left(t_0, t_0, -M \cdot M\right)}\right) \end{array} \end{array} \]
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (* (/ c0 (* w h)) (/ (* d d) (* D D)))))
   (* (/ c0 (* 2.0 w)) (+ t_0 (sqrt (fma t_0 t_0 (- (* M M))))))))
D = abs(D);
d = abs(d);
M = abs(M);
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 / (w * h)) * ((d * d) / (D * D));
	return (c0 / (2.0 * w)) * (t_0 + sqrt(fma(t_0, t_0, -(M * M))));
}
D = abs(D)
d = abs(d)
M = abs(M)
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(c0 / Float64(w * h)) * Float64(Float64(d * d) / Float64(D * D)))
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(t_0 + sqrt(fma(t_0, t_0, Float64(-Float64(M * M))))))
end
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(c0 / N[(w * h), $MachinePrecision]), $MachinePrecision] * N[(N[(d * d), $MachinePrecision] / N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$0 + N[Sqrt[N[(t$95$0 * t$95$0 + (-N[(M * M), $MachinePrecision])), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
D = |D|\\
d = |d|\\
M = |M|\\
\\
\begin{array}{l}
t_0 := \frac{c0}{w \cdot h} \cdot \frac{d \cdot d}{D \cdot D}\\
\frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{\mathsf{fma}\left(t_0, t_0, -M \cdot M\right)}\right)
\end{array}
\end{array}
Derivation
  1. Initial program 24.1%

    \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0}{w \cdot h} \cdot \frac{d \cdot d}{D \cdot D} + \sqrt{\mathsf{fma}\left(\frac{c0}{w \cdot h} \cdot \frac{d \cdot d}{D \cdot D}, \frac{c0}{w \cdot h} \cdot \frac{d \cdot d}{D \cdot D}, -M \cdot M\right)}\right) \]

Alternative 3: 31.3% accurate, 0.4× speedup?

\[\begin{array}{l} D = |D|\\ d = |d|\\ M = |M|\\ \\ \begin{array}{l} t_0 := \frac{d}{D} \cdot \frac{d}{D}\\ t_1 := \frac{\frac{c0}{h}}{w}\\ \frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(t_1, t_0, \sqrt{\mathsf{fma}\left(t_1, t_0, M\right) \cdot \left(\frac{c0 \cdot d}{D \cdot D} \cdot \frac{\frac{d}{h}}{w} - M\right)}\right) \end{array} \end{array} \]
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (* (/ d D) (/ d D))) (t_1 (/ (/ c0 h) w)))
   (*
    (/ c0 (* 2.0 w))
    (fma
     t_1
     t_0
     (sqrt
      (* (fma t_1 t_0 M) (- (* (/ (* c0 d) (* D D)) (/ (/ d h) w)) M)))))))
D = abs(D);
d = abs(d);
M = abs(M);
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (d / D) * (d / D);
	double t_1 = (c0 / h) / w;
	return (c0 / (2.0 * w)) * fma(t_1, t_0, sqrt((fma(t_1, t_0, M) * ((((c0 * d) / (D * D)) * ((d / h) / w)) - M))));
}
D = abs(D)
d = abs(d)
M = abs(M)
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(d / D) * Float64(d / D))
	t_1 = Float64(Float64(c0 / h) / w)
	return Float64(Float64(c0 / Float64(2.0 * w)) * fma(t_1, t_0, sqrt(Float64(fma(t_1, t_0, M) * Float64(Float64(Float64(Float64(c0 * d) / Float64(D * D)) * Float64(Float64(d / h) / w)) - M)))))
end
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(d / D), $MachinePrecision] * N[(d / D), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(c0 / h), $MachinePrecision] / w), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$1 * t$95$0 + N[Sqrt[N[(N[(t$95$1 * t$95$0 + M), $MachinePrecision] * N[(N[(N[(N[(c0 * d), $MachinePrecision] / N[(D * D), $MachinePrecision]), $MachinePrecision] * N[(N[(d / h), $MachinePrecision] / w), $MachinePrecision]), $MachinePrecision] - M), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
D = |D|\\
d = |d|\\
M = |M|\\
\\
\begin{array}{l}
t_0 := \frac{d}{D} \cdot \frac{d}{D}\\
t_1 := \frac{\frac{c0}{h}}{w}\\
\frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(t_1, t_0, \sqrt{\mathsf{fma}\left(t_1, t_0, M\right) \cdot \left(\frac{c0 \cdot d}{D \cdot D} \cdot \frac{\frac{d}{h}}{w} - M\right)}\right)
\end{array}
\end{array}
Derivation
  1. Initial program 29.9%

    \[\frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(\frac{\frac{c0}{h}}{w}, \frac{d}{D} \cdot \frac{d}{D}, \sqrt{\mathsf{fma}\left(\frac{\frac{c0}{h}}{w}, \frac{d}{D} \cdot \frac{d}{D}, M\right) \cdot \left(\frac{c0 \cdot d}{D \cdot D} \cdot \frac{\frac{d}{h}}{w} - M\right)}\right) \]

Alternative 4: 31.3% accurate, 0.3× speedup?

\[\begin{array}{l} D = |D|\\ d = |d|\\ M = |M|\\ \\ \begin{array}{l} t_0 := \frac{d}{D} \cdot \frac{d}{D}\\ t_1 := \frac{\frac{c0}{h}}{w}\\ \frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(t_1, t_0, \sqrt{\mathsf{fma}\left(t_0, t_1, M\right) \cdot \mathsf{fma}\left(\frac{c0 \cdot d}{D \cdot D}, \frac{\frac{d}{h}}{w}, -M\right)}\right) \end{array} \end{array} \]
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (* (/ d D) (/ d D))) (t_1 (/ (/ c0 h) w)))
   (*
    (/ c0 (* 2.0 w))
    (fma
     t_1
     t_0
     (sqrt
      (* (fma t_0 t_1 M) (fma (/ (* c0 d) (* D D)) (/ (/ d h) w) (- M))))))))
D = abs(D);
d = abs(d);
M = abs(M);
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (d / D) * (d / D);
	double t_1 = (c0 / h) / w;
	return (c0 / (2.0 * w)) * fma(t_1, t_0, sqrt((fma(t_0, t_1, M) * fma(((c0 * d) / (D * D)), ((d / h) / w), -M))));
}
D = abs(D)
d = abs(d)
M = abs(M)
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(d / D) * Float64(d / D))
	t_1 = Float64(Float64(c0 / h) / w)
	return Float64(Float64(c0 / Float64(2.0 * w)) * fma(t_1, t_0, sqrt(Float64(fma(t_0, t_1, M) * fma(Float64(Float64(c0 * d) / Float64(D * D)), Float64(Float64(d / h) / w), Float64(-M))))))
end
NOTE: D should be positive before calling this function
NOTE: d should be positive before calling this function
NOTE: M should be positive before calling this function
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(d / D), $MachinePrecision] * N[(d / D), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(c0 / h), $MachinePrecision] / w), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$1 * t$95$0 + N[Sqrt[N[(N[(t$95$0 * t$95$1 + M), $MachinePrecision] * N[(N[(N[(c0 * d), $MachinePrecision] / N[(D * D), $MachinePrecision]), $MachinePrecision] * N[(N[(d / h), $MachinePrecision] / w), $MachinePrecision] + (-M)), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
D = |D|\\
d = |d|\\
M = |M|\\
\\
\begin{array}{l}
t_0 := \frac{d}{D} \cdot \frac{d}{D}\\
t_1 := \frac{\frac{c0}{h}}{w}\\
\frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(t_1, t_0, \sqrt{\mathsf{fma}\left(t_0, t_1, M\right) \cdot \mathsf{fma}\left(\frac{c0 \cdot d}{D \cdot D}, \frac{\frac{d}{h}}{w}, -M\right)}\right)
\end{array}
\end{array}
Derivation
  1. Initial program 29.9%

    \[\frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(\frac{\frac{c0}{h}}{w}, \frac{d}{D} \cdot \frac{d}{D}, \sqrt{\mathsf{fma}\left(\frac{d}{D} \cdot \frac{d}{D}, \frac{\frac{c0}{h}}{w}, M\right) \cdot \mathsf{fma}\left(\frac{c0 \cdot d}{D \cdot D}, \frac{\frac{d}{h}}{w}, -M\right)}\right) \]

Reproduce

?
herbie shell --seed 2023276 
(FPCore (c0 w h D d M)
  :name "Henrywood and Agarwal, Equation (13)"
  :precision binary64
  (* (/ c0 (* 2.0 w)) (+ (/ (* c0 (* d d)) (* (* w h) (* D D))) (sqrt (- (* (/ (* c0 (* d d)) (* (* w h) (* D D))) (/ (* c0 (* d d)) (* (* w h) (* D D)))) (* M M))))))