def stable_softminus(x):
    """
    Calculating log(1 - exp(x)) (x < 0) in a numerically stable way
    """
    # negative normal range
    y1 = (1 - x.exp()).log()
    # negative zero
    y2 = (- x - x.pow(2)/2 - x.pow(3)/6).log()
    # negative infinity
    y3 = - x.exp() - (2*x).exp()/2 - (3*x).exp()/3
    switch2 = (1 - x.exp()) == 0
    switch3 = (1 - x.exp()) == 1
    switch1 = ~switch2 ^ switch3

    y1[~switch1] = 0
    y2[~switch2] = 0
    y3[~switch3] = 0
    y = y1 + y2 + y3

    return y