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

Average Error: 12.7 → 0.7

Time: 29.0s

Precision: binary32

\[0.0001 \leq alphax \land alphax \leq 1 \land 0.0001 \leq alphay \land alphay \leq 1 \land 2.328306437 \cdot 10^{-10} \leq u0 \land u0 \leq 1 \land 0 \leq cos2phi \land cos2phi \leq 1 \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 := \frac{cos2phi}{alphax \cdot alphax}\\ \mathbf{if}\;1 - u0 \leq 0.9546453952789307:\\ \;\;\;\;\frac{-\log \left(1 - u0\right)}{t_0 + sin2phi \cdot \frac{1}{alphay \cdot alphay}}\\ \mathbf{else}:\\ \;\;\;\;\frac{u0 + \left(0.25 \cdot {u0}^{4} + u0 \cdot \left(u0 \cdot \left(0.5 + u0 \cdot 0.3333333333333333\right)\right)\right)}{t_0 + \frac{sin2phi}{alphay \cdot alphay}}\\ \end{array} \]

FPCore

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

↓

(FPCore (alphax alphay u0 cos2phi sin2phi)
 :precision binary32
 (let* ((t_0 (/ cos2phi (* alphax alphax))))
   (if (<= (- 1.0 u0) 0.9546453952789307)
     (/ (- (log (- 1.0 u0))) (+ t_0 (* sin2phi (/ 1.0 (* alphay alphay)))))
     (/
      (+
       u0
       (+
        (* 0.25 (pow u0 4.0))
        (* u0 (* u0 (+ 0.5 (* u0 0.3333333333333333))))))
      (+ t_0 (/ sin2phi (* alphay alphay)))))))

C

float code(float alphax, float alphay, float u0, float cos2phi, float sin2phi) {
	return -logf(1.0f - u0) / ((cos2phi / (alphax * alphax)) + (sin2phi / (alphay * alphay)));
}

↓

float code(float alphax, float alphay, float u0, float cos2phi, float sin2phi) {
	float t_0 = cos2phi / (alphax * alphax);
	float tmp;
	if ((1.0f - u0) <= 0.9546453952789307f) {
		tmp = -logf(1.0f - u0) / (t_0 + (sin2phi * (1.0f / (alphay * alphay))));
	} else {
		tmp = (u0 + ((0.25f * powf(u0, 4.0f)) + (u0 * (u0 * (0.5f + (u0 * 0.3333333333333333f)))))) / (t_0 + (sin2phi / (alphay * alphay)));
	}
	return tmp;
}

Error

Bits error versus alphax

Bits error versus alphay

Bits error versus u0

Bits error versus cos2phi

Bits error versus sin2phi

Derivation

1. Split input into 2 regimes

2. if (-.f32 1 u0) < 0.954645395

   1. Initial program 1.5

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

   3. Applied div-inv_binary32 1.6

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

   if 0.954645395 < (-.f32 1 u0)

   1. Initial program 14.6

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

   2. Taylor expanded around 0 0.5

      \[\leadsto \frac{-\color{blue}{\left(-\left(u0 + \left(0.3333333333333333 \cdot {u0}^{3} + \left(0.25 \cdot {u0}^{4} + 0.5 \cdot {u0}^{2}\right)\right)\right)\right)}}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]

   3. Simplified 0.5

      \[\leadsto \frac{-\color{blue}{\left(\left(-u0\right) - \left(\left(u0 \cdot u0\right) \cdot \left(0.5 + 0.3333333333333333 \cdot u0\right) + 0.25 \cdot {u0}^{4}\right)\right)}}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]
   4. Using strategy rm

   5. Applied *-un-lft-identity_binary32 0.5

      \[\leadsto \frac{-\left(\left(-u0\right) - \left(\left(u0 \cdot u0\right) \cdot \left(0.5 + 0.3333333333333333 \cdot u0\right) + 0.25 \cdot {u0}^{4}\right)\right)}{\color{blue}{1 \cdot \left(\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}\right)}} \]

   6. Applied associate-/r*_binary32 0.5

      \[\leadsto \color{blue}{\frac{\frac{-\left(\left(-u0\right) - \left(\left(u0 \cdot u0\right) \cdot \left(0.5 + 0.3333333333333333 \cdot u0\right) + 0.25 \cdot {u0}^{4}\right)\right)}{1}}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}}} \]

   7. Simplified 0.5

      \[\leadsto \frac{\color{blue}{-\left(\left(-u0\right) - \left(0.25 \cdot {u0}^{4} + u0 \cdot \left(u0 \cdot \left(0.5 + u0 \cdot 0.3333333333333333\right)\right)\right)\right)}}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}} \]
3. Recombined 2 regimes into one program.

4. Final simplification 0.7

   \[\leadsto \begin{array}{l} \mathbf{if}\;1 - u0 \leq 0.9546453952789307:\\ \;\;\;\;\frac{-\log \left(1 - u0\right)}{\frac{cos2phi}{alphax \cdot alphax} + sin2phi \cdot \frac{1}{alphay \cdot alphay}}\\ \mathbf{else}:\\ \;\;\;\;\frac{u0 + \left(0.25 \cdot {u0}^{4} + u0 \cdot \left(u0 \cdot \left(0.5 + u0 \cdot 0.3333333333333333\right)\right)\right)}{\frac{cos2phi}{alphax \cdot alphax} + \frac{sin2phi}{alphay \cdot alphay}}\\ \end{array} \]

Reproduce

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