There are infinitely many solutions.

solve(z^(1+I) = 1, z, allsolutions):
evalc(%);

So, you don't need evalf.

BTW, seq(print(...),...)  is strange (and the output is wrong; what version are you using?). Why not:

interface(rtablesize=infinity):
Vector([seq(evalf(evalf[d](u)), d = 10 .. 3010, 300)]);

I think that readstat being  an I/O command, it may need time to complete before calling it again.
Inserting some Sleeps makes it work.
dt depends on the speed of the computer.

dt:=0.5:
fun:=proc(n)  local a, b, i, x;
a:=readstat("insert a");
Threads:-Sleep(dt);
b:=readstat("insert b");
for i to n do
Threads:-Sleep(dt);
x[i]:=readstat("insert x");
od;
end proc:
fun(3);

sqrt (the principal branch) is af course well defined, the branch cut is arbitrary (-oo,0); it could be any "arc" from 0 to oo.

But it is possible to define a custom sqrt. Here is one with a branch cut along negative imaginary axis.

restart;
mysqrt:=z -> piecewise(argument(z)>-Pi/2,sqrt(z),-sqrt(z)):
 

limit(mysqrt(-1-t*I), t=0, left), limit(mysqrt(-1-t*I), t=0, right); #continuous
                              I, I
limit(sqrt(-1-t*I), t=0, left), limit(sqrt(-1-t*I), t=0, right);     #discontinuous
                             I, -I
limit(mysqrt(-I-t), t=0, left), limit(mysqrt(-I-t), t=0, right);     #discontinuous

   
limit(sqrt(-I-t), t=0, left), limit(sqrt(-I-t), t=0, right);        #continuous        

(1/2)*sqrt(2)-I*sqrt(2)*(1/2), (1/2)*sqrt(2)-I*sqrt(2)*(1/2)

simplify(convert(z, arcsin)) assuming x>-a,x<a;

works in Maple 2017.

Download loss.mw

Compilable procedures are very restrictive. See ?compile.

- Arrays can only appear as parameters of a procedure and must have datatype=float or integer
- No symbolics; for example (in your case) int(...), orthopoly[...](...)  are out of the question. Also W(...) is not allowed.

The symbol O has nothing special, except it is protected,  to prevent it from being modified.

letters:=map(convert,[seq("A".."Z"), seq("a".."z")], symbol);
    letters := [A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R,  S, T, U, V, W, X, Y, Z, a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z]
select(type, letters, protected);
    [D, O]

Edit. The symbol I (the imaginary unit) is a special one, not detected above.

It is strange that convert(...,tanh) does not work as expected.
For example, convert(sinh(x), tanh)  works but convert(exp(x), tanh) does not.

`convert/exp2tanh`:=proc(ex)
  simplify(eval(ex, exp = (t -> (tanh(t/2)+1)/(1-tanh(t/2)))))
end;

convert(exp(3*x),exp2tanh);

       

Expr:=-(exp(-a*s)-1)*kw/((exp(-a*s)+1)*s^2):

convert(Expr,exp2tanh);

PerfectRandom := proc()
  return 4  # Chosen by fair dice roll, guaranteed to be random
end proc:

Example:

PerfectRandom();
        4

Output:=Array(-10..10,-10..10, (i,j)->i^2+j^2):
plots:-surfdata(Output);

Let f be a real function defined on [a,b].

If f has an antiderivative F (on [a,b]) then F is by definition differentiable, hence continuous.
But f could have infinitely many discontinuities and could be non-integrable (Riemann or even Lebesgue).

Note that there exist other (generalized) integrals such as Denjoy; in this context f becomes integrable
and Newton-Leibniz applies. However the Denjoy integral is not usually considered even at a university level.

The simplest example of a function not integrable (Riemann or Lebesgue) but having an antiderivative
is  f(x) = F'(x), F(x) = x^2*sin(1/x^2), x∈[-1,1]∖{0}, F(0) = 0. Note that f is discontinuous at 0.

(Using Rouben's suggestion.)

fn:=cat(kernelopts(homedir),"/Elsymfun.txt");
save(`convert/elsymfun`,fn):  close(fn):
sf:=readbytes(fn,infinity,TEXT):
sf:=StringTools:-Subs([") local vars," = ", vars::set(name)) local ", "elsy"="Elsy", "vars := indets(p);"=""], sf): 
close(fn):writebytes(fn,sf):close(fn):read(fn):

u := a*(x+y) + b*(x^2+y^2) + a*b*(x^3+y^3):
convert(u, elsymfun);
      Error, (in convert/elsymfun) polynomial is not symmetric
convert(u, Elsymfun, {x,y});
       a*b*(x+y)^3-3*a*b*(x+y)*x*y+b*(x+y)^2+a*(x+y)-2*b*x*y

Maybe something like this.

imagefile := cat(kernelopts(datadir), "/images/fjords.jpg"):
plot3d(1-x^2-y^2,x=-1..1,y=-1..1, image = imagefile);

Note that the image file can be obtained from any plot/plot3d.

restart;
select(u->type(convert(u,symbol),procedure), convert([seq(33..127)],compose,bytes,list));

        ["$", "*", "+", "-", ".", "/", "<", "=", ">", "@", "D", "^", "~"]