(f)
return result
Calulate exact roots of a solvable quintic
Calulate exact roots of a solvable quintic
Calulate exact roots of a solvable quintic
def roots_quintic(f):
    """
    Calulate exact roots of a solvable quintic
    """
    result = []
    coeff_5, coeff_4, p, q, r, s = f.all_coeffs()

    # Eqn must be of the form x^5 + px^3 + qx^2 + rx + s
    if coeff_4:
        return result

    if coeff_5 != 1:
        l = [p/coeff_5, q/coeff_5, r/coeff_5, s/coeff_5]
        if not all(coeff.is_Rational for coeff in l):
            return result
        f = Poly(f/coeff_5)
    quintic = PolyQuintic(f)

    # Eqn standardized. Algo for solving starts here
    if not f.is_irreducible:
        return result

    f20 = quintic.f20
    # Check if f20 has linear factors over domain Z
    if f20.is_irreducible:
        return result

    # Now, we know that f is solvable
    for _factor in f20.factor_list()[1]:
        if _factor[0].is_linear:
            theta = _factor[0].root(0)
            break
    d = discriminant(f)
    delta = sqrt(d)
    # zeta = a fifth root of unity
    zeta1, zeta2, zeta3, zeta4 = quintic.zeta
    T = quintic.T(theta, d)
    tol = S(1e-10)
    alpha = T[1] + T[2]*delta
    alpha_bar = T[1] - T[2]*delta
    beta = T[3] + T[4]*delta
    beta_bar = T[3] - T[4]*delta

    disc = alpha**2 - 4*beta
    disc_bar = alpha_bar**2 - 4*beta_bar

    l0 = quintic.l0(theta)

    l1 = _quintic_simplify((-alpha + sqrt(disc)) / S(2))
    l4 = _quintic_simplify((-alpha - sqrt(disc)) / S(2))

    l2 = _quintic_simplify((-alpha_bar + sqrt(disc_bar)) / S(2))
    l3 = _quintic_simplify((-alpha_bar - sqrt(disc_bar)) / S(2))

    order = quintic.order(theta, d)
    test = (order*delta.n()) - ( (l1.n() - l4.n())*(l2.n() - l3.n()) )
    # Comparing floats
    # Problems importing on top
    from sympy.utilities.randtest import comp
    if not comp(test, 0, tol):
        l2, l3 = l3, l2

    # Now we have correct order of l's
    R1 = l0 + l1*zeta1 + l2*zeta2 + l3*zeta3 + l4*zeta4
    R2 = l0 + l3*zeta1 + l1*zeta2 + l4*zeta3 + l2*zeta4
    R3 = l0 + l2*zeta1 + l4*zeta2 + l1*zeta3 + l3*zeta4
    R4 = l0 + l4*zeta1 + l3*zeta2 + l2*zeta3 + l1*zeta4

    Res = [None, [None]*5, [None]*5, [None]*5, [None]*5]
    Res_n = [None, [None]*5, [None]*5, [None]*5, [None]*5]
    sol = Symbol('sol')

    # Simplifying improves performace a lot for exact expressions
    R1 = _quintic_simplify(R1)
    R2 = _quintic_simplify(R2)
    R3 = _quintic_simplify(R3)
    R4 = _quintic_simplify(R4)

    # Solve imported here. Causing problems if imported as 'solve'
    # and hence the changed name
    from sympy.solvers.solvers import solve as _solve
    a, b = symbols('a b', cls=Dummy)
    _sol = _solve( sol**5 - a - I*b, sol)
    for i in range(5):
        _sol[i] = factor(_sol[i])
    R1 = R1.as_real_imag()
    R2 = R2.as_real_imag()
    R3 = R3.as_real_imag()
    R4 = R4.as_real_imag()

    for i, root in enumerate(_sol):
        Res[1][i] = _quintic_simplify(root.subs({ a: R1[0], b: R1[1] }))
        Res[2][i] = _quintic_simplify(root.subs({ a: R2[0], b: R2[1] }))
        Res[3][i] = _quintic_simplify(root.subs({ a: R3[0], b: R3[1] }))
        Res[4][i] = _quintic_simplify(root.subs({ a: R4[0], b: R4[1] }))

    for i in range(1, 5):
        for j in range(5):
            Res_n[i][j] = Res[i][j].n()
            Res[i][j] = _quintic_simplify(Res[i][j])
    r1 = Res[1][0]
    r1_n = Res_n[1][0]

    for i in range(5):
        if comp(im(r1_n*Res_n[4][i]), 0, tol):
            r4 = Res[4][i]
            break

    u, v = quintic.uv(theta, d)
    sqrt5 = math.sqrt(5)

    # Now we have various Res values. Each will be a list of five
    # values. We have to pick one r value from those five for each Res
    u, v = quintic.uv(theta, d)
    testplus = (u + v*delta*sqrt(5)).n()
    testminus = (u - v*delta*sqrt(5)).n()

    # Evaluated numbers suffixed with _n
    # We will use evaluated numbers for calculation. Much faster.
    r4_n = r4.n()
    r2 = r3 = None

    for i in range(5):
        r2temp_n = Res_n[2][i]
        for j in range(5):
            # Again storing away the exact number and using
            # evaluated numbers in computations
            r3temp_n = Res_n[3][j]

            if( comp( r1_n*r2temp_n**2 + r4_n*r3temp_n**2 - testplus, 0, tol) and
                comp( r3temp_n*r1_n**2 + r2temp_n*r4_n**2 - testminus, 0, tol ) ):
                r2 = Res[2][i]
                r3 = Res[3][j]
                break
        if r2:
            break

    # Now, we have r's so we can get roots
    x1 = (r1 + r2 + r3 + r4)/5
    x2 = (r1*zeta4 + r2*zeta3 + r3*zeta2 + r4*zeta1)/5
    x3 = (r1*zeta3 + r2*zeta1 + r3*zeta4 + r4*zeta2)/5
    x4 = (r1*zeta2 + r2*zeta4 + r3*zeta1 + r4*zeta3)/5
    x5 = (r1*zeta1 + r2*zeta2 + r3*zeta3 + r4*zeta4)/5
    result = [x1, x2, x3, x4, x5]

    # Now check if solutions are distinct

    result_n = []
    for root in result:
        result_n.append(root.n(5))
    result_n = sorted(result_n, key=default_sort_key)

    prev_entry = None
    for r in result_n:
        if r == prev_entry:
            # Roots are identical. Abort. Return []
            # and fall back to usual solve
            return []
        prev_entry = r

    return result
def roots_quintic ( f ) : result = [ ] coeff_5 , coeff_4 , p , q , r , s = f . all_coeffs ( ) # Eqn must be of the form x^5 + px^3 + qx^2 + rx + s if coeff_4 : return result if coeff_5 != 1 : l = [ p / coeff_5 , q / coeff_5 , r / coeff_5 , s / coeff_5 ] if not all ( coeff . is_Rational for coeff in l ) : return result f = Poly ( f / coeff_5 ) quintic = PolyQuintic ( f ) # Eqn standardized. Algo for solving starts here if not f . is_irreducible : return result f20 = quintic . f20 # Check if f20 has linear factors over domain Z if f20 . is_irreducible : return result # Now, we know that f is solvable for _factor in f20 . factor_list ( ) [ 1 ] : if _factor [ 0 ] . is_linear : theta = _factor [ 0 ] . root ( 0 ) break d = discriminant ( f ) delta = sqrt ( d ) # zeta = a fifth root of unity zeta1 , zeta2 , zeta3 , zeta4 = quintic . zeta T = quintic . T ( theta , d ) tol = S ( 1e-10 ) alpha = T [ 1 ] + T [ 2 ] * delta alpha_bar = T [ 1 ] - T [ 2 ] * delta beta = T [ 3 ] + T [ 4 ] * delta beta_bar = T [ 3 ] - T [ 4 ] * delta disc = alpha ** 2 - 4 * beta disc_bar = alpha_bar ** 2 - 4 * beta_bar l0 = quintic . l0 ( theta ) l1 = _quintic_simplify ( ( - alpha + sqrt ( disc ) ) / S ( 2 ) ) l4 = _quintic_simplify ( ( - alpha - sqrt ( disc ) ) / S ( 2 ) ) l2 = _quintic_simplify ( ( - alpha_bar + sqrt ( disc_bar ) ) / S ( 2 ) ) l3 = _quintic_simplify ( ( - alpha_bar - sqrt ( disc_bar ) ) / S ( 2 ) ) order = quintic . order ( theta , d ) test = ( order * delta . n ( ) ) - ( ( l1 . n ( ) - l4 . n ( ) ) * ( l2 . n ( ) - l3 . n ( ) ) ) # Comparing floats # Problems importing on top from sympy . utilities . randtest import comp if not comp ( test , 0 , tol ) : l2 , l3 = l3 , l2 # Now we have correct order of l's R1 = l0 + l1 * zeta1 + l2 * zeta2 + l3 * zeta3 + l4 * zeta4 R2 = l0 + l3 * zeta1 + l1 * zeta2 + l4 * zeta3 + l2 * zeta4 R3 = l0 + l2 * zeta1 + l4 * zeta2 + l1 * zeta3 + l3 * zeta4 R4 = l0 + l4 * zeta1 + l3 * zeta2 + l2 * zeta3 + l1 * zeta4 Res = [ None , [ None ] * 5 , [ None ] * 5 , [ None ] * 5 , [ None ] * 5 ] Res_n = [ None , [ None ] * 5 , [ None ] * 5 , [ None ] * 5 , [ None ] * 5 ] sol = Symbol ( 'sol' ) # Simplifying improves performace a lot for exact expressions R1 = _quintic_simplify ( R1 ) R2 = _quintic_simplify ( R2 ) R3 = _quintic_simplify ( R3 ) R4 = _quintic_simplify ( R4 ) # Solve imported here. Causing problems if imported as 'solve' # and hence the changed name from sympy . solvers . solvers import solve as _solve a , b = symbols ( 'a b' , cls = Dummy ) _sol = _solve ( sol ** 5 - a - I * b , sol ) for i in range ( 5 ) : _sol [ i ] = factor ( _sol [ i ] ) R1 = R1 . as_real_imag ( ) R2 = R2 . as_real_imag ( ) R3 = R3 . as_real_imag ( ) R4 = R4 . as_real_imag ( ) for i , root in enumerate ( _sol ) : Res [ 1 ] [ i ] = _quintic_simplify ( root . subs ( { a : R1 [ 0 ] , b : R1 [ 1 ] } ) ) Res [ 2 ] [ i ] = _quintic_simplify ( root . subs ( { a : R2 [ 0 ] , b : R2 [ 1 ] } ) ) Res [ 3 ] [ i ] = _quintic_simplify ( root . subs ( { a : R3 [ 0 ] , b : R3 [ 1 ] } ) ) Res [ 4 ] [ i ] = _quintic_simplify ( root . subs ( { a : R4 [ 0 ] , b : R4 [ 1 ] } ) ) for i in range ( 1 , 5 ) : for j in range ( 5 ) : Res_n [ i ] [ j ] = Res [ i ] [ j ] . n ( ) Res [ i ] [ j ] = _quintic_simplify ( Res [ i ] [ j ] ) r1 = Res [ 1 ] [ 0 ] r1_n = Res_n [ 1 ] [ 0 ] for i in range ( 5 ) : if comp ( im ( r1_n * Res_n [ 4 ] [ i ] ) , 0 , tol ) : r4 = Res [ 4 ] [ i ] break u , v = quintic . uv ( theta , d ) sqrt5 = math . sqrt ( 5 ) # Now we have various Res values. Each will be a list of five # values. We have to pick one r value from those five for each Res u , v = quintic . uv ( theta , d ) testplus = ( u + v * delta * sqrt ( 5 ) ) . n ( ) testminus = ( u - v * delta * sqrt ( 5 ) ) . n ( ) # Evaluated numbers suffixed with _n # We will use evaluated numbers for calculation. Much faster. r4_n = r4 . n ( ) r2 = r3 = None for i in range ( 5 ) : r2temp_n = Res_n [ 2 ] [ i ] for j in range ( 5 ) : # Again storing away the exact number and using # evaluated numbers in computations r3temp_n = Res_n [ 3 ] [ j ] if ( comp ( r1_n * r2temp_n ** 2 + r4_n * r3temp_n ** 2 - testplus , 0 , tol ) and comp ( r3temp_n * r1_n ** 2 + r2temp_n * r4_n ** 2 - testminus , 0 , tol ) ) : r2 = Res [ 2 ] [ i ] r3 = Res [ 3 ] [ j ] break if r2 : break # Now, we have r's so we can get roots x1 = ( r1 + r2 + r3 + r4 ) / 5 x2 = ( r1 * zeta4 + r2 * zeta3 + r3 * zeta2 + r4 * zeta1 ) / 5 x3 = ( r1 * zeta3 + r2 * zeta1 + r3 * zeta4 + r4 * zeta2 ) / 5 x4 = ( r1 * zeta2 + r2 * zeta4 + r3 * zeta1 + r4 * zeta3 ) / 5 x5 = ( r1 * zeta1 + r2 * zeta2 + r3 * zeta3 + r4 * zeta4 ) / 5 result = [ x1 , x2 , x3 , x4 , x5 ] # Now check if solutions are distinct result_n = [ ] for root in result : result_n . append ( root . n ( 5 ) ) result_n = sorted ( result_n , key = default_sort_key ) prev_entry = None for r in result_n : if r == prev_entry : # Roots are identical. Abort. Return [] # and fall back to usual solve return [ ] prev_entry = r return result