def commutator(op_a, op_b, op_c=None, threshold=1e-12):
    """
    Compute commutator of `op_a` and `op_b` or
    the symmetric double commutator of `op_a`, `op_b` and `op_c`.

    See McWeeny chapter 13.6 Equation of motion methods (page 479)

    | If only `op_a` and `op_b` are provided:
    |     result = A\\*B - B\\*A;
    |
    | If `op_a`, `op_b` and `op_c` are provided:
    |     result = 0.5 \\* (2\\*A\\*B\\*C + 2\\*C\\*B\\*A - B\\*A\\*C - C\\*A\\*B - A\\*C\\*B - B\\*C\\*A)

    Args:
        op_a (WeightedPauliOperator): operator a
        op_b (WeightedPauliOperator): operator b
        op_c (Optional(WeightedPauliOperator)): operator c
        threshold (float): the truncation threshold

    Returns:
        WeightedPauliOperator: the commutator

    Note:
        For the final chop, the original codes only contain the paulis with real coefficient.
    """
    op_ab = op_a * op_b
    op_ba = op_b * op_a
    if op_c is None:
        res = op_ab - op_ba
    else:
        op_ac = op_a * op_c
        op_ca = op_c * op_a
        op_abc = op_ab * op_c
        op_cba = op_c * op_ba
        op_bac = op_ba * op_c
        op_cab = op_c * op_ab
        op_acb = op_ac * op_b
        op_bca = op_b * op_ca
        tmp = op_bac + op_cab + op_acb + op_bca
        tmp = 0.5 * tmp
        res = op_abc + op_cba - tmp
    res.simplify()
    res.chop(threshold)
    return res
