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

Average Error: 13.3 → 0.5

Time: 38.4s

Precision: binary64

Cost: 7876

\[\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}} \]

↓

\[\begin{array}{l} t_0 := \mathsf{log1p}\left(-u0\right)\\ \mathbf{if}\;alphay \cdot alphay \ne 0:\\ \;\;\;\;\frac{-t_0}{\frac{cos2phi \cdot alphay}{alphax} \cdot \frac{alphay}{alphax} + sin2phi} \cdot \left(alphay \cdot alphay\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{t_0}{\left(-\frac{cos2phi}{alphax \cdot alphax}\right) - \frac{sin2phi}{alphay \cdot alphay}}\\ \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
 (let* ((t_0 (log1p (- u0))))
   (if (!= (* alphay alphay) 0.0)
     (*
      (/
       (- t_0)
       (+ (* (/ (* cos2phi alphay) alphax) (/ alphay alphax)) sin2phi))
      (* alphay alphay))
     (/
      t_0
      (- (- (/ cos2phi (* alphax alphax))) (/ sin2phi (* alphay alphay)))))))

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 t_0 = log1p(-u0);
	double tmp;
	if ((alphay * alphay) != 0.0) {
		tmp = (-t_0 / ((((cos2phi * alphay) / alphax) * (alphay / alphax)) + sin2phi)) * (alphay * alphay);
	} else {
		tmp = t_0 / (-(cos2phi / (alphax * alphax)) - (sin2phi / (alphay * alphay)));
	}
	return 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 t_0 = Math.log1p(-u0);
	double tmp;
	if ((alphay * alphay) != 0.0) {
		tmp = (-t_0 / ((((cos2phi * alphay) / alphax) * (alphay / alphax)) + sin2phi)) * (alphay * alphay);
	} else {
		tmp = t_0 / (-(cos2phi / (alphax * alphax)) - (sin2phi / (alphay * alphay)));
	}
	return 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):
	t_0 = math.log1p(-u0)
	tmp = 0
	if (alphay * alphay) != 0.0:
		tmp = (-t_0 / ((((cos2phi * alphay) / alphax) * (alphay / alphax)) + sin2phi)) * (alphay * alphay)
	else:
		tmp = t_0 / (-(cos2phi / (alphax * alphax)) - (sin2phi / (alphay * alphay)))
	return 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)
	t_0 = log1p(Float64(-u0))
	tmp = 0.0
	if (Float64(alphay * alphay) != 0.0)
		tmp = Float64(Float64(Float64(-t_0) / Float64(Float64(Float64(Float64(cos2phi * alphay) / alphax) * Float64(alphay / alphax)) + sin2phi)) * Float64(alphay * alphay));
	else
		tmp = Float64(t_0 / Float64(Float64(-Float64(cos2phi / Float64(alphax * alphax))) - Float64(sin2phi / Float64(alphay * alphay))));
	end
	return 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_] := Block[{t$95$0 = N[Log[1 + (-u0)], $MachinePrecision]}, If[Unequal[N[(alphay * alphay), $MachinePrecision], 0.0], N[(N[((-t$95$0) / N[(N[(N[(N[(cos2phi * alphay), $MachinePrecision] / alphax), $MachinePrecision] * N[(alphay / alphax), $MachinePrecision]), $MachinePrecision] + sin2phi), $MachinePrecision]), $MachinePrecision] * N[(alphay * alphay), $MachinePrecision]), $MachinePrecision], N[(t$95$0 / N[((-N[(cos2phi / N[(alphax * alphax), $MachinePrecision]), $MachinePrecision]) - N[(sin2phi / N[(alphay * alphay), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

1. Initial program 13.3

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

2. Applied egg-rr 0.5

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

3. Taylor expanded in cos2phi around 0 0.4

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

4. Simplified 0.5

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

5. Applied egg-rr 0.5

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

Alternatives

Alternative 1
Error	8.1
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^{-238}:\\ \;\;\;\;-\frac{\left(alphax \cdot alphax\right) \cdot t_0}{cos2phi}\\ \mathbf{elif}\;t_1 \leq 5 \cdot 10^{-45}:\\ \;\;\;\;\frac{\left(-u0\right) + -0.5 \cdot \left(u0 \cdot u0\right)}{\left(-\frac{cos2phi}{alphax} \cdot \frac{1}{alphax}\right) - t_1}\\ \mathbf{else}:\\ \;\;\;\;-\frac{\left(alphay \cdot alphay\right) \cdot t_0}{sin2phi}\\ \end{array} \]

Alternative 2
Error	0.6
Cost	7488

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

Alternative 3
Error	0.6
Cost	7488

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

Alternative 4
Error	19.0
Cost	7364

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

Alternative 5
Error	0.6
Cost	7360

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

Alternative 6
Error	0.6
Cost	7360

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

Alternative 7
Error	20.3
Cost	2688

\[\begin{array}{l} t_0 := \left(-\frac{cos2phi \cdot alphay}{alphax}\right) - \frac{sin2phi \cdot alphax}{alphay}\\ \left(-\frac{\left(alphay \cdot alphax\right) \cdot u0}{t_0}\right) + \frac{-0.5 \cdot \left(alphay \cdot \left(\left(u0 \cdot u0\right) \cdot alphax\right)\right)}{t_0} \end{array} \]

Alternative 8
Error	20.4
Cost	1536

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

Alternative 9
Error	20.4
Cost	1536

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

Alternative 10
Error	34.3
Cost	1348

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

Alternative 11
Error	20.4
Cost	1344

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

Alternative 12
Error	20.4
Cost	1344

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

Alternative 13
Error	26.0
Cost	1092

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

Alternative 14
Error	35.4
Cost	960

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

Alternative 15
Error	38.6
Cost	836

\[\begin{array}{l} \mathbf{if}\;sin2phi \leq 2.5 \cdot 10^{-115}:\\ \;\;\;\;\frac{u0}{cos2phi} \cdot \frac{1}{\frac{\frac{1}{alphax}}{alphax}}\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{u0}{\frac{\frac{1}{alphay}}{alphay}}}{sin2phi}\\ \end{array} \]

Alternative 16
Error	35.4
Cost	832

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

Alternative 17
Error	35.4
Cost	832

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

Alternative 18
Error	38.6
Cost	708

\[\begin{array}{l} \mathbf{if}\;sin2phi \leq 3 \cdot 10^{-115}:\\ \;\;\;\;\frac{u0}{\frac{\frac{cos2phi}{alphax}}{alphax}}\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{u0}{sin2phi}}{\frac{\frac{1}{alphay}}{alphay}}\\ \end{array} \]

Alternative 19
Error	38.6
Cost	708

\[\begin{array}{l} \mathbf{if}\;sin2phi \leq 3 \cdot 10^{-115}:\\ \;\;\;\;\frac{u0}{\frac{\frac{cos2phi}{alphax}}{alphax}}\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{u0}{\frac{\frac{1}{alphay}}{alphay}}}{sin2phi}\\ \end{array} \]

Alternative 20
Error	38.6
Cost	580

\[\begin{array}{l} \mathbf{if}\;sin2phi \leq 4.8 \cdot 10^{-112}:\\ \;\;\;\;\frac{alphax \cdot alphax}{cos2phi} \cdot u0\\ \mathbf{else}:\\ \;\;\;\;\frac{alphay \cdot alphay}{sin2phi} \cdot u0\\ \end{array} \]

Alternative 21
Error	38.6
Cost	580

\[\begin{array}{l} \mathbf{if}\;sin2phi \leq 8.5 \cdot 10^{-115}:\\ \;\;\;\;\frac{u0}{\frac{\frac{cos2phi}{alphax}}{alphax}}\\ \mathbf{else}:\\ \;\;\;\;\frac{alphay \cdot alphay}{sin2phi} \cdot u0\\ \end{array} \]

Alternative 22
Error	38.6
Cost	580

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

Alternative 23
Error	54.8
Cost	448

\[\frac{alphax \cdot alphax}{cos2phi} \cdot u0 \]

Reproduce

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