The cure is simple: define ff1 etc via unapply, e,g.

ff1 := unapply( display([line([0, 0, 0], [d[1][1], d[1][2], d[1][3]], colour = red, thickness = 4)]),  t);

Your code works. The only problem is that a,b,c,x3  are complex (in general) so in printf you must replace for example  %2.8f  with %2.8Zf  (and without spaces!).

You must assign ranges for all the variables (or for none).
You may set of course e.g.  A0 = -infinity .. infinity  if you have not a better range.

It is unlikely that the print occurs only at the end, unless a procedure consumes almost all of the running time.
Maybe you must scroll the page while running in order to see the messages.

Anyway, you can set e.g.
printlevel:=10;

before running the procedure and see exactly what it is doing. If necessary increase slightly the printlevel value.

SmithForm does not accept polynomials with complex coefficients.
But there exists LinearAlgebra:-Generic:-SmithForm which works in any Euclidean domain. In order to use it in C[z], some simple procedures must be provided by the user (see the help page). For example:

C[EuclideanNorm] := a -> degree(a,z):
C[Gcdex] := (a,b,s,t)->gcdex(a,b,z,s,t):
...

After that, you can use:
LinearAlgebra:-Generic:-SmithForm[C](A, output=['S','U','V']):

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).