@_Maxim_ 

It works, but by chance, for Digits=30 the coefficient of the 1/s term is 0.0  instead of  x.y*10^(-...).

evalf[10]( evalf[40](series(f(1.*s), s = 0, 1)) );
        (1.*10^(-40)*I)/s+O(s)

@Ian Jones 

You must have typos; it is not validated numerically.

@Carl Love 

I mean this:

Download digits.mw

@Carl Love 

But in "most" situations evalf[d1](e)  is enough because Digits is automatically increased by evalf if necessary.

1. The third argument of Cubes must be an m-dimensional list. So, use:
Cubes(3, 3, [a, b, c], h);

2. For an enumeration of N^m  see ?vectoint, ?inttovec

@_Maxim_ 

Nice analysis!

Have you used Maple for the convolution, or you did it by hand?

Cubes:=proc(m::posint,n::posint,a::list,h:=1) # for the cube [a, a+~h] in R^m
local k,u,T,C;
T:=combinat:-cartprod([[seq(k/n*h, k=0..n-1)] $ m ]):
C:=('T[nextvalue]()')$(n^m):
seq( [a+u, a+u+~h/n], u=[C])
end:

Cubes(2,2,[a,b],h);

This is probably because you like the symmetry.
But  [x=3,y=2], [x=2, y=10] etc, are also correct.

@Christian Wolinski 

@arjangash If you really have no problem about optimization, it should be obvious.
Here is a simpler example:

restart;
Z:=(a,b,g)->a+b+g:
constr:=(a,b,g) -> (a>=0,b>=0, a^2+b^2-2*a+1<=g):
Zmax:= g -> Optimization:-Maximize(Z(a,b,g), [constr(a,b,g)])[1]:
plot(Zmax, 0.1 .. 2.2);

@acer 

You are right, I simply forgot about having reported it, sorry! Should I delete the question?

In my opinion the problem is not to find workarounds (which are not complicated) but to understand where exactly is Maple's fault.

Maple can be used to compute the limit of the sequence

n -> sqrt(n) * (sin@@n)(1);

Can MMA compute it?

@Markiyan Hirnyk 

It would have been nice if MMA had given the simplified result sqrt(2)/4.
Probably it simply uses the fact that the sequence exp(I*n) is dense in the unit circle.

Maybe some day we will get something like:
f := sin(n)/(3 + cos(n)):
limsup[discrete](f,n)
     sqrt(2)/4;

For the moment:
maximize(f, n=0..2*Pi);
    
 

Just a short remark.
If the sequence is given by a procedure f, then its limit is L iff
limit( f(floor(x)), x=infinity) = L.

This way, the assumption is not necessary.

Unfortunately this does not work because floor has a very poor integration in Maple [I don't know why].
A simple (and embarrassing) example is:

simplify( floor(x) ) assuming x>3,x<4;
    floor(x)

@Markiyan Hirnyk 

It is OK for your surface defined by x^2+y^2+z^2=1, x^2-z^2<=1
but not for Mariusz' which was the boundary of the body  x^2+y^2+z^2<=1, x^2-z^2<=1.
In the previous comment both situations were answered.
But indeed in Mathematica code I see an "<"; probably this means the exclusion of this part of the boundary [I do not know the MMA syntax for this].