(A, u, DE)
return (A, u)
Generate a system for the constant solutions.

    Given a differential field (K, D) with constant field C = Const(K), a Matrix
    A, and a vector (Matrix) u with coefficients in K, returns the tuple
    (B, v, s), where B is a Matrix with coefficients in C and v is a vector
    (Matrix) such that either v has coefficients in C, in which case s is True
    and the solutions in C of Ax == u are exactly all the solutions of Bx == v,
    or v has a non-constant coefficient, in which case s is False Ax == u has no
    constant solution.

    This algorithm is used both in solving parametric problems and in
    determining if an element a of K is a derivative of an element of K or the
    logarithmic derivative of a K-radical using the structure theorem approach.

    Because Poly does not play well with Matrix yet, this algorithm assumes that
    all matrix entries are Basic expressions.
Generate a system for the constant solutions.
Generate a system for the constant solutions .
def constant_system(A, u, DE):
    """
    Generate a system for the constant solutions.

    Given a differential field (K, D) with constant field C = Const(K), a Matrix
    A, and a vector (Matrix) u with coefficients in K, returns the tuple
    (B, v, s), where B is a Matrix with coefficients in C and v is a vector
    (Matrix) such that either v has coefficients in C, in which case s is True
    and the solutions in C of Ax == u are exactly all the solutions of Bx == v,
    or v has a non-constant coefficient, in which case s is False Ax == u has no
    constant solution.

    This algorithm is used both in solving parametric problems and in
    determining if an element a of K is a derivative of an element of K or the
    logarithmic derivative of a K-radical using the structure theorem approach.

    Because Poly does not play well with Matrix yet, this algorithm assumes that
    all matrix entries are Basic expressions.
    """
    if not A:
        return A, u
    Au = A.row_join(u)
    Au = Au.rref(simplify=cancel)[0]
    # Warning: This will NOT return correct results if cancel() cannot reduce
    # an identically zero expression to 0.  The danger is that we might
    # incorrectly prove that an integral is nonelementary (such as
    # risch_integrate(exp((sin(x)**2 + cos(x)**2 - 1)*x**2), x).
    # But this is a limitation in computer algebra in general, and implicit
    # in the correctness of the Risch Algorithm is the computability of the
    # constant field (actually, this same correctness problem exists in any
    # algorithm that uses rref()).
    #
    # We therefore limit ourselves to constant fields that are computable
    # via the cancel() function, in order to prevent a speed bottleneck from
    # calling some more complex simplification function (rational function
    # coefficients will fall into this class).  Furthermore, (I believe) this
    # problem will only crop up if the integral explicitly contains an
    # expression in the constant field that is identically zero, but cannot
    # be reduced to such by cancel().  Therefore, a careful user can avoid this
    # problem entirely by being careful with the sorts of expressions that
    # appear in his integrand in the variables other than the integration
    # variable (the structure theorems should be able to completely decide these
    # problems in the integration variable).

    Au = Au.applyfunc(cancel)
    A, u = Au[:, :-1], Au[:, -1]

    for j in range(A.cols):
        for i in range(A.rows):
            if A[i, j].has(*DE.T):
                # This assumes that const(F(t0, ..., tn) == const(K) == F
                Ri = A[i, :]
                # Rm+1; m = A.rows
                Rm1 = Ri.applyfunc(lambda x: derivation(x, DE, basic=True)/
                    derivation(A[i, j], DE, basic=True))
                Rm1 = Rm1.applyfunc(cancel)
                um1 = cancel(derivation(u[i], DE, basic=True)/
                    derivation(A[i, j], DE, basic=True))

                for s in range(A.rows):
                    # A[s, :] = A[s, :] - A[s, i]*A[:, m+1]
                    Asj = A[s, j]
                    A.row_op(s, lambda r, jj: cancel(r - Asj*Rm1[jj]))
                    # u[s] = u[s] - A[s, j]*u[m+1
                    u.row_op(s, lambda r, jj: cancel(r - Asj*um1))

                A = A.col_join(Rm1)
                u = u.col_join(Matrix([um1]))

    return (A, u)
def constant_system ( A , u , DE ) : if not A : return A , u Au = A . row_join ( u ) Au = Au . rref ( simplify = cancel ) [ 0 ] # Warning: This will NOT return correct results if cancel() cannot reduce # an identically zero expression to 0.  The danger is that we might # incorrectly prove that an integral is nonelementary (such as # risch_integrate(exp((sin(x)**2 + cos(x)**2 - 1)*x**2), x). # But this is a limitation in computer algebra in general, and implicit # in the correctness of the Risch Algorithm is the computability of the # constant field (actually, this same correctness problem exists in any # algorithm that uses rref()). # # We therefore limit ourselves to constant fields that are computable # via the cancel() function, in order to prevent a speed bottleneck from # calling some more complex simplification function (rational function # coefficients will fall into this class).  Furthermore, (I believe) this # problem will only crop up if the integral explicitly contains an # expression in the constant field that is identically zero, but cannot # be reduced to such by cancel().  Therefore, a careful user can avoid this # problem entirely by being careful with the sorts of expressions that # appear in his integrand in the variables other than the integration # variable (the structure theorems should be able to completely decide these # problems in the integration variable). Au = Au . applyfunc ( cancel ) A , u = Au [ : , : - 1 ] , Au [ : , - 1 ] for j in range ( A . cols ) : for i in range ( A . rows ) : if A [ i , j ] . has ( * DE . T ) : # This assumes that const(F(t0, ..., tn) == const(K) == F Ri = A [ i , : ] # Rm+1; m = A.rows Rm1 = Ri . applyfunc ( lambda x : derivation ( x , DE , basic = True ) / derivation ( A [ i , j ] , DE , basic = True ) ) Rm1 = Rm1 . applyfunc ( cancel ) um1 = cancel ( derivation ( u [ i ] , DE , basic = True ) / derivation ( A [ i , j ] , DE , basic = True ) ) for s in range ( A . rows ) : # A[s, :] = A[s, :] - A[s, i]*A[:, m+1] Asj = A [ s , j ] A . row_op ( s , lambda r , jj : cancel ( r - Asj * Rm1 [ jj ] ) ) # u[s] = u[s] - A[s, j]*u[m+1 u . row_op ( s , lambda r , jj : cancel ( r - Asj * um1 ) ) A = A . col_join ( Rm1 ) u = u . col_join ( Matrix ( [ um1 ] ) ) return ( A , u )