Maple can handle some answers with countable, discrete infinities of solutions. It can also handle very simple cases of finite numbers of solutions which are real ranges. As far as I know, there is no code to deal with solutions which consist of a countable union of ranges. In this particular case, something could probably be done and return something like Range(Open(2*_Z1*Pi),Open((2*_Z1+1)*Pi)). But only extremely simple cases could be done.

Whenever it matters what the name of a variable is in a symbolic computation which should be name-independent, then it's a bug. For summation, i is often a problem. For integration, it is X and _X (even though the latter is well-documented as being reserved for Maple's use).

Whenever it matters what the name of a variable is in a symbolic computation which should be name-independent, then it's a bug. For summation, i is often a problem. For integration, it is X and _X (even though the latter is well-documented as being reserved for Maple's use).

degree with one argument should be deprecated - it was known to be a really bad idea 15 years ago, but no one ever had the guts to do something about it (it might break some people's code which accidentally worked). In any case, if you fix that, then the resulting code works: foldl((y, v)-> map((x)-> `if`(degree(x,indets(lhs(v))) = degree(lhs(v),indets(lhs(v))), subs(lhs(v) = rhs(v), x), x), y), Poly, a[1] = a[1]^2, a[2]^2 = a[2], a[4]^3 = a[4]^2); gives what you want. (Sorry for the long delay in answering, I was away at a conference).

degree with one argument should be deprecated - it was known to be a really bad idea 15 years ago, but no one ever had the guts to do something about it (it might break some people's code which accidentally worked). In any case, if you fix that, then the resulting code works: foldl((y, v)-> map((x)-> `if`(degree(x,indets(lhs(v))) = degree(lhs(v),indets(lhs(v))), subs(lhs(v) = rhs(v), x), x), y), Poly, a[1] = a[1]^2, a[2]^2 = a[2], a[4]^3 = a[4]^2); gives what you want. (Sorry for the long delay in answering, I was away at a conference).

Safer yet would be f:=proc(n::name, x,y); point(n, x,y); end proc; otherwise you'll get rather strange error messages if someone forgets to quote their input.

Granted - that was just something I whipped up. My version does not require the matrices to be square by the way (or did I goof?). And I am sure I could make it work for all 2-dimensional rtables without much additional effort - but the code would be less readable. The code I posted was meant to be somewhat pedagogical.

Granted - that was just something I whipped up. My version does not require the matrices to be square by the way (or did I goof?). And I am sure I could make it work for all 2-dimensional rtables without much additional effort - but the code would be less readable. The code I posted was meant to be somewhat pedagogical.

This code which I believe is more idiomatic is asymptotically faster: Kron := proc(A::Matrix, B::Matrix) local rA, cA, rB, cB; (rA, cA) := LinearAlgebra:-Dimensions(A); (rB, cB) := LinearAlgebra:-Dimensions(B); Matrix(rA*rB, cA*cB, proc(i,j) local iA,jA,iB,jB; iB := irem(i-1,rA,'iA'); jB := irem(j-1,cA,'jA'); A[iA+1,jA+1]*B[iB+1,jB+1] end); end; It is too bad that one has to play indexing games (irem is 0-based, Matrix is 1-based) as it obscures the result. Also I really wish there was a side-effect free version of commands like irem!

This code which I believe is more idiomatic is asymptotically faster: Kron := proc(A::Matrix, B::Matrix) local rA, cA, rB, cB; (rA, cA) := LinearAlgebra:-Dimensions(A); (rB, cB) := LinearAlgebra:-Dimensions(B); Matrix(rA*rB, cA*cB, proc(i,j) local iA,jA,iB,jB; iB := irem(i-1,rA,'iA'); jB := irem(j-1,cA,'jA'); A[iA+1,jA+1]*B[iB+1,jB+1] end); end; It is too bad that one has to play indexing games (irem is 0-based, Matrix is 1-based) as it obscures the result. Also I really wish there was a side-effect free version of commands like irem!

You have been given a lot of help already and some useful hints. You'll probably get a lot more help here if you do parts of the work yourself, and then come back when you get stuck. If you explained a bit more why you want this (it's not part of an assignment, is it?), you may be able to get more people to help you. Remember, such problems are quite difficult, so they take time. If I am intrigued by a problem, I'll spend time on it (as will many others here).

You have been given a lot of help already and some useful hints. You'll probably get a lot more help here if you do parts of the work yourself, and then come back when you get stuck. If you explained a bit more why you want this (it's not part of an assignment, is it?), you may be able to get more people to help you. Remember, such problems are quite difficult, so they take time. If I am intrigued by a problem, I'll spend time on it (as will many others here).

It is currently on the front page of YouTube [but that is unlikely to remain so for much longer].

I don't have a full proof, but it is easy to see that the convergence will be slow by plotting (sin(x)/x)^n for various n for x from 0 to 1. While the area under the curve does go down, it does so very, very slowly (while the rest of the curve goes to zero amazingly fast). My proof sketch revolves around showing that there is a breakpoint function alpha(n)>0 such that the integral 0..alpha(n) and alpha(n)..infinity both go to 0.

I don't have a full proof, but it is easy to see that the convergence will be slow by plotting (sin(x)/x)^n for various n for x from 0 to 1. While the area under the curve does go down, it does so very, very slowly (while the rest of the curve goes to zero amazingly fast). My proof sketch revolves around showing that there is a breakpoint function alpha(n)>0 such that the integral 0..alpha(n) and alpha(n)..infinity both go to 0.