What is p ? It's not hard to code parts of your algorithm, but I can't run it unless I know what p is. Also what is n in gamma(n-p) ? It is important to know whether or not these values change as you run the algorithm. I'm going to assume that they don't.
Here are a couple of hints for how I would do this. Make procedures for alpha, epsilon, and gamma. You can't actually use the word gamma however since Maple thinks that is the constant gamma=0.577... I will use "gam" instead. Use the remember option in the procedures so that you will store computed results (and won't recompute them). This will automatically give you good asymptotic complexity.
alpha := proc(i)
option remember;
`if`( alpha(i-1) + epsilon(i-1) + epsilon(i-2) <= 1, 0, 1)
end proc:
alpha(0) := 1: # store initial value for alpha
gam := proc(i)
option remember;
`if`(i <= p, 0, epsilon(n-p)); # should this be epsilon(i-p) ?
end proc:
Now I'm not sure about epsilon because the following facts seem contradictory:
epsi(i+1)=epsi(i) # this would be easy to compute :)
epsi(i) = gamma(i) - epsi(i-1) - alpha(i) (mod 2) # we could use this
# this seems to define epsilon for a whole bunch of values
# however it contradicts epsi(i+1)=epsi(i)
eps(i) = 0 if i=2k where k is a natural number less than or equal to p
eps(i) = 1 if i=2k-1 where k is a natiural number and i is less than or equal to p
Anyways, it shouldn't be too hard to write this program. Even with a bad algorithm you should be able to get up to i=1000. If you need very high values however then you will have to look closely at how epsilon is computed and try to predict what certain intermediate values will be.

What is p ? It's not hard to code parts of your algorithm, but I can't run it unless I know what p is. Also what is n in gamma(n-p) ? It is important to know whether or not these values change as you run the algorithm. I'm going to assume that they don't. Here are a couple of hints for how I would do this. Make procedures for alpha, epsilon, and gamma. You can't actually use the word gamma however since Maple thinks that is the constant gamma=0.577... I will use "gam" instead. Use the remember option in the procedures so that you will store computed results (and won't recompute them). This will automatically give you good asymptotic complexity. alpha := proc(i) option remember; `if`( alpha(i-1) + epsilon(i-1) + epsilon(i-2) <= 1, 0, 1) end proc: alpha(0) := 1: # store initial value for alpha gam := proc(i) option remember; `if`(i <= p, 0, epsilon(n-p)); # should this be epsilon(i-p) ? end proc: Now I'm not sure about epsilon because the following facts seem contradictory: epsi(i+1)=epsi(i) # this would be easy to compute :) epsi(i) = gamma(i) - epsi(i-1) - alpha(i) (mod 2) # we could use this # this seems to define epsilon for a whole bunch of values # however it contradicts epsi(i+1)=epsi(i) eps(i) = 0 if i=2k where k is a natural number less than or equal to p eps(i) = 1 if i=2k-1 where k is a natiural number and i is less than or equal to p Anyways, it shouldn't be too hard to write this program. Even with a bad algorithm you should be able to get up to i=1000. If you need very high values however then you will have to look closely at how epsilon is computed and try to predict what certain intermediate values will be.

Try holding shift and pressing up or down to select the line, then hit delete. Alternatively you could select it with the mouse.

I do this all the time. First you should write g as a polynomial, ie:
g := 5*x^2+8*y^2+5*z^2-4*y*z+8*x*z+4*x*y;
If you have it as a function then just use g(x,y,z) in place of g.
C := [coeffs(g, {x,y,z}, 'M')]:
C := table(zip(`=`, [M], C), 'sparse');
Now C[x^2] will give you 5, etc. and any monomial which does not appear will give you zero because of the 'sparse' option. The first line puts the coefficients of g into a list and puts the corresponding monomials into M (as a sequence). The zip function merges two lists, in this case using equals, so the result of zip(`=`, [M], C); is the list [x^2=5, y^2=8, ...]; If you just want that list, then stop there. If you want to access coefficients use a table.

I like using Maple because it has a nice programming language. The system itself is pretty easy to use, and you can get up to speed in a couple of hours. The manual used to have a very nice introduction to the system, showing you how to write programs and plot functions and the like, however I couldn't find it in the current version of the manual. In any case, you can download the manuals here (you have to register): http://maplesoft.com/products/maple/manuals/
Unfortunately there is no trial version of Maple 10 available, but they do have some animations on the site.
Mathematica has a trial version, I suggest you download it and try it out. Personally I find it a bit harder to use than Maple, with less mathematical functionality in the areas I am interested in. It really depends on what you're doing however, since Mathematica is stronger than Maple in some areas, weaker in others.
I've never used Matlab, however it is the standard tool for numerical work.

Yes, the result of any logical expression containing FAIL is FAIL.

It is strange to see this problem, but you should force Maple to do a garbage collection using the gc() command at the end of the loop. This will clear any intermediate objects created during your computations and trash any remembered results.

make the polynomial a function:
p := (x,y) -> x^2+y;
Now try p(1,2), and for that matter, p(a,b);

You need to make each point a list, ie: points := { [1,1], [3,4], [1,5], [0,-4], [-3,1] }; # set of points hull := simplex[convexhull](points); P := plots[pointplot](points): # draw the points Q := plots[pointplot]([op(hull), hull[1]], connect=true): # draw the convex hull plots[display](P,Q); # display them both together Note that I had to repeat the first point in the convex hull in order to close the polygon. Try [op(hull), hull[1]]; to see what I mean.

Maple gets a little fussy when you have variables as exponents. Try using the expand and factor commands. f := x^(n+r)*(n+r-1)*c[n]+(n+r)*x^(1+n+r)*c[n]; factor(expand(f)); # try this Interestingly enough, factor(f) won't give you the same answer.

I would suggest turning on antialiasing, but I'm not sure if Maple 9.5 supports it. It would be in the Options dialog under the Display tab.

I noticed this too and reported the bugs to Maplesoft. One problem is that when you create a new section it does not automatically make a new "Document Block". You can use Format->Create Document Block to correct this manually whenever you create a new section.

You can't really do this. The closest thing is MapleNet, which will let people run your Maple worksheets in a web browser over the internet.

Well, you can fill the matrices with symbolic entries (a11, a12, etc.) and use LinearAlgebra:-LinearSolve, but this will be slow.
What you really want is

Cramer's Rule. I'm not sure if there's a Maple command which implements it.

Just to clarify, do you need a primitive element to generate GF(2^32) ? Or do you want to test whether polynomials with coefficients in GF(2^32) are irreducible ?