symmetry log of sum of exp

Percentage Accurate: 53.3% → 98.4%
Time: 7.7s
Alternatives: 14
Speedup: 1.5×

Specification

?
\[\begin{array}{l} \\ \log \left(e^{a} + e^{b}\right) \end{array} \]
(FPCore (a b) :precision binary64 (log (+ (exp a) (exp b))))
double code(double a, double b) {
	return log((exp(a) + exp(b)));
}
real(8) function code(a, b)
    real(8), intent (in) :: a
    real(8), intent (in) :: b
    code = log((exp(a) + exp(b)))
end function
public static double code(double a, double b) {
	return Math.log((Math.exp(a) + Math.exp(b)));
}
def code(a, b):
	return math.log((math.exp(a) + math.exp(b)))
function code(a, b)
	return log(Float64(exp(a) + exp(b)))
end
function tmp = code(a, b)
	tmp = log((exp(a) + exp(b)));
end
code[a_, b_] := N[Log[N[(N[Exp[a], $MachinePrecision] + N[Exp[b], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]
\begin{array}{l}

\\
\log \left(e^{a} + e^{b}\right)
\end{array}

Sampling outcomes in binary64 precision:

Local Percentage Accuracy vs ?

The average percentage accuracy by input value. Horizontal axis shows value of an input variable; the variable is choosen in the title. Vertical axis is accuracy; higher is better. Red represent the original program, while blue represents Herbie's suggestion. These can be toggled with buttons below the plot. The line is an average while dots represent individual samples.

Accuracy vs Speed?

Herbie found 14 alternatives:

AlternativeAccuracySpeedup
The accuracy (vertical axis) and speed (horizontal axis) of each alternatives. Up and to the right is better. The red square shows the initial program, and each blue circle shows an alternative.The line shows the best available speed-accuracy tradeoffs.

Initial Program: 53.3% accurate, 1.0× speedup?

\[\begin{array}{l} \\ \log \left(e^{a} + e^{b}\right) \end{array} \]
(FPCore (a b) :precision binary64 (log (+ (exp a) (exp b))))
double code(double a, double b) {
	return log((exp(a) + exp(b)));
}
real(8) function code(a, b)
    real(8), intent (in) :: a
    real(8), intent (in) :: b
    code = log((exp(a) + exp(b)))
end function
public static double code(double a, double b) {
	return Math.log((Math.exp(a) + Math.exp(b)));
}
def code(a, b):
	return math.log((math.exp(a) + math.exp(b)))
function code(a, b)
	return log(Float64(exp(a) + exp(b)))
end
function tmp = code(a, b)
	tmp = log((exp(a) + exp(b)));
end
code[a_, b_] := N[Log[N[(N[Exp[a], $MachinePrecision] + N[Exp[b], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]
\begin{array}{l}

\\
\log \left(e^{a} + e^{b}\right)
\end{array}

Alternative 1: 98.4% accurate, 1.0× speedup?

\[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \mathsf{log1p}\left(e^{a}\right) + \frac{b}{1 + e^{a}} \end{array} \]
NOTE: a and b should be sorted in increasing order before calling this function.
(FPCore (a b) :precision binary64 (+ (log1p (exp a)) (/ b (+ 1.0 (exp a)))))
assert(a < b);
double code(double a, double b) {
	return log1p(exp(a)) + (b / (1.0 + exp(a)));
}
assert a < b;
public static double code(double a, double b) {
	return Math.log1p(Math.exp(a)) + (b / (1.0 + Math.exp(a)));
}
[a, b] = sort([a, b])
def code(a, b):
	return math.log1p(math.exp(a)) + (b / (1.0 + math.exp(a)))
a, b = sort([a, b])
function code(a, b)
	return Float64(log1p(exp(a)) + Float64(b / Float64(1.0 + exp(a))))
end
NOTE: a and b should be sorted in increasing order before calling this function.
code[a_, b_] := N[(N[Log[1 + N[Exp[a], $MachinePrecision]], $MachinePrecision] + N[(b / N[(1.0 + N[Exp[a], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\mathsf{log1p}\left(e^{a}\right) + \frac{b}{1 + e^{a}}
\end{array}
Derivation
  1. Initial program 52.2%

    \[\log \left(e^{a} + e^{b}\right) \]
  2. Add Preprocessing
  3. Taylor expanded in b around 0

    \[\leadsto \color{blue}{\log \left(1 + e^{a}\right) + \frac{b}{1 + e^{a}}} \]
  4. Step-by-step derivation
    1. +-commutativeN/A

      \[\leadsto \color{blue}{\frac{b}{1 + e^{a}} + \log \left(1 + e^{a}\right)} \]
    2. *-rgt-identityN/A

      \[\leadsto \frac{\color{blue}{b \cdot 1}}{1 + e^{a}} + \log \left(1 + e^{a}\right) \]
    3. associate-*r/N/A

      \[\leadsto \color{blue}{b \cdot \frac{1}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
    4. lower-+.f64N/A

      \[\leadsto \color{blue}{b \cdot \frac{1}{1 + e^{a}} + \log \left(1 + e^{a}\right)} \]
    5. associate-*r/N/A

      \[\leadsto \color{blue}{\frac{b \cdot 1}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
    6. *-rgt-identityN/A

      \[\leadsto \frac{\color{blue}{b}}{1 + e^{a}} + \log \left(1 + e^{a}\right) \]
    7. lower-/.f64N/A

      \[\leadsto \color{blue}{\frac{b}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
    8. +-commutativeN/A

      \[\leadsto \frac{b}{\color{blue}{e^{a} + 1}} + \log \left(1 + e^{a}\right) \]
    9. lower-+.f64N/A

      \[\leadsto \frac{b}{\color{blue}{e^{a} + 1}} + \log \left(1 + e^{a}\right) \]
    10. lower-exp.f64N/A

      \[\leadsto \frac{b}{\color{blue}{e^{a}} + 1} + \log \left(1 + e^{a}\right) \]
    11. lower-log1p.f64N/A

      \[\leadsto \frac{b}{e^{a} + 1} + \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
    12. lower-exp.f6471.7

      \[\leadsto \frac{b}{e^{a} + 1} + \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
  5. Applied rewrites71.7%

    \[\leadsto \color{blue}{\frac{b}{e^{a} + 1} + \mathsf{log1p}\left(e^{a}\right)} \]
  6. Final simplification71.7%

    \[\leadsto \mathsf{log1p}\left(e^{a}\right) + \frac{b}{1 + e^{a}} \]
  7. Add Preprocessing

Alternative 2: 57.3% accurate, 0.7× speedup?

\[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{b} + e^{a} \leq 1:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\log \left(\left(1 + b\right) + e^{a}\right)\\ \end{array} \end{array} \]
NOTE: a and b should be sorted in increasing order before calling this function.
(FPCore (a b)
 :precision binary64
 (if (<= (+ (exp b) (exp a)) 1.0)
   (* (fma 0.125 b 0.5) b)
   (log (+ (+ 1.0 b) (exp a)))))
assert(a < b);
double code(double a, double b) {
	double tmp;
	if ((exp(b) + exp(a)) <= 1.0) {
		tmp = fma(0.125, b, 0.5) * b;
	} else {
		tmp = log(((1.0 + b) + exp(a)));
	}
	return tmp;
}
a, b = sort([a, b])
function code(a, b)
	tmp = 0.0
	if (Float64(exp(b) + exp(a)) <= 1.0)
		tmp = Float64(fma(0.125, b, 0.5) * b);
	else
		tmp = log(Float64(Float64(1.0 + b) + exp(a)));
	end
	return tmp
end
NOTE: a and b should be sorted in increasing order before calling this function.
code[a_, b_] := If[LessEqual[N[(N[Exp[b], $MachinePrecision] + N[Exp[a], $MachinePrecision]), $MachinePrecision], 1.0], N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision], N[Log[N[(N[(1.0 + b), $MachinePrecision] + N[Exp[a], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]
\begin{array}{l}
[a, b] = \mathsf{sort}([a, b])\\
\\
\begin{array}{l}
\mathbf{if}\;e^{b} + e^{a} \leq 1:\\
\;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\

\mathbf{else}:\\
\;\;\;\;\log \left(\left(1 + b\right) + e^{a}\right)\\


\end{array}
\end{array}
Derivation
  1. Split input into 2 regimes
  2. if (+.f64 (exp.f64 a) (exp.f64 b)) < 1

    1. Initial program 7.3%

      \[\log \left(e^{a} + e^{b}\right) \]
    2. Add Preprocessing
    3. Taylor expanded in a around 0

      \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
    4. Step-by-step derivation
      1. lower-log1p.f64N/A

        \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
      2. lower-exp.f645.6

        \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
    5. Applied rewrites5.6%

      \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
    6. Taylor expanded in b around 0

      \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
    7. Step-by-step derivation
      1. Applied rewrites2.7%

        \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
      2. Taylor expanded in b around inf

        \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
      3. Step-by-step derivation
        1. Applied rewrites9.8%

          \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]

        if 1 < (+.f64 (exp.f64 a) (exp.f64 b))

        1. Initial program 96.4%

          \[\log \left(e^{a} + e^{b}\right) \]
        2. Add Preprocessing
        3. Taylor expanded in b around 0

          \[\leadsto \log \left(e^{a} + \color{blue}{\left(1 + b\right)}\right) \]
        4. Step-by-step derivation
          1. lower-+.f6493.6

            \[\leadsto \log \left(e^{a} + \color{blue}{\left(1 + b\right)}\right) \]
        5. Applied rewrites93.6%

          \[\leadsto \log \left(e^{a} + \color{blue}{\left(1 + b\right)}\right) \]
      4. Recombined 2 regimes into one program.
      5. Final simplification52.0%

        \[\leadsto \begin{array}{l} \mathbf{if}\;e^{b} + e^{a} \leq 1:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\log \left(\left(1 + b\right) + e^{a}\right)\\ \end{array} \]
      6. Add Preprocessing

      Alternative 3: 58.8% accurate, 0.7× speedup?

      \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{b} \leq 1.0000000000000004:\\ \;\;\;\;0.5 \cdot b + \mathsf{log1p}\left(e^{a}\right)\\ \mathbf{else}:\\ \;\;\;\;\log \left(e^{b} + e^{a}\right)\\ \end{array} \end{array} \]
      NOTE: a and b should be sorted in increasing order before calling this function.
      (FPCore (a b)
       :precision binary64
       (if (<= (exp b) 1.0000000000000004)
         (+ (* 0.5 b) (log1p (exp a)))
         (log (+ (exp b) (exp a)))))
      assert(a < b);
      double code(double a, double b) {
      	double tmp;
      	if (exp(b) <= 1.0000000000000004) {
      		tmp = (0.5 * b) + log1p(exp(a));
      	} else {
      		tmp = log((exp(b) + exp(a)));
      	}
      	return tmp;
      }
      
      assert a < b;
      public static double code(double a, double b) {
      	double tmp;
      	if (Math.exp(b) <= 1.0000000000000004) {
      		tmp = (0.5 * b) + Math.log1p(Math.exp(a));
      	} else {
      		tmp = Math.log((Math.exp(b) + Math.exp(a)));
      	}
      	return tmp;
      }
      
      [a, b] = sort([a, b])
      def code(a, b):
      	tmp = 0
      	if math.exp(b) <= 1.0000000000000004:
      		tmp = (0.5 * b) + math.log1p(math.exp(a))
      	else:
      		tmp = math.log((math.exp(b) + math.exp(a)))
      	return tmp
      
      a, b = sort([a, b])
      function code(a, b)
      	tmp = 0.0
      	if (exp(b) <= 1.0000000000000004)
      		tmp = Float64(Float64(0.5 * b) + log1p(exp(a)));
      	else
      		tmp = log(Float64(exp(b) + exp(a)));
      	end
      	return tmp
      end
      
      NOTE: a and b should be sorted in increasing order before calling this function.
      code[a_, b_] := If[LessEqual[N[Exp[b], $MachinePrecision], 1.0000000000000004], N[(N[(0.5 * b), $MachinePrecision] + N[Log[1 + N[Exp[a], $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Log[N[(N[Exp[b], $MachinePrecision] + N[Exp[a], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]
      
      \begin{array}{l}
      [a, b] = \mathsf{sort}([a, b])\\
      \\
      \begin{array}{l}
      \mathbf{if}\;e^{b} \leq 1.0000000000000004:\\
      \;\;\;\;0.5 \cdot b + \mathsf{log1p}\left(e^{a}\right)\\
      
      \mathbf{else}:\\
      \;\;\;\;\log \left(e^{b} + e^{a}\right)\\
      
      
      \end{array}
      \end{array}
      
      Derivation
      1. Split input into 2 regimes
      2. if (exp.f64 b) < 1.0000000000000004

        1. Initial program 52.2%

          \[\log \left(e^{a} + e^{b}\right) \]
        2. Add Preprocessing
        3. Taylor expanded in b around 0

          \[\leadsto \color{blue}{\log \left(1 + e^{a}\right) + \frac{b}{1 + e^{a}}} \]
        4. Step-by-step derivation
          1. +-commutativeN/A

            \[\leadsto \color{blue}{\frac{b}{1 + e^{a}} + \log \left(1 + e^{a}\right)} \]
          2. *-rgt-identityN/A

            \[\leadsto \frac{\color{blue}{b \cdot 1}}{1 + e^{a}} + \log \left(1 + e^{a}\right) \]
          3. associate-*r/N/A

            \[\leadsto \color{blue}{b \cdot \frac{1}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
          4. lower-+.f64N/A

            \[\leadsto \color{blue}{b \cdot \frac{1}{1 + e^{a}} + \log \left(1 + e^{a}\right)} \]
          5. associate-*r/N/A

            \[\leadsto \color{blue}{\frac{b \cdot 1}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
          6. *-rgt-identityN/A

            \[\leadsto \frac{\color{blue}{b}}{1 + e^{a}} + \log \left(1 + e^{a}\right) \]
          7. lower-/.f64N/A

            \[\leadsto \color{blue}{\frac{b}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
          8. +-commutativeN/A

            \[\leadsto \frac{b}{\color{blue}{e^{a} + 1}} + \log \left(1 + e^{a}\right) \]
          9. lower-+.f64N/A

            \[\leadsto \frac{b}{\color{blue}{e^{a} + 1}} + \log \left(1 + e^{a}\right) \]
          10. lower-exp.f64N/A

            \[\leadsto \frac{b}{\color{blue}{e^{a}} + 1} + \log \left(1 + e^{a}\right) \]
          11. lower-log1p.f64N/A

            \[\leadsto \frac{b}{e^{a} + 1} + \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
          12. lower-exp.f6472.0

            \[\leadsto \frac{b}{e^{a} + 1} + \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
        5. Applied rewrites72.0%

          \[\leadsto \color{blue}{\frac{b}{e^{a} + 1} + \mathsf{log1p}\left(e^{a}\right)} \]
        6. Taylor expanded in a around 0

          \[\leadsto \frac{1}{2} \cdot b + \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
        7. Step-by-step derivation
          1. Applied rewrites52.8%

            \[\leadsto 0.5 \cdot b + \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]

          if 1.0000000000000004 < (exp.f64 b)

          1. Initial program 51.7%

            \[\log \left(e^{a} + e^{b}\right) \]
          2. Add Preprocessing
        8. Recombined 2 regimes into one program.
        9. Final simplification52.8%

          \[\leadsto \begin{array}{l} \mathbf{if}\;e^{b} \leq 1.0000000000000004:\\ \;\;\;\;0.5 \cdot b + \mathsf{log1p}\left(e^{a}\right)\\ \mathbf{else}:\\ \;\;\;\;\log \left(e^{b} + e^{a}\right)\\ \end{array} \]
        10. Add Preprocessing

        Alternative 4: 56.7% accurate, 1.0× speedup?

        \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{a} \leq 0:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\mathsf{log1p}\left(e^{a}\right)\\ \end{array} \end{array} \]
        NOTE: a and b should be sorted in increasing order before calling this function.
        (FPCore (a b)
         :precision binary64
         (if (<= (exp a) 0.0) (* (fma 0.125 b 0.5) b) (log1p (exp a))))
        assert(a < b);
        double code(double a, double b) {
        	double tmp;
        	if (exp(a) <= 0.0) {
        		tmp = fma(0.125, b, 0.5) * b;
        	} else {
        		tmp = log1p(exp(a));
        	}
        	return tmp;
        }
        
        a, b = sort([a, b])
        function code(a, b)
        	tmp = 0.0
        	if (exp(a) <= 0.0)
        		tmp = Float64(fma(0.125, b, 0.5) * b);
        	else
        		tmp = log1p(exp(a));
        	end
        	return tmp
        end
        
        NOTE: a and b should be sorted in increasing order before calling this function.
        code[a_, b_] := If[LessEqual[N[Exp[a], $MachinePrecision], 0.0], N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision], N[Log[1 + N[Exp[a], $MachinePrecision]], $MachinePrecision]]
        
        \begin{array}{l}
        [a, b] = \mathsf{sort}([a, b])\\
        \\
        \begin{array}{l}
        \mathbf{if}\;e^{a} \leq 0:\\
        \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\
        
        \mathbf{else}:\\
        \;\;\;\;\mathsf{log1p}\left(e^{a}\right)\\
        
        
        \end{array}
        \end{array}
        
        Derivation
        1. Split input into 2 regimes
        2. if (exp.f64 a) < 0.0

          1. Initial program 8.8%

            \[\log \left(e^{a} + e^{b}\right) \]
          2. Add Preprocessing
          3. Taylor expanded in a around 0

            \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
          4. Step-by-step derivation
            1. lower-log1p.f64N/A

              \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
            2. lower-exp.f643.8

              \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
          5. Applied rewrites3.8%

            \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
          6. Taylor expanded in b around 0

            \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
          7. Step-by-step derivation
            1. Applied rewrites3.9%

              \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
            2. Taylor expanded in b around inf

              \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
            3. Step-by-step derivation
              1. Applied rewrites18.4%

                \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]

              if 0.0 < (exp.f64 a)

              1. Initial program 66.3%

                \[\log \left(e^{a} + e^{b}\right) \]
              2. Add Preprocessing
              3. Taylor expanded in b around 0

                \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
              4. Step-by-step derivation
                1. lower-log1p.f64N/A

                  \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                2. lower-exp.f6463.3

                  \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
              5. Applied rewrites63.3%

                \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
            4. Recombined 2 regimes into one program.
            5. Add Preprocessing

            Alternative 5: 56.2% accurate, 1.3× speedup?

            \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{a} \leq 2 \cdot 10^{-39}:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(a \cdot a, -0.005208333333333333, 0.125\right), a, 0.5\right), a, \log 2\right)\\ \end{array} \end{array} \]
            NOTE: a and b should be sorted in increasing order before calling this function.
            (FPCore (a b)
             :precision binary64
             (if (<= (exp a) 2e-39)
               (* (fma 0.125 b 0.5) b)
               (fma (fma (fma (* a a) -0.005208333333333333 0.125) a 0.5) a (log 2.0))))
            assert(a < b);
            double code(double a, double b) {
            	double tmp;
            	if (exp(a) <= 2e-39) {
            		tmp = fma(0.125, b, 0.5) * b;
            	} else {
            		tmp = fma(fma(fma((a * a), -0.005208333333333333, 0.125), a, 0.5), a, log(2.0));
            	}
            	return tmp;
            }
            
            a, b = sort([a, b])
            function code(a, b)
            	tmp = 0.0
            	if (exp(a) <= 2e-39)
            		tmp = Float64(fma(0.125, b, 0.5) * b);
            	else
            		tmp = fma(fma(fma(Float64(a * a), -0.005208333333333333, 0.125), a, 0.5), a, log(2.0));
            	end
            	return tmp
            end
            
            NOTE: a and b should be sorted in increasing order before calling this function.
            code[a_, b_] := If[LessEqual[N[Exp[a], $MachinePrecision], 2e-39], N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision], N[(N[(N[(N[(a * a), $MachinePrecision] * -0.005208333333333333 + 0.125), $MachinePrecision] * a + 0.5), $MachinePrecision] * a + N[Log[2.0], $MachinePrecision]), $MachinePrecision]]
            
            \begin{array}{l}
            [a, b] = \mathsf{sort}([a, b])\\
            \\
            \begin{array}{l}
            \mathbf{if}\;e^{a} \leq 2 \cdot 10^{-39}:\\
            \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\
            
            \mathbf{else}:\\
            \;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(a \cdot a, -0.005208333333333333, 0.125\right), a, 0.5\right), a, \log 2\right)\\
            
            
            \end{array}
            \end{array}
            
            Derivation
            1. Split input into 2 regimes
            2. if (exp.f64 a) < 1.99999999999999986e-39

              1. Initial program 8.8%

                \[\log \left(e^{a} + e^{b}\right) \]
              2. Add Preprocessing
              3. Taylor expanded in a around 0

                \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
              4. Step-by-step derivation
                1. lower-log1p.f64N/A

                  \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                2. lower-exp.f643.8

                  \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
              5. Applied rewrites3.8%

                \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
              6. Taylor expanded in b around 0

                \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
              7. Step-by-step derivation
                1. Applied rewrites3.9%

                  \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
                2. Taylor expanded in b around inf

                  \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
                3. Step-by-step derivation
                  1. Applied rewrites18.4%

                    \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]

                  if 1.99999999999999986e-39 < (exp.f64 a)

                  1. Initial program 66.3%

                    \[\log \left(e^{a} + e^{b}\right) \]
                  2. Add Preprocessing
                  3. Taylor expanded in b around 0

                    \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
                  4. Step-by-step derivation
                    1. lower-log1p.f64N/A

                      \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                    2. lower-exp.f6463.3

                      \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                  5. Applied rewrites63.3%

                    \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                  6. Taylor expanded in a around 0

                    \[\leadsto \log 2 + \color{blue}{a \cdot \left(\frac{1}{2} + a \cdot \left(\frac{1}{8} + \frac{-1}{192} \cdot {a}^{2}\right)\right)} \]
                  7. Step-by-step derivation
                    1. Applied rewrites62.7%

                      \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(a \cdot a, -0.005208333333333333, 0.125\right), a, 0.5\right), \color{blue}{a}, \log 2\right) \]
                  8. Recombined 2 regimes into one program.
                  9. Add Preprocessing

                  Alternative 6: 56.1% accurate, 1.4× speedup?

                  \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{a} \leq 0:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(0.125, a, 0.5\right), a, \log 2\right)\\ \end{array} \end{array} \]
                  NOTE: a and b should be sorted in increasing order before calling this function.
                  (FPCore (a b)
                   :precision binary64
                   (if (<= (exp a) 0.0)
                     (* (fma 0.125 b 0.5) b)
                     (fma (fma 0.125 a 0.5) a (log 2.0))))
                  assert(a < b);
                  double code(double a, double b) {
                  	double tmp;
                  	if (exp(a) <= 0.0) {
                  		tmp = fma(0.125, b, 0.5) * b;
                  	} else {
                  		tmp = fma(fma(0.125, a, 0.5), a, log(2.0));
                  	}
                  	return tmp;
                  }
                  
                  a, b = sort([a, b])
                  function code(a, b)
                  	tmp = 0.0
                  	if (exp(a) <= 0.0)
                  		tmp = Float64(fma(0.125, b, 0.5) * b);
                  	else
                  		tmp = fma(fma(0.125, a, 0.5), a, log(2.0));
                  	end
                  	return tmp
                  end
                  
                  NOTE: a and b should be sorted in increasing order before calling this function.
                  code[a_, b_] := If[LessEqual[N[Exp[a], $MachinePrecision], 0.0], N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision], N[(N[(0.125 * a + 0.5), $MachinePrecision] * a + N[Log[2.0], $MachinePrecision]), $MachinePrecision]]
                  
                  \begin{array}{l}
                  [a, b] = \mathsf{sort}([a, b])\\
                  \\
                  \begin{array}{l}
                  \mathbf{if}\;e^{a} \leq 0:\\
                  \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\
                  
                  \mathbf{else}:\\
                  \;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(0.125, a, 0.5\right), a, \log 2\right)\\
                  
                  
                  \end{array}
                  \end{array}
                  
                  Derivation
                  1. Split input into 2 regimes
                  2. if (exp.f64 a) < 0.0

                    1. Initial program 8.8%

                      \[\log \left(e^{a} + e^{b}\right) \]
                    2. Add Preprocessing
                    3. Taylor expanded in a around 0

                      \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
                    4. Step-by-step derivation
                      1. lower-log1p.f64N/A

                        \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                      2. lower-exp.f643.8

                        \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
                    5. Applied rewrites3.8%

                      \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                    6. Taylor expanded in b around 0

                      \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
                    7. Step-by-step derivation
                      1. Applied rewrites3.9%

                        \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
                      2. Taylor expanded in b around inf

                        \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
                      3. Step-by-step derivation
                        1. Applied rewrites18.4%

                          \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]

                        if 0.0 < (exp.f64 a)

                        1. Initial program 66.3%

                          \[\log \left(e^{a} + e^{b}\right) \]
                        2. Add Preprocessing
                        3. Taylor expanded in b around 0

                          \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
                        4. Step-by-step derivation
                          1. lower-log1p.f64N/A

                            \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                          2. lower-exp.f6463.3

                            \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                        5. Applied rewrites63.3%

                          \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                        6. Taylor expanded in a around 0

                          \[\leadsto \log 2 + \color{blue}{a \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot a\right)} \]
                        7. Step-by-step derivation
                          1. Applied rewrites62.7%

                            \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, a, 0.5\right), \color{blue}{a}, \log 2\right) \]
                        8. Recombined 2 regimes into one program.
                        9. Add Preprocessing

                        Alternative 7: 55.8% accurate, 1.4× speedup?

                        \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{a} \leq 2 \cdot 10^{-39}:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(0.5, a, \log 2\right)\\ \end{array} \end{array} \]
                        NOTE: a and b should be sorted in increasing order before calling this function.
                        (FPCore (a b)
                         :precision binary64
                         (if (<= (exp a) 2e-39) (* (fma 0.125 b 0.5) b) (fma 0.5 a (log 2.0))))
                        assert(a < b);
                        double code(double a, double b) {
                        	double tmp;
                        	if (exp(a) <= 2e-39) {
                        		tmp = fma(0.125, b, 0.5) * b;
                        	} else {
                        		tmp = fma(0.5, a, log(2.0));
                        	}
                        	return tmp;
                        }
                        
                        a, b = sort([a, b])
                        function code(a, b)
                        	tmp = 0.0
                        	if (exp(a) <= 2e-39)
                        		tmp = Float64(fma(0.125, b, 0.5) * b);
                        	else
                        		tmp = fma(0.5, a, log(2.0));
                        	end
                        	return tmp
                        end
                        
                        NOTE: a and b should be sorted in increasing order before calling this function.
                        code[a_, b_] := If[LessEqual[N[Exp[a], $MachinePrecision], 2e-39], N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision], N[(0.5 * a + N[Log[2.0], $MachinePrecision]), $MachinePrecision]]
                        
                        \begin{array}{l}
                        [a, b] = \mathsf{sort}([a, b])\\
                        \\
                        \begin{array}{l}
                        \mathbf{if}\;e^{a} \leq 2 \cdot 10^{-39}:\\
                        \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\
                        
                        \mathbf{else}:\\
                        \;\;\;\;\mathsf{fma}\left(0.5, a, \log 2\right)\\
                        
                        
                        \end{array}
                        \end{array}
                        
                        Derivation
                        1. Split input into 2 regimes
                        2. if (exp.f64 a) < 1.99999999999999986e-39

                          1. Initial program 8.8%

                            \[\log \left(e^{a} + e^{b}\right) \]
                          2. Add Preprocessing
                          3. Taylor expanded in a around 0

                            \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
                          4. Step-by-step derivation
                            1. lower-log1p.f64N/A

                              \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                            2. lower-exp.f643.8

                              \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
                          5. Applied rewrites3.8%

                            \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                          6. Taylor expanded in b around 0

                            \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
                          7. Step-by-step derivation
                            1. Applied rewrites3.9%

                              \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
                            2. Taylor expanded in b around inf

                              \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
                            3. Step-by-step derivation
                              1. Applied rewrites18.4%

                                \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]

                              if 1.99999999999999986e-39 < (exp.f64 a)

                              1. Initial program 66.3%

                                \[\log \left(e^{a} + e^{b}\right) \]
                              2. Add Preprocessing
                              3. Taylor expanded in b around 0

                                \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
                              4. Step-by-step derivation
                                1. lower-log1p.f64N/A

                                  \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                2. lower-exp.f6463.3

                                  \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                              5. Applied rewrites63.3%

                                \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                              6. Taylor expanded in a around 0

                                \[\leadsto \log 2 + \color{blue}{\frac{1}{2} \cdot a} \]
                              7. Step-by-step derivation
                                1. Applied rewrites62.5%

                                  \[\leadsto \mathsf{fma}\left(0.5, \color{blue}{a}, \log 2\right) \]
                              8. Recombined 2 regimes into one program.
                              9. Add Preprocessing

                              Alternative 8: 55.8% accurate, 1.4× speedup?

                              \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{a} \leq 2 \cdot 10^{-39}:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\mathsf{log1p}\left(1 + a\right)\\ \end{array} \end{array} \]
                              NOTE: a and b should be sorted in increasing order before calling this function.
                              (FPCore (a b)
                               :precision binary64
                               (if (<= (exp a) 2e-39) (* (fma 0.125 b 0.5) b) (log1p (+ 1.0 a))))
                              assert(a < b);
                              double code(double a, double b) {
                              	double tmp;
                              	if (exp(a) <= 2e-39) {
                              		tmp = fma(0.125, b, 0.5) * b;
                              	} else {
                              		tmp = log1p((1.0 + a));
                              	}
                              	return tmp;
                              }
                              
                              a, b = sort([a, b])
                              function code(a, b)
                              	tmp = 0.0
                              	if (exp(a) <= 2e-39)
                              		tmp = Float64(fma(0.125, b, 0.5) * b);
                              	else
                              		tmp = log1p(Float64(1.0 + a));
                              	end
                              	return tmp
                              end
                              
                              NOTE: a and b should be sorted in increasing order before calling this function.
                              code[a_, b_] := If[LessEqual[N[Exp[a], $MachinePrecision], 2e-39], N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision], N[Log[1 + N[(1.0 + a), $MachinePrecision]], $MachinePrecision]]
                              
                              \begin{array}{l}
                              [a, b] = \mathsf{sort}([a, b])\\
                              \\
                              \begin{array}{l}
                              \mathbf{if}\;e^{a} \leq 2 \cdot 10^{-39}:\\
                              \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\
                              
                              \mathbf{else}:\\
                              \;\;\;\;\mathsf{log1p}\left(1 + a\right)\\
                              
                              
                              \end{array}
                              \end{array}
                              
                              Derivation
                              1. Split input into 2 regimes
                              2. if (exp.f64 a) < 1.99999999999999986e-39

                                1. Initial program 8.8%

                                  \[\log \left(e^{a} + e^{b}\right) \]
                                2. Add Preprocessing
                                3. Taylor expanded in a around 0

                                  \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
                                4. Step-by-step derivation
                                  1. lower-log1p.f64N/A

                                    \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                  2. lower-exp.f643.8

                                    \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
                                5. Applied rewrites3.8%

                                  \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                6. Taylor expanded in b around 0

                                  \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
                                7. Step-by-step derivation
                                  1. Applied rewrites3.9%

                                    \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
                                  2. Taylor expanded in b around inf

                                    \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
                                  3. Step-by-step derivation
                                    1. Applied rewrites18.4%

                                      \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]

                                    if 1.99999999999999986e-39 < (exp.f64 a)

                                    1. Initial program 66.3%

                                      \[\log \left(e^{a} + e^{b}\right) \]
                                    2. Add Preprocessing
                                    3. Taylor expanded in b around 0

                                      \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
                                    4. Step-by-step derivation
                                      1. lower-log1p.f64N/A

                                        \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                      2. lower-exp.f6463.3

                                        \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                                    5. Applied rewrites63.3%

                                      \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                    6. Taylor expanded in a around 0

                                      \[\leadsto \mathsf{log1p}\left(1 + a\right) \]
                                    7. Step-by-step derivation
                                      1. Applied rewrites62.4%

                                        \[\leadsto \mathsf{log1p}\left(1 + a\right) \]
                                    8. Recombined 2 regimes into one program.
                                    9. Add Preprocessing

                                    Alternative 9: 57.6% accurate, 1.5× speedup?

                                    \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ 0.5 \cdot b + \mathsf{log1p}\left(e^{a}\right) \end{array} \]
                                    NOTE: a and b should be sorted in increasing order before calling this function.
                                    (FPCore (a b) :precision binary64 (+ (* 0.5 b) (log1p (exp a))))
                                    assert(a < b);
                                    double code(double a, double b) {
                                    	return (0.5 * b) + log1p(exp(a));
                                    }
                                    
                                    assert a < b;
                                    public static double code(double a, double b) {
                                    	return (0.5 * b) + Math.log1p(Math.exp(a));
                                    }
                                    
                                    [a, b] = sort([a, b])
                                    def code(a, b):
                                    	return (0.5 * b) + math.log1p(math.exp(a))
                                    
                                    a, b = sort([a, b])
                                    function code(a, b)
                                    	return Float64(Float64(0.5 * b) + log1p(exp(a)))
                                    end
                                    
                                    NOTE: a and b should be sorted in increasing order before calling this function.
                                    code[a_, b_] := N[(N[(0.5 * b), $MachinePrecision] + N[Log[1 + N[Exp[a], $MachinePrecision]], $MachinePrecision]), $MachinePrecision]
                                    
                                    \begin{array}{l}
                                    [a, b] = \mathsf{sort}([a, b])\\
                                    \\
                                    0.5 \cdot b + \mathsf{log1p}\left(e^{a}\right)
                                    \end{array}
                                    
                                    Derivation
                                    1. Initial program 52.2%

                                      \[\log \left(e^{a} + e^{b}\right) \]
                                    2. Add Preprocessing
                                    3. Taylor expanded in b around 0

                                      \[\leadsto \color{blue}{\log \left(1 + e^{a}\right) + \frac{b}{1 + e^{a}}} \]
                                    4. Step-by-step derivation
                                      1. +-commutativeN/A

                                        \[\leadsto \color{blue}{\frac{b}{1 + e^{a}} + \log \left(1 + e^{a}\right)} \]
                                      2. *-rgt-identityN/A

                                        \[\leadsto \frac{\color{blue}{b \cdot 1}}{1 + e^{a}} + \log \left(1 + e^{a}\right) \]
                                      3. associate-*r/N/A

                                        \[\leadsto \color{blue}{b \cdot \frac{1}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
                                      4. lower-+.f64N/A

                                        \[\leadsto \color{blue}{b \cdot \frac{1}{1 + e^{a}} + \log \left(1 + e^{a}\right)} \]
                                      5. associate-*r/N/A

                                        \[\leadsto \color{blue}{\frac{b \cdot 1}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
                                      6. *-rgt-identityN/A

                                        \[\leadsto \frac{\color{blue}{b}}{1 + e^{a}} + \log \left(1 + e^{a}\right) \]
                                      7. lower-/.f64N/A

                                        \[\leadsto \color{blue}{\frac{b}{1 + e^{a}}} + \log \left(1 + e^{a}\right) \]
                                      8. +-commutativeN/A

                                        \[\leadsto \frac{b}{\color{blue}{e^{a} + 1}} + \log \left(1 + e^{a}\right) \]
                                      9. lower-+.f64N/A

                                        \[\leadsto \frac{b}{\color{blue}{e^{a} + 1}} + \log \left(1 + e^{a}\right) \]
                                      10. lower-exp.f64N/A

                                        \[\leadsto \frac{b}{\color{blue}{e^{a}} + 1} + \log \left(1 + e^{a}\right) \]
                                      11. lower-log1p.f64N/A

                                        \[\leadsto \frac{b}{e^{a} + 1} + \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                      12. lower-exp.f6471.7

                                        \[\leadsto \frac{b}{e^{a} + 1} + \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                                    5. Applied rewrites71.7%

                                      \[\leadsto \color{blue}{\frac{b}{e^{a} + 1} + \mathsf{log1p}\left(e^{a}\right)} \]
                                    6. Taylor expanded in a around 0

                                      \[\leadsto \frac{1}{2} \cdot b + \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                                    7. Step-by-step derivation
                                      1. Applied rewrites52.0%

                                        \[\leadsto 0.5 \cdot b + \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                                      2. Add Preprocessing

                                      Alternative 10: 55.4% accurate, 1.5× speedup?

                                      \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \begin{array}{l} \mathbf{if}\;e^{a} \leq 0:\\ \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\ \mathbf{else}:\\ \;\;\;\;\mathsf{log1p}\left(1\right)\\ \end{array} \end{array} \]
                                      NOTE: a and b should be sorted in increasing order before calling this function.
                                      (FPCore (a b)
                                       :precision binary64
                                       (if (<= (exp a) 0.0) (* (fma 0.125 b 0.5) b) (log1p 1.0)))
                                      assert(a < b);
                                      double code(double a, double b) {
                                      	double tmp;
                                      	if (exp(a) <= 0.0) {
                                      		tmp = fma(0.125, b, 0.5) * b;
                                      	} else {
                                      		tmp = log1p(1.0);
                                      	}
                                      	return tmp;
                                      }
                                      
                                      a, b = sort([a, b])
                                      function code(a, b)
                                      	tmp = 0.0
                                      	if (exp(a) <= 0.0)
                                      		tmp = Float64(fma(0.125, b, 0.5) * b);
                                      	else
                                      		tmp = log1p(1.0);
                                      	end
                                      	return tmp
                                      end
                                      
                                      NOTE: a and b should be sorted in increasing order before calling this function.
                                      code[a_, b_] := If[LessEqual[N[Exp[a], $MachinePrecision], 0.0], N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision], N[Log[1 + 1.0], $MachinePrecision]]
                                      
                                      \begin{array}{l}
                                      [a, b] = \mathsf{sort}([a, b])\\
                                      \\
                                      \begin{array}{l}
                                      \mathbf{if}\;e^{a} \leq 0:\\
                                      \;\;\;\;\mathsf{fma}\left(0.125, b, 0.5\right) \cdot b\\
                                      
                                      \mathbf{else}:\\
                                      \;\;\;\;\mathsf{log1p}\left(1\right)\\
                                      
                                      
                                      \end{array}
                                      \end{array}
                                      
                                      Derivation
                                      1. Split input into 2 regimes
                                      2. if (exp.f64 a) < 0.0

                                        1. Initial program 8.8%

                                          \[\log \left(e^{a} + e^{b}\right) \]
                                        2. Add Preprocessing
                                        3. Taylor expanded in a around 0

                                          \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
                                        4. Step-by-step derivation
                                          1. lower-log1p.f64N/A

                                            \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                          2. lower-exp.f643.8

                                            \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
                                        5. Applied rewrites3.8%

                                          \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                        6. Taylor expanded in b around 0

                                          \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
                                        7. Step-by-step derivation
                                          1. Applied rewrites3.9%

                                            \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
                                          2. Taylor expanded in b around inf

                                            \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
                                          3. Step-by-step derivation
                                            1. Applied rewrites18.4%

                                              \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]

                                            if 0.0 < (exp.f64 a)

                                            1. Initial program 66.3%

                                              \[\log \left(e^{a} + e^{b}\right) \]
                                            2. Add Preprocessing
                                            3. Taylor expanded in b around 0

                                              \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
                                            4. Step-by-step derivation
                                              1. lower-log1p.f64N/A

                                                \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                              2. lower-exp.f6463.3

                                                \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                                            5. Applied rewrites63.3%

                                              \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                            6. Taylor expanded in a around 0

                                              \[\leadsto \mathsf{log1p}\left(1\right) \]
                                            7. Step-by-step derivation
                                              1. Applied rewrites62.0%

                                                \[\leadsto \mathsf{log1p}\left(1\right) \]
                                            8. Recombined 2 regimes into one program.
                                            9. Add Preprocessing

                                            Alternative 11: 11.9% accurate, 25.3× speedup?

                                            \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \end{array} \]
                                            NOTE: a and b should be sorted in increasing order before calling this function.
                                            (FPCore (a b) :precision binary64 (* (fma 0.125 b 0.5) b))
                                            assert(a < b);
                                            double code(double a, double b) {
                                            	return fma(0.125, b, 0.5) * b;
                                            }
                                            
                                            a, b = sort([a, b])
                                            function code(a, b)
                                            	return Float64(fma(0.125, b, 0.5) * b)
                                            end
                                            
                                            NOTE: a and b should be sorted in increasing order before calling this function.
                                            code[a_, b_] := N[(N[(0.125 * b + 0.5), $MachinePrecision] * b), $MachinePrecision]
                                            
                                            \begin{array}{l}
                                            [a, b] = \mathsf{sort}([a, b])\\
                                            \\
                                            \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b
                                            \end{array}
                                            
                                            Derivation
                                            1. Initial program 52.2%

                                              \[\log \left(e^{a} + e^{b}\right) \]
                                            2. Add Preprocessing
                                            3. Taylor expanded in a around 0

                                              \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
                                            4. Step-by-step derivation
                                              1. lower-log1p.f64N/A

                                                \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                              2. lower-exp.f6449.2

                                                \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
                                            5. Applied rewrites49.2%

                                              \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                            6. Taylor expanded in b around 0

                                              \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
                                            7. Step-by-step derivation
                                              1. Applied rewrites47.3%

                                                \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
                                              2. Taylor expanded in b around inf

                                                \[\leadsto {b}^{2} \cdot \left(\frac{1}{8} + \color{blue}{\frac{1}{2} \cdot \frac{1}{b}}\right) \]
                                              3. Step-by-step derivation
                                                1. Applied rewrites6.9%

                                                  \[\leadsto \mathsf{fma}\left(0.125, b, 0.5\right) \cdot b \]
                                                2. Add Preprocessing

                                                Alternative 12: 5.3% accurate, 27.6× speedup?

                                                \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \left(b \cdot b\right) \cdot 0.125 \end{array} \]
                                                NOTE: a and b should be sorted in increasing order before calling this function.
                                                (FPCore (a b) :precision binary64 (* (* b b) 0.125))
                                                assert(a < b);
                                                double code(double a, double b) {
                                                	return (b * b) * 0.125;
                                                }
                                                
                                                NOTE: a and b should be sorted in increasing order before calling this function.
                                                real(8) function code(a, b)
                                                    real(8), intent (in) :: a
                                                    real(8), intent (in) :: b
                                                    code = (b * b) * 0.125d0
                                                end function
                                                
                                                assert a < b;
                                                public static double code(double a, double b) {
                                                	return (b * b) * 0.125;
                                                }
                                                
                                                [a, b] = sort([a, b])
                                                def code(a, b):
                                                	return (b * b) * 0.125
                                                
                                                a, b = sort([a, b])
                                                function code(a, b)
                                                	return Float64(Float64(b * b) * 0.125)
                                                end
                                                
                                                a, b = num2cell(sort([a, b])){:}
                                                function tmp = code(a, b)
                                                	tmp = (b * b) * 0.125;
                                                end
                                                
                                                NOTE: a and b should be sorted in increasing order before calling this function.
                                                code[a_, b_] := N[(N[(b * b), $MachinePrecision] * 0.125), $MachinePrecision]
                                                
                                                \begin{array}{l}
                                                [a, b] = \mathsf{sort}([a, b])\\
                                                \\
                                                \left(b \cdot b\right) \cdot 0.125
                                                \end{array}
                                                
                                                Derivation
                                                1. Initial program 52.2%

                                                  \[\log \left(e^{a} + e^{b}\right) \]
                                                2. Add Preprocessing
                                                3. Taylor expanded in a around 0

                                                  \[\leadsto \color{blue}{\log \left(1 + e^{b}\right)} \]
                                                4. Step-by-step derivation
                                                  1. lower-log1p.f64N/A

                                                    \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                                  2. lower-exp.f6449.2

                                                    \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{b}}\right) \]
                                                5. Applied rewrites49.2%

                                                  \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{b}\right)} \]
                                                6. Taylor expanded in b around 0

                                                  \[\leadsto \log 2 + \color{blue}{b \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot b\right)} \]
                                                7. Step-by-step derivation
                                                  1. Applied rewrites47.3%

                                                    \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, b, 0.5\right), \color{blue}{b}, \log 2\right) \]
                                                  2. Taylor expanded in b around inf

                                                    \[\leadsto \frac{1}{8} \cdot {b}^{\color{blue}{2}} \]
                                                  3. Step-by-step derivation
                                                    1. Applied rewrites4.0%

                                                      \[\leadsto \left(b \cdot b\right) \cdot 0.125 \]
                                                    2. Add Preprocessing

                                                    Alternative 13: 3.2% accurate, 27.6× speedup?

                                                    \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ \left(0.125 \cdot a\right) \cdot a \end{array} \]
                                                    NOTE: a and b should be sorted in increasing order before calling this function.
                                                    (FPCore (a b) :precision binary64 (* (* 0.125 a) a))
                                                    assert(a < b);
                                                    double code(double a, double b) {
                                                    	return (0.125 * a) * a;
                                                    }
                                                    
                                                    NOTE: a and b should be sorted in increasing order before calling this function.
                                                    real(8) function code(a, b)
                                                        real(8), intent (in) :: a
                                                        real(8), intent (in) :: b
                                                        code = (0.125d0 * a) * a
                                                    end function
                                                    
                                                    assert a < b;
                                                    public static double code(double a, double b) {
                                                    	return (0.125 * a) * a;
                                                    }
                                                    
                                                    [a, b] = sort([a, b])
                                                    def code(a, b):
                                                    	return (0.125 * a) * a
                                                    
                                                    a, b = sort([a, b])
                                                    function code(a, b)
                                                    	return Float64(Float64(0.125 * a) * a)
                                                    end
                                                    
                                                    a, b = num2cell(sort([a, b])){:}
                                                    function tmp = code(a, b)
                                                    	tmp = (0.125 * a) * a;
                                                    end
                                                    
                                                    NOTE: a and b should be sorted in increasing order before calling this function.
                                                    code[a_, b_] := N[(N[(0.125 * a), $MachinePrecision] * a), $MachinePrecision]
                                                    
                                                    \begin{array}{l}
                                                    [a, b] = \mathsf{sort}([a, b])\\
                                                    \\
                                                    \left(0.125 \cdot a\right) \cdot a
                                                    \end{array}
                                                    
                                                    Derivation
                                                    1. Initial program 52.2%

                                                      \[\log \left(e^{a} + e^{b}\right) \]
                                                    2. Add Preprocessing
                                                    3. Taylor expanded in b around 0

                                                      \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
                                                    4. Step-by-step derivation
                                                      1. lower-log1p.f64N/A

                                                        \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                                      2. lower-exp.f6449.0

                                                        \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                                                    5. Applied rewrites49.0%

                                                      \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                                    6. Taylor expanded in a around 0

                                                      \[\leadsto \log 2 + \color{blue}{a \cdot \left(\frac{1}{2} + \frac{1}{8} \cdot a\right)} \]
                                                    7. Step-by-step derivation
                                                      1. Applied rewrites47.8%

                                                        \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(0.125, a, 0.5\right), \color{blue}{a}, \log 2\right) \]
                                                      2. Taylor expanded in a around inf

                                                        \[\leadsto \frac{1}{8} \cdot {a}^{\color{blue}{2}} \]
                                                      3. Step-by-step derivation
                                                        1. Applied rewrites4.3%

                                                          \[\leadsto \left(0.125 \cdot a\right) \cdot a \]
                                                        2. Add Preprocessing

                                                        Alternative 14: 2.6% accurate, 50.7× speedup?

                                                        \[\begin{array}{l} [a, b] = \mathsf{sort}([a, b])\\ \\ 0.5 \cdot a \end{array} \]
                                                        NOTE: a and b should be sorted in increasing order before calling this function.
                                                        (FPCore (a b) :precision binary64 (* 0.5 a))
                                                        assert(a < b);
                                                        double code(double a, double b) {
                                                        	return 0.5 * a;
                                                        }
                                                        
                                                        NOTE: a and b should be sorted in increasing order before calling this function.
                                                        real(8) function code(a, b)
                                                            real(8), intent (in) :: a
                                                            real(8), intent (in) :: b
                                                            code = 0.5d0 * a
                                                        end function
                                                        
                                                        assert a < b;
                                                        public static double code(double a, double b) {
                                                        	return 0.5 * a;
                                                        }
                                                        
                                                        [a, b] = sort([a, b])
                                                        def code(a, b):
                                                        	return 0.5 * a
                                                        
                                                        a, b = sort([a, b])
                                                        function code(a, b)
                                                        	return Float64(0.5 * a)
                                                        end
                                                        
                                                        a, b = num2cell(sort([a, b])){:}
                                                        function tmp = code(a, b)
                                                        	tmp = 0.5 * a;
                                                        end
                                                        
                                                        NOTE: a and b should be sorted in increasing order before calling this function.
                                                        code[a_, b_] := N[(0.5 * a), $MachinePrecision]
                                                        
                                                        \begin{array}{l}
                                                        [a, b] = \mathsf{sort}([a, b])\\
                                                        \\
                                                        0.5 \cdot a
                                                        \end{array}
                                                        
                                                        Derivation
                                                        1. Initial program 52.2%

                                                          \[\log \left(e^{a} + e^{b}\right) \]
                                                        2. Add Preprocessing
                                                        3. Taylor expanded in b around 0

                                                          \[\leadsto \color{blue}{\log \left(1 + e^{a}\right)} \]
                                                        4. Step-by-step derivation
                                                          1. lower-log1p.f64N/A

                                                            \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                                          2. lower-exp.f6449.0

                                                            \[\leadsto \mathsf{log1p}\left(\color{blue}{e^{a}}\right) \]
                                                        5. Applied rewrites49.0%

                                                          \[\leadsto \color{blue}{\mathsf{log1p}\left(e^{a}\right)} \]
                                                        6. Taylor expanded in a around 0

                                                          \[\leadsto \log 2 + \color{blue}{\frac{1}{2} \cdot a} \]
                                                        7. Step-by-step derivation
                                                          1. Applied rewrites47.8%

                                                            \[\leadsto \mathsf{fma}\left(0.5, \color{blue}{a}, \log 2\right) \]
                                                          2. Taylor expanded in a around inf

                                                            \[\leadsto \frac{1}{2} \cdot a \]
                                                          3. Step-by-step derivation
                                                            1. Applied rewrites7.6%

                                                              \[\leadsto 0.5 \cdot a \]
                                                            2. Add Preprocessing

                                                            Reproduce

                                                            ?
                                                            herbie shell --seed 2024332 
                                                            (FPCore (a b)
                                                              :name "symmetry log of sum of exp"
                                                              :precision binary64
                                                              (log (+ (exp a) (exp b))))