You can't prove such a theorem using a CAS.
But you can check it for any value or range of values. For example:

seq( [n, isprime(n), mods( (n-1)!, n ) = -1],  n=2..100);

It is not uncommon in mathematics to solve a problem using two metods and obtain very different forms as answers. And sometimes it is very difficult to prove directly that the results are the same.

In our case this can be done using Maple itself.

a:=21*arccos(RootOf(448*_Z^7+192*_Z^6-784*_Z^5-288*_Z^4+392*_Z^3+108*_Z^2-49*_Z-6, index = 2))/Pi:
b:=-(21*I)*ln(RootOf(7*_Z^14+6*_Z^13+6*_Z+7, index = 2))/Pi:
evalf(a-b);
        0.+1.060497132*10^(-10)*I
#      so, a=b is plausible
c:=simplify(Pi/21*(a-b)):
expand(  cos(c) - 1):
convert(%, RootOf):
evala(%);
          0

At this moment we are almost done (provided that Maple's computations are correct).
We know that c is a multiple of 2*Pi, so it would remain e.g. to show that |c| < 1 to conclude that c=0; this is not difficult.

I always recommend in such situations to try a parametrization of the surface.
In your example this needs some work but (almost) nothing is free in our world.

Xmax:=1.195918936:
Ymax:= X -> -X^2+X+1/80+(1/80)*sqrt(-160*X^2+320*X+161):
Ymin:= X -> piecewise(X<.2263095011,0,X<.7082018428, -X^2+X+1/80-(1/80)*sqrt(-160*X^2+320*X+161),0):
ZXY:=(X,Y)-> -(40/19)*X^4+(80/19)*X^3-(80/19)*X^2*Y-(40/19)*X^2+(80/19)*X*Y-(40/19)*Y^2+(1/19)*X+(1/19)*Y+1/19:
P:=NULL:
for i in [-1,1] do for j in [-1,1] do for k in [-1,1] do
P:=P, plot3d( [i*sqrt(abs(X)), j*sqrt(abs(Y)), k*sqrt(abs(ZXY(X,Y)))], X=0..Xmax, Y=Ymin(X)..Ymax(X)):
od:od:od:
plots:-display(P);

The advantage is that the intersection points are more clearly separated ==>  better visibility for higher number of lines.

UniformCuts:=proc(N::posint)
local NN:=floor(N/2)*2+1,
      R:=1.25/sin(Pi/2./NN), d:=arccos(1/R),k,p,a,b;
uses plots,plottools;
for k from 0 to N-1 do
  a[k]:=[R*cos(2*k/NN*Pi+d), R*sin(2*k/NN*Pi+d)];
  b[k]:=[R*cos(2*k/NN*Pi-d), R*sin(2*k/NN*Pi-d)];
od:
p:=seq(line(a[k],b[k]),k=0..N-1):
display(disk([0,0],R,color=yellow),p, axes=none, color=red);
end:

UniformCuts(9);

# animation
A:=seq( plots:-display(seq(op(i,%),i=1..k)), k=1..9+1 ):
plots:-display(A, axes=none, insequence, color=red);

Note first that the problem is very ill-conditioned (see previous answer), so unless your coefficients are (almost) exact, the problem is (almost) nonsense.
So, let us assume the coefficients are exact and work with 200 digits!

Unfortunately fsolve refuses to work if the coefficients are complex (but it should!).
Fortunately, using a change of variable it is possible to get real coeffs (in this case) and we can find all the roots via fsolve:

restart;
Digits:=200:
aq := (2.626145*10^(-111)*I)*beta^41+(2.723460372*10^(-55)*I)*beta^25-1.125718*10^(-103)*beta^38-4.42696*10^(-96)*beta^36+(4.038976*10^(-119)*I)*beta^43+(1.897840*10^(-135)*I)*beta^47-1.4537*10^(-128)*beta^44-2.393897*10^(-75)*beta^30-5.345113*10^(-69)*beta^28-(8.88232*10^(-14)*I)*beta^7+(3.78162*10^(-127)*I)*beta^45+(1.87236321*10^(-39)*I)*beta^19-4.22943*10^(-111)*beta^40-5.98764*10^(-82)*beta^32+(1.73215*10^(-27)*I)*beta^13+(3.14576*10^(-144)*I)*beta^49+1.0000002*10^(-150)*beta^50+1.000483*10^(-142)*beta^48+(0.5707492e-4*I)*beta-1.860356732*10^(-56)*beta^26-1.202764308*10^(-50)*beta^24+0.1078870970e-3*beta^4-.1337634356*beta^2+(4.558807*10^(-82)*I)*beta^33+11.99907662+(1.552482870*10^(-49)*I)*beta^23+1.50289073*10^(-33)*beta^16+(4.738072*10^(-89)*I)*beta^35-2.560992731*10^(-45)*beta^22-1.821396189*10^(-40)*beta^20+1.*10^(-159)*beta^52-1.025045*10^(-118)*beta^42-2.708307310*10^(-27)*beta^14-2.880*10^(-137)*beta^46+3.616579272*10^(-22)*beta^12-2.775313578*10^(-17)*beta^10+(2.247453*10^(-75)*I)*beta^31+(4.317773*10^(-69)*I)*beta^29-(8.4742656*10^(-63)*I)*beta^27+(1.086756*10^(-103)*I)*beta^39+(2.875650*10^(-96)*I)*beta^37+(2.927689932*10^(-44)*I)*beta^21+(5.19084*10^(-18)*I)*beta^9+(3.3203077*10^(-35)*I)*beta^17-1.867365177*10^(-8)*beta^6+1.091287414*10^(-12)*beta^8-3.549248092*10^(-36)*beta^18-6.89128*10^(-89)*beta^34-(1.32011*10^(-22)*I)*beta^11+(8.67973*10^(-32)*I)*beta^15+(5.131768*10^(-10)*I)*beta^5-(6.362604*10^(-7)*I)*beta^3-(5.75387*10^(-153)*I)*beta^51:
f:=subs(beta=I*x,aq):
X:=fsolve(f,complex):
B:=X*~I:  # the solutions
seq( evalf(abs(eval(aq, beta=u))), u=[B]): #check
evalf[10]([%]);

[4.183842838*10^(-147), 4.183842838*10^(-147), 6.717604704*10^(-151), 6.717604704*10^(-151), 1.671175541*10^(-152), 1.671175541*10^(-152), 1.531062653*10^(-155), 1.531062653*10^(-155), 2.747698848*10^(-160), 2.747698848*10^(-160), 2.743673043*10^(-162), 2.743673043*10^(-162), 4.127466023*10^(-185), 4.127466023*10^(-185), 8.756351146*10^(-191), 8.756351146*10^(-191), 1.390759103*10^(-192), 1.390759103*10^(-192), 6.530643371*10^(-123), 6.530643371*10^(-123), 1.269161412*10^(-194), 1.269161412*10^(-194), 2.849455036*10^(-195), 2.849455036*10^(-195), 6.407673697*10^(-198), 6.407673697*10^(-198), 7.870992411*10^(-199), 7.870992411*10^(-199), 2.957280721*10^(-198), 2.957280721*10^(-198), 7.563118735*10^(-195), 7.563118735*10^(-195), 3.655881493*10^(-194), 3.655881493*10^(-194), 2.802759752*10^(-192), 2.802759752*10^(-192), 2.265640916*10^(-190), 2.265640916*10^(-190), 1.052966328*10^(-184), 1.052966328*10^(-184), 7.249928057*10^(-162), 7.249928057*10^(-162), 1.960460423*10^(-158), 1.960460423*10^(-158), 1.691321460*10^(-153), 1.691321460*10^(-153), 3.687720329*10^(-151), 3.687720329*10^(-151), 7.453401398*10^(-149), 7.453401398*10^(-149), 1.172077955*10^(-146), 4.492662674*10^(-8)]

The sequence of all the roots is in B.

vx:=[1, 3, 5, 10, 15];
                       [1, 3, 5, 10, 15]
vy:=[2.58, 8.57, 15, 30.10, 45];
                  [2.58, 8.57, 15, 30.10, 45]
f:=Statistics:-LinearFit([1,x], vx, vy, x);
      -0.3817701863353998 + 3.034083850931677 x
p1:=plot(vx,vy, style=point,symbolsize=20,color=red):
p2:=plot(f,x=0..16):
plots:-display(p1,p2,gridlines);

So, you want to generate an n x n random unitary matrix, k (<n) columns being fixed.
You could do the following:
1. Generate an n x k matrix A with the fixed columns.
2. Generate an  n x (n-k) random matrix B.
3. C := <A | B>   (=square matrix)
4. Check that det(C)<>0; if not, goto 2.
5. Apply the Gram-Schmidt process to C resulting a matrix G; it will be unitary.
6. Move the first k columns of G to their original positions. This is the desired matrix.
 

with(plots):  with(plottools):
 H:=(20/3)*sqrt(3.):
 Sph := sphere([0, 0, 0], 10,   color=gold, transparency=0.4):
 Cyls := display([seq(cylinder([0,0,h], sqrt(100-h^2),2*abs(h), color=red, strips=50),
     h=-10 .. -H/2.0, 0.2)], insequence):
 display([Sph, Cyls], style=surface, lightmodel=light1, axes=none);

dsolve({deqv, v(0) = v__0}, v(s)) assuming k>0, v__0>0;

A space filling curve very easy to implement (and understand) is due to Lebesgue and simplified by Schoenberg:

F:=proc(t)
   local u:=abs(2*round(t/2)-t);
   piecewise(u<=1/3,0, u<2/3, 3*u-1,1)
end:
PN:=4:
PX:=proc(t)  add(F(9^k*t)/2^(k+1),  k=0..PN) end:
PY:=proc(t)  add(F(9^k*3*t)/2^(k+1),k=0..PN) end:

plot(PX,0..1,numpoints=1000);

plot([PX,PY,0..1], numpoints=10000);

(Of course the real space filling curve is for PN=infinity).

If you want the numeric result, use:

int(f,-1..1,numeric);
   1.177418042
For a symbolici result, Markiyan suggestion should work but it gives a wrong answer (Maple's bug!):

g := proc(x) local r; r := sin(x)+x*cos(x); if is(abs(r) < 1/2) then sin(x) else cos(x) end if; end proc;

int(g, -1..1);evalf(%);
    2*sin(1)
   1.682941970

A piecewise attempt also fails for symbolics!

Workaround for a symbolic solution:

x0:=RootOf(2*_Z*cos(_Z)+2*sin(_Z)-1,0..1):
int(cos(x),x=-1..-x0)+int(sin(x),x=-x0..x0)+int(cos(x),x=x0..1);
     2*sin(1)-2*sin(RootOf(2*_Z*cos(_Z)+2*sin(_Z)-1, 0 .. 1))
evalf(%);
    1.177416231

EDIT. The value of the minimum given by minimize is of course wrong (bug!).

loc:=%[2][][1];
               {x = 3405.798701, y = 2169.889609}
eval([eq1,eq2,eq3],loc);
                       [0.01, 0., -0.01]

P.S. You may be interested in  ?geometry,RadicalCenter

EQ:=convert(lhs(eq1),set) union convert(lhs(eq2),set):
sol:=solve(EQ):
eval([fn1,gn1],sol);

N.B. You should avoid evalm and use "+", ".".

Use:

[A$2] .~ LM;

or

(A$2) .~ LM;

GetCoeff:=proc(f,u)
  local x:=indets(f,name);
  if u=0 or u=1 then return eval(f,x=~0) fi; 
  eval(foldl(coeff,f,op(u)), x=~0);
end proc: