Here are a pair of procedures for it. The procedure Sumof2Squares computes the first square (the 259 in your example) then loops with procedure FermatDescent to lower the coefficient (the 34 in your case). By using trace, we can observe the iterations.

Notice that I use mods(..., m) rather than the more-usual (... mod m). This guarantees that the residue will be between -m/2 and m/2. This is needed to guarantee that the coefficient is reduced at each iteration and thus ends at 1.
 

Download FermatDescent.mw

There are several algorithms for this sum-of-two-squares problem. Fermat Descent is not the fastest, but it's not horribly slow either.

convert(exp(x^2), compose, FormalPowerSeries, GAMMA);

The issue is operator precedence, as VV pointed out. As usual, the precedence can be overridden by using parentheses:

(3 = 0) mod 3;

When doing pure integer arithmetic, it's better to use irem rather than mod:

irem(3,3);
irem(0,3);

A simple linear transformation (x-> 2*x-1) of 0..1 makes it -1..1.

RM:= Matrix((10$2), 2*rand(0..1) - 1);

Here I've performed arithmetic directly on the procedure rand(0..1), which produces a new procedure. This is different than

RM:= Matrix((10$2), ()-> 2*rand(1) - 1);

although both produce a matrix of uniformly distributed -1s and 1s.

It's not hosted by Maplesoft, but I highly recommend the question & answer site math.stackexchange.com

Some general mathematical discussions are welcome here on MaplePrimes. If your questions haven't been answered, it's because no one knows what to say, not because we think that your questions are irrelevant to this site. So, no offense is intended by a failure to answer. If material is truly irrelevant, it gets deleted or its irrelevance is quickly commented on.

In Maple 12, you can sort numeric (i.e., decimal) values by magnitude by

sort([Y], (x,y)-> abs(x) < abs(y));

In recent versions of Maple, this can be done more efficiently by

sort([Y], key= abs);

Like this

Pts:= [[0, 438.912], [0.0295, 402.336], [0.125, 365.76], [0.2009, 341.376]]:
plot(
   [Pts$2], style= [point, line], symbol= diamond, symbolsize= 12, 
   color= blue, gridlines, view= [0..0.25, 0..500]
);

Using printlevel is a crude tool for debugging only. It's impossible to get it to display only the information that you want. Thus it is never intended for the display the output of a finished procedure.

If you have a procedure that generates plots which are not the official return value of the procedure, but you want to see them anyway, the use print in the procedure, as in

print(plots[display]( ...your exact same display command here...));

The point-data matrices can be extracted from a plot like this:

P:= implicitplot(...your exact same command...);
(M1,M2):= map2(op, -1, [plottools:-getdata(P)])[];

Now M1 and M2 will both be two-column matrices, the first for disf and the second for dispbeamf.

How about

convert(simplify(eval(diff(f(x),x)*cos(f(x)), x= RootOf(f(x),x))), D);

Use eliminate instead of solve. The second set returned is the conditions (which must be understood to be equated to 0).

You're asking for a 3D plot, yet nowhere in your command do you specify which functions are to be represented on each of the three axes.

If you want all decimal numbers to be displayed with 3 digits after the decimal point (or decimal comma, if that's what your culture uses), give the command

interface(displayprecision= 3);

This can also be set from the Tools -> Options menu.

Factoring polynomials over the rationals is a completely solved problem, whose solution has been completely implemented in Maple for a great many years, at least 17. If you mean to express the given polynomial as a sum of some small number (but greater than one) of factored polynomials, that's a different problem that's not called "factoring" strictly speaking.

If you have an example of a rational-coefficient polynomial that Maple can't completely factor, please report that as a bug, or post it here.

You can see a plaintext form of your most-recent output by the command lprint(%).