Neither the command-line interface (CLI) nor the print command is needed to realize this bug. Apparently, the bug is in the kernel. Any theories about how this arises? It seems that there's two distinct copies of t23.

@Muhammad Usman 

(vv: Thanks for explaining the raising to the 1/p power.)

Muhammad:

Please try this corrected translation of your algorithm. Note that it is divided into two procedures.

`mod/ReduceExponents`:= (f::polynom, p::posint)->
     #This procedure should only be used when all exponents of f are multiples of p.
     evalindets(f, `^`, t-> op(1,t)^iquo(op(2,t), p))
:

`mod/SFF`:= proc(f::polynom, p::posint)
local R:= 1, x:= indets(f, name)[], i, c, w, g, y, z;
     if x=NULL then return f mod p end if; #f is degree 0.
     g:= diff(f,x);
     if g=0 then return thisproc(ReduceExponents(f) mod p, p)^p end if;
     c:= Gcd(f,g) mod p;
     w:= Quo(f,c,x) mod p;
     for i while w <> 1 do
          y:= Gcd(w,c) mod p;
          z:= Quo(w,y,x) mod p;
          R:= R*z^i;
          w:= y;
          c:= Quo(c,y,x) mod p
     end do;
     R*thisproc(ReduceExponents(c) mod p, p)^p #returns R if c=1
end proc:

f:= x^11 + 2*x^9 + 2*x^8 + x^6 + x^5 + 2*x^3 + 2*x^2 + 1:
SFF(f) mod 3;

@Muhammad Usman 

Ah, yes, I didn't include the f^(1/p) and c^(1/p) steps. I don't know how to handle that. How can this be done without fully factoring the polynomial? I see that in the case you mention, you can just replace x^3 with x, but what about the general case?

@Just_Some_Guy 

In your worksheet, you have a vector of independent-variable values, X, which doesn't appear in your original Question. Also, what you called E---the independent variable in the original---has become dependent variable f in the above. Which is correct?

@Maple4evergr8 Before upgrading, you may want to confirm that that bug has indeed been fixed. I'm not sure. I'm sure that someone else reading this who has Maple 2015 can answer that question.

@Bendesarts You need to experiment with different initial conditions and different t-ranges until you find what you want. In the command below, the intial conditions are [1,1] and the t-range is 0..2.

DEtools[DEplot](EqSys, [x,z](t), t= 0..2, [[x,z](0)=~[1,1]]);

You need initial conditions: x(0) = ? y(0) = ?

@Preben Alsholm Thanks for the explanation. Although the work is done "by hand," the triangle explanation is sufficiently detailed to get a vote up from me. I know that Maple completely refuses to reverse the order of definite double integrals. I wonder if Mathematica or any other CAS has any capability at it for regions more complicated than rectangles.

@Preben Alsholm I have no doubt that your Answer is correct---the results can be verified with intsolve---but I wonder What's the rule of calculus that lets you change the limits of integration from x= 0..y to x= t..y? 

A variable p is used in the algorithm, in exponents only. I don't see where its value is set. So, what is p?

@J F Ogilvie Your reputation is next to your name on every post, next to your name at the top of the page (if you're logged in), and next to your name in many other places. Yours is currently 63.

@Preben Alsholm Thanks. Noted and corrected.

@Axel Vogt 

The links to a Mathematica marketing page and to a blog posted by Acer in a Reply to the Answer below suggest that Mathematica has implemented a general algorithm for highly oscillatory integrals rather than an ad hoc method to address the contest problem. I'm not a Mathematica user, so I can't test this.

@michjoh 

My primary purpose in the Answer above was to explain how to do a double break because sometimes it can't reasonably be avoided, and it also seemed to be the way that your thinking was headed. However, your code above uses an explicit loop to search through a table. This is not an efficient way to search. The time complexity is O(M*N^2) or possibly even O(M*N^2*ln(N)) (with the extra ln(N) because of the possibility of table indices not being stored linearly---I'm not sure how they're stored). Tom Leslie came up with a better way to do it: He uses an array for storage and the built-in member for the searches. Search by member still uses time proportional to the size of the object being searched, but is significantly faster than a for loop because it's built-in. I think that I can reduce the search time to logarithmic by using a table to mark the visited spots (kind of like marking your trail when you're lost in the woods so that you know when you've crossed your own trail). I'll probably have a chance to code it later today. It'll have time complexity O(M*N*ln(N)).

@michjoh I'm sorry that I didn't see this "Follow-up question" at the time that you posted it. On a quick reading, I think that this question has subsequently been answered in your subsequent thread "Break to form new lines". So, before I delve into this, please let me know if an answer is needed here.