OK, that's a different question entirely.  If alpha is a root of a polynomial p(z) of degree d over a field F, the sequence alpha^n satisfies a linear recurrence of order d corresponding to the polynomial.
If X(k) is the vector with entries alpha^(j+k), j = 0 to d-1, we have X(k) = M^k . X(0) where M is a certain matrix.   In this case you want to look at the first row of M^k mod 2 (mod 2 because the field has characteristic 2) where k = 12345678987654321.
You can get that using LinearAlgebra[Modular][MatrixPower].  Thus:

> p:= x^32+x^26+x^19+x^15+x^13+x^11+x^9+x^8+x^4+x+1;
  M:=LinearAlgebra[CompanionMatrix](p,x)^%T;
  Mp:= LinearAlgebra[Modular][MatrixPower](2, M, 12345678987654321);
  add(Mp[1,i]*alpha^(i-1),i=1..32);

1+alpha^2+alpha^3+alpha^9+alpha^13+alpha^14+alpha^15+alpha^16+alpha^17+alpha^18+alpha^19+alpha^20+alpha^24+alpha^25+alpha^28+alpha^30+alpha^31

Yes, sorry, that's what I meant.  I copied from a place where I had already extracted the first operand.

Yes, sorry, that's what I meant.  I copied from a place where I had already extracted the first operand.

This is a rather interesting equation.  It can be written as

> isolate(eq, sin(b^2));

sin(b^2) = 25*tan((1/18)*Pi)*(2273/125+289*cos(2*b)*(1/25))/(578*cot(b))+25/289

The right side is periodic in b with period Pi.  Except for narrow intervals around Pi/2 and its translates by Pi, its value will be between -1 and 1.  Thus for large values of n, there will be
many solutions in the interval [n*Pi, (n+1)*Pi); it is not hard to approximate how many such solutions (but it's late and I have to go to bed...)

This is a rather interesting equation.  It can be written as

> isolate(eq, sin(b^2));

sin(b^2) = 25*tan((1/18)*Pi)*(2273/125+289*cos(2*b)*(1/25))/(578*cot(b))+25/289

The right side is periodic in b with period Pi.  Except for narrow intervals around Pi/2 and its translates by Pi, its value will be between -1 and 1.  Thus for large values of n, there will be
many solutions in the interval [n*Pi, (n+1)*Pi); it is not hard to approximate how many such solutions (but it's late and I have to go to bed...)

These expressions are not equivalent. 

exp(-I*Pi) = -1, so  (exp(-I*Pi))^(1/2-(1/100)*(I*8)) = (-1)^(1/2 - 2*I/25)

Now Maple uses the principal branch for complex powers: z^w = exp(w*Log(z)) where Log is the principal  branch of ln(z), i.e. the one with imaginary part in (-Pi, Pi].  Log(-1) = I*Pi, and so
(-1)^(1/2-2*I/25) = exp(I*Pi*(1/2 - 2*I/25)) = exp(2*Pi/25 + Pi*I/2) = I*exp(2*Pi/25).  This is indeed approximately 1.285730980*I.

On the other hand,

exp(-I*Pi*(1/2)) = -I

exp(-8*Pi*(1/100)) = exp(-2*Pi/25)

exp(-I*Pi*(1/2))*exp(-8*Pi*(1/100)) = -I*exp(-2*Pi/25), which is indeed approximately
-0.7777676792*I.
 

The point is that you have to be careful with familiar "identities" such as (a^b)^c = a^(b*c) and (a*b)^c = a^c * b^c, which are not always true when complex numbers are involved. 

These expressions are not equivalent. 

exp(-I*Pi) = -1, so  (exp(-I*Pi))^(1/2-(1/100)*(I*8)) = (-1)^(1/2 - 2*I/25)

Now Maple uses the principal branch for complex powers: z^w = exp(w*Log(z)) where Log is the principal  branch of ln(z), i.e. the one with imaginary part in (-Pi, Pi].  Log(-1) = I*Pi, and so
(-1)^(1/2-2*I/25) = exp(I*Pi*(1/2 - 2*I/25)) = exp(2*Pi/25 + Pi*I/2) = I*exp(2*Pi/25).  This is indeed approximately 1.285730980*I.

On the other hand,

exp(-I*Pi*(1/2)) = -I

exp(-8*Pi*(1/100)) = exp(-2*Pi/25)

exp(-I*Pi*(1/2))*exp(-8*Pi*(1/100)) = -I*exp(-2*Pi/25), which is indeed approximately
-0.7777676792*I.
 

The point is that you have to be careful with familiar "identities" such as (a^b)^c = a^(b*c) and (a*b)^c = a^c * b^c, which are not always true when complex numbers are involved. 

Here it an example of  Euler's equations with odeplot.

>  i:= [1,2,3]; M:= [3,4,2];
   des:= seq(i[j]*diff(omega[j](t),t) + 
    (i[1+(j-2 mod 3)]-i[1+(j mod 3)])*omega[1+(j-2 mod 3)](t)
    *omega[1+(j mod 3)](t) = M[j], j=1..3);
   sol:= dsolve({des,omega[1](0) = 0, omega[2](0)=0, omega[3](0)=0},
     numeric);
   plots[odeplot](sol,[omega[1](t),omega[2](t),omega[3](t)],t=0..10,
     numpoints=1000,axes=box);

Here it an example of  Euler's equations with odeplot.

>  i:= [1,2,3]; M:= [3,4,2];
   des:= seq(i[j]*diff(omega[j](t),t) + 
    (i[1+(j-2 mod 3)]-i[1+(j mod 3)])*omega[1+(j-2 mod 3)](t)
    *omega[1+(j mod 3)](t) = M[j], j=1..3);
   sol:= dsolve({des,omega[1](0) = 0, omega[2](0)=0, omega[3](0)=0},
     numeric);
   plots[odeplot](sol,[omega[1](t),omega[2](t),omega[3](t)],t=0..10,
     numpoints=1000,axes=box);

What version of Maple are you using?  My code should work in recent versions.
There really isn't such a thing in Maple as "plotting points on the same graph sequentially".  Maple plots a graph all at once.  If you want to create the illusion of plotting sequentially, you can use animation as I did.

What version of Maple are you using?  My code should work in recent versions.
There really isn't such a thing in Maple as "plotting points on the same graph sequentially".  Maple plots a graph all at once.  If you want to create the illusion of plotting sequentially, you can use animation as I did.

For some reason, when solution curves are to be included, even though the constraint should not have any effect, Maple thinks there's an imaginary component.  You might combine a dfieldplot (with the constraint) and a DEplot of solutions (with no arrows) using display in the plots package.

For some reason, when solution curves are to be included, even though the constraint should not have any effect, Maple thinks there's an imaginary component.  You might combine a dfieldplot (with the constraint) and a DEplot of solutions (with no arrows) using display in the plots package.

This does not diagonalize to Diag(2,1,0).  The characteristic polynomial mod 3 is t^3 + t^2, so 1 is not an eigenvalue. 
0 is an eigenvalue of algebraic multiplicity 2, but geometric multiplicity 1: the eigenvectors for eigenvalue 0 are all multiples of <1,1,1>, and v = <1,1,0> is a vector such that M^2 . v = 0 but M . v <> 0.  Thus this is an example that in nonzero characteristic, symmetric matrices may have nontrivial Jordan forms.

This does not diagonalize to Diag(2,1,0).  The characteristic polynomial mod 3 is t^3 + t^2, so 1 is not an eigenvalue. 
0 is an eigenvalue of algebraic multiplicity 2, but geometric multiplicity 1: the eigenvectors for eigenvalue 0 are all multiples of <1,1,1>, and v = <1,1,0> is a vector such that M^2 . v = 0 but M . v <> 0.  Thus this is an example that in nonzero characteristic, symmetric matrices may have nontrivial Jordan forms.