length(tekssifer);

All operations can be done with a one-line procedure if you're willing to input the operations in prefix form. The previous two Answers also require prefix form anyway, but my prefix form is identical to the original operations.

`&B`:= (ex::uneval)-> convert(eval(map(convert, ex, decimal, 2)), binary):

Your examples:

&B `+`(1101, 111);

     10100

&B `-`(11000, 1011);

     1101

&B `*`(1011, 1101);

     10001111

Integer division is done in Maple with the command iquo.

&B iquo(10010011, 1011);

     1101

Here is your program with corrected syntax and plot of U:

restart:
H:= 0.5:
q__0:= 10^5:
Es:= 4*10^9:
p__0:= 10^6:
q__s:= 10^10:
C:= q-> q^(H+1.5):
G:= q-> int(qs^3*C(q), qs= q__0..q):
P:= q-> 1/sqrt(G(q)):
p:= xi-> p__0/P(xi*q__0):
w:= (q,xi)-> 1/(Pi*(int(C(qs)*qs^3, qs= xi*q__0..q)))^(1/2):

U:= xi-> int(q*C(q)*w(q,xi)*Int(exp(-(w(q,xi)*ps/Es)^2)/ps, ps= p(xi)..infinity), q= xi*q__0..q__s):

plot(U, 1.5..5);

How do you expect the plotting routines to deal with the imaginary parts? If you use Re to extract the real part of your function, then you can get a plot. Like this:

ex1:= (x,t,z)-> Re(-1.132e11*exp(9.9e6*x + I*(1.95e6*z-2.98e15*t)));
plots:-implicitplot3d(
     ex1(x,t,z), x= -10..0, t= 0..10, z= 0..10,
     axes= boxed, style= patchcontour, scaling= constrained, shading= z,
     grid= [20$3]
);

PS: Because your function can never be 0, as Preben pointed out, the plot produced by the code above is just random noise.

The following recursive procedure solves your dice-sum problem. It allows for any number of dice, all possibly different, each with any number of sides, each with anything printed on the sides. The only requirement is that each die be "fair", having the same probability of landing on each of its sides. Each die is represented as a list of the values on its sides.

DiceSum:= (S, D::list(list))-> DiceSumRec(S,D) / `*`(nops~(D)[]):

DiceSumRec:= proc(S,D)
option remember;
local k;
     `if`(nops(D)=1, `if`(S in D[1], 1, 0), add(thisproc(S-k, D[2..]), k= D[1]))
end proc:

For example, the probability of getting a sum of 10 when rolling three standard six-sided dice is

DiceSum(10, [[$1..6] $ 3]);

Several points to help you along:

1. The output of Pollard's p-1 algorithm, when it works, is simply a proper factor of n. It may be composite or prime.
2. The algorithm only works if n has a prime factor p such that p-1 has only relatively small prime factors. No practical cryptography system uses such n. "Relatively small" means less than B as used in my procedure below. If B is set to n^(1/6), the probability that the algorithm will work is 1/27 (see the Wikipedia article).
3. Detecting primality is much easier than finding a factor. The algorithm shouldn't be used to test primality. The algorithm should reject input that is prime.
4. Since you use 2 as the base of your exponentations, your procedure needs to reject even input.
5. To find the gcd of integers, use igcd rather than gcd. The latter will work, but it is slower as it is intended for polynomials.
6. To do a modular exponentiation with a large exponent, use a &^ b mod n rather than a^b mod n. The former is much faster.

Here's my Pollard procedure, based on the Wikipedia article Pollard's p-1 algorithm:

Pollard:= proc(
     n::And(posint, odd, Not({identical(1), prime})),
     {B::posint:= max(trunc(n^(1/4)), 10^4)}
)
description "Pollard's p-1 algorithm";
local q:= 1, lnn:= evalf(ln(n)), g, a:= 2;
     while q <= B do
          q:= nextprime(q);
          a:= a &^ (q^trunc(lnn/ln(q))) mod n;
          g:= igcd(a-1,n);
          if g > 1 then return g end if
     end do;
     FAIL
end proc;

Do you mean this?

plots:-display(
     plots:-intersectplot(x+y=1, x^2+y^2+z^2=1, x= -1..1, y= -1..1, z= -1..1),
     plot3d(1, theta= 0..2*Pi, phi= 0..Pi, coords= spherical),
     scaling= constrained, axes= boxed
);

Simply use coeff. There's nothing special about it being a Lie derivative; it's treated as any other sum of terms.

coeff(L1fh, D_x);

You should use explicit multiplication signs and not rely on juxtaposition implying multiplication.

int(1/(x*(1-x))*(1+2*x), x);

1/2*u^2*(1-u)^2*(1+2*u)^2;

Now the integration is performed, and the second expression is interpretted correctly.

teksbiasa:= `Hello! Bob`:
ListTools:-Flatten(
     ListTools:-Pad[1](
          (w-> w+~nops(w))~(convert~(StringTools:-Split(teksbiasa, ` `), bytes)),
          32
     )[..-2]
);

The following plots the first 5122 zeros (the plot in the paper shows 10,000) in about 18 minutes. To get more zeros, expand the imaginary range, which I set at -9999..9999. The zeros are rather evenly distributed throughout the imaginary range.

R:= [RootFinding:-Analytic(add(1/n^s, n= 1..5), s, re= -3..1, im= -9999..9999)]:
nops(R);

     5122

To normalize the roots modulo 2*Pi*I/ln(5), as is done in the paper, use frem.

M:= evalf(2*Pi/ln(5)):  
Frem:= proc(x,y) local r:= frem(x,y); `if`(r<0, r+y, r) end proc:

One way to plot complex points (the easiest way as far as I'm concerned) is to separate them into real and imaginary parts with Re and Im.

P:= map(z-> [Re(z), Frem(Im(z),M)], R):
plot(P, style= point, symbol= point, symbolsize= 1);

Sure, it's easy:

sumif:= (L::{list,set}, B)-> `+`(select(B, L, _rest)[]):

Now let's say that I have a list:

R:= rand(9):  L:= ['R()' $ 9];

And I want to sum all those elements that are greater than 5:

sumif(L, `>`, 5);

 30

Another way: If L is the list of digits, do

(parse@cat@op)(L);

ListTools:-Reverse(convert(mylist, base, 10));

The reverse of [2,7,5] is [5,7,2]. You have it as [5,7,5].