###############################################################################
#                                                                             #
#   RationalSum( f , x );                                                     #
#   Indefinite summation of rational functions                                #
#   Input:                                                                    #
#     f - rational function in x                                              #
#     x - name                                                                #
#   Output:                                                                   #
#     [ g , r ]                                                               #
#   where g and r are rational functions in x such that                       #
#     g(x+1) - g(x) + r = f                                                   #
#   and r has minimal degree of denominator                                   #
#   Optional arguments:                                                       #
#     'factorization' = false                                                 #
#       - full factorization of denominator of f should not be used           #
#     'method' = <integer>                                                    #
#       - specifies which summation algorithm should be used:                 #
#        0: factorization-based [Polyakov, 2011] (default one)                #
#        1: recursive [Abramov, 1975]                                         #
#        2: linear-algebraic [Abramov, 1995]                                  #
#        3: GGSZ [Gerhard et al, 2003]                                        #
#     'normalize'= true                                                       #
#       - for method 0, the output should be in normalized form               #
#     'minimize' = "sum denominator"                                          #
#       - g should have the denominator of minimal degree                     #
#     'minimize' = "numerator"                                                #
#       - r should have the numerator of minimal degree                       #
#         (not guaranteed)                                                    #
#     'minimized' = <name>                                                    #
#       - if 'minimize' = "numerator" option is used, <name> is               #
#         assigned true if the degree of numerator is guaranteed              #
#         to be minimal and false otherwise                                   #
#     'alignment' = "left"                                                    #
#       - for methods 0 and 1, remainder should be aligned to the             #
#         left side of f(x) (i.e. if p(x) is an irreducible factor of         #
#         the remainder's denominator, then it is also divides the            #
#         denominator of f(x), and p(x-k) does not divide the denominator     #
#         of f(x) for any positive integer k)                                 #
#     'alignment' = "right"                                                   #
#       - for methods 0 and 1, remainder should be aligned to the             #
#         right side of f(x) (i.e. if p(x) is an irreducible factor of        #
#         the remainder's denominator, then it is also divides the            #
#         denominator of f(x), and p(x+k) does not divide the denominator     #
#         of f(x) for any positive integer k)                                 #
#     'alignment' = "peak degree"                                             #
#       - for methods 0 and 1, remainder should be aligned to the             #
#         peak degrees of the denominator of f(x) [Pirastu, 1995]             #
#     'dispersion' = <nonnegative integer>                                    #
#       - for methods 1, 2, and 3, the given integer should be used           #
#         as a dispersion value (this may result in mistakes if the actual    #
#         dispersion is greater)                                              #
#                                                                             #
#                                          Stanislav Polyakov, April 2012     #
#                                                                             #
###############################################################################

RationalSum := module()
    local ModuleApply, Rational,          
          SumProperFraction, ShiftClasses, NonzeroDispersion, SortShifts,
          ArrangeNumerators, PartialFractionNumerators,
          SumComponent, DegreeL, ShiftResult, Rewind,
          ResultToFunction,
          Dispersion, DispersionF,
          Abramov75, APart,
          Abramov95, denominators,
          GGSZ,
          AdditionalMinimization, 
          MinimizeNum, NumSystem, NumSystemSolve, unsolvable,
          OptimizeSum, HandleApart, OptSFSum, GetRelevantPart, AreBalanced,
          factor_allowed, sum_method, disp, align, normal_result,
          add_min, degree_is_minimal, dim_report; 

ModuleApply := proc( F , x::name )
    local ARGS, S, opt, fact, meth, do_norm, rationalMinType,
          rem_align, min_option, DIM;

    if not type( F , ratpoly( anything , x ) ) then
        error "input must be a rational function in %1", x;
    end if;

    ( opt, ARGS ) := selectremove( z->type( z, name=anything ) , [args] );

# Factorization of the denominators
    factor_allowed := true;
    if hasoption( opt, 'factorization'=boolean, 'fact', 'opt' ) then
        factor_allowed := fact;
    end if;

# Summation algorithm: 0 - factorization-based straightforward,
#                      1 - Abramov75, 
#                      2 - Abramov95,
#                      3 - GGSZ
    if factor_allowed then
        sum_method := 0;
    else
        sum_method := 1;
    end if;
    if hasoption( opt, 'method'=integer, 'meth', 'opt' ) then
       if meth >= 0 and meth <= 3 then
            sum_method := meth;
        else
            error "no such method";
        end if;
    end if;
    if not factor_allowed and sum_method = 0 then
        error "method 0 (straightforward) requires factorization";
    end if;

# Normalization of the result
    normal_result := false;
    if hasoption( opt, 'normalize'=boolean, 'do_norm', 'opt' ) then
        normal_result := do_norm;
    end if;

# Dispersion
    disp := -1;
    if hasoption( opt, 'dispersion'=nonnegint, 'disp', 'opt' ) then
        if method = 0 then
            WARNING("dispersion is not used by method 0");
        else
            WARNING("%1 will be used as the dispersion of %2", disp , F);
        end if
    end if;

# Additional minimization
    rationalMinType := "-";
    if hasoption( opt, 'minimize'=string, 'min_option', 'opt' ) then
        rationalMinType := min_option;
    end if;
    if hasoption(opt,'minimize'=name,'min_option','opt') then
        rationalMinType := convert(min_option, string);
    end if;

    add_min := -1;
    if rationalMinType = "numerator" then
         add_min := 2;
    elif rationalMinType = "sum denominator" then
         add_min := 1;
    end if;

# Variable that should be assigned true if the numerator degree
# of the remainder is minimal, false otherwise
    dim_report := false;
    if hasoption( opt , 'minimized'=name, 'DIM', 'opt' ) then
        dim_report := true;
        degree_is_minimal := DIM;
    end if;

# Remainder alignment
    align := -1;
    if hasoption( opt, 'alignment'=string, 'rem_align', 'opt' ) then
        if add_min > -1 then
            WARNING("align option does not work with additional minimization");
        elif sum_method >= 2 then
            WARNING("this summation method does not support align option");
        elif rem_align = "left" then
            align := -1;
        elif rem_align = "peak degree" then
            align := 0;
        elif rem_align = "right" then
            align := 1;
        end if;
    end if;

    if nops(opt)>0 then
        error "wrong arguments"
    end if;

    Rational( F , x );

end proc:



    ############################################################################
    # Main summation procedure. Extracts polynomial part, calculates           #
    # disoersion (unless straightforward method is used for summation),        #
    # and calls the summation procedures.                                      #
    ############################################################################

    Rational := proc( f, x )
        local num, den, cont, p, c, s, r, ds, i, y;

        ( num, den ) := numer(f) , denom(f);
        # Handle the polynomial part
        p := 0; 
        if degree( expand(num) , x ) >= degree( expand(den) , x ) then
            num := rem( num , den , x , 'p' );
            p := `sum/polynom`( p , x );
        end if;
        num := primpart( num , x , 'cont' );

        if num = 0 or degree( den , x ) = 1 then 
            return [ p , cont*num/den ];
        end if;

        if sum_method = 0 then  #  apply straightforward method
            ( s , r ) := SumProperFraction( num , den , x );
            if normal_result then
                return [ Normalizer(p + cont*s) , cont*r ];
            end if; 
            return [ p + cont*s , cont*r ];
        end if;

        # Compute dispersion (unless it is provided) and dispersion set
        # (in the list form, values shifted by -1)
        if disp = -1 then
            if factor_allowed then
                ds := DispersionF( den , x );
            else
                ds := Dispersion( den , x );
            end if;
            disp := max( 0 , ds[] );
            ds := map( y -> y-1 , ds );
        elif sum_method > 1 then
            ds := [ seq( i , i = 0..disp-1 ) ];
        end if;

        if disp = 0 then
            if add_min =2 then
                ( s , r ) := AdditionalMinimization( 0 , num/den , x );
                return [ p + cont*s , cont*r ];
            else
                return [ p , cont*num/den ];
            end if;
        end if;

       
        if   sum_method = 1 then   # apply recursive Abramov's algorithm
            ( s , r ) := Abramov75( [ 0 , num/den ] , x , disp )[];
        elif sum_method = 2 then   # apply linear algebraic Abramov's algorithm 
            ( s , r ) := Abramov95( collect( num , x) , den , x , ds )[];
        elif sum_method = 3 then   # apply algorithm GGSZ
            ( s , r ) := GGSZ( num/den , x , ds )[];
        end if;
        if r <> 0 and add_min > -1 then # additional minimization
            ( s , r ) := AdditionalMinimization( s , r , x );
        end if;

        [ p + cont*s , cont*r ];

    end proc:


    SumProperFraction := proc( num , den , x::name )
        local f, c, C, S, d1, d0, n0, n1, c1, t, N, R, M, s, r, y;

    # factorizing
        f := factors( den ); 

    # removing the coefficient (anything not depending on x)
        C := f[1];            
        ( f , c ) := selectremove( y -> evalb( degree( y[1] , x ) > 0 ) , f[2] );
        C := C * convert( map( y -> y[1]^y[2] , c ) , `*` );

    # dividing the irreducible polynomials into shift-equivalence classes
    # (q_i(x+k) = q_j(x) for some integer k if and only if
    #  q_i and q_j belong to the same class)
        S := ShiftClasses( f , x );

    # stop if dispersion is zero
        if add_min <> 2 and { map( op , map2( op , 2 , S ))[] } = { 0 } then
            return 0 , num/den;
        end if;

    # extracting the part with zero dispersion
        if add_min <> 2 then
            S := NonzeroDispersion( S , x , 'd1' , 'd0' );
            gcdex( d1 , d0 , num , x , 'n0' , 'n1' );
            ( c1 , n1 ) := denom( n1 ) , numer ( n1 );
            C := C * c1;
            n0 := n0 * c1;
        else
            d1 := `*`( map( y -> `*`( zip(`^`,y[4],y[3])[] ) , S)[] );
            ( d0 , n0 , n1 ) := 1 , 0 , num;
        end if;

    # sorting the polynomials within shift classes according to their shifts
        S := map( SortShifts , S );

    # finding partial fraction decomposition
        N := PartialFractionNumerators( S , n1 , d1 , x ); # gcdex-based procedure

    # summation
        R := map( y -> SumComponent( y , ArrangeNumerators(y, N) , x ) , S );

        if add_min = 2 then   # Minimize remainder numerator
            M := MinimizeNum( map2( op , 2 , R ) , x );
            if dim_report then
                assign( degree_is_minimal , M[2] );
            end if;
            R := ShiftResult( R , M[1] , x );
        end if;

    # converting result to function format
        ( s , r ) := ResultToFunction( map2( op , 1 , R ) , map2( op , 2 , R ) , x );
        ( s , r ) := s/C , r/C + n0/d0/C;
        if normal_result then
            r := Normalizer(r);
        end if;
        s , r;

    end proc:


    ############################################################################
    #   Procedure splits the factorization ff of a polynomial into             #
    #   shift-equivalence classes.                                             #
    #   The shift-equivalence class is written in the form of a list with the  #
    #   following elements: irreducible (base) polynomial, list of shifts,     #
    #   list of corresponding degrees, list of shifted polynomials, and the    #
    #   second coefficient of the base polynomial.                             #
    #   Example: Let the factorization of                                      #
    #     p = x^8-x^2 = x^2*(x-1)*(x+1)*(x^2-x+1)*(x^2+x+1)                    #
    #   be written in the form                                                 #
    #     [ [x-1,1], [x^2-x+1,1],  [x,2], [x+1,1], [x^2+x+1,1] ].              #
    #   Then the procedure will return                                         #
    #     [ [ x-1 , [0,1,2] , [1,2,1] , [x-1,x,x+1] , -1 ] ,                   #
    #       [ x^2-x+1 , [0,1] , [1,1] , [x^2-x+1,x^2+x+1] , -1 ] ].            #
    ############################################################################

    ShiftClasses := proc( ff::list , x::name )
        local f , F , i , j , c , d , searching , alpha;

        f := map( y -> [ op(y) , coeff( y[1], x, degree(y[1],x)-1 ) / lcoeff( y[1], x ) ] , ff );
        F := [];
        for i from 1 to nops(f) do
            searching := true;
            d := degree( f[i][1] , x );
            for j from 1 to nops(F) while searching do
                if degree( F[j][1] , x ) = d then
                    alpha := radnormal( ( f[i][3] - F[j][5] ) / d );
                    if type( alpha , integer ) then
                        if Testzero( f[i][1] - eval( F[j][1] , x = x+alpha ) ) then
                            F[j][2] := [ op(F[j][2]) , alpha ];
                            F[j][3] := [ op(F[j][3]) , f[i][2] ];
                            F[j][4] := [ op(F[j][4]) , f[i][1] ];
                            searching := false;
                        end if;
                    end if;                
                end if;
            end do;
            if searching then
                F := [ F[] , [ f[i][1] , [ 0 ] , [ f[i][2] ] , [ f[i][1] ] , f[i][3] ] ] ;
            end if;
        end do;
        F;

    end proc:


    NonzeroDispersion := proc( L::list , x::name , d1 , d0 )
        local S, T, p;
        ( T , S ) := selectremove( y -> evalb( y[2] = [0] ) , L );
        d1 := `*`( map( y -> `*`( zip(`^`,y[4],y[3])[] ) , S)[] );
        d0 := `*`( map( y -> y[1]^y[3][1] , T )[] );
        S;
    end proc:



    SortShifts := proc( S::list )
        local s, m, k, y, z, t;

        s := S[1..4];
        m := min( s[2][] );
        if m <> 0 then
            member( m , s[2] , 'k' );
            s[1] := s[4][k];
            s[2] := map( y -> y-m , s[2] );
        end if;
        if sort( s[2] ) <> s[2] then
            t := [ seq( [ s[2][i] , s[3][i] , s[4][i] ] , i = 1..nops(s[2]) ) ];
            t := sort( t , ( y , z ) -> evalb( y[1] <= z[1] ) );
            s := [ s[1] , map2( op, 1, t ), map2( op, 2, t ), map2( op, 3, t ) ];
        end if;

        s;

    end proc:


    PartialFractionNumerators := proc( S , num , den , x )
        local dt , nt , c, N , i , j , u , v , ns , k , r;

        ( nt , dt ) := num , den;
        c := 1;
        N := [];
        for i from 1 to nops(S) do 
            for j from 1 to nops(S[i][4]) do
                dt := quo( dt , S[i][4][j]^S[i][3][j] , x);
                gcdex( S[i][4][j]^S[i][3][j] , dt , nt , x , 'u' , 'v' );
                v := v/c;
                u := u/c;
                c := denom(u); 
                nt := numer(u);
                ns := [];
                for k from 2 to S[i][3][j] do
                    r := rem( v , S[i][4][j] , x , 'v' );
                    ns := [ r , ns[] ];
                end do;
                N := [ N[] ,  [ S[i][4][j] , v , ns[] ] ];        
            end do;
        end do;
        N;
    end proc:


    ############################################################################
    #   Given a single shift-equivalence class of the denominator (in the form #
    #   of ShiftClasses output) and the list of numerators of partial          #
    #   fraction decomposition, the procedure chooses the numerators           #
    #   corresponding to the denominators of the shift-equivalence class       #
    #   and shifts them accordingly to denominator shifts.                     #
    ############################################################################

    ArrangeNumerators := proc( f::list , N::list )
        local q,L,i,k,d,j,y;

        q := map2(op,1,N);
        L := [];
        for i from 1 to nops(f[4]) do
            member( f[4][i] , q , 'k' );
            L := [ L[] , expand( eval( N[k][2..-1] , x = x - f[2][i] )) ];
        end do:
        d := max( f[3][] );
        map( y -> [ y[] , seq( 0 , j = nops(y)+1..d ) ] , L );   # even the lengths

    end proc:


    ############################################################################
    #   Sum decomposition of a rational function with a denominator belonging  #
    #   to a single shift-equivalence class (i.e. for any irreducible p,q that #
    #   divide this denominator p(x)=q(x+k) for some integer k).               #
    #   Minimization of the denominator of summed part ("rewind") is called    #
    #   from this procedure.                                                   #
    #   Input: f - list containing the basic polynomial, list of its shifts,   #
    #              list of corresponding degrees, list of shifted polynomials  #
    #          N - list of corresponding numerators                            #
    #          x - summation variable                                          #
    ############################################################################

    SumComponent := proc( f::list , N::list , x::name )
        local s , d , l , r , deg_l , deg_r , num_l , num_r , h , i , k , j , u;

        ( s , d ) := [] , nops( N[1] );
        ( l , r ) := 1 , nops( f[2] );
        if l = r then
            return [ [] , add( N[1][-i]*f[1]^(i-1) , i = 1..d ) / f[1]^d ];
        end if;
        ( deg_l , deg_r ) := f[3][l] , f[3][r];
        ( num_l , num_r ) := N[l] , N[r];

        while r > l do
            if align < 0 or ( align = 0 and deg_l >= deg_r ) then   # move left
                s := [ [ num_r , f[2][r-1]..f[2][r]-1 ] , s[] ];
                if r-1 = l then
                    h := [ [ seq( num_r[i] + num_l[i] , i = 1..d ) ] , f[2][l] , f[4][l] ];
                else
                    num_r := [ seq( num_r[i] + N[r-1][i] , i = 1..d ) ];
                    deg_r := DegreeL(num_r);
                end if;
                r := r-1;
            else    # move right
                s := [ [ -num_l , f[2][l]..f[2][l+1]-1 ] , s[] ];
                if l+1 = r then
                    h := [ [ seq( num_r[i] + num_l[i] , i = 1..d ) ] , f[2][r] , f[4][r] ];
                else
                    num_l := [ seq( num_l[i] + N[l+1][i] , i = 1..d ) ];
                    deg_l := DegreeL(num_l);
                end if;
                l := l+1;
            end if;
        end do;

        if add_min = 1 then
            ( s , h ) := Rewind( h , s , x );
        end if;

        k := map( y -> DegreeL( y[1] ) , s );
        u := [ seq( add( s[j][1][-i+k[j]-d]*f[1]^(i-1) , i = 1..k[j] ) , j = 1..nops(k) ) ];
        # [ numerator, denominator, denominator degree, shifts ]
        s := [ seq( [ u[j] , f[1] , k[j] , s[j][2] ] , j = 1..nops(s) ) ];
        k := DegreeL( h[1] );
        h[1] := add( h[1][-i+k-d]*f[1]^(i-1) , i = 1..k );
        [ s , expand( eval( h[1] , x = x+h[2] )) / h[3]^k ];

    end proc:


    DegreeL := proc( L::list )
        local n;
        n := nops( L );
        while n > 0 and L[n] = 0 do
            n := n-1;
        end do;
        n;
    end proc:


    ############################################################################
    #   Procedure finds the sum decomposition of a rational function with a    #
    #   minimal degree of the denominator of summed part.                      #
    #   The function denominator factors must belong to a single               #
    #   shift-equivalence class.                                               #
    #   Input: R - remainder in the left-aligned decomposition                 #
    #          s - list of the sum components in the form                      #
    #              [ [ g1 , a1..b1 ] , ... , [ gk , ak..bk ] ]                 #
    #              ( where G = g1(x+a1)+g(x+a1+1)+...+g1(x+b1)+...+gk(x+bk) )  #
    ############################################################################
    
    Rewind := proc( R , s::list , x::name )
        local L, m, v, p, w, k, g, d;

        ( L, m, v, p, w ) := [] , 0 , 0 , 0 , 0;
        for k from 1 to nops(s) do
            g := [ seq( s[k][1][i] - R[1][i] , i = 1..nops( R[1] ) ) ];
            d := DegreeL( g ) - DegreeL( s[k][1] );
            w := w + d * ( rhs(s[k][2]) - lhs(s[k][2]) + 1 );
            if w < v then
                v := w;
                m := k;
                p := rhs(s[k][2])+1;
            end if;
            L := [ L[] , [ g , s[k][2] ] ];
        end do;

        [ L[1..m][] , s[m+1..nops(s)][] ] , expand( eval( R , x = x+p ) );

    end proc:


    ResultToFunction := proc( S , R , x )
        local s , j , SSRF , SumSeries , f , p , q , u , v;

        s := map( op , S );
        if not normal_result then
            s := map( y -> add( eval( y[1]/y[2]^y[3] , x=x+j ) , j=y[4] ) , s );
            return convert( s , `+` ) , convert( R , `+` );
        end if;

        # normalizing result
        if add_min = 2 then
            s := map( y -> add( eval( y[1]/y[2]^y[3] , x=x+j ) , j=y[4] ) , s );
            return Normalizer( convert( s , `+` ) ) , Normalizer( convert( R , `+` ) );
        end if;

        SSRF := proc( P , Q )
            local Ps,Qs,Ns;
            if nops(P) > 3 then
                Ps := [ ListTools:-LengthSplit( P , 2 ) ];
                Qs := [ ListTools:-LengthSplit( Q , 2 ) ];
                Ns := zip( ( y , z ) -> [ SSRF(y,z) ] , Ps , Qs );
                return SSRF( map2( op , 1 , Ns ) , map2( op , 2 , Ns ) );
            elif nops(P) = 3 then
                return expand( P[1]*Q[2]*Q[3] + Q[1]*P[2]*Q[3] + Q[1]*Q[2]*P[3] ) , Q[1]*Q[2]*Q[3];
            elif nops(P) = 2 then
                return expand( P[1]*Q[2] + P[2]*Q[1] ) , Q[1]*Q[2];
            end if;
            P[1] , Q[1];
        end proc;

        SumSeries := proc( L , x )
            local dens, nums, k;
            dens := [ seq( expand( eval( L[2] , x = x+k ) )^L[3] , k=L[4] ) ];
            nums := [ seq( eval( L[1] , x = x+k ) , k=L[4] ) ];
            [ SSRF( nums , dens ) ];
        end proc;

        f := map( SumSeries , s , x );
        ( p , q ) := SSRF( map2( op , 1 , f ) , map2( op , 2 , f ) );
        ( u , v ) := SSRF( map( numer , R ) , map( denom , R ) );
        p/q , u/v;

    end proc:


    AdditionalMinimization := proc( S , R , x::name )
        local s, r, f, rl, t, M, i, j, sl, y;

        if add_min = 1 then
            ( s , r ) := OptimizeSum( S , R , x )[];
            return s , r;
        end if; # combined minimization can be added here
 
        r := radnormal(R);
        f := factors( denom(r) );
        rl := convert( [ numer(r)/f[1] , f[2][] ] , parfrac , x )[2..-1];
        rl := map( y -> add( y[-t]*y[1]^(t-1) , t=1..nops(y)-1 )/y[1]^(nops(y)-1) , rl );
        M := MinimizeNum( rl , x );
        if dim_report then
            assign( degree_is_minimal , M[2] );
        end if;
        rl := [ seq( [ rl[i] , M[1][i] ] , i = 1..nops(rl) ) ];
        sl := map( y -> -sign(y[2])*add(eval(y[1], x=x+t) , t = min(0,y[2])..max(-1,y[2]-1)), rl );
        r := convert( map( y -> eval( y[1] , x = x+y[2] ) , rl ) , `+` );

        S + convert( sl , `+` ) , r;
         
    end proc:


    #################################################################
    # Finds the remainder with the minimal degree of the numerator  #
    # and the necessary shift values                                #
    #################################################################

    MinimizeNum := proc( fractions::list , x::name )
        local n, degs, min_dd, lcs, i, p, k, vars, t, y;
 
        n := nops(fractions);
        if n = 1 or select( type , indets(fractions) , name ) minus {x} <> {} then  
            return [ [ seq( 0 , p = 1..n ) ] , evalb(n=1) ];
        end if;

        degs := map( y -> degree( numer(y), x ) - degree( denom(y), x ) , fractions );
        min_dd := max( degs[] );
        lcs := map( y -> lcoeff( numer(y), x ) / lcoeff( denom(y), x ) , fractions );
        lcs := [ seq( [ min_dd - degs[i] , lcs[i] ] , i = 1..n ) ];
        lcs := map( y -> if y[1] = 0 then y[2] else 0 end if , lcs );
        if `+`( lcs[] ) <> 0 then
            return [ [ seq( 0 , p = 1..n ) ] , true ];
        end if;

        vars := [ seq( k[p] , p = 1..n ) ];
        member( min_dd , degs , 't' );
        vars[t] := 0;

        NumSystem( fractions , vars , n , x );

    end proc:


    NumSystem := proc( fractions::list , vars::list , n::posint , x::name )
        local s, num, d, i, j, p, S, sol, w, m, z, d_i_m, v, y;

        s := [ seq( eval( fractions[p] , x = x+vars[p] ) , p = 1..n ) ];
        num := expand( numer( normal( convert( s , `+`))));
        d := degree( num , x );
        S := [ seq( coeff( num, x, d-i ) , i = 0..d-1 ) ];

        if nops(S) > 0 and degree( S[1] ) = 1 then
            sol := isolve( {S[1]} );
            if sol = NULL or sol = {} then
                return [ [ seq( 0 , p = 1..n ) ] , true ];
            end if;
            if type( sol[1] , set ) then
                sol := sol[1];
            end if;
            w := [ op( indets( sol ) minus { vars[] } ) ];
            w := [ seq( w[j] = m[j] , j = 1..nops(w) ) ];
            sol := eval( sol , w );
            S := map( expand , eval( S , sol ) );
            S := remove( y -> evalb( y=0 ) , S );
            z := NumSystemSolve( S , 'd_i_m' );
            sol := [ z[] ,  eval( sol , z )[] ];
        else
            sol := NumSystemSolve( S , 'd_i_m' );
        end if;

        v := eval( vars , sol );
        [ eval( v , map( y -> y=0 , indets(v) )) , evalb( d_i_m or n<=3 ) ];

    end proc:


    NumSystemSolve := proc( S::list , d_i_m )
        local vars , i , p , previous , lsol , sol;

       if nops(S) = 0 then
            return {};
        end if;
        vars := indets(S);
        ( i , previous , lsol ) := 1 , true, [ seq( vars[p]=0 , p = 1..nops(vars) ) ];
        while i <= nops(S) do
            sol := isolve( { S[1..i][] } );
            if sol = NULL or sol = {} then
                if unsolvable( { S[1..i][] } ) then 
                    i := nops(S);
                else 
                    previous := false;
                end if;
            else
                lsol := sol;
                previous := true;
            end if;
            i := i+1;
        end do;
        if nops(lsol) > 0 and type( lsol[1] , set ) then
            lsol := lsol[1];
        end if;

        d_i_m := previous;
        lsol;
        
    end proc:


    #####################################################################
    # Procedure tries to define if the system has no integer solutions. #
    # Rewritten (Aug 2008)                                              #
    #####################################################################

    unsolvable := proc(system)
        local S , L , sol , monom , i , alpha , sols;

        S := { map( expand , system )[] };
        monom := map( x -> if type(x,`+`) then op(x); else x; end if , S );
        monom := map( x -> if degree(x)>1 then x / lcoeff(x); end if , monom );
        monom := sort( [ monom[] ] , (x,y) -> evalb( degree(x) > degree(y) ) );

        for i from 1 to nops(monom) do
            S := eval(S,monom[i]=alpha[i]);
        end do;
        sols := isolve(S);

        evalb(sols = NULL or sols = {});

    end proc:


    ShiftResult := proc( R::list , s::list , x::name )
        local ShiftG, L, i, y;
        ShiftG := proc( g , r , k )
            [ g[] , [ -sign(k) * r , 1 , 1 , min(0,k)..max(-1,k-1) ] ];
        end proc;
        L := [ seq( [ R[i][] , s[i] ] , i = 1..nops(s) ) ];
        return map( y -> [ ShiftG( y[] ) , eval( y[2] , x=x+y[3] ) ] , L );
    end proc:


    ##########################  Other summation algorithms  ##########################


    Dispersion := proc( p , x )
        local f , k , r;
 
        f := sqrfree( p , x )[2];
        f := convert( map2( op , 1 , f ) , `*` );
        r := roots( resultant( eval( f , x = x+k ) , f , x ) , k );
        sort( select( type , map2( op , 1 , r ) , posint ) );

    end proc:


    DispersionF := proc( p , x )
        local f, F, s, i, j, y;

        f := factors( p ); 
        f := select( y -> evalb( degree( y[1] , x ) > 0 ) , f[2] );
        F := ShiftClasses( f , x );
        s := map2( op , 2 , F );
        s := map( y -> seq( seq( abs(y[j]-y[i]) , j=i+1..nops(y) ) , i=1..nops(y)-1 ) , s );
        { s[] };        

    end proc:


    Abramov75 := proc( F::list , x::name , dis::nonnegint )
        local g, f, p, q, cp, vp, wp, vm, wm, a, b, u, r;

        ( g , f ) := op(F);   
        if f = 0 or dis = 0 then
            return F;
        end if;

        ( p , q ) := numer(f), denom(f);
        cp := gcd( q , eval( q , x = x+dis ) );
        if cp = 1 then
            return Abramov75( F , x , dis-1 );
        end if;

        APart( q , cp , 'vp' , 'wp' , x );
        APart( q , eval( cp , x = x-dis ) , 'vm' , 'wm' , x );
        if align = -1 or ( align = 0 and degree( vm , x) > degree( vp , x) ) then
            gcdex( vp, wp, p, x, 'b', 'a' );
            u := eval( a/vp , x = x-1 );
            r := radnormal( b/wp + u );
        else
            gcdex( vm, wm, p, x, 'b', 'a' );
            u := -a/vm;
            r := radnormal( b/wm + eval( a/vm , x = x+1 ) );
        end if;

        Abramov75( [ g+u , r ] , x , dis-1 );

    end proc:

    APart := proc( p , h , v , w , x )
        local v1, w1, g, q;
        g := gcd( p , h );
        ( v1 , w1 ) := 1 , p;
        while g <> 1 do
            q := radnormal( w1 / g );
            if type( q , polynom( anything , x ) ) then
                w1 := q;
            end if;
            v1 := v1 * g;
            g := gcd( w1 , v1 );
        end do;
        ( v , w ) := v1 , w1;
    end proc:


    Abramov95 := proc( num , den , x::name , ds )
        local p1, p2, i, j, d1, d2, n1, n2, deg1, deg2, f, f1, f2,
              constnumb, a, lmatr, lrside, work , res;

        ( d1 , d2 ) := denominators( den, x, ds );
        ( deg1, deg2 ) := degree( d1 , x ) - 1 , degree( d2 , x ) - 1;
        p1 := `+`( seq( a[i+1]*x^i , i = 0..deg1 ) );
        p2 := `+`( seq( a[deg1+i+2]*x^i , i = 0..deg2 ) );
        f := normal( eval( p1/d1 , x = x+1 ) - p1/d1 + p2/d2 );
        ( f1 , f2 ) := numer(f) , denom(f);
        n2 := quo( f2 , den , x );
        n1 := collect( num*n2 , x );
        deg1 := max( degree( f1, x ) , degree( d1, x ) );

        constnumb := deg1 + deg2 + 2;
        lmatr := matrix( 1+deg1 , constnumb );
        lrside := vector( 1+deg1 );
        for i from 0 to deg1 do
            lrside[i+1] := coeff( n1 , x , i );
            work := coeff( f1 , x , i );
            for j to constnumb do
                lmatr[ i+1 , j ] := coeff( work , a[j] , 1 );
            end do;
        end do;
        a := linalg['linsolve']( lmatr , lrside );
        [ normal( eval( p1/d1 )) , normal( eval( p2/d2 )) ];

    end proc:


    denominators := proc( den_f , x::name , ds )
        local j, i, b, u, g, h, lh, d, s, k, m, a;

        a := eval( den_f , x = x-1 );
        ( u , b , g ) := 1 , den_f , den_f;
        lh := sort( [ ds[] ] );
        m := nops(lh);
        lh := [ seq( lh[-i] , i = 1 .. m ) ];
        h := Array(1 .. m, lh);
        d := Array(1 .. m, [ seq( gcd( a , eval(b, x = x+i) ), i = lh ) ]);
        s := Array(1 .. m);
        k := 0;
        while k < m do 
            k := k+1; 
            s[k] := mul( eval( d[k] , x = x-i ) , i = 0 .. h[k] ); 
            u := u*s[k]; 
            g := quo( g , gcd( g, eval( d[k] , x = x+1 )*s[k] ) , x ); 
            j := k; 
            for i from k+1 to m do 
                d[i] := quo( d[i] , gcd( d[i], s[k] ) , x);
                if degree( d[i] , x ) > 0 then 
                    j := j+1;
                    d[j] := d[i];
                    h[j] := h[i];
                end if 
            end do;
            m := j;
        end do; 
        for k from m by -1 to 1 do
            g := `sum/ratsum_nonrat`( g, d[k], h[k], s[k], x );
        end do;
        u, g;
    end proc:



    ##########################################################
    # Input : rational functions G, R in x s.t. Dis(R(x))=0  #
    # Output : rational functions G1, R1, satifsying         #
    #  (i) G1(x+1)-G1(x)+R1(x) = G(x+1)-G(x)+R(x),           #
    #  (ii) degree(denom(R1)) = degree(denom(R)), and        #
    #  (iii) G1 has the minimal degree of denominator.       #
    ##########################################################

    OptimizeSum := proc( G , R , x::name )
        local N, sf, S;

        if R = 0 or G = 0 then
          return [ G , R ];
        end if;

        N := radnormal(R); 
        sf := sqrfree(denom(N));
        S := HandleApart( G, numer(N)/sf[1] , sf[2] , x, 0, 0, 0, 0);

        return [ S[1] + G , S[2] ];

    end proc:


    ############################################################################
    # Input : P - numerator of the remainder                                   #
    #         L - denominator in the [ [q1,d1] , ... , [qn,dn] ] form          #
    ############################################################################

    HandleApart := proc( G , P , L::list , x::name , pos , v , Pos , V )
        local H, g, r, p, n, i, q, Q, s, t, opt, y;

        ( H , g , r , p , n ) := G , 0 , 0 , P , nops(L);
        for i from 1 to n do
            q := L[i][1]^L[i][2];
            Q := convert( map( y -> y[1]^y[2] , L[i+1..n] ) , `*` );
            gcdex( Q , q , p , x , 's' , 't' );
            p := t;
            opt := OptSFSum( H, s, L[i][1], L[i][2], x, pos, v, Pos, V );
            ( g , r , H ) := g + opt[1] , r + opt[2] , opt[3];
        end do;

        return [ g , r , H ];

    end proc:


    OptSFSum := proc( G , N , q , d , x , pos , v , Pos , V )
        local R, H, k, K, c, b, dir, dd, g, h, f, i, y;

        if N = 0 or degree( denom(G) , x ) = 0 then
            return [ 0 , N/q^d , G ];
        end if;
 
        R := N/q^d;
        H := G;
        ( k , K , c , b ) := pos , Pos , v , V;
        dir := 1;
        dd := 1;
        while dd > 0 and abs(k) <= disp do
            if dir < 0 then
                k := k - 1;
            end if;

            ( g , h ) := GetRelevantPart( H , eval( q^d , x = x+k ) , x );
            if g = 0 then
                c := c + d;
            else
                g := eval( g , x = x-k );
                if AreBalanced( g , radnormal( g - dir*R ) , q , x , 'f' ) then
                    c := c + f;
                else
                    f := map( y -> [ y , d ] , f );
                    return HandleApart( H, N, f, x, k, c, K, b );
                end if;
            end if;
            H := h;

            if dir > 0 then
                k := k + 1;
            end if;
            if c < b then
                ( K , b ) := k , c;
            end if;
            dd := degree( denom(H) , x ) - degree( q , x ) * ( c - b );
            if dir > 0 and ( dd <= 0 or k > disp ) and align > -1 then
                ( dir , k , c , dd ) := -1 , 0 , 0 , 1;                         
            end if;
        end do;

        if K >= 0 then
            [ -add( eval( R , x = x+i ), i = 0..K-1 ) , eval( R , x = x+K ) , H ];
        else
            [ add( eval( R , x = x+i ), i = K..-1 )  ,  eval( R , x = x+K ) , H ];
        end if;

    end proc:


    GetRelevantPart := proc( R , p , x::name )
        local P, Q, q, r, s, t;

        if type( R , `+` ) and not type( R , polynom( anything , x ) ) then
            return map( GetRelevantPart , R , p , x );             
        end if; 
        ( P , Q ) := numer(R) , denom(R);
        ( q , r , t ) := Q , 1 , 0;
        while degree(q) > 0 and expand( t - p ) <> 0 do
            r := r * gcd( q , p , 's' , 't' );
            q := s;
        end do:
        Q := quo( Q , r , x );
        gcdex( Q , r , P , x , 's' , 't' );

        s/r , t/Q;  

    end proc:



    ############################################################################
    # Returns true if the ratio of the denominators of B and A is an integer   #
    # power of P. In this case, res is assigned that power, otherwise it is    # 
    # assigned the decomposition of P w.r.t. degrees of its factors in B/A.    #
    ############################################################################

    AreBalanced := proc( A , B , P , x::name , res )
        local R, r, d, p;

        R := radnormal( denom(B) / denom(A) );
        r := sqrfree( R , x )[2];
        d := map2( op , 2 , r );
        r := map2( op , 1 , r );
        p := convert( r , `*` );
        if p <> P then
            if degree( p , x ) < degree( P , x ) then
                r := [ r[] , quo( P , p , x ) ];
            elif nops(r) > 1 then
                r[1] := r[1] * lcoeff(P,x)/lcoeff(p,x);
            end if;
        end if;

        if nops(r) = 1 then
            res := [ d[] , 0 ][1];
        else
            res := r;
        end if;
        evalb( nops(r) = 1 );

    end proc:


    ############################################################################
    # GGSZ algorithm                                                           #
    ############################################################################

GGSZ := module()
export ModuleApply;
local GGSZnormal, ShiftlessDecomposition, GcdFreeBasisImplementation, PolyRefine;
option `Copyright (c) 2003 by Maplesoft. All rights reserved.`;
   
   
#
# Compute an additive decomposition of the rational function F
# w.r.t. x, i.e., rational functions R,H satisfying
#       F = (E-1)@R + H and autodispersion(denom(H))=0,
# where E is the shift operator w.r.t. x and
# @ denotes operator application
# returns the pair [R,H]
#
# The advantage of this algorithm over the one implemented in
# SumTools:-Hypergeometric:-SumDecomposition for rational functions
# is that it is much faster when the dispersion is large but
# the output is small.
#
# Reference: Gerhard, Giesbrecht, Storjohann & Zima (2003),
#            Shiftless factorization and polynomial-time rational summation,
#            Proceedings ISSAC 2003, 119-126.
#
# Created:   November 2003
# Authors: 	 Juergen Gerhard, Mark Giesbrecht, Arne Storjohann & Eugene Zima
# Contact:   jgerhard@maplesoft.com
#
# Modified by Stanislav Polyakov, January 2012
# (now uses previously computed dispersion set)
#

ModuleApply := proc(F,x,ds)
local f,g,c,SLD,gg,z,r,R,H,offset,i,gi,hi,n,ri,k,t,j,S;
option `Copyright (c) 2003 by Maplesoft. All rights reserved.`;
   
   ASSERT(x::'name');
   ASSERT(F::'ratpoly'('anything',x),
          "expect a rational function as input"):
   f,g := (numer,denom)(F);
   
   # compute shiftless decomposition g = c*prod_ij gi(x+hij)^eij
   c,SLD := op(ShiftlessDecomposition(g,x,ds)):
   
   # compute partial fraction decomposition
   # F = h + sum_ijk rijk/gi(x+hij)^k
   # gg := linear list of pairs [gi(x+hij),eij]
   gg := map(z->seq([eval(z[1],x=x+op([2,j,1],z)),op([2,j,2],z)],
                    j=1..nops(z[2])),SLD):
   # This is Allan's cool new programmer level interface :-)
   r := convert([f/c,op(gg)],'parfrac',x):
   
   # initialize rational/nonrational part of the antidifference
   R,H := SumTools:-IndefiniteSum:-Polynomial(r[1],x),0:

   # remove polynomial part and denominators from r
   r := map2(subsop,1=NULL,subsop(1=NULL,r)):
   offset := 1: # offset into r
   
   # for each shift class ...
   for i from 1 to nops(SLD) do
      # rewrite PFD in terms of linear difference operators Mik:
      # sum_jk rijk(x)/gi(x+hij)^k = sum_k Mik @ 1/gi(x)^k,
      # with Mik = sum_j rijk E^hij corresponding to column k
      # of the matrix ri = (rijk)_jk
      gi,hi := op(SLD[i]):
      hi := map2(op,1,hi): # list of all shifts
      n := nops(hi): # number of shifts = number of rows of ri
      ri := Matrix(r[offset..offset+n-1]):
      offset := offset+n:
      
      # partial summation: antidifference(Mik @ 1/gi^k)
      # = (LQ(Mik,E-1) @ 1/gi^k) + antidifference(sum_j rijk(x-hij)/gi^k)
      for k from 1 to LinearAlgebra:-ColumnDimension(ri) do
         # Left division of Mik = sum_j E^sj tj by E-1,
         # where sj = hij and tj = rijk(x-hij), yields:
         #
         # Mik = (tn + t(n-1) + ... + t1 + t0) +
         #       (E-1)[(E^(sn-1)    +...+E^s(n-1))  tn +
         #             (E^(s(n-1)-1)+...+E^s(n-2)) (tn+t(n-1)) + 
         #             ... +
         #             (E^(s1-1)    +...+E^s0)     (tn+t(n-1)+...+t1)]
         t := Normalizer(eval(ri[n,k],x=x-hi[n])):
         for j from n-1 to 1 by -1 do
            if not Testzero(t) then
               # since the gi are pairwise shift coprime, all
               # denominators are coprime and we can use GGSZnormal
               # below which is much faster than Normalizer
               S := GGSZnormal(add(eval(t/gi^k,x=x+z), z=hi[j]..hi[j+1]-1), x):
               R := R + S[1]/S[3];
            end if:
            t := Normalizer(t + eval(ri[j,k],x=x-hi[j])):
            # invariant: t = tn + ... + tj
         end do:
         H := Normalizer(H + t/gi^k):
      end do:
   end do:
   R := Normalizer(R):
   
   ASSERT([R,H]::['ratpoly'('anything',x),'ratpoly'('anything',x)]);
   [R,H]
end proc;
   

# like Normalizer, but assumes that all denominators are
# pairwise coprime and returns three operands
# Normalizer(numer(l)),Normalizer(denom(l)),denom(l)
# (the first two operands are in expanded form, the
#  third operand is in factored form)
# This is much more efficient than Normalizer for many summands
GGSZnormal := proc(l,x)
local n,r,F,G;
option `Copyright (c) 2003 by Maplesoft. All rights reserved.`;
   n := nops(l):
   # "10" is a magic number coming from some experiments
   if not type(l,`+`) or n<10 then
      r := Normalizer(l):
      Normalizer(collect(numer(r),x)),
      Normalizer(collect(denom(r),x)),denom(r)
   else
      # divide and conquer
      F := GGSZnormal(`+`(op(1..floor(n/2),l)),x):
      G := GGSZnormal(`+`(op(1+floor(n/2)..n,l)),x):
      Normalizer(collect(F[1]*G[2]+G[1]*F[2],x)),
      Normalizer(collect(F[2]*G[2],x)),F[3]*G[3]
   end if
end proc:


# modified version of PolynomialTools:-ShiftlessDecomposition
# with dispersion set (ds) provided as input

ShiftlessDecomposition := proc (f, x::name , ds )::list; 
    local B, B1, C, g, n, i, k, ff, c, A, h, z, f_i, e, gg, he_pairs, classes, multiplicity, shifts;

    if not type(f,('polynom')('anything',x)) then 
        error "expecting a polynomial as 1st argument" 
    end if; 
    ff := collect(Normalizer(f),x); 
    if degree(ff,x) <= 0 then 
        return [ff, []] 
    end if; 
    c, g := op(sqrfree(ff,x)); 
    g := sort(g,(a, b) -> a[2] <= b[2]); 
    e := map2(op,2,g); 
    g := map2(op,1,g); 
    gg := convert(g,`*`); 
    k := nops(g); 
    A := map(`+`,ds,1); 
    B := convert(g,'set'); 
    for h in A while true do
        B := `union`(`union`(B,{seq(gcd(gg,PolynomialTools:-Translate(g[i],x,h)),i = 1 .. k)}),
                     {seq(gcd(gg,PolynomialTools:-Translate(g[i],x,-h)),i = 1 .. k)});
        B := GcdFreeBasisImplementation(B);
    end do;
    B := convert(map(collect,B,x),list); 
    i := 0; 
    classes := table(); 
    B := map(z -> [FAIL, z],B);  
    while 0 < nops(B) do 
        B1 := B[1][2]; 
        B := map(z -> [PolynomialTools:-ShiftEquivalent(B1,z[2],x,'integer'), z[2]],B); 
        C, B := selectremove(z -> type(z[1],'integer'),B); 
        shifts := sort(map2(op,1,C)); 
        h := shifts[1]; 
        f_i := op([1, 2],select(z -> z[1] = h,C)); 
        i := i+1; 
        classes[i] := [f_i, map(`-`,shifts,h)];
    end do;
    n := i; 
    g := Array(1 .. nops(g),g);

    multiplicity := proc (f_i, h) 
        local j, `fi(x+h)`; 
        `fi(x+h)` := PolynomialTools:-Translate(f_i,x,h);
        for j to k do
           if divide(g[j],`fi(x+h)`,('g')[j]) then 
               return [h, e[j]];
           end if;
        end do;
        ASSERT(false,"should not happen: fi(x+h) does not divide the input"); 
    end proc;

    c := lcoeff(ff,x); 
    for i to n do 
        C := classes[i]; 
        he_pairs := map2(multiplicity,op(C));
        c := Normalizer(c/(lcoeff(C[1],x)^convert(map2(op,2,he_pairs),`+`))); 
        classes[i] := [C[1], he_pairs];
    end do; 
    return [c, sort(convert(classes,list),(a, b) -> degree(a[1],x) <= degree(b[1],x))];

end proc:


GcdFreeBasisImplementation := proc(S)
   local NonConstant, B, l, M, r, A, i;

   NonConstant := proc(a)
      evalb(0 < degree(a))
   end proc:


   ASSERT(S::('set(polynom)'));
   if nops(S) < 2 then
      B := S
   elif nops(S) = 2 then
      l, M, r := PolyRefine(op(S));
      B := `union`(M,{r, l})
   else
      A := GcdFreeBasisImplementation(subsop(1 = NULL,S));
      l := S[1];
      B := {};
      for i to nops(A) do
         l, M, r := PolyRefine(l,A[i]);
         B := `union`(`union`(B,M),{r})
      end do;
      B := `union`(B,{l})
   end if;
   select(NonConstant,B);
end proc:


PolyRefine := proc(a, b)
   local g, ag, bg, l1, l2, r1, r2, M1, M2;
   ASSERT(a::('polynom') and b::('polynom'));
   g := gcd(a,b,ag,bg);
   if degree(g) = 0 then
      a, {}, b
   else
      l1, M1, r1 := PolyRefine(ag,g);
      l2, M2, r2 := PolyRefine(r1,bg);
      l1, `union`(`union`(M1,M2),{l2}), r2
   end if
end proc:

end module:


end module:

####################################################################################################