To compute A + B
within float32
where both A and B is ~1e100, we can do:
log(A + B) = log(exp(a) + exp(b)) = a + softplus(b - a)
where a = log(A)
and b = log(B)
.
/** | |
* This Source Code Form is subject to the terms of the Mozilla Public | |
* License, v. 2.0. If a copy of the MPL was not distributed with this | |
* file, You can obtain one at https://mozilla.org/MPL/2.0/. | |
*/ | |
#pragma once | |
#include <complex> | |
#include <algorithm> | |
// Real softplus. | |
template <typename T> | |
T softplus(T x) | |
{ | |
using namespace std; | |
return log(1 + exp(-abs(x))) + max(x, T(0)); | |
} | |
// Complex softplus. | |
template <typename T> | |
std::complex<T> softplus(std::complex<T> z) | |
{ | |
using namespace std; | |
T r = real(z); | |
T i = imag(z); | |
T logval0 = softplus(2*r); | |
T logval2; | |
T ci = cos(i); | |
if (ci > 0) { | |
T logval1 = r + log(2*ci); | |
logval2 = 0.5*(logval0 + softplus(logval1 - logval0)); | |
} else { | |
T logval1 = r + log(2*abs(ci)); | |
logval2 = 0.5*(logval0 + log(1 - exp(logval1 - logval0))); | |
} | |
T arg = asin(exp(r - logval2) * sin(i)); | |
if (ci < 0 && r > log(-1/ci)) { | |
arg = M_PI - arg; | |
} | |
return std::complex<T>(logval2, arg); | |
} | |
To compute A + B
within float32
where both A and B is ~1e100, we can do:
log(A + B) = log(exp(a) + exp(b)) = a + softplus(b - a)
where a = log(A)
and b = log(B)
.