Disney BSSRDF, PDF of scattering profile

Percentage Accurate: 99.6% → 99.6%

Time: 15.0s

Alternatives: 22

Speedup: 0.7×

Specification

\[\left(0 \leq s \land s \leq 256\right) \land \left(10^{-6} < r \land r < 1000000\right)\]

Math

\[\begin{array}{l} \\ \frac{0.25 \cdot e^{\frac{-r}{s}}}{\left(\left(2 \cdot \pi\right) \cdot s\right) \cdot r} + \frac{0.75 \cdot e^{\frac{-r}{3 \cdot s}}}{\left(\left(6 \cdot \pi\right) \cdot s\right) \cdot r} \end{array} \]

FPCore

(FPCore (s r)
 :precision binary32
 (+
  (/ (* 0.25 (exp (/ (- r) s))) (* (* (* 2.0 PI) s) r))
  (/ (* 0.75 (exp (/ (- r) (* 3.0 s)))) (* (* (* 6.0 PI) s) r))))

C

float code(float s, float r) {
	return ((0.25f * expf((-r / s))) / (((2.0f * ((float) M_PI)) * s) * r)) + ((0.75f * expf((-r / (3.0f * s)))) / (((6.0f * ((float) M_PI)) * s) * r));
}

Julia

function code(s, r)
	return Float32(Float32(Float32(Float32(0.25) * exp(Float32(Float32(-r) / s))) / Float32(Float32(Float32(Float32(2.0) * Float32(pi)) * s) * r)) + Float32(Float32(Float32(0.75) * exp(Float32(Float32(-r) / Float32(Float32(3.0) * s)))) / Float32(Float32(Float32(Float32(6.0) * Float32(pi)) * s) * r)))
end

MATLAB

function tmp = code(s, r)
	tmp = ((single(0.25) * exp((-r / s))) / (((single(2.0) * single(pi)) * s) * r)) + ((single(0.75) * exp((-r / (single(3.0) * s)))) / (((single(6.0) * single(pi)) * s) * r));
end

Initial Program: 99.6% accurate, 1.0× speedup

Math

\[\begin{array}{l} \\ \frac{0.25 \cdot e^{\frac{-r}{s}}}{\left(\left(2 \cdot \pi\right) \cdot s\right) \cdot r} + \frac{0.75 \cdot e^{\frac{-r}{3 \cdot s}}}{\left(\left(6 \cdot \pi\right) \cdot s\right) \cdot r} \end{array} \]

FPCore

(FPCore (s r)
 :precision binary32
 (+
  (/ (* 0.25 (exp (/ (- r) s))) (* (* (* 2.0 PI) s) r))
  (/ (* 0.75 (exp (/ (- r) (* 3.0 s)))) (* (* (* 6.0 PI) s) r))))

C

float code(float s, float r) {
	return ((0.25f * expf((-r / s))) / (((2.0f * ((float) M_PI)) * s) * r)) + ((0.75f * expf((-r / (3.0f * s)))) / (((6.0f * ((float) M_PI)) * s) * r));
}

Julia

function code(s, r)
	return Float32(Float32(Float32(Float32(0.25) * exp(Float32(Float32(-r) / s))) / Float32(Float32(Float32(Float32(2.0) * Float32(pi)) * s) * r)) + Float32(Float32(Float32(0.75) * exp(Float32(Float32(-r) / Float32(Float32(3.0) * s)))) / Float32(Float32(Float32(Float32(6.0) * Float32(pi)) * s) * r)))
end

MATLAB

function tmp = code(s, r)
	tmp = ((single(0.25) * exp((-r / s))) / (((single(2.0) * single(pi)) * s) * r)) + ((single(0.75) * exp((-r / (single(3.0) * s)))) / (((single(6.0) * single(pi)) * s) * r));
end

Alternative 1: 99.6% accurate, 0.7× speedup

Math

\[\begin{array}{l} \\ \frac{0.125}{r} \cdot \mathsf{fma}\left({e}^{\left(\frac{-0.16666666666666666 \cdot r}{s}\right)}, \frac{{e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}}{\pi \cdot s}, \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \end{array} \]

FPCore

(FPCore (s r)
 :precision binary32
 (*
  (/ 0.125 r)
  (fma
   (pow E (/ (* -0.16666666666666666 r) s))
   (/ (pow E (* (/ r s) -0.16666666666666666)) (* PI s))
   (/ (exp (/ (- r) s)) (* PI s)))))

C

float code(float s, float r) {
	return (0.125f / r) * fmaf(powf(((float) M_E), ((-0.16666666666666666f * r) / s)), (powf(((float) M_E), ((r / s) * -0.16666666666666666f)) / (((float) M_PI) * s)), (expf((-r / s)) / (((float) M_PI) * s)));
}

Julia

function code(s, r)
	return Float32(Float32(Float32(0.125) / r) * fma((Float32(exp(1)) ^ Float32(Float32(Float32(-0.16666666666666666) * r) / s)), Float32((Float32(exp(1)) ^ Float32(Float32(r / s) * Float32(-0.16666666666666666))) / Float32(Float32(pi) * s)), Float32(exp(Float32(Float32(-r) / s)) / Float32(Float32(pi) * s))))
end

Derivation

1. Initial program 99.6%

   \[\frac{0.25 \cdot e^{\frac{-r}{s}}}{\left(\left(2 \cdot \pi\right) \cdot s\right) \cdot r} + \frac{0.75 \cdot e^{\frac{-r}{3 \cdot s}}}{\left(\left(6 \cdot \pi\right) \cdot s\right) \cdot r} \]
2. Add Preprocessing

3. Taylor expanded in r around inf

   \[\leadsto \color{blue}{\frac{\frac{1}{8} \cdot \frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{1}{8} \cdot \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}}{r}} \]
4. Step-by-step derivation

   1. distribute-lft-out N/A

      \[\leadsto \frac{\color{blue}{\frac{1}{8} \cdot \left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right)}}{r} \]

   2. *-commutative N/A

      \[\leadsto \frac{\color{blue}{\left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right) \cdot \frac{1}{8}}}{r} \]

   3. associate-/l* N/A

      \[\leadsto \color{blue}{\left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right) \cdot \frac{\frac{1}{8}}{r}} \]

   4. lower-*.f32 N/A

      \[\leadsto \color{blue}{\left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right) \cdot \frac{\frac{1}{8}}{r}} \]

5. Applied rewrites 99.7%

   \[\leadsto \color{blue}{\left(\frac{e^{-0.3333333333333333 \cdot \frac{r}{s}}}{\pi \cdot s} + \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{0.125}{r}} \]
6. Step-by-step derivation

7. Applied rewrites 99.7%

   \[\leadsto \left(\frac{{e}^{\left(\frac{-0.3333333333333333}{s} \cdot r\right)}}{\pi \cdot s} + \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{0.125}{r} \]
8. Step-by-step derivation

9. Applied rewrites 99.7%

   \[\leadsto \mathsf{fma}\left({e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}, \frac{{e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}}{\pi \cdot s}, \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{\color{blue}{0.125}}{r} \]
10. Step-by-step derivation

11. Applied rewrites 99.8%

    \[\leadsto \mathsf{fma}\left({e}^{\left(\frac{-0.16666666666666666 \cdot r}{s}\right)}, \frac{{e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}}{\pi \cdot s}, \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{0.125}{r} \]

12. Final simplification 99.8%

    \[\leadsto \frac{0.125}{r} \cdot \mathsf{fma}\left({e}^{\left(\frac{-0.16666666666666666 \cdot r}{s}\right)}, \frac{{e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}}{\pi \cdot s}, \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \]
13. Add Preprocessing

Alternative 2: 99.6% accurate, 0.7× speedup

Math

\[\begin{array}{l} \\ \begin{array}{l} t_0 := {e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}\\ \mathsf{fma}\left(t\_0, \frac{t\_0}{\pi \cdot s}, \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{0.125}{r} \end{array} \end{array} \]

FPCore

(FPCore (s r)
 :precision binary32
 (let* ((t_0 (pow E (* (/ r s) -0.16666666666666666))))
   (* (fma t_0 (/ t_0 (* PI s)) (/ (exp (/ (- r) s)) (* PI s))) (/ 0.125 r))))

C

float code(float s, float r) {
	float t_0 = powf(((float) M_E), ((r / s) * -0.16666666666666666f));
	return fmaf(t_0, (t_0 / (((float) M_PI) * s)), (expf((-r / s)) / (((float) M_PI) * s))) * (0.125f / r);
}

Julia

function code(s, r)
	t_0 = Float32(exp(1)) ^ Float32(Float32(r / s) * Float32(-0.16666666666666666))
	return Float32(fma(t_0, Float32(t_0 / Float32(Float32(pi) * s)), Float32(exp(Float32(Float32(-r) / s)) / Float32(Float32(pi) * s))) * Float32(Float32(0.125) / r))
end

Derivation

1. Initial program 99.6%

   \[\frac{0.25 \cdot e^{\frac{-r}{s}}}{\left(\left(2 \cdot \pi\right) \cdot s\right) \cdot r} + \frac{0.75 \cdot e^{\frac{-r}{3 \cdot s}}}{\left(\left(6 \cdot \pi\right) \cdot s\right) \cdot r} \]
2. Add Preprocessing

3. Taylor expanded in r around inf

   \[\leadsto \color{blue}{\frac{\frac{1}{8} \cdot \frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{1}{8} \cdot \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}}{r}} \]
4. Step-by-step derivation

   1. distribute-lft-out N/A

      \[\leadsto \frac{\color{blue}{\frac{1}{8} \cdot \left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right)}}{r} \]

   2. *-commutative N/A

      \[\leadsto \frac{\color{blue}{\left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right) \cdot \frac{1}{8}}}{r} \]

   3. associate-/l* N/A

      \[\leadsto \color{blue}{\left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right) \cdot \frac{\frac{1}{8}}{r}} \]

   4. lower-*.f32 N/A

      \[\leadsto \color{blue}{\left(\frac{e^{-1 \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)} + \frac{e^{\frac{-1}{3} \cdot \frac{r}{s}}}{s \cdot \mathsf{PI}\left(\right)}\right) \cdot \frac{\frac{1}{8}}{r}} \]

5. Applied rewrites 99.5%

   \[\leadsto \color{blue}{\left(\frac{e^{-0.3333333333333333 \cdot \frac{r}{s}}}{\pi \cdot s} + \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{0.125}{r}} \]
6. Step-by-step derivation

7. Applied rewrites 99.5%

   \[\leadsto \left(\frac{{e}^{\left(\frac{-0.3333333333333333}{s} \cdot r\right)}}{\pi \cdot s} + \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{0.125}{r} \]
8. Step-by-step derivation

9. Applied rewrites 99.6%

   \[\leadsto \mathsf{fma}\left({e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}, \frac{{e}^{\left(\frac{r}{s} \cdot -0.16666666666666666\right)}}{\pi \cdot s}, \frac{e^{\frac{-r}{s}}}{\pi \cdot s}\right) \cdot \frac{\color{blue}{0.125}}{r} \]
10. Add Preprocessing

Reproduce

herbie shell --seed 2024230 
(FPCore (s r)
  :name "Disney BSSRDF, PDF of scattering profile"
  :precision binary32
  :pre (and (and (<= 0.0 s) (<= s 256.0)) (and (< 1e-6 r) (< r 1000000.0)))
  (+ (/ (* 0.25 (exp (/ (- r) s))) (* (* (* 2.0 PI) s) r)) (/ (* 0.75 (exp (/ (- r) (* 3.0 s)))) (* (* (* 6.0 PI) s) r))))