@Markiyan Hirnyk 

I will have a try again this evening.

I am not aware of that parts of Geometry/Topology. But if it is a priori known to
be an integer then you can do it numerically - as soon as one has some sound
error estimates for the numerical result. I think Maple can do that.

Do you have the non-parametrized version of the integral as well?

Hm - then why you do not use the geometric reasoning?

Currently I only see a symmetry in t = Pi (allowing to reduce it to 2*...) and another
in s = Pi/2 (allowing a split in s=Pi, where the first up to Pi/2 is the same as 2 times
up to Pi and the second split is s = Pi .. 2*Pi - having another symmetry in Pi/2 + Pi).

May be that helps to pass from trig to algebraic?

May be the old Maple 9 would do that as default (I do not have it installed) as the comment 'regression' by acer may indicate.

Here is a way to do it using results due to Koepf in the book "Hypergeometric Summation", 1998.

It handles the case n=0 and 1 <= n correctly.

Thus a correct and short answer to the inital task is 1/3 of the term
sinh(1/2)*cos(sqrt(3)/2)+sqrt(3)*cosh(1/2)*sin(sqrt(3)/2)+sinh(1);

factorial_sum.mws
hsum13.mpl

factorial_sum.pdf

Edited: newer code version at http://www.mathematik.uni-kassel.de/~koepf/hsum15.mpl

May be the old Maple 9 would do that as default (I do not have it installed) as the comment 'regression' by acer may indicate.

Here is a way to do it using results due to Koepf in the book "Hypergeometric Summation", 1998.

It handles the case n=0 and 1 <= n correctly.

Thus a correct and short answer to the inital task is 1/3 of the term
sinh(1/2)*cos(sqrt(3)/2)+sqrt(3)*cosh(1/2)*sin(sqrt(3)/2)+sinh(1);

factorial_sum.mws
hsum13.mpl

factorial_sum.pdf

Edited: newer code version at http://www.mathematik.uni-kassel.de/~koepf/hsum15.mpl

Using the given answer of Maple 9 for n=0 there seems to be a bug
in converting to StandardFunctions:

hypergeom([], [1/3, 1/2, 2/3, 5/6, 7/6], 1/46656):
% = convert(%, StandardFunctions);

  hypergeom([],[1/3, 1/2, 2/3, 5/6, 7/6],1/46656) =
    1/6*3^(1/2)*exp(-1/2)*sin(1/2*3^(1/2))+
    1/6*3^(1/2)*exp(1/2)*sin(1/2*3^(1/2))+
    1/6*exp(1/2)*cos(1/2*3^(1/2))

evalf(%);
                 1.00019841285900 = 0.673955774283463

Using the given answer of Maple 9 for n=0 there seems to be a bug
in converting to StandardFunctions:

hypergeom([], [1/3, 1/2, 2/3, 5/6, 7/6], 1/46656):
% = convert(%, StandardFunctions);

  hypergeom([],[1/3, 1/2, 2/3, 5/6, 7/6],1/46656) =
    1/6*3^(1/2)*exp(-1/2)*sin(1/2*3^(1/2))+
    1/6*3^(1/2)*exp(1/2)*sin(1/2*3^(1/2))+
    1/6*exp(1/2)*cos(1/2*3^(1/2))

evalf(%);
                 1.00019841285900 = 0.673955774283463

For my answer I would say: a kind of "at least not a full homework solution" according to
your initial remark ( ~ 'z < 0 ? and then ?'). Yes, I used 0 < z.

Stokes will be difficult to be handled beyond curves as boundraries, I guess.

Admitted, that this task (and the other in dim=3) is about change of coordinates. But:

S:=solve({x^2+y^2 <= 2*y,  y<=x, 0<x, 0<y}, {x,y}, parametric) returns 3 solutions, where
two of them are curves and can be ignored for integration, being meager.

  S[2];
                                           2       1/2
              [0 < y, y < 1, y < x, x < (-y  + 2 y)   ]

For better understanding plot it (x becomes the ordinata)

  plot([sqrt(2*y - y^2),y], y=0..1, color=[red,blue]);

and we are seeking within the domain between the red (=circle) and blue (=line) bounds.

And that is what Preben did: one integrates between y = 0 .. 1 and depending on that
for such y one needs those x which satisfy y < x < sqrt(2*y - y^2). Which results in
the double integral (Cartesian coordinates) he wrote down and then solved it.

Stokes will be difficult to be handled beyond curves as boundraries, I guess.

Admitted, that this task (and the other in dim=3) is about change of coordinates. But:

S:=solve({x^2+y^2 <= 2*y,  y<=x, 0<x, 0<y}, {x,y}, parametric) returns 3 solutions, where
two of them are curves and can be ignored for integration, being meager.

  S[2];
                                           2       1/2
              [0 < y, y < 1, y < x, x < (-y  + 2 y)   ]

For better understanding plot it (x becomes the ordinata)

  plot([sqrt(2*y - y^2),y], y=0..1, color=[red,blue]);

and we are seeking within the domain between the red (=circle) and blue (=line) bounds.

And that is what Preben did: one integrates between y = 0 .. 1 and depending on that
for such y one needs those x which satisfy y < x < sqrt(2*y - y^2). Which results in
the double integral (Cartesian coordinates) he wrote down and then solved it.

I think it silently uses the integration bounds, " ... assuming 0 < y, y < 1".

One can do it in 2 steps: the inner integral, then simplify assuming about y.
After that it comes out as well

I think it silently uses the integration bounds, " ... assuming 0 < y, y < 1".

One can do it in 2 steps: the inner integral, then simplify assuming about y.
After that it comes out as well

Yes, 'solve' for the condition gives 2 curves - which are meager, the 3rd solution is a circle segment

Yes, 'solve' for the condition gives 2 curves - which are meager, the 3rd solution is a circle segment

I give up, sorry.

I give up, sorry.