Series expansion of the square root gives wrong result

[astholkohtz]

Whenever the constant term on the square root is a sum of 2 or more symbols, the series() function gets the expansion wrong

In [1]: from sympy import *

In [2]: a,b,c,x = symbols('a b c x')

In [3]: sqrt(a+b + c*x).series(x,0,1)
Out[3]: 5*sqrt(a + b)/16 + 15*b/(16*sqrt(a + b)) - 5*b**2*sqrt(a + b)/(16*(a**2 + 2*a*b + b**2)) + b**3*sqrt(a + b)/(16*(a**3 + 3*a**2*b + 3*a*b**2 + b**3)) + 15*a/(16*sqrt(a + b)) - 5*a*b*sqrt(a + b)/(8*(a**2 + 2*a*b + b**2)) + 3*a*b**2*sqrt(a + b)/(16*(a**3 + 3*a**2*b + 3*a*b**2 + b**3)) - 5*a**2*sqrt(a + b)/(16*(a**2 + 2*a*b + b**2)) + 3*a**2*b*sqrt(a + b)/(16*(a**3 + 3*a**2*b + 3*a*b**2 + b**3)) + a**3*sqrt(a + b)/(16*(a**3 + 3*a**2*b + 3*a*b**2 + b**3)) + O(x)

In [4]: sqrt(a+b+c + c*x).series(x,0,1)
Out[4]: 5*sqrt(a + b + c)/16 + 15*c/(16*sqrt(a + b + c)) - 5*c**2*sqrt(a + b + c)/(16*(a**2 + 2*a*b + 2*a*c + b**2 + 2*b*c + c**2)) + c**3*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) + 15*b/(16*sqrt(a + b + c)) - 5*b*c*sqrt(a + b + c)/(8*(a**2 + 2*a*b + 2*a*c + b**2 + 2*b*c + c**2)) + 3*b*c**2*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) - 5*b**2*sqrt(a + b + c)/(16*(a**2 + 2*a*b + 2*a*c + b**2 + 2*b*c + c**2)) + 3*b**2*c*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) + b**3*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) + 15*a/(16*sqrt(a + b + c)) - 5*a*c*sqrt(a + b + c)/(8*(a**2 + 2*a*b + 2*a*c + b**2 + 2*b*c + c**2)) + 3*a*c**2*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) - 5*a*b*sqrt(a + b + c)/(8*(a**2 + 2*a*b + 2*a*c + b**2 + 2*b*c + c**2)) + 3*a*b*c*sqrt(a + b + c)/(8*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) + 3*a*b**2*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) - 5*a**2*sqrt(a + b + c)/(16*(a**2 + 2*a*b + 2*a*c + b**2 + 2*b*c + c**2)) + 3*a**2*c*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) + 3*a**2*b*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) + a**3*sqrt(a + b + c)/(16*(a**3 + 3*a**2*b + 3*a**2*c + 3*a*b**2 + 6*a*b*c + 3*a*c**2 + b**3 + 3*b**2*c + 3*b*c**2 + c**3)) + O(x)

As long as the constant term is only one, series() gives the correct result

In [5]: sqrt(a + c*x).series(x,0,1)
Out[5]: sqrt(a) + O(x)

In [6]: sqrt(a + (a+b)*x).series(x,0,1)
Out[6]: sqrt(a) + O(x)

In [7]: sqrt(a + (a+b+c)*x).series(x,0,1)
Out[7]: sqrt(a) + O(x)

In [8]: sqrt(a + (a+b+c)*x).series(x,0,2)
Out[8]: x*(sqrt(a)/2 + b/(2*sqrt(a)) + c/(2*sqrt(a))) + sqrt(a) + O(x**2)

[oscarbenjamin]

It seems that SymPy comes out with an unnecessarily complicated form:

In [56]: sqrt(a+b+c*x).series(x, 0, 1)                                                                                                                        
Out[56]: 
    _______                        2   _______                3   _______                                     _______                 2   _______           
5⋅╲╱ a + b        15⋅b          5⋅b ⋅╲╱ a + b                b ⋅╲╱ a + b                 15⋅a         5⋅a⋅b⋅╲╱ a + b             3⋅a⋅b ⋅╲╱ a + b            
─────────── + ──────────── - ──────────────────── + ────────────────────────────── + ──────────── - ─────────────────── + ────────────────────────────── - ─
     16            _______      ⎛ 2            2⎞      ⎛ 3      2          2    3⎞        _______     ⎛ 2            2⎞      ⎛ 3      2          2    3⎞    
              16⋅╲╱ a + b    16⋅⎝a  + 2⋅a⋅b + b ⎠   16⋅⎝a  + 3⋅a ⋅b + 3⋅a⋅b  + b ⎠   16⋅╲╱ a + b    8⋅⎝a  + 2⋅a⋅b + b ⎠   16⋅⎝a  + 3⋅a ⋅b + 3⋅a⋅b  + b ⎠   1

     2   _______                2     _______                    3   _______                
  5⋅a ⋅╲╱ a + b              3⋅a ⋅b⋅╲╱ a + b                    a ⋅╲╱ a + b                 
─────────────────── + ────────────────────────────── + ────────────────────────────── + O(x)
  ⎛ 2            2⎞      ⎛ 3      2          2    3⎞      ⎛ 3      2          2    3⎞       
6⋅⎝a  + 2⋅a⋅b + b ⎠   16⋅⎝a  + 3⋅a ⋅b + 3⋅a⋅b  + b ⎠   16⋅⎝a  + 3⋅a ⋅b + 3⋅a⋅b  + b ⎠       

In [57]: simplify(_)                                                                                                                                          
Out[57]: 
 3        2      2      3       
b  + 3⋅a⋅b  + 3⋅a ⋅b + a  + O(x)
────────────────────────────────
    _______ ⎛ 2            2⎞   
  ╲╱ a + b ⋅⎝a  + 2⋅a⋅b + b ⎠   

In [58]: factor(_)                                                                                                                                            
Out[58]: 
 3        2      2      3       
b  + 3⋅a⋅b  + 3⋅a ⋅b + a  + O(x)
────────────────────────────────
                  5/2           
           (a + b)  

The numerator there is (a+b)**3 so it should somehow simplify to sqrt(a+b).

[sachin-4099]

On the current master:

In [1]: a,b,c,x = symbols('a b c x')

In [2]: sqrt(a+b + c*x).series(x,0,1)
Out[2]:
  _______
\/ a + b  + O(x)

In [3]: sqrt(a+b+c + c*x).series(x,0,1)
Out[3]:
  ___________
\/ a + b + c  + O(x)

This issue seems resolved and thus, only requires a test in order to close.