Symbolics.coeff silently returns incorrect results for products that need expansion, unbounded recursion on fractions.

[snoeyink]

This is related to comments in the discussion between @t-bltg and @bowenszhu on the implementation of Symbolics.coeff in #677, #685, #549. There are four related issues in the current implementation:

1. The implementation of Symbolics.coeff for a product (the ismul(p) case) works correctly if a single term of the product contains the sym that the coeffient is being computed for. It is incorrect if more than one does: e.g., it returns the coefficient of $x^4$ instead of the requested $x^2$ on an unexpanded product:

using Symbolics
@variables x y 
@testbroken isequal(Symbolics.coeff((3x^2+2)*(2x^2+1),x^2), 6)
@test isequal(Symbolics.coeff(expand((3x^2+2)*(2x^2+1)),x^2), 7)
@testbroken isequal(Symbolics.coeff((3x^2+2)*y*(2x^2+y),x^2), 6y)
@test isequal(Symbolics.coeff(expand((3x^2+2)*y*(2x^2+y)),x^2), 4y + 3(y^2))

2. One cannot ask for the constant coefficient. Maple & Mathematica gives a syntax that allows the variable and the power to be given separately so that x^0 coefficient can be correctly computed.

@testbroken isequal(Symbolics.coeff((3x^2+2)*(2x^2+1),x^0), 0)

3. Since no isdiv(p) case is implemented an expression Symbolics.coeff( x/y, x) renders $p = x/y$ isa Symbolic == true, so coeff(x/y, x) is called again, leading to unbounded recursion (until the stack overflows).

4. $x^{-1}$ is translated to $1/x$, so an expression Symbolics.coeff( 3/x, x^-1) produces sym = 1/x This should be caught unless negative powers are to be allowed. (Right now this expression gives the stack overflow due to the first division, but even if that is fixed, p will never match sym = 1/x.

As the comments in the above PRs suggest that users should call expand themselves (Maple and Mathematic don't), so the minimal fix is to throw an error when the sym variable appears in more than one term of the product. I would then extend this to division: if the variable appears only in the numerator, then we get its coefficient and multiply by 1/p.den, otherwise throw an error. Coefficients with negative powers could also produce an error and coeff(p, sym, expon) could be used to special case the case for expon=0.

A more permissive implementation could use a keyword to allow automatic expansion of products and another to support handling negative powers. I can do the minimal fix, but I wanted to get feedback, esp. from @t-bltg and @bowenszhu if possible, before making a PR.

[t-bltg]

Thanks for pinging @snoeyink.
As you have noticed, there were several limitations in the implementation of coeff in #677.
I think you have interesting ideas here (especially comparing to other symbolic computing softwares), unfortunately, I haven't the time to work on these atm, but I believe any improvement to the current state would be very much welcome.