are `expand_mul(expand_multinomial(expr))` sufficient for the task? Ideally we should use whatever Poly is using, perhaps even sending expr to Poly to get the flattened result: `if f.as_expr() == Poly(expr).as_expr() and exp == 1`. And by all means put the `exp == 1` test first since that is a cheap test which, if False, will circumvent the more expensice step of expanding.