Data.Array.Repa.Algorithms.Pixel:doubleRmsOfRGB8 from repa-algorithms-3.4.0.1

Average Error: 38.6 → 27.1

Time: 1.9min

Precision: binary64

Cost: 7496

Math

\[\sqrt{\frac{\left(x \cdot x + y \cdot y\right) + z \cdot z}{3}} \]

↓

\[\begin{array}{l} \mathbf{if}\;z \leq -2.1 \cdot 10^{+43}:\\ \;\;\;\;-\frac{z}{\sqrt{3}}\\ \mathbf{elif}\;z \leq 3.15 \cdot 10^{+81}:\\ \;\;\;\;\sqrt{\frac{\left(x \cdot x + y \cdot y\right) + z \cdot z}{3}}\\ \mathbf{else}:\\ \;\;\;\;z \cdot \sqrt{0.3333333333333333}\\ \end{array} \]

FPCore

(FPCore (x y z)
 :precision binary64
 (sqrt (/ (+ (+ (* x x) (* y y)) (* z z)) 3.0)))

↓

(FPCore (x y z)
 :precision binary64
 (if (<= z -2.1e+43)
   (- (/ z (sqrt 3.0)))
   (if (<= z 3.15e+81)
     (sqrt (/ (+ (+ (* x x) (* y y)) (* z z)) 3.0))
     (* z (sqrt 0.3333333333333333)))))

C

double code(double x, double y, double z) {
	return sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0));
}

↓

double code(double x, double y, double z) {
	double tmp;
	if (z <= -2.1e+43) {
		tmp = -(z / sqrt(3.0));
	} else if (z <= 3.15e+81) {
		tmp = sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0));
	} else {
		tmp = z * sqrt(0.3333333333333333);
	}
	return tmp;
}

Fortran

real(8) function code(x, y, z)
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    real(8), intent (in) :: z
    code = sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0d0))
end function

↓

real(8) function code(x, y, z)
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    real(8), intent (in) :: z
    real(8) :: tmp
    if (z <= (-2.1d+43)) then
        tmp = -(z / sqrt(3.0d0))
    else if (z <= 3.15d+81) then
        tmp = sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0d0))
    else
        tmp = z * sqrt(0.3333333333333333d0)
    end if
    code = tmp
end function

Java

public static double code(double x, double y, double z) {
	return Math.sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0));
}

↓

public static double code(double x, double y, double z) {
	double tmp;
	if (z <= -2.1e+43) {
		tmp = -(z / Math.sqrt(3.0));
	} else if (z <= 3.15e+81) {
		tmp = Math.sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0));
	} else {
		tmp = z * Math.sqrt(0.3333333333333333);
	}
	return tmp;
}

Python

def code(x, y, z):
	return math.sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0))

↓

def code(x, y, z):
	tmp = 0
	if z <= -2.1e+43:
		tmp = -(z / math.sqrt(3.0))
	elif z <= 3.15e+81:
		tmp = math.sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0))
	else:
		tmp = z * math.sqrt(0.3333333333333333)
	return tmp

Julia

function code(x, y, z)
	return sqrt(Float64(Float64(Float64(Float64(x * x) + Float64(y * y)) + Float64(z * z)) / 3.0))
end

↓

function code(x, y, z)
	tmp = 0.0
	if (z <= -2.1e+43)
		tmp = Float64(-Float64(z / sqrt(3.0)));
	elseif (z <= 3.15e+81)
		tmp = sqrt(Float64(Float64(Float64(Float64(x * x) + Float64(y * y)) + Float64(z * z)) / 3.0));
	else
		tmp = Float64(z * sqrt(0.3333333333333333));
	end
	return tmp
end

MATLAB

function tmp = code(x, y, z)
	tmp = sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0));
end

↓

function tmp_2 = code(x, y, z)
	tmp = 0.0;
	if (z <= -2.1e+43)
		tmp = -(z / sqrt(3.0));
	elseif (z <= 3.15e+81)
		tmp = sqrt(((((x * x) + (y * y)) + (z * z)) / 3.0));
	else
		tmp = z * sqrt(0.3333333333333333);
	end
	tmp_2 = tmp;
end

Wolfram

code[x_, y_, z_] := N[Sqrt[N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y * y), $MachinePrecision]), $MachinePrecision] + N[(z * z), $MachinePrecision]), $MachinePrecision] / 3.0), $MachinePrecision]], $MachinePrecision]

↓

code[x_, y_, z_] := If[LessEqual[z, -2.1e+43], (-N[(z / N[Sqrt[3.0], $MachinePrecision]), $MachinePrecision]), If[LessEqual[z, 3.15e+81], N[Sqrt[N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y * y), $MachinePrecision]), $MachinePrecision] + N[(z * z), $MachinePrecision]), $MachinePrecision] / 3.0), $MachinePrecision]], $MachinePrecision], N[(z * N[Sqrt[0.3333333333333333], $MachinePrecision]), $MachinePrecision]]]

Target

Original	38.6
Target	26.4
Herbie	27.1

\[\begin{array}{l} \mathbf{if}\;z < -6.396479394109776 \cdot 10^{+136}:\\ \;\;\;\;\frac{-z}{\sqrt{3}}\\ \mathbf{elif}\;z < 7.320293694404182 \cdot 10^{+117}:\\ \;\;\;\;\frac{\sqrt{\left(z \cdot z + x \cdot x\right) + y \cdot y}}{\sqrt{3}}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{0.3333333333333333} \cdot z\\ \end{array} \]

Derivation

1. Split input into 3 regimes

2. if z < -2.10000000000000002e43

   1. Initial program 49.1

      \[\sqrt{\frac{\left(x \cdot x + y \cdot y\right) + z \cdot z}{3}} \]

   2. Simplified 49.1

      \[\leadsto \color{blue}{\sqrt{\mathsf{fma}\left(z, z, \mathsf{fma}\left(y, y, {x}^{2}\right)\right) \cdot 0.3333333333333333}} \]
      Proof

   3. Taylor expanded in z around -inf 22.7

      \[\leadsto \color{blue}{-1 \cdot \left(z \cdot \sqrt{0.3333333333333333}\right)} \]

   4. Simplified 22.7

      \[\leadsto \color{blue}{-z \cdot \sqrt{0.3333333333333333}} \]
      Proof

   5. Applied egg-rr 22.7

      \[\leadsto -\color{blue}{z \cdot \sqrt{0.3333333333333333}} \]

   6. Simplified 22.7

      \[\leadsto -\color{blue}{\frac{z}{\sqrt{3}}} \]
      Proof

   if -2.10000000000000002e43 < z < 3.1500000000000002e81

   1. Initial program 30.8

      \[\sqrt{\frac{\left(x \cdot x + y \cdot y\right) + z \cdot z}{3}} \]

   if 3.1500000000000002e81 < z

   1. Initial program 51.5

      \[\sqrt{\frac{\left(x \cdot x + y \cdot y\right) + z \cdot z}{3}} \]

   2. Simplified 51.6

      \[\leadsto \color{blue}{\sqrt{\mathsf{fma}\left(z, z, \mathsf{fma}\left(y, y, {x}^{2}\right)\right) \cdot 0.3333333333333333}} \]
      Proof

   3. Taylor expanded in z around inf 20.2

      \[\leadsto \color{blue}{z \cdot \sqrt{0.3333333333333333}} \]
3. Recombined 3 regimes into one program.

Alternatives

Alternative 1
Error	36.2
Cost	7900

\[\begin{array}{l} t_0 := \sqrt{0.3333333333333333} \cdot x\\ t_1 := -t_0\\ \mathbf{if}\;z \leq -8 \cdot 10^{-16}:\\ \;\;\;\;-\frac{z}{\sqrt{3}}\\ \mathbf{elif}\;z \leq -7.5 \cdot 10^{-204}:\\ \;\;\;\;y \cdot \sqrt{0.3333333333333333}\\ \mathbf{elif}\;z \leq -1 \cdot 10^{-252}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;z \leq -1.18 \cdot 10^{-288}:\\ \;\;\;\;\frac{y}{\sqrt{3}}\\ \mathbf{elif}\;z \leq 3.1 \cdot 10^{-161}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;z \leq 2.4 \cdot 10^{-107}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;z \leq 2.6 \cdot 10^{+29}:\\ \;\;\;\;\sqrt{\frac{x \cdot x + z \cdot z}{3}}\\ \mathbf{else}:\\ \;\;\;\;z \cdot \sqrt{0.3333333333333333}\\ \end{array} \]

Alternative 2
Error	36.9
Cost	7448

\[\begin{array}{l} t_0 := \sqrt{0.3333333333333333} \cdot x\\ t_1 := -t_0\\ t_2 := z \cdot \sqrt{0.3333333333333333}\\ \mathbf{if}\;z \leq -2.9 \cdot 10^{-16}:\\ \;\;\;\;-t_2\\ \mathbf{elif}\;z \leq -4.4 \cdot 10^{-205}:\\ \;\;\;\;y \cdot \sqrt{0.3333333333333333}\\ \mathbf{elif}\;z \leq -2.6 \cdot 10^{-249}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;z \leq -5.9 \cdot 10^{-289}:\\ \;\;\;\;\frac{y}{\sqrt{3}}\\ \mathbf{elif}\;z \leq 8.2 \cdot 10^{-167}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;z \leq 6.1 \cdot 10^{-33}:\\ \;\;\;\;t_1\\ \mathbf{else}:\\ \;\;\;\;t_2\\ \end{array} \]

Alternative 3
Error	36.9
Cost	7448

\[\begin{array}{l} t_0 := \sqrt{0.3333333333333333} \cdot x\\ t_1 := -t_0\\ \mathbf{if}\;z \leq -4.2 \cdot 10^{-16}:\\ \;\;\;\;-\frac{z}{\sqrt{3}}\\ \mathbf{elif}\;z \leq -4.2 \cdot 10^{-205}:\\ \;\;\;\;y \cdot \sqrt{0.3333333333333333}\\ \mathbf{elif}\;z \leq -3.3 \cdot 10^{-252}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;z \leq -1.36 \cdot 10^{-288}:\\ \;\;\;\;\frac{y}{\sqrt{3}}\\ \mathbf{elif}\;z \leq 5.2 \cdot 10^{-164}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;z \leq 5.5 \cdot 10^{-32}:\\ \;\;\;\;t_1\\ \mathbf{else}:\\ \;\;\;\;z \cdot \sqrt{0.3333333333333333}\\ \end{array} \]

Alternative 4
Error	44.5
Cost	6988

\[\begin{array}{l} t_0 := \sqrt{0.3333333333333333} \cdot x\\ \mathbf{if}\;z \leq -4 \cdot 10^{+73}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;z \leq -1.32 \cdot 10^{-288}:\\ \;\;\;\;y \cdot \sqrt{0.3333333333333333}\\ \mathbf{elif}\;z \leq 1.08 \cdot 10^{-50}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;z \cdot \sqrt{0.3333333333333333}\\ \end{array} \]

Alternative 5
Error	44.5
Cost	6988

\[\begin{array}{l} \mathbf{if}\;z \leq -1.2 \cdot 10^{+73}:\\ \;\;\;\;\frac{x}{\sqrt{3}}\\ \mathbf{elif}\;z \leq -1.7 \cdot 10^{-288}:\\ \;\;\;\;y \cdot \sqrt{0.3333333333333333}\\ \mathbf{elif}\;z \leq 1.25 \cdot 10^{-50}:\\ \;\;\;\;\sqrt{0.3333333333333333} \cdot x\\ \mathbf{else}:\\ \;\;\;\;z \cdot \sqrt{0.3333333333333333}\\ \end{array} \]

Alternative 6
Error	36.7
Cost	6988

\[\begin{array}{l} t_0 := z \cdot \sqrt{0.3333333333333333}\\ \mathbf{if}\;z \leq -3.6 \cdot 10^{-16}:\\ \;\;\;\;-t_0\\ \mathbf{elif}\;z \leq -1.35 \cdot 10^{-288}:\\ \;\;\;\;y \cdot \sqrt{0.3333333333333333}\\ \mathbf{elif}\;z \leq 9.3 \cdot 10^{-51}:\\ \;\;\;\;\sqrt{0.3333333333333333} \cdot x\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

Alternative 7
Error	37.7
Cost	6856

\[\begin{array}{l} t_0 := y \cdot \sqrt{0.3333333333333333}\\ \mathbf{if}\;y \leq -4.8 \cdot 10^{+89}:\\ \;\;\;\;-t_0\\ \mathbf{elif}\;y \leq 4.5 \cdot 10^{-40}:\\ \;\;\;\;z \cdot \sqrt{0.3333333333333333}\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

Alternative 8
Error	44.9
Cost	6724

\[\begin{array}{l} \mathbf{if}\;z \leq 10500:\\ \;\;\;\;y \cdot \sqrt{0.3333333333333333}\\ \mathbf{else}:\\ \;\;\;\;z \cdot \sqrt{0.3333333333333333}\\ \end{array} \]

Alternative 9
Error	52.4
Cost	6592

\[y \cdot \sqrt{0.3333333333333333} \]

Reproduce

herbie shell --seed 2023033 
(FPCore (x y z)
  :name "Data.Array.Repa.Algorithms.Pixel:doubleRmsOfRGB8 from repa-algorithms-3.4.0.1"
  :precision binary64

  :herbie-target
  (if (< z -6.396479394109776e+136) (/ (- z) (sqrt 3.0)) (if (< z 7.320293694404182e+117) (/ (sqrt (+ (+ (* z z) (* x x)) (* y y))) (sqrt 3.0)) (* (sqrt 0.3333333333333333) z)))

  (sqrt (/ (+ (+ (* x x) (* y y)) (* z z)) 3.0)))