text: from sympy.core.function import Add, ArgumentIndexError, Function
  from sympy.core.power import Pow
  from sympy.core.singleton import S
  from sympy.core.sorting import default_sort_key
  from sympy.functions.elementary.exponential import exp, log
  
  
  def _logaddexp(x1, x2, *, evaluate=True):
      return log(Add(exp(x1, evaluate=evaluate), exp(x2, evaluate=evaluate), evaluate=evaluate))
  
  
  _two = S.One*2
  _ln2 = log(_two)
  
  
  def _lb(x, *, evaluate=True):
      return log(x, evaluate=evaluate)/_ln2
  
  
  def _exp2(x, *, evaluate=True):
      return Pow(_two, x, evaluate=evaluate)
  
  
  def _logaddexp2(x1, x2, *, evaluate=True):
      return _lb(Add(_exp2(x1, evaluate=evaluate),
                     _exp2(x2, evaluate=evaluate), evaluate=evaluate))
  
  
  class logaddexp(Function):
      """ Logarithm of the sum of exponentiations of the inputs.
  
      Helper class for use with e.g. numpy.logaddexp
  
      See Also
      ========
  
      https://numpy.org/doc/stable/reference/generated/numpy.logaddexp.html
      """
      nargs = 2
  
      def __new__(cls, *args):
          return Function.__new__(cls, *sorted(args, key=default_sort_key))
  
      def fdiff(self, argindex=1):
          """
          Returns the first derivative of this function.
          """
          if argindex == 1:
              wrt, other = self.args
          elif argindex == 2:
              other, wrt = self.args
          else:
              raise ArgumentIndexError(self, argindex)
          return S.One/(S.One + exp(other-wrt))
  
      def _eval_rewrite_as_log(self, x1, x2, **kwargs):
          return _logaddexp(x1, x2)
  
      def _eval_evalf(self, *args, **kwargs):
          return self.rewrite(log).evalf(*args, **kwargs)
  
      def _eval_simplify(self, *args, **kwargs):
          a, b = map(lambda x: x.simplify(**kwargs), self.args)
          candidate = _logaddexp(a, b)
          if candidate != _logaddexp(a, b, evaluate=False):
              return candidate
          else:
              return logaddexp(a, b)
  
  
  class logaddexp2(Function):
      """ Logarithm of the sum of exponentiations of the inputs in base-2.
  
      Helper class for use with e.g. numpy.logaddexp2
  
      See Also
      ========
  
      https://numpy.org/doc/stable/reference/generated/numpy.logaddexp2.html
      """
      nargs = 2
  
      def __new__(cls, *args):
          return Function.__new__(cls, *sorted(args, key=default_sort_key))
  
      def fdiff(self, argindex=1):
          """
          Returns the first derivative of this function.
          """
          if argindex == 1:
              wrt, other = self.args
          elif argindex == 2:
              other, wrt = self.args
          else:
              raise ArgumentIndexError(self, argindex)
          return S.One/(S.One + _exp2(other-wrt))
  
      def _eval_rewrite_as_log(self, x1, x2, **kwargs):
          return _logaddexp2(x1, x2)
  
      def _eval_evalf(self, *args, **kwargs):
          return self.rewrite(log).evalf(*args, **kwargs)
  
      def _eval_simplify(self, *args, **kwargs):
          a, b = map(lambda x: x.simplify(**kwargs).factor(), self.args)
          candidate = _logaddexp2(a, b)
          if candidate != _logaddexp2(a, b, evaluate=False):
              return candidate
          else:
              return logaddexp2(a, b)
  