
(FPCore (B x) :precision binary64 (+ (- (* x (/ 1.0 (tan B)))) (/ 1.0 (sin B))))
double code(double B, double x) {
return -(x * (1.0 / tan(B))) + (1.0 / sin(B));
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = -(x * (1.0d0 / tan(b))) + (1.0d0 / sin(b))
end function
public static double code(double B, double x) {
return -(x * (1.0 / Math.tan(B))) + (1.0 / Math.sin(B));
}
def code(B, x): return -(x * (1.0 / math.tan(B))) + (1.0 / math.sin(B))
function code(B, x) return Float64(Float64(-Float64(x * Float64(1.0 / tan(B)))) + Float64(1.0 / sin(B))) end
function tmp = code(B, x) tmp = -(x * (1.0 / tan(B))) + (1.0 / sin(B)); end
code[B_, x_] := N[((-N[(x * N[(1.0 / N[Tan[B], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]) + N[(1.0 / N[Sin[B], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\left(-x \cdot \frac{1}{\tan B}\right) + \frac{1}{\sin B}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 10 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (B x) :precision binary64 (+ (- (* x (/ 1.0 (tan B)))) (/ 1.0 (sin B))))
double code(double B, double x) {
return -(x * (1.0 / tan(B))) + (1.0 / sin(B));
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = -(x * (1.0d0 / tan(b))) + (1.0d0 / sin(b))
end function
public static double code(double B, double x) {
return -(x * (1.0 / Math.tan(B))) + (1.0 / Math.sin(B));
}
def code(B, x): return -(x * (1.0 / math.tan(B))) + (1.0 / math.sin(B))
function code(B, x) return Float64(Float64(-Float64(x * Float64(1.0 / tan(B)))) + Float64(1.0 / sin(B))) end
function tmp = code(B, x) tmp = -(x * (1.0 / tan(B))) + (1.0 / sin(B)); end
code[B_, x_] := N[((-N[(x * N[(1.0 / N[Tan[B], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]) + N[(1.0 / N[Sin[B], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\left(-x \cdot \frac{1}{\tan B}\right) + \frac{1}{\sin B}
\end{array}
(FPCore (B x) :precision binary64 (- (/ 1.0 (sin B)) (/ x (tan B))))
double code(double B, double x) {
return (1.0 / sin(B)) - (x / tan(B));
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = (1.0d0 / sin(b)) - (x / tan(b))
end function
public static double code(double B, double x) {
return (1.0 / Math.sin(B)) - (x / Math.tan(B));
}
def code(B, x): return (1.0 / math.sin(B)) - (x / math.tan(B))
function code(B, x) return Float64(Float64(1.0 / sin(B)) - Float64(x / tan(B))) end
function tmp = code(B, x) tmp = (1.0 / sin(B)) - (x / tan(B)); end
code[B_, x_] := N[(N[(1.0 / N[Sin[B], $MachinePrecision]), $MachinePrecision] - N[(x / N[Tan[B], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{1}{\sin B} - \frac{x}{\tan B}
\end{array}
Initial program 99.7%
distribute-lft-neg-in99.7%
+-commutative99.7%
*-commutative99.7%
remove-double-neg99.7%
distribute-frac-neg299.7%
tan-neg99.7%
cancel-sign-sub-inv99.7%
*-commutative99.7%
associate-*r/99.8%
*-rgt-identity99.8%
tan-neg99.8%
distribute-neg-frac299.8%
distribute-neg-frac99.8%
remove-double-neg99.8%
Simplified99.8%
Final simplification99.8%
(FPCore (B x) :precision binary64 (if (or (<= x -3.0) (not (<= x 3500000000.0))) (* x (/ (cos B) (- (sin B)))) (+ (/ 1.0 (sin B)) (* x (/ -1.0 B)))))
double code(double B, double x) {
double tmp;
if ((x <= -3.0) || !(x <= 3500000000.0)) {
tmp = x * (cos(B) / -sin(B));
} else {
tmp = (1.0 / sin(B)) + (x * (-1.0 / B));
}
return tmp;
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
real(8) :: tmp
if ((x <= (-3.0d0)) .or. (.not. (x <= 3500000000.0d0))) then
tmp = x * (cos(b) / -sin(b))
else
tmp = (1.0d0 / sin(b)) + (x * ((-1.0d0) / b))
end if
code = tmp
end function
public static double code(double B, double x) {
double tmp;
if ((x <= -3.0) || !(x <= 3500000000.0)) {
tmp = x * (Math.cos(B) / -Math.sin(B));
} else {
tmp = (1.0 / Math.sin(B)) + (x * (-1.0 / B));
}
return tmp;
}
def code(B, x): tmp = 0 if (x <= -3.0) or not (x <= 3500000000.0): tmp = x * (math.cos(B) / -math.sin(B)) else: tmp = (1.0 / math.sin(B)) + (x * (-1.0 / B)) return tmp
function code(B, x) tmp = 0.0 if ((x <= -3.0) || !(x <= 3500000000.0)) tmp = Float64(x * Float64(cos(B) / Float64(-sin(B)))); else tmp = Float64(Float64(1.0 / sin(B)) + Float64(x * Float64(-1.0 / B))); end return tmp end
function tmp_2 = code(B, x) tmp = 0.0; if ((x <= -3.0) || ~((x <= 3500000000.0))) tmp = x * (cos(B) / -sin(B)); else tmp = (1.0 / sin(B)) + (x * (-1.0 / B)); end tmp_2 = tmp; end
code[B_, x_] := If[Or[LessEqual[x, -3.0], N[Not[LessEqual[x, 3500000000.0]], $MachinePrecision]], N[(x * N[(N[Cos[B], $MachinePrecision] / (-N[Sin[B], $MachinePrecision])), $MachinePrecision]), $MachinePrecision], N[(N[(1.0 / N[Sin[B], $MachinePrecision]), $MachinePrecision] + N[(x * N[(-1.0 / B), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;x \leq -3 \lor \neg \left(x \leq 3500000000\right):\\
\;\;\;\;x \cdot \frac{\cos B}{-\sin B}\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{\sin B} + x \cdot \frac{-1}{B}\\
\end{array}
\end{array}
if x < -3 or 3.5e9 < x Initial program 99.5%
Taylor expanded in x around inf 97.5%
mul-1-neg97.5%
associate-/l*97.5%
distribute-rgt-neg-in97.5%
distribute-neg-frac97.5%
Simplified97.5%
if -3 < x < 3.5e9Initial program 99.9%
Taylor expanded in B around 0 99.4%
Final simplification98.5%
(FPCore (B x) :precision binary64 (if (or (<= x -3.0) (not (<= x 3500000000.0))) (* (cos B) (/ x (- (sin B)))) (+ (/ 1.0 (sin B)) (* x (/ -1.0 B)))))
double code(double B, double x) {
double tmp;
if ((x <= -3.0) || !(x <= 3500000000.0)) {
tmp = cos(B) * (x / -sin(B));
} else {
tmp = (1.0 / sin(B)) + (x * (-1.0 / B));
}
return tmp;
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
real(8) :: tmp
if ((x <= (-3.0d0)) .or. (.not. (x <= 3500000000.0d0))) then
tmp = cos(b) * (x / -sin(b))
else
tmp = (1.0d0 / sin(b)) + (x * ((-1.0d0) / b))
end if
code = tmp
end function
public static double code(double B, double x) {
double tmp;
if ((x <= -3.0) || !(x <= 3500000000.0)) {
tmp = Math.cos(B) * (x / -Math.sin(B));
} else {
tmp = (1.0 / Math.sin(B)) + (x * (-1.0 / B));
}
return tmp;
}
def code(B, x): tmp = 0 if (x <= -3.0) or not (x <= 3500000000.0): tmp = math.cos(B) * (x / -math.sin(B)) else: tmp = (1.0 / math.sin(B)) + (x * (-1.0 / B)) return tmp
function code(B, x) tmp = 0.0 if ((x <= -3.0) || !(x <= 3500000000.0)) tmp = Float64(cos(B) * Float64(x / Float64(-sin(B)))); else tmp = Float64(Float64(1.0 / sin(B)) + Float64(x * Float64(-1.0 / B))); end return tmp end
function tmp_2 = code(B, x) tmp = 0.0; if ((x <= -3.0) || ~((x <= 3500000000.0))) tmp = cos(B) * (x / -sin(B)); else tmp = (1.0 / sin(B)) + (x * (-1.0 / B)); end tmp_2 = tmp; end
code[B_, x_] := If[Or[LessEqual[x, -3.0], N[Not[LessEqual[x, 3500000000.0]], $MachinePrecision]], N[(N[Cos[B], $MachinePrecision] * N[(x / (-N[Sin[B], $MachinePrecision])), $MachinePrecision]), $MachinePrecision], N[(N[(1.0 / N[Sin[B], $MachinePrecision]), $MachinePrecision] + N[(x * N[(-1.0 / B), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;x \leq -3 \lor \neg \left(x \leq 3500000000\right):\\
\;\;\;\;\cos B \cdot \frac{x}{-\sin B}\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{\sin B} + x \cdot \frac{-1}{B}\\
\end{array}
\end{array}
if x < -3 or 3.5e9 < x Initial program 99.5%
distribute-lft-neg-in99.5%
+-commutative99.5%
*-commutative99.5%
remove-double-neg99.5%
distribute-frac-neg299.5%
tan-neg99.5%
cancel-sign-sub-inv99.5%
*-commutative99.5%
associate-*r/99.7%
*-rgt-identity99.7%
tan-neg99.7%
distribute-neg-frac299.7%
distribute-neg-frac99.7%
remove-double-neg99.7%
Simplified99.7%
tan-quot99.6%
associate-/r/99.7%
Applied egg-rr99.7%
Taylor expanded in x around inf 97.5%
mul-1-neg97.5%
associate-*l/97.6%
distribute-lft-neg-in97.6%
*-commutative97.6%
distribute-neg-frac297.6%
Simplified97.6%
if -3 < x < 3.5e9Initial program 99.9%
Taylor expanded in B around 0 99.4%
Final simplification98.6%
(FPCore (B x)
:precision binary64
(if (<= B 0.43)
(+
(* B 0.16666666666666666)
(+ (/ 1.0 B) (+ (* x (/ -1.0 B)) (* x (* B 0.3333333333333333)))))
(/ 1.0 (sin B))))
double code(double B, double x) {
double tmp;
if (B <= 0.43) {
tmp = (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333))));
} else {
tmp = 1.0 / sin(B);
}
return tmp;
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
real(8) :: tmp
if (b <= 0.43d0) then
tmp = (b * 0.16666666666666666d0) + ((1.0d0 / b) + ((x * ((-1.0d0) / b)) + (x * (b * 0.3333333333333333d0))))
else
tmp = 1.0d0 / sin(b)
end if
code = tmp
end function
public static double code(double B, double x) {
double tmp;
if (B <= 0.43) {
tmp = (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333))));
} else {
tmp = 1.0 / Math.sin(B);
}
return tmp;
}
def code(B, x): tmp = 0 if B <= 0.43: tmp = (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333)))) else: tmp = 1.0 / math.sin(B) return tmp
function code(B, x) tmp = 0.0 if (B <= 0.43) tmp = Float64(Float64(B * 0.16666666666666666) + Float64(Float64(1.0 / B) + Float64(Float64(x * Float64(-1.0 / B)) + Float64(x * Float64(B * 0.3333333333333333))))); else tmp = Float64(1.0 / sin(B)); end return tmp end
function tmp_2 = code(B, x) tmp = 0.0; if (B <= 0.43) tmp = (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333)))); else tmp = 1.0 / sin(B); end tmp_2 = tmp; end
code[B_, x_] := If[LessEqual[B, 0.43], N[(N[(B * 0.16666666666666666), $MachinePrecision] + N[(N[(1.0 / B), $MachinePrecision] + N[(N[(x * N[(-1.0 / B), $MachinePrecision]), $MachinePrecision] + N[(x * N[(B * 0.3333333333333333), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(1.0 / N[Sin[B], $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;B \leq 0.43:\\
\;\;\;\;B \cdot 0.16666666666666666 + \left(\frac{1}{B} + \left(x \cdot \frac{-1}{B} + x \cdot \left(B \cdot 0.3333333333333333\right)\right)\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{\sin B}\\
\end{array}
\end{array}
if B < 0.429999999999999993Initial program 99.8%
Taylor expanded in B around 0 74.1%
Taylor expanded in x around 0 74.1%
sub-neg74.1%
distribute-rgt-in74.1%
*-commutative74.1%
distribute-neg-frac74.1%
metadata-eval74.1%
Applied egg-rr74.1%
if 0.429999999999999993 < B Initial program 99.6%
Taylor expanded in x around 0 56.3%
Final simplification69.5%
(FPCore (B x) :precision binary64 (/ (- 1.0 x) (sin B)))
double code(double B, double x) {
return (1.0 - x) / sin(B);
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = (1.0d0 - x) / sin(b)
end function
public static double code(double B, double x) {
return (1.0 - x) / Math.sin(B);
}
def code(B, x): return (1.0 - x) / math.sin(B)
function code(B, x) return Float64(Float64(1.0 - x) / sin(B)) end
function tmp = code(B, x) tmp = (1.0 - x) / sin(B); end
code[B_, x_] := N[(N[(1.0 - x), $MachinePrecision] / N[Sin[B], $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{1 - x}{\sin B}
\end{array}
Initial program 99.7%
distribute-lft-neg-in99.7%
+-commutative99.7%
*-commutative99.7%
remove-double-neg99.7%
distribute-frac-neg299.7%
tan-neg99.7%
cancel-sign-sub-inv99.7%
*-commutative99.7%
associate-*r/99.8%
*-rgt-identity99.8%
tan-neg99.8%
distribute-neg-frac299.8%
distribute-neg-frac99.8%
remove-double-neg99.8%
Simplified99.8%
tan-quot99.7%
associate-/r/99.8%
Applied egg-rr99.8%
Taylor expanded in B around inf 99.7%
div-sub99.8%
Simplified99.8%
Taylor expanded in B around 0 80.6%
Final simplification80.6%
(FPCore (B x) :precision binary64 (+ (* B 0.16666666666666666) (+ (/ 1.0 B) (+ (* x (/ -1.0 B)) (* x (* B 0.3333333333333333))))))
double code(double B, double x) {
return (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333))));
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = (b * 0.16666666666666666d0) + ((1.0d0 / b) + ((x * ((-1.0d0) / b)) + (x * (b * 0.3333333333333333d0))))
end function
public static double code(double B, double x) {
return (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333))));
}
def code(B, x): return (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333))))
function code(B, x) return Float64(Float64(B * 0.16666666666666666) + Float64(Float64(1.0 / B) + Float64(Float64(x * Float64(-1.0 / B)) + Float64(x * Float64(B * 0.3333333333333333))))) end
function tmp = code(B, x) tmp = (B * 0.16666666666666666) + ((1.0 / B) + ((x * (-1.0 / B)) + (x * (B * 0.3333333333333333)))); end
code[B_, x_] := N[(N[(B * 0.16666666666666666), $MachinePrecision] + N[(N[(1.0 / B), $MachinePrecision] + N[(N[(x * N[(-1.0 / B), $MachinePrecision]), $MachinePrecision] + N[(x * N[(B * 0.3333333333333333), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
B \cdot 0.16666666666666666 + \left(\frac{1}{B} + \left(x \cdot \frac{-1}{B} + x \cdot \left(B \cdot 0.3333333333333333\right)\right)\right)
\end{array}
Initial program 99.7%
Taylor expanded in B around 0 56.0%
Taylor expanded in x around 0 56.0%
sub-neg56.0%
distribute-rgt-in56.0%
*-commutative56.0%
distribute-neg-frac56.0%
metadata-eval56.0%
Applied egg-rr56.0%
Final simplification56.0%
(FPCore (B x) :precision binary64 (+ (* B 0.16666666666666666) (+ (/ 1.0 B) (* x (+ (* B 0.3333333333333333) (/ -1.0 B))))))
double code(double B, double x) {
return (B * 0.16666666666666666) + ((1.0 / B) + (x * ((B * 0.3333333333333333) + (-1.0 / B))));
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = (b * 0.16666666666666666d0) + ((1.0d0 / b) + (x * ((b * 0.3333333333333333d0) + ((-1.0d0) / b))))
end function
public static double code(double B, double x) {
return (B * 0.16666666666666666) + ((1.0 / B) + (x * ((B * 0.3333333333333333) + (-1.0 / B))));
}
def code(B, x): return (B * 0.16666666666666666) + ((1.0 / B) + (x * ((B * 0.3333333333333333) + (-1.0 / B))))
function code(B, x) return Float64(Float64(B * 0.16666666666666666) + Float64(Float64(1.0 / B) + Float64(x * Float64(Float64(B * 0.3333333333333333) + Float64(-1.0 / B))))) end
function tmp = code(B, x) tmp = (B * 0.16666666666666666) + ((1.0 / B) + (x * ((B * 0.3333333333333333) + (-1.0 / B)))); end
code[B_, x_] := N[(N[(B * 0.16666666666666666), $MachinePrecision] + N[(N[(1.0 / B), $MachinePrecision] + N[(x * N[(N[(B * 0.3333333333333333), $MachinePrecision] + N[(-1.0 / B), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
B \cdot 0.16666666666666666 + \left(\frac{1}{B} + x \cdot \left(B \cdot 0.3333333333333333 + \frac{-1}{B}\right)\right)
\end{array}
Initial program 99.7%
Taylor expanded in B around 0 56.0%
Taylor expanded in x around 0 56.0%
Final simplification56.0%
(FPCore (B x) :precision binary64 (if (or (<= x -34000000.0) (not (<= x 1.0))) (/ x (- B)) (/ 1.0 B)))
double code(double B, double x) {
double tmp;
if ((x <= -34000000.0) || !(x <= 1.0)) {
tmp = x / -B;
} else {
tmp = 1.0 / B;
}
return tmp;
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
real(8) :: tmp
if ((x <= (-34000000.0d0)) .or. (.not. (x <= 1.0d0))) then
tmp = x / -b
else
tmp = 1.0d0 / b
end if
code = tmp
end function
public static double code(double B, double x) {
double tmp;
if ((x <= -34000000.0) || !(x <= 1.0)) {
tmp = x / -B;
} else {
tmp = 1.0 / B;
}
return tmp;
}
def code(B, x): tmp = 0 if (x <= -34000000.0) or not (x <= 1.0): tmp = x / -B else: tmp = 1.0 / B return tmp
function code(B, x) tmp = 0.0 if ((x <= -34000000.0) || !(x <= 1.0)) tmp = Float64(x / Float64(-B)); else tmp = Float64(1.0 / B); end return tmp end
function tmp_2 = code(B, x) tmp = 0.0; if ((x <= -34000000.0) || ~((x <= 1.0))) tmp = x / -B; else tmp = 1.0 / B; end tmp_2 = tmp; end
code[B_, x_] := If[Or[LessEqual[x, -34000000.0], N[Not[LessEqual[x, 1.0]], $MachinePrecision]], N[(x / (-B)), $MachinePrecision], N[(1.0 / B), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;x \leq -34000000 \lor \neg \left(x \leq 1\right):\\
\;\;\;\;\frac{x}{-B}\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{B}\\
\end{array}
\end{array}
if x < -3.4e7 or 1 < x Initial program 99.6%
Taylor expanded in B around 0 55.6%
Taylor expanded in x around inf 54.7%
neg-mul-154.7%
distribute-neg-frac54.7%
Simplified54.7%
if -3.4e7 < x < 1Initial program 99.8%
Taylor expanded in B around 0 55.8%
Taylor expanded in x around 0 55.1%
Final simplification55.0%
(FPCore (B x) :precision binary64 (/ (- 1.0 x) B))
double code(double B, double x) {
return (1.0 - x) / B;
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = (1.0d0 - x) / b
end function
public static double code(double B, double x) {
return (1.0 - x) / B;
}
def code(B, x): return (1.0 - x) / B
function code(B, x) return Float64(Float64(1.0 - x) / B) end
function tmp = code(B, x) tmp = (1.0 - x) / B; end
code[B_, x_] := N[(N[(1.0 - x), $MachinePrecision] / B), $MachinePrecision]
\begin{array}{l}
\\
\frac{1 - x}{B}
\end{array}
Initial program 99.7%
Taylor expanded in B around 0 55.7%
Final simplification55.7%
(FPCore (B x) :precision binary64 (/ 1.0 B))
double code(double B, double x) {
return 1.0 / B;
}
real(8) function code(b, x)
real(8), intent (in) :: b
real(8), intent (in) :: x
code = 1.0d0 / b
end function
public static double code(double B, double x) {
return 1.0 / B;
}
def code(B, x): return 1.0 / B
function code(B, x) return Float64(1.0 / B) end
function tmp = code(B, x) tmp = 1.0 / B; end
code[B_, x_] := N[(1.0 / B), $MachinePrecision]
\begin{array}{l}
\\
\frac{1}{B}
\end{array}
Initial program 99.7%
Taylor expanded in B around 0 55.7%
Taylor expanded in x around 0 31.6%
Final simplification31.6%
herbie shell --seed 2024059
(FPCore (B x)
:name "VandenBroeck and Keller, Equation (24)"
:precision binary64
(+ (- (* x (/ 1.0 (tan B)))) (/ 1.0 (sin B))))