Beckmann Distribution sample, tan2theta, alphax != alphay, u1 <= 0.5

Average Error: 13.2 → 0.6

Time: 1.8min

Precision: binary64

Cost: 7620

\[\left(\left(\left(\left(0.0001 \leq alphax \land alphax \leq 1\right) \land \left(0.0001 \leq alphay \land alphay \leq 1\right)\right) \land \left(2.328306437 \cdot 10^{-10} \leq u0 \land u0 \leq 1\right)\right) \land \left(0 \leq cos2phi \land cos2phi \leq 1\right)\right) \land 0 \leq sin2phi\]

Math

\[\frac{-\log \left(1 - u0\right)}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]

↓

\[\frac{\mathsf{log1p}\left(-u0\right)}{\frac{sin2phi}{-alphay \cdot alphay} - \begin{array}{l} \mathbf{if}\;cos2phi \ne 0:\\ \;\;\;\;\frac{1}{\frac{alphax \cdot alphax}{cos2phi}}\\ \mathbf{else}:\\ \;\;\;\;\frac{cos2phi}{alphax \cdot alphax}\\ \end{array}} \]

FPCore

(FPCore (alphax alphay u0 cos2phi sin2phi)
 :precision binary64
 (/
  (- (log (- 1.0 u0)))
  (+ (/ cos2phi (* alphax alphax)) (/ sin2phi (* alphay alphay)))))

↓

(FPCore (alphax alphay u0 cos2phi sin2phi)
 :precision binary64
 (/
  (log1p (- u0))
  (-
   (/ sin2phi (- (* alphay alphay)))
   (if (!= cos2phi 0.0)
     (/ 1.0 (/ (* alphax alphax) cos2phi))
     (/ cos2phi (* alphax alphax))))))

C

double code(double alphax, double alphay, double u0, double cos2phi, double sin2phi) {
	return -log((1.0 - u0)) / ((cos2phi / (alphax * alphax)) + (sin2phi / (alphay * alphay)));
}

↓

double code(double alphax, double alphay, double u0, double cos2phi, double sin2phi) {
	double tmp;
	if (cos2phi != 0.0) {
		tmp = 1.0 / ((alphax * alphax) / cos2phi);
	} else {
		tmp = cos2phi / (alphax * alphax);
	}
	return log1p(-u0) / ((sin2phi / -(alphay * alphay)) - tmp);
}

Java

public static double code(double alphax, double alphay, double u0, double cos2phi, double sin2phi) {
	return -Math.log((1.0 - u0)) / ((cos2phi / (alphax * alphax)) + (sin2phi / (alphay * alphay)));
}

↓

public static double code(double alphax, double alphay, double u0, double cos2phi, double sin2phi) {
	double tmp;
	if (cos2phi != 0.0) {
		tmp = 1.0 / ((alphax * alphax) / cos2phi);
	} else {
		tmp = cos2phi / (alphax * alphax);
	}
	return Math.log1p(-u0) / ((sin2phi / -(alphay * alphay)) - tmp);
}

Python

def code(alphax, alphay, u0, cos2phi, sin2phi):
	return -math.log((1.0 - u0)) / ((cos2phi / (alphax * alphax)) + (sin2phi / (alphay * alphay)))

↓

def code(alphax, alphay, u0, cos2phi, sin2phi):
	tmp = 0
	if cos2phi != 0.0:
		tmp = 1.0 / ((alphax * alphax) / cos2phi)
	else:
		tmp = cos2phi / (alphax * alphax)
	return math.log1p(-u0) / ((sin2phi / -(alphay * alphay)) - tmp)

Julia

function code(alphax, alphay, u0, cos2phi, sin2phi)
	return Float64(Float64(-log(Float64(1.0 - u0))) / Float64(Float64(cos2phi / Float64(alphax * alphax)) + Float64(sin2phi / Float64(alphay * alphay))))
end

↓

function code(alphax, alphay, u0, cos2phi, sin2phi)
	tmp = 0.0
	if (cos2phi != 0.0)
		tmp = Float64(1.0 / Float64(Float64(alphax * alphax) / cos2phi));
	else
		tmp = Float64(cos2phi / Float64(alphax * alphax));
	end
	return Float64(log1p(Float64(-u0)) / Float64(Float64(sin2phi / Float64(-Float64(alphay * alphay))) - tmp))
end

Wolfram

code[alphax_, alphay_, u0_, cos2phi_, sin2phi_] := N[((-N[Log[N[(1.0 - u0), $MachinePrecision]], $MachinePrecision]) / N[(N[(cos2phi / N[(alphax * alphax), $MachinePrecision]), $MachinePrecision] + N[(sin2phi / N[(alphay * alphay), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

↓

code[alphax_, alphay_, u0_, cos2phi_, sin2phi_] := N[(N[Log[1 + (-u0)], $MachinePrecision] / N[(N[(sin2phi / (-N[(alphay * alphay), $MachinePrecision])), $MachinePrecision] - If[Unequal[cos2phi, 0.0], N[(1.0 / N[(N[(alphax * alphax), $MachinePrecision] / cos2phi), $MachinePrecision]), $MachinePrecision], N[(cos2phi / N[(alphax * alphax), $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]), $MachinePrecision]

Derivation

1. Initial program 13.2

   \[\frac{-\log \left(1 - u0\right)}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]

2. Applied egg-rr 0.6

   \[\leadsto \frac{-\color{blue}{\mathsf{log1p}\left(-u0\right)}}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]

3. Applied egg-rr 0.6

   \[\leadsto \frac{-\mathsf{log1p}\left(-u0\right)}{\color{blue}{\begin{array}{l} \color{blue}{\mathbf{if}\;cos2phi \ne 0:\\ \;\;\;\;\frac{1}{\frac{{alphax}^{2}}{cos2phi}}\\ \mathbf{else}:\\ \;\;\;\;\frac{cos2phi}{{alphax}^{2}}\\ } \end{array}} + \frac{sin2phi}{alphay \cdot alphay}} \]

4. Applied egg-rr 23.3

   \[\leadsto \color{blue}{\sqrt[3]{\frac{-\mathsf{log1p}\left(-u0\right)}{\begin{array}{l} \mathbf{if}\;cos2phi \ne 0:\\ \;\;\;\;\frac{1}{\frac{{alphax}^{2}}{cos2phi}}\\ \mathbf{else}:\\ \;\;\;\;\frac{cos2phi}{{alphax}^{2}}\\ \end{array} + \frac{sin2phi}{{alphay}^{2}}}} \cdot \sqrt[3]{{\left(\frac{-\mathsf{log1p}\left(-u0\right)}{\begin{array}{l} \mathbf{if}\;cos2phi \ne 0:\\ \;\;\;\;\frac{1}{\frac{{alphax}^{2}}{cos2phi}}\\ \mathbf{else}:\\ \;\;\;\;\frac{cos2phi}{{alphax}^{2}}\\ \end{array} + \frac{sin2phi}{{alphay}^{2}}}\right)}^{2}}} \]

5. Simplified 0.6

   \[\leadsto \color{blue}{\frac{\mathsf{log1p}\left(-u0\right)}{\frac{sin2phi}{-alphay \cdot alphay} - \begin{array}{l} \mathbf{if}\;cos2phi \ne 0:\\ \;\;\;\;\frac{1}{\frac{alphax \cdot alphax}{cos2phi}}\\ \mathbf{else}:\\ \;\;\;\;\frac{cos2phi}{alphax \cdot alphax}\\ \end{array}}} \]
   Proof

Alternatives

Alternative 1
Error	7.8
Cost	7752

\[\begin{array}{l} t_0 := \mathsf{log1p}\left(-u0\right)\\ t_1 := \frac{sin2phi}{alphay \cdot alphay}\\ \mathbf{if}\;t_1 \leq 10^{-188}:\\ \;\;\;\;\left(\left(-alphax\right) \cdot \frac{alphax}{cos2phi}\right) \cdot t_0\\ \mathbf{elif}\;t_1 \leq 2 \cdot 10^{-67}:\\ \;\;\;\;\frac{0.5 \cdot \left(u0 \cdot u0\right) + u0}{\frac{cos2phi}{alphax \cdot alphax} + t_1}\\ \mathbf{else}:\\ \;\;\;\;\left(-alphay\right) \cdot \left(t_0 \cdot \frac{alphay}{sin2phi}\right)\\ \end{array} \]

Alternative 2
Error	7.8
Cost	7752

\[\begin{array}{l} t_0 := \mathsf{log1p}\left(-u0\right)\\ t_1 := \frac{sin2phi}{alphay \cdot alphay}\\ \mathbf{if}\;t_1 \leq 10^{-188}:\\ \;\;\;\;\frac{alphax \cdot alphax}{cos2phi} \cdot \left(-t_0\right)\\ \mathbf{elif}\;t_1 \leq 2 \cdot 10^{-67}:\\ \;\;\;\;\frac{0.5 \cdot \left(u0 \cdot u0\right) + u0}{\frac{cos2phi}{alphax \cdot alphax} + t_1}\\ \mathbf{else}:\\ \;\;\;\;\left(-alphay\right) \cdot \left(t_0 \cdot \frac{alphay}{sin2phi}\right)\\ \end{array} \]

Alternative 3
Error	9.1
Cost	7364

\[\begin{array}{l} t_0 := \frac{sin2phi}{alphay \cdot alphay}\\ \mathbf{if}\;t_0 \leq 2 \cdot 10^{-67}:\\ \;\;\;\;\frac{0.5 \cdot \left(u0 \cdot u0\right) + u0}{\frac{\frac{cos2phi}{alphax}}{alphax} + t_0}\\ \mathbf{else}:\\ \;\;\;\;\left(-alphay\right) \cdot \left(\mathsf{log1p}\left(-u0\right) \cdot \frac{alphay}{sin2phi}\right)\\ \end{array} \]

Alternative 4
Error	0.6
Cost	7360

\[\frac{-\mathsf{log1p}\left(-u0\right)}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]

Alternative 5
Error	34.2
Cost	1476

\[\begin{array}{l} \mathbf{if}\;\frac{sin2phi}{alphay \cdot alphay} \leq 10^{-188}:\\ \;\;\;\;\frac{\left(alphax \cdot alphax\right) \cdot u0 + 0.5 \cdot \left(\left(u0 \cdot alphax\right) \cdot \left(u0 \cdot alphax\right)\right)}{cos2phi}\\ \mathbf{else}:\\ \;\;\;\;\frac{u0}{\frac{cos2phi}{alphax \cdot alphax} + \frac{\frac{sin2phi}{alphay}}{alphay}}\\ \end{array} \]

Alternative 6
Error	20.6
Cost	1216

\[\frac{\left(0.5 \cdot \left(u0 + 2\right)\right) \cdot u0}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]

Alternative 7
Error	35.4
Cost	832

\[\frac{u0}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]

Alternative 8
Error	35.4
Cost	832

\[\frac{u0}{\frac{cos2phi}{alphax \cdot alphax} + \frac{\frac{sin2phi}{alphay}}{alphay}} \]

Alternative 9
Error	38.4
Cost	772

\[\begin{array}{l} \mathbf{if}\;sin2phi \leq 5.5 \cdot 10^{-156}:\\ \;\;\;\;\frac{\left(alphax \cdot alphax\right) \cdot u0}{cos2phi}\\ \mathbf{else}:\\ \;\;\;\;\frac{-1}{sin2phi} \cdot \left(alphay \cdot \left(\left(-u0\right) \cdot alphay\right)\right)\\ \end{array} \]

Alternative 10
Error	38.4
Cost	580

\[\begin{array}{l} \mathbf{if}\;sin2phi \leq 2.15 \cdot 10^{-157}:\\ \;\;\;\;\frac{\left(alphax \cdot alphax\right) \cdot u0}{cos2phi}\\ \mathbf{else}:\\ \;\;\;\;\frac{u0 \cdot \left(alphay \cdot alphay\right)}{sin2phi}\\ \end{array} \]

Alternative 11
Error	41.4
Cost	448

\[\frac{alphay}{sin2phi} \cdot \left(u0 \cdot alphay\right) \]

Alternative 12
Error	41.4
Cost	448

\[\frac{u0 \cdot \left(alphay \cdot alphay\right)}{sin2phi} \]

Reproduce

herbie shell --seed 2023033 
(FPCore (alphax alphay u0 cos2phi sin2phi)
  :name "Beckmann Distribution sample, tan2theta, alphax != alphay, u1 <= 0.5"
  :precision binary64
  :pre (and (and (and (and (and (<= 0.0001 alphax) (<= alphax 1.0)) (and (<= 0.0001 alphay) (<= alphay 1.0))) (and (<= 2.328306437e-10 u0) (<= u0 1.0))) (and (<= 0.0 cos2phi) (<= cos2phi 1.0))) (<= 0.0 sin2phi))
  (/ (- (log (- 1.0 u0))) (+ (/ cos2phi (* alphax alphax)) (/ sin2phi (* alphay alphay)))))