If you have a finite sum such as  s(x) = a[0] + a[1]*x + ... + a[9]*x^9
it is not possible to find a (unique) function whose Taylor (or other) series starts  with s(x).
[A similar problem is to find "the" irrational number for which the first decimals are 1.234]

Using your random generators,  the equality fails in about 2% of the tests.

restart;
with(RandomTools):
NN := 1000:
err := 0: 
to NN do 
d := Generate(integer(range = 1 .. 120)); 
T := Generate(integer(range = 1 .. d)); 
B := Generate(integer(range = 1 .. d)); 
Y := Generate(integer(range = 0 .. B)); 
Pr := Generate(integer(range = 0 .. d-T)); 
py := Generate(integer(range = 0 .. min(Y, d-T, Pr))); 
e1 := binomial(Y, py)*(binomial(d-Y, Pr-py)*binomial(d-Pr, T)-binomial(d-B, T)*binomial(d-T-Y, Pr-py))/(binomial(d-T, Pr)*(binomial(d, T)-binomial(d-B, T))); 
e2 := add((binomial(Y, n)*binomial(d-Y, T-n)/(binomial(d, T)-binomial(d-B, T))+n-max(n, binomial(d-B, T)/(binomial(d, T)-binomial(d-B, T))))*binomial(Y-n, py)*binomial(d-T-Y+n, Pr-py)/binomial(d-T, Pr), n = 0 .. Y); 

if e1 <> e2 then err := err+1; print(e1 <> e2, 'd'=d, 'T'=T, 'B'=B, 'Y'=Y, 'Pr'=Pr, 'pr'=py) end if end do: 
fails=err,'NN'=NN;

It is not a good idea to use floats in symbolic computations (if you can avoid them).

Replacing 0.2 by 1/5 (or convert( ...,rational) ), everything is OK.

Alternatively, increase Digits.

with(plots):
with(plottools):
Body := display([
point([1, -.2, 0], colour = red, symbolsize = 50, symbol = solidsphere), 
point([-1, -.2, 0], colour = red, symbolsize = 50, symbol = solidsphere), 
point([0, 1, 0], colour = purple, symbolsize = 50, symbol = solidsphere), 
line([0, 0, 0], [1, -.2, 0]), 
line([0, 0, 0], [-1, -.2, 0]), 
line([0, 0, 0], [0, 1, 0]), line([0, 0, 0], [.5, 0, 0], colour = red, thickness = 4), 
line([0, 0, 0], [0, .5, 0], colour = yellow, thickness = 4), 
line([0, 0, 0], [0, 0, .5], colour = blue, thickness = 4)], scaling = constrained):
w1 := .35: w2 := .45: w3 := .15:
animate( rotate, [Body,w1*t,w2*t,w3*t], t=0..8*Pi,frames=100 );

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