Advanced Search

Carl Love

Carl Love

28100 Reputation

25 Badges

13 years, 105 days
Himself
Wayland, Massachusetts, United States
Contact Carl Love
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

Not the right "catch string"...

Carl Love 28100
April 20 2015
0 0

@matmxhu 

Your "catch string" should not include the first part, "Error, (in fsolve/polyill) "; it should be simply "time expired". So try this:

g:= proc()  [solve(eq)] end proc:

try
     sols:= timelimit(300, g())
catch "time expired":
     sols:= SolveEquations(eq)
end try;

Doesn't check...

Carl Love 28100
April 20 2015
0 0

@testht06 

No, it doesn't check out. You can see for yourself simply by computing R^4 using the proposed fourth root as the R. 

R:= <0,1,0,0;
          0,0,1,0;
          0,0,0,1;
          1,1,x,x^2>:

Expand~(R^4) mod p;

R converted back to GF(2^8)...

Carl Love 28100
April 20 2015
0 0

@testht06 

You wrote:

Please help me convert the entries of R (and R^4) from GF(2^16) back to GF(2^8).

Okay, it's done, amazingly enough. And R has a very simple representation in the original field. This conversion back to the original field will not be possible in the general case. I had to construct an explicit field isomorphism element by element to do this. I could not come up with an algebraic trick for it. I hope that someone here can help me with that.

I think, GF(2^8) is the sub set of GF(2^16).

Yes, but that subset is not a subfield, in particular it is not the subfield of GF(2^16) which is isomorphic to our original GF(2^8). To see that such a situation is not even possible, consider this: y^7 is a member of that subset, but its square, y^14, is not. So the subset is not closed under multiplication and is thus not a subfield.

And while convert R and R^4 back to GF(2^8), the entries ≤ 255 will be not changed

No, they must be changed, for the reason in the preceding paragraph. The only elements that can stay the same are the members of the prime subfield---in this case just 0 and 1.

Furthermore, for the calculation of the roots of the factor (t2 +231t +30), can You choose an irreducible polynomial ext2 over GF(2^8)? (ie ext2 = t2 + at + b, where a, b in GF(2^8). GF((28)2 is the composition field)

Well t^2+231T+30 itself is irreducible, and that is what I chose, essentially. But the only way that Maple will work with a composite field is by "flattening" the representation from (2^8)^2 to 2^16 via Primfield.

Here is the new worksheet, which is mostly the same as the previous worksheet, but it has isomorphism at the bottom and the conversion of R back to the original field.



Fractional powers of matrices over finite fields

Author: Carl J Love, 2015-Apr-18.

> 

restart:

Define a field extension for GF(2^8) using an irreducible polynomial of degree 8.

> 

p:= 2:  ext1:= Z^8+Z^4+Z^3+Z+1:  d:= degree(ext1, Z):

Verify irreducibility:

> 

Factors(ext1) mod p;

[1, [[Z^8+Z^4+Z^3+Z+1, 1]]]

Coefficients in the new field are expressed as polynomials in abstract roots of the extending polynomial expressed via the RootOf function. We wish to write and to see this root as a simple variable rather than as a complicated RootOf function, so we use alias.

> 

alias(x= RootOf(ext1)):

We will extract the fourth root of the matrix below. First we will extract its eigenvalues in this field, and then we'll construct a splitting field to get the remaining eigenvalues.

> 

ROOT:= 4:
> 

M:= < 1,     x,             1,           x^2;
      x^2,   x^3+1,         x^2+x,       x^4+1;
      x^4+1, x^5+x^2+x,     x^4+x^3,     x^6+x;
      x^6+x, x^7+x^4+x^2+1, x^6+x^5+x^2, x^3+x+1 >:
> 

C:= Charpoly(M,t) mod p;

t^4+(x^3+x+1+x^4)*t^3+t^2+x^4*t+1

The field's elements can be compactly represented as integers by a natural base-p encoding scheme. This is for our viewing purposes only; it is not used for computation.

> 

encode:= (ex,x)-> eval(ex, x= p):
> 

encode(C,x);

t^4+27*t^3+t^2+16*t+1

> 

CF:= Factors(C) mod p;

[1, [[x^7*t+x^6*t+x^5*t+x^4+x^3+x^2*t+x^2+x*t+t^2+x+t, 1], [x^5+x^2+t+1, 1], [x^7+x^6+x^4+x^3+t+1, 1]]]

> 

encode(CF,x);

[1, [[t^2+231*t+30, 1], [37+t, 1], [217+t, 1]]]

So we have two eigenvalues in our field (represented as 217 and 37 in the encoding scheme), and two more that are the roots of the irreducible quadratic.

 

Extract the nonlinear irreducible parts (this command works even if there is more than one irreducible factor):

> 

irr:= map(collect, select(q-> degree(q,t) > 1, map2(op,1,CF[2])), t);

[t^2+(x^7+x^6+x^5+x^2+x+1)*t+x^4+x^3+x^2+x]

> 

encode(irr,x);

[t^2+231*t+30]

> 

q:= irr[1]:

Now construct the splitting field that contains q's roots.

> 

alias(z= RootOf(q)):
> 

PF:= Primfield({x,z}) mod p:

Extract the new degree 16 field extension:

> 

ext2:= (indets(PF, 'RootOf'(anything)) minus {x,z})[];

RootOf(_Z^16+_Z^13+_Z^10+_Z^9+_Z^2+_Z+1)

Note that ext2 is a randomly chosen irreducible (mod 2) polynomial of degree 16. If you execute the Primfield command again, you may get a different polynomial. The important thing is how the old field's entries are expressed in the new field. All these relationships are stored in PF. We compactify the visual representation of the relationships with a new alias.

> 

alias(y= ext2):
> 

PF;

[[y = z^15*x+z^14+z^11+z^8+z^4+z^3+z], [z = y^15+y^14+y^8+y^5+y^4+y^3+y^2, x = y^15+y^13+y^11+y^10+y^8+y^7+y^6+y^4+y^2+y]]

The first equation is not of much use to us here. One of the others expresses one of the roots of q,which we named z, in the new field. The other expresses the old field extension x in terms of the new, y.

> 

newx:= eval(x, PF[2]);

y^15+y^13+y^11+y^10+y^8+y^7+y^6+y^4+y^2+y

Now re-express q in the new field.

> 

newq:= collect(Expand(eval(q, x= newx)) mod p, t);

t^2+(y^13+y^12+y^11+y^9+y^8+y^7+y^5+y^4+y^2+y)*t+y^15+y^13+y^11+y^7+y^6+y^5+y^2+y+1

> 

encode(%,y);

t^2+15286*t+43239

Now get both roots:

> 

Rq:= Roots(newq) mod p;

[[y^15+y^14+y^8+y^5+y^4+y^3+y^2, 1], [y^15+y^14+y^13+y^12+y^11+y^9+y^7+y^3+y, 1]]

> 

encode(%,y);

[[49468, 1], [64138, 1]]

We can see that one of the roots is, of course, the z. Now verify that the second root is the first one raised to the 2^8 power:

> 

Expand(eval(z, PF[2])^(p^d)) mod p;

y^15+y^14+y^13+y^12+y^11+y^9+y^7+y^3+y

> 

encode(%,y);

64138

Re-express the original matrix in the newly extended field.

> 

newM:= Expand~(eval(M, x= newx)) mod p:

Now I need the Eigenvalues and Eigenvector procedure that I presented for an earlier Question.

> 

`mod/Eigenvalues`:= proc(M::Matrix(square), p::posint)
     local x;
     Roots(Charpoly(M, x) mod p) mod p
end proc:
> 

`mod/Eigenvector`:= proc(M::Matrix(square), lambda, p::posint)
#Returns an eigenvector of M associated with the eigenvalue lambda.
uses LA= LinearAlgebra;
local
     n:= op([1,1], M),
     rank, t,
     sol:= Linsolve(M - lambda*LA:-IdentityMatrix(n), Vector([0$n]), rank, t) mod p
;
     if rank=n then
          error "Second argument must be an eigenvalue of the first."
     end if;
     eval(sol, indets(sol, specindex(posint,t))[]= 1)
end proc:
> 

eigs:= Eigenvalues(newM) mod p;

[[y^15+y^14+y^8+y^5+y^4+y^3+y^2, 1], [y^15+y^14+y^13+y^12+y^11+y^9+y^7+y^3+y, 1], [y^15+y^14+y^11+y^10+y^9+y^7+y^2+y, 1], [y^13+y^12+y^11+y^9+y^8+y^7+y^2+y+1, 1]]

We know that this matrix will pass the following check, but the check will be needed in general.

> 

if nops(eigs) <> op([1,1],M) then
    error "Case of repeated or missing eigenvalues not handled"
end if;

Construct the matrix P of eigenvectors:

> 

P:= `<|>`(seq(Eigenvector(newM, lambda[1]) mod p, lambda= eigs)):

Compute the roots of the eigenvalues.

> 

ROOTs:= map(lambda-> Roots(t^ROOT - lambda[1]) mod p, eigs);

[[[y^15+y^14+y^10+y^9+y^6+y^5+y^4+y^2+y, 4]], [[y^15+y^14+y^11+y^10+y^5+y^4+y^3+y^2+y+1, 4]], [[y^10+y^9+y^8+y^4+y, 4]], [[y^14+y^13+y^11+y^9+y^8+y^3+y, 4]]]

Check that each eigenvalue has an appropriate root in the current field.

> 

if [] in {ROOTs[]} then
     error "Case of root of eigenvalue not existing in the current field not yet implemented"
end if:

Construct the diagonal matrix of the roots of the eigenvalues

> 

Diag:= LinearAlgebra:-DiagonalMatrix(map(x-> op(1,op(x)), ROOTs)):

Construct the "root" matrix

> 

R:= Expand~(P.Diag.(Inverse(P) mod p)) mod p;

R := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 1, (4, 1) = 1, (4, 2) = y^15+y^13+y^11+y^10+y^8+y^7+y^6+y^4+y^2+y, (4, 3) = 1, (4, 4) = y^14+y^13+y^10+y^9+y^6+y^4+1})

Verify that R is indeed the root of the starting matrix.

> 

Expand~(R^ROOT - newM) mod p;

Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0})

We wish to convert---to whatever extent possible---the entries of these matrices from their representations in GF(2^16) back to GF(2^8). To this end, we explicitly construct a field isomorphism between GF(2^8) and the subfield of GF(2^16) that it is isomorphic to. Unfortunately, this is done by a direct element-by-element mapping: Every element of GF(2^8) is converted to its GF(2^16) form and then this mapping is inverted. I wish that I knew an algebraic trick to do this, but I can't come up with anything.

> 

X:= <seq(x^k, k= 0..d-1)>:
Y:= Expand~(<seq(newx^k, k= 0..d-1)>) mod p:
> 

assign(((C-> Iso(Expand(C.Y) mod p) = C.X)@`<,>`)~(combinat:-permute([seq(k$d, k= 0..p-1)], d))):  

Apply the isomorphism to R:

> 

R:= Iso~(R);

R := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 1, (4, 1) = 1, (4, 2) = x, (4, 3) = 1, (4, 4) = x^2})

Once again, verify that R^4 is the starting matrix, this time in its original form.

> 

Expand~(R^ROOT - M) mod p;

Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0})

> 

 



Download Matrix_powers_finite_field.mw

 

new, unrelated Question...

Carl Love 28100
April 19 2015
0 0

@matmxhu 

This is a new question completely unrelated to the theme of this thread: limiting the time of processes! Please post as a new Question, a new thread!

It is surprising...

Carl Love 28100
April 19 2015
0 0

@tomleslie 

Yes, I find it very surprising that frem has no symbolic capability and that it cannot be evaluated by evalhf. My solution seems, so far, to correct all problems.

You wrote:

[A]ll "solutions" I have seen are based on the premise: please don't check the arguments to frem() until floats are actually presented.

My (second) solution does not rely on that.

[I]f frem(x,y) when presented with one or two symbolic arguments. simply returned frem(x,y), then the issue would be avoided.

It is trivial to write an frem variant that does that, but it wouldn't solve all problems. For example, you wouldn't be able to symbolically integrate or differentiate it. My solution relies on the fact that floor is fully symbolic.

Can integrate backwards...

Carl Love 28100
April 19 2015
0 0

@patient For what it's worth (which I don't know at this point), the FindSingularity procedure that I gave will work for backwards integration, say using 0 as the initial point and integrating towards -1. No modification is needed---simply call FindSingularity(sol, -1).

A better frem...

Carl Love 28100
April 19 2015
0 0

@Preben Alsholm 

Try a time comparison of your procedural dsolve with my symbolic version of frem:

Frem:= (x,y)-> x - y*floor(x/y+1/2);

Ironically, Frem can be evaluated by evalhf, but frem can't be. For that reason, I think that it'll be faster.

Numeric integration digits...

Carl Love 28100
April 19 2015
0 0

@Axel Vogt 

What you called "acer's trick" was my trick. It seems that the setting of Digits determines whether declaring a method is needed. At Digits = 10, no method is needed. My guess is that the native Maple methods can't guarantee the accuracy at the higher setting.

Regarding the symbolic integration, see my new Answer about a symbolic equivalent to frem. It integrates easily.

So what?...

Carl Love 28100
April 19 2015
0 0

@Markiyan Hirnyk What does it matter that neither you nor I understand the first part of the error message? All that matters is that the program stops because the time limit has expired, that the error is trapable, and that the user can continue from there.

The problem with that...

Carl Love 28100
April 19 2015
0 0

@Markiyan Hirnyk The problem with using cpulimit is How do you continue the worksheet after the limit has expired? The timelimit command is more effective.

Please supply code...

Carl Love 28100
April 19 2015
0 0

Please supply the complete code that leads to the error message, preferably in plaintext form. If you have difficulty with the plaintext form, or the code is too lengthy, then post a worksheet using the green uparrow tool that is the last item on the second row of the toolbar in the MaplePrimes editor.

Please post worksheet...

Carl Love 28100
April 19 2015
0 0

Hello Raquel,

Please post a complete worksheet showing, at the least, the construction of Matrices K and M and the generation of the error message. To attach a worksheet to a post, use the green uparrow tool, which is the last item on the second row of the toolbar in the MaplePrimes editor.

What are Loading RealDomain; and solve; ? Are those commands that you gave? If so, they are meaningless. I think that you mean

with(RealDomain);
solve(LinearAlgebra:-Determinant(M), omega);

Personally, I would only use RealDomain after first trying to use solve without RealDomain.

Works for me...

Carl Love 28100
April 18 2015
0 4

It works for me. Perhaps your x has a value and you need to do a restart.

True...

Carl Love 28100
April 18 2015
0 0

@testht06 

Yes, that is true. But the computation R^4 = newM is only done to verify for the reader that R was computed correctly; it's not a necessary part of the computation of R. If you want, I could explicitly show that R^4 = M (rather than newM) by converting the entries of R^4 from GF(2^16) back to GF(2^8). (This may be computationally intensive---I'm not sure.) None of this would change the computation of R itself, which is (I believe) the main goal.

Different parameter levels...

Carl Love 28100
April 18 2015
0 0

@zmq 

The parameter of recurse is w, not v. If it was Departures(G,v), you'd be computing the departures of the root vertex repeatedly. Does that answer your question?

First 493 494 495 496 497 498 499 Last Page 495 of 709