Fix for issue 3891

Pow._eval_nseries: do expand_mul for output if b0!=0.  Now:
```
In [1]: x,a,b = symbols('x,a,b')

In [2]: (sqrt(a+b*x+x**2)).series(x,0,3)
Out[2]: 
             ⎛             2  ⎞                
  b⋅x      2 ⎜   1        b   ⎟     ___    ⎛ 3⎞
─────── + x ⋅⎜─────── - ──────⎟ + ╲╱ a  + O⎝x ⎠
    ___      ⎜    ___      3/2⎟                
2⋅╲╱ a       ⎝2⋅╲╱ a    8⋅a   ⎠                


```