Beckmann Sample, normalization factor

Average Accuracy: 99.4% → 99.6%

Time: 1.6min

Precision: binary64

\[\left(0 < cosTheta \land cosTheta < 0.9999\right) \land \left(-1 < c \land c < 1\right)\]

Math

\[\frac{1}{\left(1 + c\right) + \left(\frac{1}{\sqrt{\pi}} \cdot \frac{\sqrt{\left(1 - cosTheta\right) - cosTheta}}{cosTheta}\right) \cdot e^{\left(-cosTheta\right) \cdot cosTheta}} \]

↓

\[\begin{array}{l} t_0 := \sqrt{1 - \left(cosTheta + cosTheta\right)}\\ t_1 := {\left(\sqrt{\pi}\right)}^{-1}\\ \frac{1}{\left(1 + c\right) + \begin{array}{l} \mathbf{if}\;t_1 \ne 0:\\ \;\;\;\;\frac{t_0}{\frac{cosTheta}{t_1}}\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{t_0}{\sqrt{\pi}}}{cosTheta}\\ \end{array} \cdot e^{\left(-cosTheta\right) \cdot cosTheta}} \end{array} \]

FPCore

(FPCore (cosTheta c)
 :precision binary64
 (/
  1.0
  (+
   (+ 1.0 c)
   (*
    (* (/ 1.0 (sqrt PI)) (/ (sqrt (- (- 1.0 cosTheta) cosTheta)) cosTheta))
    (exp (* (- cosTheta) cosTheta))))))

↓

(FPCore (cosTheta c)
 :precision binary64
 (let* ((t_0 (sqrt (- 1.0 (+ cosTheta cosTheta)))) (t_1 (pow (sqrt PI) -1.0)))
   (/
    1.0
    (+
     (+ 1.0 c)
     (*
      (if (!= t_1 0.0) (/ t_0 (/ cosTheta t_1)) (/ (/ t_0 (sqrt PI)) cosTheta))
      (exp (* (- cosTheta) cosTheta)))))))

C

double code(double cosTheta, double c) {
	return 1.0 / ((1.0 + c) + (((1.0 / sqrt(((double) M_PI))) * (sqrt(((1.0 - cosTheta) - cosTheta)) / cosTheta)) * exp((-cosTheta * cosTheta))));
}

↓

double code(double cosTheta, double c) {
	double t_0 = sqrt((1.0 - (cosTheta + cosTheta)));
	double t_1 = pow(sqrt(((double) M_PI)), -1.0);
	double tmp;
	if (t_1 != 0.0) {
		tmp = t_0 / (cosTheta / t_1);
	} else {
		tmp = (t_0 / sqrt(((double) M_PI))) / cosTheta;
	}
	return 1.0 / ((1.0 + c) + (tmp * exp((-cosTheta * cosTheta))));
}

Java

public static double code(double cosTheta, double c) {
	return 1.0 / ((1.0 + c) + (((1.0 / Math.sqrt(Math.PI)) * (Math.sqrt(((1.0 - cosTheta) - cosTheta)) / cosTheta)) * Math.exp((-cosTheta * cosTheta))));
}

↓

public static double code(double cosTheta, double c) {
	double t_0 = Math.sqrt((1.0 - (cosTheta + cosTheta)));
	double t_1 = Math.pow(Math.sqrt(Math.PI), -1.0);
	double tmp;
	if (t_1 != 0.0) {
		tmp = t_0 / (cosTheta / t_1);
	} else {
		tmp = (t_0 / Math.sqrt(Math.PI)) / cosTheta;
	}
	return 1.0 / ((1.0 + c) + (tmp * Math.exp((-cosTheta * cosTheta))));
}

Python

def code(cosTheta, c):
	return 1.0 / ((1.0 + c) + (((1.0 / math.sqrt(math.pi)) * (math.sqrt(((1.0 - cosTheta) - cosTheta)) / cosTheta)) * math.exp((-cosTheta * cosTheta))))

↓

def code(cosTheta, c):
	t_0 = math.sqrt((1.0 - (cosTheta + cosTheta)))
	t_1 = math.pow(math.sqrt(math.pi), -1.0)
	tmp = 0
	if t_1 != 0.0:
		tmp = t_0 / (cosTheta / t_1)
	else:
		tmp = (t_0 / math.sqrt(math.pi)) / cosTheta
	return 1.0 / ((1.0 + c) + (tmp * math.exp((-cosTheta * cosTheta))))

Julia

function code(cosTheta, c)
	return Float64(1.0 / Float64(Float64(1.0 + c) + Float64(Float64(Float64(1.0 / sqrt(pi)) * Float64(sqrt(Float64(Float64(1.0 - cosTheta) - cosTheta)) / cosTheta)) * exp(Float64(Float64(-cosTheta) * cosTheta)))))
end

↓

function code(cosTheta, c)
	t_0 = sqrt(Float64(1.0 - Float64(cosTheta + cosTheta)))
	t_1 = sqrt(pi) ^ -1.0
	tmp = 0.0
	if (t_1 != 0.0)
		tmp = Float64(t_0 / Float64(cosTheta / t_1));
	else
		tmp = Float64(Float64(t_0 / sqrt(pi)) / cosTheta);
	end
	return Float64(1.0 / Float64(Float64(1.0 + c) + Float64(tmp * exp(Float64(Float64(-cosTheta) * cosTheta)))))
end

MATLAB

function tmp = code(cosTheta, c)
	tmp = 1.0 / ((1.0 + c) + (((1.0 / sqrt(pi)) * (sqrt(((1.0 - cosTheta) - cosTheta)) / cosTheta)) * exp((-cosTheta * cosTheta))));
end

↓

function tmp_2 = code(cosTheta, c)
	t_0 = sqrt((1.0 - (cosTheta + cosTheta)));
	t_1 = sqrt(pi) ^ -1.0;
	tmp = 0.0;
	if (t_1 ~= 0.0)
		tmp = t_0 / (cosTheta / t_1);
	else
		tmp = (t_0 / sqrt(pi)) / cosTheta;
	end
	tmp_2 = 1.0 / ((1.0 + c) + (tmp * exp((-cosTheta * cosTheta))));
end

Wolfram

code[cosTheta_, c_] := N[(1.0 / N[(N[(1.0 + c), $MachinePrecision] + N[(N[(N[(1.0 / N[Sqrt[Pi], $MachinePrecision]), $MachinePrecision] * N[(N[Sqrt[N[(N[(1.0 - cosTheta), $MachinePrecision] - cosTheta), $MachinePrecision]], $MachinePrecision] / cosTheta), $MachinePrecision]), $MachinePrecision] * N[Exp[N[((-cosTheta) * cosTheta), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

↓

code[cosTheta_, c_] := Block[{t$95$0 = N[Sqrt[N[(1.0 - N[(cosTheta + cosTheta), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = N[Power[N[Sqrt[Pi], $MachinePrecision], -1.0], $MachinePrecision]}, N[(1.0 / N[(N[(1.0 + c), $MachinePrecision] + N[(If[Unequal[t$95$1, 0.0], N[(t$95$0 / N[(cosTheta / t$95$1), $MachinePrecision]), $MachinePrecision], N[(N[(t$95$0 / N[Sqrt[Pi], $MachinePrecision]), $MachinePrecision] / cosTheta), $MachinePrecision]] * N[Exp[N[((-cosTheta) * cosTheta), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

1. Initial program 99.4%

   \[\frac{1}{\left(1 + c\right) + \left(\frac{1}{\sqrt{\pi}} \cdot \frac{\sqrt{\left(1 - cosTheta\right) - cosTheta}}{cosTheta}\right) \cdot e^{\left(-cosTheta\right) \cdot cosTheta}} \]

2. Applied egg-rr 99.4%

   \[\leadsto \frac{1}{\left(1 + c\right) + \color{blue}{\left(\frac{\sqrt{1 - \left(cosTheta + cosTheta\right)}}{cosTheta} \cdot {\left(\sqrt{\pi}\right)}^{-1}\right)} \cdot e^{\left(-cosTheta\right) \cdot cosTheta}} \]

3. Simplified 99.5%

   \[\leadsto \frac{1}{\left(1 + c\right) + \color{blue}{\frac{\sqrt{\left(1 - cosTheta\right) - cosTheta} \cdot {\left(\sqrt{\pi}\right)}^{-1}}{cosTheta}} \cdot e^{\left(-cosTheta\right) \cdot cosTheta}} \]
   Proof

4. Applied egg-rr 99.6%

   \[\leadsto \frac{1}{\left(1 + c\right) + \color{blue}{\begin{array}{l} \color{blue}{\mathbf{if}\;{\left(\sqrt{\pi}\right)}^{-1} \ne 0:\\ \;\;\;\;\frac{\sqrt{1 - \left(cosTheta + cosTheta\right)}}{\frac{cosTheta}{{\left(\sqrt{\pi}\right)}^{-1}}}\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{\sqrt{1 - \left(cosTheta + cosTheta\right)}}{\sqrt{\pi}}}{cosTheta}\\ } \end{array}} \cdot e^{\left(-cosTheta\right) \cdot cosTheta}} \]

Reproduce

herbie shell --seed 2023187 
(FPCore (cosTheta c)
  :name "Beckmann Sample, normalization factor"
  :precision binary64
  :pre (and (and (< 0.0 cosTheta) (< cosTheta 0.9999)) (and (< -1.0 c) (< c 1.0)))
  (/ 1.0 (+ (+ 1.0 c) (* (* (/ 1.0 (sqrt PI)) (/ (sqrt (- (- 1.0 cosTheta) cosTheta)) cosTheta)) (exp (* (- cosTheta) cosTheta))))))