@Declan 

f := (i,j) -> `if`(j=1, (i+4)^2, (i+4)^3):

@Declan 

The function (Maple procedure) f defines the (i,j) element of the matrix (i = row, j = column).

If j=1 then f(i,j) = i^2, otherwise (j=2 in our case), f(i,j) = i^3.

@Declan 

Why don't you use the clear and short preudocode algorithm at

https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

Note also that you should not use FactorInteger (see the pseudocode above).
(BTW, FactorInteger is in the MmaTranslator; otherwise ifactors does the job).

Also note that your fea(x, e, n)  exists in Maple as  x &^ e mod n.

@JessyOw 

Have you tried the provided hint? Then the answer should be clear: you can't.

Yes, each letter is a polygon.
Unfortunately, the code is not mine and AFAIK such vector fonts are copyright protected (I do not know why, but that's it), so I cannot post them because I do not want to have problems.

@John Fredsted 

Groebner[Basis]( [x - a*b*y,y - a*b*x, a^2-1,b^2-1],plex(y,x,a,b)):
select(has,%,[x,y]);

@John Fredsted 

For polynomial systems you could use Groebner (which is what eliminate does, but you have more options).

Groebner[Basis](eqs,plex(x,a,y));

@Axel Vogt 

I asked the question because I was hoping that Maple 2016 will compute the integral directly. But both answers confirmed that this is not the case. They also contain the easy solution (integration by parts, any order) common for Maple 2015 and 2016.

@Christian Wolinski 

OP seems to be a beginner. Don't you think that a more useful procedure would have been

Cubic:=proc(p,x,a,b)
  local s,f,horizontals;
  f:=unapply(p,x);
  s:=solve(diff(p,x),x);
  horizontals:=map(f,[s])[];
  plot([p, horizontals], x=a..b, color=[red,blue,blue]);
end;

Cubic(1/3 *x^ 3 - 6/10 *x^ 2 - 187/100* x + 1, x, -3, 3);

@I_Mariusz 

Only if generalized solutions are allowed, because the solution is not a classical one
(it is not differentiable for x=t).

@I_Mariusz 

What solution do you have in mind for:

pde1 := diff(u(t, x), t)+diff(u(t, x), x) = 0; bc1 := u(t, 0) = 0, u(0, x) = sin(x);

The series example series(harmonic(n-1), n) does not seem to be very inspired because

series(harmonic(n)-1/n, n);

already works.

The announced asymptotics for hypergeometric functions will be very useful I think.

@Kitonum 

When looking for N = N(eps) we are generaly non interested in finding the minimum (being impossible in most cases).

@Markiyan Hirnyk 

Unfortunately, the solution found by Maple is wrong.

(N < 0 for any epsilon>0)!

(I also started this way, but seeing that Maple has problems, I switched to asympt).

@Christian Wolinski 

You seem to be an adept of obfuscation. Why not like this (S=something)?

u:=-Omega*a*sqrt(2)*sqrt(-Omega^2*a^2-2*k*m+sqrt(Omega^2*a^2*(Omega^2*a^2+4*k*m)))/(-Omega^2*a^2+sqrt(Omega^2*a^2*(Omega^2*a^2+4*k*m))):

Q := -4*k*m+2*Omega^2*a^2*S+Omega^2*a^2*S^2:  
# S is the denominator*positive,  Q=0.

normal(simplify(u,[Q])) assuming a>0,k>0,Omega>0,m>0,S>0;
     -I

(you took S<-2 the negative root of Q just to be "more interesting").