@Carl Love   Your code is compact and elegant, but conceptually more complex and less effective:

t := time():  SpM(1000):  time()-t;

t := time():  Am(1000):  time()-t;

Of course, this can be done in different ways. Here are 2 other ways:

restart;

f:=x+y+z:

V:=<seq(diff(f, var), var=[x,y,z])>;

restart;

f:=x+y+z:  V:=Vector():  var:=[x,y,z]:

for i to nops(var) do

V(i):=diff(f, var[i]):

od:

V;

@Carl Love   It is sufficient to find the exact roots of the polynomial (the first command in my answer)

Another example:

factor(3*x^2-7);

factor(3*x^2-7, {sqrt(3), sqrt(7)});

 -op(1,%)*3*``(expand(op(2,%)/sqrt(3)))*``(-expand(op(3,%)/sqrt(3)));

@Carl Love   In your example we can do exact symbolic factorization:

A := factor(x^4+1, sqrt(2));

B := ``(op(1, A)*op(2, A))*op(3, A);

@Markiyan Hirnyk  My program has the following differences from these:

1) The program works for systems of all sizes, and these tutors only if the every dimension of the matrix  no greater than 5.

2) My program in the during of Gauss - Jordan elimination , simulating manual solution produces some simplifying transformations. For example, solving the system  {19*x+20*y=39, x+y=2}, to get rid from fractions, useful to swap positions of equations.

@Axel Vogt   Yes, if we increase Digits, the error disappears:

Digits := 30:

KsT := 0.10e11: KaT := 4.6*0.10e6: KsE2 := 3.14*0.10e10: KaE2 := 4.21*0.10e6: KsDHT := 3*0.10e6: KaDHT := 3.5*0.10e6:

eq1 := E2 = fE2*(1+(KaE2+Ca)/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsE2*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq2 := T = fT*(1+KaT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq3 := DHT = fDHT*(1+KaDHT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsDHT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

E2, T, DHT, Ca, Cshbg := 1, 2, 3, 4, 5:

solve({eq1, eq2, eq3}, useassumptions) assuming fE2 > 0, fT > 0, fDHT > 0;

     {fDHT = 0.300792608884171123405264656623e-5, fE2 = 3.02872370161763976538808187091*10^(-10), fT    = 1.90199449982478958012781172005*10^(-10)}

@velimir74  I do not know the cause of this extra solution. In fact, this is a bug, because easy to verify that this is not a solution:

KsT := 0.10e11: KaT := 4.6*0.10e6: KsE2 := 3.14*0.10e10: KaE2 := 4.21*0.10e6: KsDHT := 3*0.10e6: KaDHT := 3.5*0.10e6:

eq1 := E2 = fE2*(1+(KaE2+Ca)/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsE2*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq2 := T = fT*(1+KaT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq3 := DHT = fDHT*(1+KaDHT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsDHT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

E2, T, DHT, Ca, Cshbg:=1, 2, 3, 4, 5:

eval({eq1, eq2, eq3}, {fE2 = 7.15539742720000*10^10, fT = 3.30345789526199, fDHT =

4.72656319086592*10^14});

                               {1 = .7155397427e11, 2 = 3.303457895, 3 = .4726563191e15}

@wingjammer

UniqueColumns1:=proc(A::Matrix)

local L;

uses LinearAlgebra, ListTools;

L:=[Categorize((x,y)->Equal(x,y), [seq(A[..,i],i=1..ColumnDimension(A))])];

seq(L[i,1], i=1..nops(L));

end;

Your example:

Q := Matrix([[1,1,1,0,0,0,-2,-1], [0,0,0,1,1,1,0,-3]]);

UniqueColumns1(Q);

Addition: another way.

UniqueColumns2:=proc(A::Matrix)

uses LinearAlgebra, ListTools;

MakeUnique([seq(A[..,i], i=1..ColumnDimension(A))], 1, (x,y)->Equal(x,y));

end; 

@Carl Love   

sqrt(139.^2+105^2);

                      174.2010333

@as5987   If you are using the sequence then you should write

with(plots):  with(plottools):

d := seq(circle([1, 1/i], 1/i, color = blue), i=0.5..20, 0.1):

Example:

display(d[1], d[10], d[50]);

The number in brackets should be positive integer.  It shows the sequence number a member of the sequence.

@as5987  plottools[disk]  command creates plot structure of filled circle centered at  (1, 0)  and the radius  0.01 .  The same can be written shorter

Point1 := disk(eval([x,y], Sol[1]), 0.01, color = red):

@Carl Love   Thanks for the helpful comments.

@casperyc   I do not know what you call efficient. Here is variant without loop:

restart;

a:=Vector([2,3,4,5]): b:=convert(a,list):

map(ListTools[SearchAll],select(x->x>=3, b), b);

                                 [2, 3, 4]

@Carl Love  Your code is the most simple and effective solution to the problem. I learned from it about the useful  foldl  command. 

@Markiyan Hirnyk  I do not understand what you mean by this example? The sum computed symbolically absolutly exactly. Try to do the same by your code.