@Kitonum Thanks. This's justified (and faster). But if the problem has twenty variables, I'll have to write 20 do loops, which seems to be somewhat inadvisable. (Yes, the original problem asks for twenty variables instead of nine variables, but as this problem appears to be too expensive for a low-end computer, I just post a modified version here. (I choose "nine" because (-(-5)^3*(9-1))^(1/3) is simply a nice number.))

@vv I second that. I once wanted to use a more uncommon function in my question (to avoid certain non-essential "tricks"), but at last I gave up, as a simple function can achieve the desired effect as well.

@Axel Vogt Thanks. But it appears that this doesn't work for other functions.

Another potential workaround is: 

> m := Matrix([[1, 1, 3], [0, 1, 0], [0, I, 2]], datatype = complex(sfloat)):
> UseHardwareFloats := false:++Digits:
> add(residue~('z'/(exp('z') - 1)*LinearAlgebra:-MatrixInver\
> se('z'*LinearAlgebra:-IdentityMatrix(LinearAlgebra:-Dimens\
> ion(m)) - m), 'z' = lambda), lambda in ListTools:-MakeUniq\
> ue(LinearAlgebra:-Eigenvalues(m, output = list), 1, testfloat, 1));
                  [0.58197670686 + 0. I    -0.33869688735 + 0.20926639797 I    -0.80682426410 + 0. I]
                  [                                                                                 ]
                  [     0. + 0. I                0.58197670686 + 0. I                0. + 0. I      ]
                  [                                                                                 ]
                  [     0. + 0. I                0. - 0.26894142136 I          0.31303528550 + 0. I ]

Unfortunately, the default value of Digits is not enough here. (And I don't think that this is efficient (though not ineffective) for other functions.)

@dharr Thanks for your detailed analysis.
I read that "many LinearAlgebra package functions use NAG routines to carry out computations", while the LinearAlgebra['MatrixFunction'] does not make use of this external library. This is strange, since such routines are available (according to F01 (Matop) : matrix functions) in fact. 

@dharr Thanks. But I do not understand why the obsolete linalg:-matfunc (whose use is discouraged) performs better even with the default 10-digit precision. Equally, can't the built-in LinearAlgebra:-MatrixFunction go beyond the basics?

@Thomas Richard Actually, the OP's purpose is to display "x^2-2*x-3*x+6,x^2+3*x-7*x-21,x^2-9*x-10*x+90" instead of "x^2-5*x+6,x^2-4*x-21,x^2-19*x+90". For such a purpose, one may make use of InertForm:-Display(… have to be processed in advance …, inert = false). 

However, this is not very easy (though not impossible) to handle. In Mathematica, one can simply use: 

In[1]:= HoldForm[x^2-a*x-b*x+a*b]/.Distribute[{{a,b}}\[Rule]{{2,3},{-3,7},{9,10}},List,Rule,List,Rule/*Thread]

          2                     2                      2
Out[1]= {x  - 2 x - 3 x + 2 3, x  - -3 x - 7 x - 3 7, x  - 9 x - 10 x + 9 10}

There are also some other Hold* in Mma, yet I cannot find similar direct commands in Maple.

@Carl Love Many thanks. Surprisingly, I find that even if I don't compile it, the execution time of your implementation is still shorter than that of mine (well, tens of seconds faster).
Does loop-based code (instead of the so-called functional operation, which is used in my initial translation) typically provide better performance in similar cases?

To obtain some information about what the compiler is doing, you can set infolevel[Compiler] to a positive value (1 or 3). What can you find?

@dharr Many thanks. I believe that if GraphTheory:-AllPairsDistance (or networks:-allpairs) is fixed, this task will be more convenient to be finished.

@acer & @tomleslie Many thanks. I have just submitted a software change request.

@tomleslie According to updates,Maple2021,GraphTheory, GraphTheory['AllPairsDistance'] has had their efficiency improved. But somebody tells me that he cannot reproduce these in a "legacy" Maple. Can anyone compare the output from Maple 2020 and that from Maple2021 respectively?

@Axel Vogt & @Preben Alsholm Thanks. I believe that the built-in `sum` (as well as `int` mentioned in the previous question Integration problem) lacks something like automatic preprocessing and postprocessing, so the users have to complete such steps manually, which makes Maple‘s learning curve steeper.

@acer The much "ideal" result may be ; unfortunately, Maple cannot confirm this: 

:-simplify(diff(-x^(6-alpha)*Ei(alpha-5, beta*x), x)-x^(5-alpha)/exp(beta*x));
 = 
             (5 - alpha)                      
         -6 x            Ei(alpha - 5, beta x)

               (5 - alpha)                            
            + x            Ei(alpha - 5, beta x) alpha

               (6 - alpha)                            
            + x            beta Ei(-6 + alpha, beta x)

               (5 - alpha)             
            - x            exp(-beta x)

However, a less ideal result () can be obtained as follows: 

:-simplify(eval(convert(int(eval(exp(-beta*x)*x^(5-alpha), 5-alpha = a), x), GAMMA), a = 5-alpha));
 = 
          1   /        alpha  (-alpha)                   
      - ----- \(beta x)      x         (-GAMMA(6 - alpha)
            6                                            
        beta                                             

                                    \
         + GAMMA(6 - alpha, beta x))/
:-simplify(diff((GAMMA(-alpha + 6) - GAMMA(-alpha + 6, beta*x))*(beta*x)^alpha/(beta^6*x^alpha), x) - exp(-beta*x)*x^(5 - alpha));
 = 
                               0

@C_R 

It has been explained in evalapply (cf. `convert/algebraic`). But one may also use: 

(_ -> content(_) %* primpart(_))(sqrt(2*b - 4*a));
 = 
                    / (1/2)            (1/2)\
                  %*\2     , (-2 a + b)     /

@C_R As for "a command that factors sqrt(2) out of sqrt(2*a+2*b)", one may use: 

> (content*primpart)(sqrt(2*a+2*b));
 = 
                       (1/2)        (1/2)
                      2      (a + b)     

However, I don't think that these names are as common as factor (or the new command Physics:-Library:-FactorSum) is in general; Mathematica christens such symbols "FactorTerms", which is much easier to find for ordinary users (for example, if one searchs for factor, one will see "Factor", "FactorTerms", "FactorSquareFree" and so on).