unprotect('BeautyClinic');

module BeautyClinic()

option
    package;

export
    AddThenBeautify,
    Beautify;

local
    compareLexicographically,
    computeResult,
    dissect,
    mergeRows,
    preprocessAssumptions,
    processOneSum,
    processOneSumType,
    processTermsMultivariate,
    processTermsUnivariate,
    roundFloat,
    roundRational,
    sortRows;


    processOneSumType := {':-`+`', ':-`*`', ':-anything'^':-rational', ':-Not'(':-`+`')^':-numeric'};

    Beautify := proc(expr,
                     digits :: posint := Digits,
                     assumptions :: list({`<`, `<=`}) := [],
                     {output :: {identical(float), identical(rational)} := ':-float'},
                     {underflow :: truefalse := true},
                     {process_exact :: truefalse := false},
                     {round :: truefalse := true},
                     {beautify_top_level :: truefalse := true},
                     $)
    local
        Options,
        result;
        
        Options := op(remove(x -> lhs(x) = ':-beautify_top_level', [_options]));
        result := evalindets(expr, ':-`+`',
                             sum -> frontend(processOneSum, [sum, 0], 
                                             [processOneSumType, {}], digits, assumptions, Options));
        
        # For floats occurring outside of any sums:
        if round and beautify_top_level then
            return subsindets(result, ':-float',
                              `if`(output = ':-float',
                                   x -> roundFloat(digits, x),
                                   x -> roundRational(digits, x)));
                              
        else
            return result;
        end if;
    end proc;

    AddThenBeautify := proc(expr_1 :: algebraic,
                            expr_2 :: algebraic,
                            digits :: posint := Digits,
                            assumptions :: list({`<`, `<=`}) := [],
                            {output :: {identical(float), identical(rational)} := ':-float'},
                            {underflow :: truefalse := true},
                            {process_exact :: truefalse := false},
                            {round :: truefalse := true},
                            $)
    local
        sum_1,
        sum_2;
        
        if type(expr_1, ':-`+`') then
            sum_1 := `+`(op(Beautify([op(expr_1)], digits, assumptions, _options, ':-beautify_top_level' = false)));
        else
            sum_1 := Beautify(expr_1, digits, assumptions, _options, ':-beautify_top_level' = false);
        end if;

        if type(expr_2, ':-`+`') then
            sum_2 := `+`(op(Beautify([op(expr_2)], digits, assumptions, _options, ':-beautify_top_level' = false)));
        else
            sum_2 := Beautify(expr_2, digits, assumptions, _options, ':-beautify_top_level' = false);
        end if;

        return frontend(processOneSum, [sum_1, sum_2],
                        [processOneSumType, {}], digits, assumptions, _options);
    end proc;

    processOneSum := proc(expr_1 :: algebraic,
                          expr_2 :: algebraic,
                          digits :: posint,
                          assumptions :: list({`<`, `<=`}),
                          {output :: {identical(float), identical(rational)} := ':-float'},
                          {underflow :: truefalse := true},
                          {process_exact :: truefalse := false},
                          {round :: truefalse := true},
                          $)
    local
        assumption_matrix,
        coefficient_logs,
        coefficient_magnitudes,
        coefficients,
        coefficients_1,
        coefficients_2,
        comparator,
        exact,
        exact_1,
        exact_2,
        exponents,
        exponents_1,
        exponents_2,
        precisions,
        variables;

        if (not (underflow or round)) or (expr_1 = 0 and expr_2 = 0) then
            return expr_1 + expr_2;
        end if;

        variables := convert(indets([expr_1, expr_2], ':-name') minus {constants}, ':-list');

        comparator := compareLexicographically(nops(variables));

        coefficients_1, exponents_1, exact_1 := dissect(expr_1, variables);
        coefficients_2, exponents_2, exact_2 := dissect(expr_2, variables);

        sortRows(coefficients_1, exponents_1, exact_1, comparator);
        sortRows(coefficients_2, exponents_2, exact_2, comparator);

        coefficients, exponents, exact, coefficient_magnitudes, precisions :=
        mergeRows(coefficients_1, exponents_1, exact_1,
                  coefficients_2, exponents_2, exact_2,
                  comparator, digits, underflow, process_exact);

        coefficient_logs := log[10]~(coefficient_magnitudes);

        assumption_matrix := preprocessAssumptions(assumptions, variables);

        if nops(variables) = 1 then
            processTermsUnivariate(coefficients, coefficient_magnitudes, coefficient_logs,
                                   exponents, exact, precisions, digits, assumption_matrix, process_exact);
        else
            processTermsMultivariate(coefficients, coefficient_magnitudes, coefficient_logs, 
                                     exponents, exact, precisions, digits, assumption_matrix, process_exact);
        end if;
        
        return computeResult(variables, coefficients, exponents, exact, precisions, 
                             `if`(output = ':-float', roundFloat, roundRational),
                              underflow, process_exact, round);
        
    end proc;

    preprocessAssumptions := proc(assumptions :: list({`<`, `<=`}),
                                  variables :: list(name),
                                  $)
    local
        allowedType,
        assumption,
        assumption_no_abs,
        bound,
        expr,
        k,
        m,
        n,
        naked_names,
        var;
        
        m := Matrix(nops(assumptions), nops(variables)+1, ':-datatype' = ':-float');
        n := 0;

        allowedType := {x -> x, '':-`*`''}({':-numeric', ':-name', ':-name'^':-numeric'});

        for assumption in assumptions do
            assumption_no_abs := eval(assumption, abs = (x -> x));
            expr := evalf(lhs(assumption_no_abs) / rhs(assumption_no_abs));
            bound := 1/eval(expr, map(`=`, variables, 1));
            if not type(bound, ':-numeric') then
                userinfo(4, ':-BeautyClinic', "ignoring assumption that depends on variables that do not occur "
                         "in the expression:", assumption);
                next;
            end if;

            expr := expr * bound;
            expr := eval(expr, map(`=`, select(`=`, indets(expr, ':-float'), 1), 1));
            if not type(expr, allowedType) or type(bound, ':-nonreal') then
                WARNING("ignoring assumption of an unsupported type: %1", assumption);
                next;
            end if;

            bound := abs(bound);
            if eval(rhs(assumption_no_abs), map(`=`, variables, 1)) < 0 then
                # The sense of the inequality has been inverted -
                # invert it back.
                bound := 1/bound;
                expr := 1/expr;
            end if;

            n := n + 1;
            for k to nops(variables) do
                if has(expr, variables[k]) then
                    m[n, k] := variables[k] * diff(expr, variables[k]) / expr;
                end if;
            end do;

            m[n, -1] := log[10](bound);

            # Check if we need to warn about the interpretation of the assumption:
            assumption_no_abs := evalf(eval(assumption, abs = (x -> 1)));
            naked_names := map(indets, [op(assumption_no_abs)], ':-name');
            if (not member({}, naked_names)) # This case is complicated and usually deserves a warning.
            or (naked_names[1] = {} and naked_names[2] <> {} and lhs(assumption_no_abs) < 0)
            or (naked_names[1] <> {} and naked_names[2] = {} and rhs(assumption_no_abs) > 0) then 
                WARNING("interpreting assumption %1 as %2", assumption, expand(abs(expr)) <= bound);
            else
                userinfo(3, ':-BeautyClinic', "interpreting assumption", assumption, 
                         "as", expand(abs(expr)) <= bound);
            end if;


        end do;

        if n = nops(assumptions) then
            return m;
        else
            return m[.. n];
        end if;
    end proc;

    # Gather terms together in order to obtain generalized
    # polynomial. Rounds every non-exact coefficient to the correct number
    # of digits; if this number is 0, omits that term.
    computeResult := proc(variables :: list,
                          coefficients :: Vector,
                          exponents :: Matrix,
                          exact :: Vector,
                          precisions :: Vector,
                          rounder :: appliable,
                          underflow :: truefalse,
                          process_exact :: truefalse,
                          round :: truefalse,
                          $)
    local
        i,
        j,
        m,
        n,
        results;
        
        n, m := op(1, exponents);
        results := Vector(n);
        
        for i to n do
            if precisions[i] > 0 then
                if exact[i] = 0 and round then 
                    results[i] := rounder(precisions[i], coefficients[i]) *
                    mul(variables[j]^exponents[i, j], j=1..m);
                else
                    results[i] := coefficients[i] * mul(variables[j]^exponents[i, j], j=1..m);
                end if;
            elif exact[i] = 0 and not underflow then
                results[i] := rounder(0, coefficients[i]) *
                mul(variables[j]^exponents[i, j], j=1..m);
            elif exact[i] = 1 and not process_exact then
                results[i] := coefficients[i] * mul(variables[j]^exponents[i, j], j=1..m);
            end if;
        end do;

        return sort(add(i, i in results));
    end proc;

    roundFloat := proc(precision :: numeric,
                       coefficient :: algebraic,
                       $)
        if precision > 0 then
            if type(coefficient, ':-float') and ilog10(op(1, coefficient)) < precision then
                # There are currently fewer than precision digits -
                # don't invent any new ones.
                return coefficient;
            else
                return evalf[ceil(precision)](coefficient);
            end if;
        else
            return evalf[1](coefficient);
        end if;
    end proc;

    roundRational := proc(precision :: numeric,
                          coefficient :: algebraic,
                          $)
    local
        cf,
        convergents,
        delta,
        i,
        n,
        result;
        
        if Im(coefficient) <> 0 then
            return procname(precision, Re(coefficient)) + I * procname(precision, Im(coefficient));
        elif coefficient = 0 then
            return coefficient;
        end if;

        if precision > 0 then
            delta := 1./10^precision;
            n := 5;
            numtheory:-cfrac(coefficient, n, 'convergents');
            for i while abs((coefficient - convergents[i])/coefficient) > delta do
                if i >= n then
                    n := 2*n;
                    numtheory:-cfrac(coefficient, n, 'convergents');
                end if;
            end do;
            return convergents[i];
        else
            if coefficient > 1 or coefficient < 0 then
                return floor(coefficient);
            else
                return 1 / floor(1/coefficient);
            end if;
        end if;
    end proc;

    # Computes precisions to account for inter-term propagation of inaccuracy.
    processTermsMultivariate := proc(coefficients :: Vector,
                                     coefficient_magnitudes :: Vector,
                                     coefficient_logs :: Vector,
                                     exponents :: Matrix,
                                     exact :: Vector,
                                     precisions :: Vector,
                                     digits :: posint,
                                     assumptions :: Matrix,
                                     process_exact :: truefalse,
                                     $)
    local
        constraint_matrix,
        constraint_matrix_left_part,
        minus_float_exponents,
        m,
        many_ones,
        n,
        n_assumptions,
        negative_coefficient_logs,
        objective,
        solution,
        term_index,
        upper_bounds;
        
        n, m := op(1, exponents);
        n_assumptions := op([1, 1], assumptions);

        # First m columns are logs of variables. Last column is the
        # objective variable. 
        # First n-1 rows of constraints correspond to term
        # ratios. Following n_assumptions rows correspond to
        # assumptions. Last two rows say that the objective variable should
        # stay between -digits-1/2 and 0: we need to bound it away
        # from +infinity to detect unbounded solutions - a value > 0
        # would lead to no change anyway; and a value < -digits does
        # not improve things beyond -digits (we choose -digits - 1/2
        # to combat roundoff). If there is no feasible point, then
        # that must be because of the last constraint, so we can set
        # precision to 0 in that case.
        constraint_matrix := Matrix(n + n_assumptions + 1, m + 1, ':-datatype' = ':-float'[8]);
        constraint_matrix[.., m + 1] := 1;
        constraint_matrix[-1, m + 1] := -1;
        constraint_matrix_left_part := ArrayTools:-Alias(constraint_matrix,
                                                         0, [n + n_assumptions + 1, m]);

        if n_assumptions > 0 then 
            ArrayTools:-BlockCopy(assumptions, 0, n_assumptions, n_assumptions, m, 
                                  constraint_matrix, n - 1, n + n_assumptions + 1, n_assumptions, m);
            constraint_matrix[n .. -3, m + 1] := 0;
        end if;

        objective := Vector(m + 1, ':-datatype' = ':-float'[8]);
        objective[-1] := -1;

        upper_bounds := Vector(n + n_assumptions + 1, ':-datatype' = ':-float'[8]);
        ArrayTools:-Copy(n_assumptions, assumptions, m * n_assumptions, 1, upper_bounds, n - 1, 1);
        upper_bounds[-2] := 1/2;
        upper_bounds[-1] := digits + 1/2;
        negative_coefficient_logs := - coefficient_logs;
        minus_float_exponents := Matrix(exponents, ':-datatype' = ':-float'[8]);
        many_ones := Matrix(n + n_assumptions + 1, 1, ':-datatype' = ':-float'[8]);
        many_ones[.. n - 1] := 1;

        for term_index to n do
            if (exact[term_index] = 0 or process_exact) 
            and coefficients[term_index] <> 0 and precisions[term_index] > 0 then

                try
                    constraint_matrix_left_part[.. term_index - 1] := minus_float_exponents[.. term_index - 1];
                    constraint_matrix_left_part[term_index .. n - 1] := minus_float_exponents[term_index + 1 ..];
                    LinearAlgebra:-LA_Main:-MatrixAdd(constraint_matrix_left_part, 
                                                      many_ones . (- minus_float_exponents[term_index .. 
                                                                                           term_index]),
                                                      1, 1, ':-inplace' = true, ':-outputoptions' = []);
                    
                    upper_bounds[.. term_index - 1] := negative_coefficient_logs[.. term_index - 1];
                    upper_bounds[term_index .. n - 1] := negative_coefficient_logs[term_index + 1 ..];
                    
                    if abs(coefficients[term_index]) = coefficient_magnitudes[term_index] then
                        upper_bounds[.. -3] := map[':-evalhf'](`+`, upper_bounds[.. -3], 
                                                               coefficient_logs[term_index]);
                    else
                        upper_bounds[.. -3] := map[':-evalhf'](`+`, upper_bounds[.. -3],
                                                               log[10](abs(evalf(coefficients[term_index]))));
                    end if;
                    
                    solution := Optimization:-LPSolve(objective, [constraint_matrix, upper_bounds])[1];
                    
                    precisions[term_index] := 
                    max(0, min(precisions[term_index], digits - solution));
                catch "no feasible solution found":
                    precisions[term_index] := 0;
                end try;
                
            end if;
        end do;
        
    end proc;
    
    processTermsUnivariate := module()
    local 
        ModuleApply, 
        assignOnePrecision,
        assignPrecisions,
        determineUpperSegments,
        interpretAssumptions;

        ModuleApply := proc(coefficients :: Vector,
                            coefficient_magnitudes :: Vector,
                            coefficient_logs :: Vector,
                            exponents :: Matrix,
                            exact :: Vector,
                            precisions :: Vector,
                            digits :: posint,
                            assumptions :: Matrix,
                            process_exact :: truefalse,
                            $)
        local
            max_x,
            min_x,
            upper_segment_indices,
            upper_segment_x_coordinates,
            n,
            n_segments;
            
            n := op(1, coefficients);
            
            if n <= 1 then
                # Can't eliminate any term.
                return;
            end if;
            
            min_x, max_x := interpretAssumptions(assumptions);

            upper_segment_indices,
            upper_segment_x_coordinates,
            n_segments := determineUpperSegments(n, min_x, max_x, coefficient_logs, exponents);
            
            assignPrecisions(n,
                             coefficient_logs,
                             exponents,
                             exact,
                             precisions,
                             digits,
                             process_exact,
                             upper_segment_indices, 
                             upper_segment_x_coordinates,
                             n_segments,
                             max_x);
        end proc;

        interpretAssumptions := proc(assumptions :: Matrix, $)
        local
            i,
            min_x := HFloat(-infinity),
            max_x := HFloat(infinity);

            for i to op([1, 1], assumptions) do
                if assumptions[i, 1] > 0 then
                    max_x := min(max_x, assumptions[i, 2]/assumptions[i, 1]);
                else
                    min_x := max(min_x, assumptions[i, 2]/assumptions[i, 1]);
                end if;
            end do;
            
            return min_x, max_x;
        end proc;
        
        assignPrecisions := proc(n :: posint,
                                 coefficient_logs :: Vector,
                                 exponents :: Matrix,
                                 exact :: Vector,
                                 precisions :: Vector,
                                 digits :: posint,
                                 process_exact :: truefalse,
                                 upper_segment_indices :: Vector,
                                 upper_segment_x_coordinates :: Vector,
                                 n_segments :: posint,
                                 max_x :: {numeric, infinity},
                                 $)
        local
            i,
            j,
            left_index,
            right_index,
            segment_y_coordinate,
            upper_segment_y_coordinate;

            # Deal with segments before first segment in upper bound:
            right_index := upper_segment_indices[1];

            upper_segment_y_coordinate := coefficient_logs[right_index] +
            upper_segment_x_coordinates[1] * exponents[right_index, 1];
            for j from n to upper_segment_indices[1] + 1 by -1 do
                assignOnePrecision(j, upper_segment_x_coordinates[1], upper_segment_y_coordinate,
                                   coefficient_logs, exponents, exact, precisions, digits, process_exact);
            end do;

            # Deal with segments in between segments in upper bounds:
            for i to n_segments - 1 do
                left_index := right_index;
                right_index := upper_segment_indices[i+1];

                upper_segment_y_coordinate := coefficient_logs[right_index] + 
                upper_segment_x_coordinates[i + 1] * exponents[right_index, 1];

                for j from left_index - 1 to right_index + 1 by -1 do
                    assignOnePrecision(j, upper_segment_x_coordinates[i+1], upper_segment_y_coordinate,
                                       coefficient_logs, exponents, exact, precisions, digits, process_exact);
                end do;
            end do;

            # Deal with segments after last segment in upper bound:
            upper_segment_y_coordinate := coefficient_logs[right_index] + max_x * exponents[right_index, 1];
            for j from right_index - 1 to 1 by -1 do
                assignOnePrecision(j, max_x, upper_segment_y_coordinate,
                                   coefficient_logs, exponents, exact, precisions, digits, process_exact);
            end do;

        end proc;

        assignOnePrecision := proc(j :: posint,
                                   x :: numeric,
                                   y :: numeric,
                                   coefficient_logs :: Vector,
                                   exponents :: Matrix,
                                   exact :: Vector,
                                   precisions :: Vector,
                                   digits :: posint,
                                   process_exact :: truefalse,
                                   $)
        local
            segment_y;
            
            if exact[j] = 0 or process_exact then
                segment_y := coefficient_logs[j] + x * exponents[j, 1];
                
                precisions[j] := min(precisions[j], digits - (y - segment_y));
            end if;
        end proc;
        
        determineUpperSegments := proc(n :: posint,
                                       min_x :: {numeric, infinity},
                                       max_x :: {numeric, infinity},
                                       coefficient_logs :: Vector,
                                       exponents :: Matrix,
                                       $)
        local
            i,
            intersection_x_coordinate,
            j,
            line_index,
            upper_segment_indices,
            upper_segment_x_coordinates, # kth entry is the x-coordinate to the *left* of the kth segment
            n_segments;
            
            upper_segment_indices := Vector(n, ':-datatype' = ':-integer'[kernelopts(wordsize)/8]);
            upper_segment_x_coordinates := Vector(n, ':-datatype' = ':-float'[8]); 
            
            upper_segment_indices[1] := n;
            upper_segment_x_coordinates[1] := min_x;
            n_segments := 1;
            
            for i from n-1 to 1 by -1 do
                # Find the rightmost segment (say the kth of the ones
                # we have so far) in the upper segments such that the
                # intersection point of the kth with the new segment
                # is to the right of the upper segment boundary
                # currently listed as the kth.

                j := n_segments;
                line_index := upper_segment_indices[n_segments];
                intersection_x_coordinate := (coefficient_logs[line_index] - coefficient_logs[i]) /
                (exponents[i, 1] - exponents[line_index, 1]);

                while intersection_x_coordinate < upper_segment_x_coordinates[j] and j > 1 do
                    j := j - 1;
                    
                    line_index := upper_segment_indices[j];
                    
                    intersection_x_coordinate := (coefficient_logs[line_index] - coefficient_logs[i]) /
                    (exponents[i, 1] - exponents[line_index, 1]);
                end do;
                
                if intersection_x_coordinate < min_x then
                    n_segments := 1;
                    upper_segment_indices[1] := i;
                    upper_segment_x_coordinates[1] := min_x;
                elif intersection_x_coordinate <= max_x then
                    # We intersect the jth upper segment.
                
                    n_segments := j + 1;
                
                    upper_segment_indices[n_segments] := i;
                    upper_segment_x_coordinates[n_segments] := intersection_x_coordinate;
                # else # in this case, the ith segment doesn't occur
                       # in the upper segments, so we just skip it.
                end if;
            end do;
            
            return upper_segment_indices[.. n_segments],
            upper_segment_x_coordinates[.. n_segments],
            n_segments;
        end proc;
        end module;
    
    # Merge the terms of two generalized polynomials, deduce precision,
    # and compute the maximum magnitude of the coefficient for a given degree.
    mergeRows := proc(coefficients_1 :: Vector, 
                      exponents_1 :: Matrix,
                      exact_1 :: Vector,
                      coefficients_2 :: Vector,
                      exponents_2 :: Matrix,
                      exact_2 :: Vector,
                      comparator :: procedure,
                      digits :: posint,
                      underflow :: truefalse,
                      process_exact :: truefalse,
                      $)
    local
        coefficients,
        exact,
        exponents,
        i_1,
        i_2,
        k,
        m,
        maxcoefficients,
        n_1,
        n_2,
        precisions,
        ratio;
        
        n_1 := op(1, coefficients_1);
        n_2 := op(1, coefficients_2);
        m := op([1, 2], exponents_1);
        
        coefficients := Vector(n_1 + n_2);
        exponents := Matrix(n_1 + n_2, m);
        exact := Vector(n_1 + n_2, ':-datatype' = ':-integer'[kernelopts(wordsize)/8]);
        maxcoefficients := Vector(n_1 + n_2, ':-datatype' = ':-float'[8]);
        precisions := Vector(n_1 + n_2, ':-datatype' = ':-float'[8],
                             ':-fill' = digits);
        
        i_1 := 1;
        i_2 := 1;

        for k while i_1 <= n_1 and i_2 <= n_2 do

            if EqualEntries(exponents_1[i_1], exponents_2[i_2]) then

                exponents[k] := exponents_1[i_1];
                exact[k] := exact_1[i_1] * exact_2[i_2];

                if exact_1[i_1] + exact_2[i_2] = 1 then
                    # Exactly one of the coefficients is exact; make
                    # the float infect the other one:
                    coefficients[k] := evalf(coefficients_1[i_1] + coefficients_2[i_2]);
                else
                    coefficients[k] := coefficients_1[i_1] + coefficients_2[i_2];
                end if;

                if exact[k] = 0 then 
                    maxcoefficients[k] := max(abs(coefficients[k]),
                                              abs(coefficients_1[i_1]), 
                                              abs(coefficients_2[i_2]));
                else
                    maxcoefficients[k] := max(abs(evalf(coefficients[k])),
                                              abs(evalf(coefficients_1[i_1])), 
                                              abs(evalf(coefficients_2[i_2])));
                end if;

                if underflow and (exact[k] = 0 or process_exact) then
                    ratio := abs(coefficients[k]) / maxcoefficients[k];
                    precisions[k] := max(log[10](ratio) + digits, 0);
                end if;
                
                i_1 := i_1 + 1;
                i_2 := i_2 + 1;

            elif comparator(exponents_1[i_1], exponents_2[i_2]) then

                coefficients[k] := coefficients_1[i_1];
                exponents[k] := exponents_1[i_1];
                exact[k] := exact_1[i_1];
                if exact[k] = 0 then 
                    maxcoefficients[k] := abs(coefficients_1[i_1]); 
                else
                    maxcoefficients[k] := abs(evalf(coefficients_1[i_1]));
                end if;

                i_1 := i_1 + 1;

            else

                coefficients[k] := coefficients_2[i_2];
                exponents[k] := exponents_2[i_2];
                exact[k] := exact_2[i_2];
                if exact[k] = 0 then
                    maxcoefficients[k] := abs(coefficients_2[i_2]); 
                else
                    maxcoefficients[k] := abs(evalf(coefficients_2[i_2]));
                end if;

                i_2 := i_2 + 1;

            end if;
        end do;

        if i_1 <= n_1 then
            # Copy rest of #1

            for k from k while i_1 <= n_1 do

                coefficients[k] := coefficients_1[i_1];
                exponents[k] := exponents_1[i_1];
                exact[k] := exact_1[i_1];
                if exact[k] = 0 then 
                    maxcoefficients[k] := abs(coefficients_1[i_1]); 
                else
                    maxcoefficients[k] := abs(evalf(coefficients_1[i_1]));
                end if;

                i_1 := i_1 + 1;
            end do;
        else
            # Copy rest of #2

            for k from k while i_2 <= n_2 do

                coefficients[k] := coefficients_2[i_2];
                exponents[k] := exponents_2[i_2];
                exact[k] := exact_2[i_2];
                if exact[k] = 0 then
                    maxcoefficients[k] := abs(coefficients_2[i_2]); 
                else
                    maxcoefficients[k] := abs(evalf(coefficients_2[i_2]));
                end if;

                i_2 := i_2 + 1;
            end do;
        end if;
        
        return coefficients[..k-1], exponents[..k-1], exact[..k-1], 
        maxcoefficients[..k-1], precisions[..k-1];
    end proc;

    # Sort exponents lexicographically (for now) and sort coefficients and
    # exact with the same permutation.
    sortRows := proc(coefficients :: Vector,
                     exponents :: Matrix,
                     exact :: Vector,
                     comparator :: procedure,
                     $)
    local
        columns,
        index,
        rows;
        
        rows, columns := op(1, exponents);
        index := sort([seq(1 .. rows)], (i, j) -> comparator(exponents[i], exponents[j]));

        coefficients[..] := coefficients[index];
        exponents[..] := exponents[index];
        exact[..] := exact[index];
    end proc;

    compareLexicographically := proc(n, $)
        return proc(row1, row2)
               local k;
                   for k to n do
                       if row1[k] > row2[k] then
                           return true;
                       elif row1[k] < row2[k] then
                           return false;
                       end if;
                   end do;
                   return true;
               end proc;
    end proc;

    # Take an expression in variables apart into vectors and matrices of
    # coefficients, exponents, and booleans indicating exactness.
    dissect := proc(expr :: algebraic,
                    variables :: list(name),
                    $)
    local
        coefficient,
        coefficients,
        exact,
        expanded_expr,
        exponents,
        i,
        power_denominator,
        term,
        terms;

        expanded_expr := collect(expr, variables, ':-distributed');

        if type(expanded_expr, ':-`+`') then
            terms := [op(expanded_expr)];
        elif expanded_expr = 0 then
            terms := [];
        else
            terms := [expanded_expr];
        end if;

        coefficients := Vector(nops(terms));
        exponents := Matrix(nops(terms), nops(variables));
        exact := Vector(nops(terms), ':-datatype' = ':-integer'[kernelopts(wordsize)/8]);

        power_denominator := ilcm(op(map(power -> denom(op(2, power)), indets(terms, ':-`^`'))));
        i := 0;
        
        for term in terms do
            i := i + 1;

            exponents[i] := Vector(degree~(term^power_denominator, variables) / power_denominator);
            coefficient := eval(term, zip(`=`, variables, 1));
            
            if type(coefficient, ':-float') then
                coefficients[i] := coefficient;
            elif hastype(coefficient, ':-float') then
                coefficients[i] := evalf(coefficient);
            else
                exact[i] := 1;
                coefficients[i] := coefficient;
            end if;

        end do;

        return coefficients, exponents, exact;
    end proc;
end module;

#savelib('BeautyClinic');