@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].

@Mariusz Iwaniuk

Download surface-int.mw

@digerdiga 

For polynomials you should use:

restart;
q:=(x-1)*(x-2)*(x-3):
sol:=c -> fsolve(q+c,x,complex):
r:=c -> max(Re~(select(t->abs(Im(t))<1e-10, [sol(c)]))):
plot(r, -8..8, discont, color=green);

Or, better:

restart;
q:=(x-1)*(x-2)*(x-3):
sol:=c -> fsolve(q+c,x):
r:=c -> max(sol(c)):
plot(r, -8..8, discont, color=green);

@Kitonum 

I am sorry that you, the proposer of the test, have trouble in getting the answer.
It seems that SolveTools:-SemiAlgebraic has a bug in the 32 bit version [or, it could be a memory management problem].
On a 32 bit computer I also get the same error for:

SolveTools:-SemiAlgebraic(
[(x-4)^2+y^2 <= 25, 9 <= x^2+(y-3)^2, 16 <= (x+w)^2+y^2],
[x, y], parameters = [w]);

which is called by procedure.

I hope that this bug will be addressed soon; then MultiIntPoly will work for the 32 bit version too.

@Avel Vogt:  are you on 32 bit too?

@max125 

plots:-animate(plot, [(x-1)*(x-2)*(x-3) + c, x=0.5 .. 3.5], c=0..1);

or

Explore(plot((x-1)*(x-2)*(x-3) + c, x=0.5 .. 3.5, view=-1..3), c=0. .. 1., animate);
(you may remove the animate option and use the slider instead).

@rlopez 

It is useful to note that if S₀ is a C^oo selection in an interval [a,b] then there is a continuous selection S in R which extends S₀; but a C^1 selection S extending S₀ may not exist. This is of course valid for eigenvalues too.

@digerdiga 

1. If we are interested only in real roots, the discontinuity cannot be avoided, as the example shows.

2. Yes, h(y)>0.384... is enough. The exact value is 2*sqrt(3)/9  but this is not important.

3. When all the complex roots are considered, the continuity can be obtained.
For example,
r1(c) = min({|z| : p(z) = 0},  r2(c) = max({|z| : p(z) = 0}
are continuous.

restart;
q:=(x-1)*(x-2)*(x-3):
sol:=c -> solve(q+c,x,explicit):
r1:=unapply( min(abs~([sol(c)])),c):
r2:=unapply( max(abs~([sol(c)])),c):
r:=c -> max(select(t->Im(t)=0, [sol(c)])):
plot(r, -8..8, discont, color=green);
plot([r1,r2], -8..8, color=[red,blue]);

@Rouben Rostamian  

I just want to mention that the piecewise method also handles non-convex domains (polygonal or not):

poly:= y<=x and y>=-x and (y>=4*x-3  or y<= -4*x+3) and x>=0 and x<=1;
F:=piecewise(poly, x^2+y^2, undefined);
G:=piecewise(poly, 3*sin(x)-y^2, undefined):
fieldplot([F,G], x =0 .. 1, y = -1 .. 1, arrows = SLIM, color = x) ;

@Carl Love 
Except that mine obtains the monomials of total degree exactly d.
To get "up to d" as OP wants it is even simpler:

mondeg:=proc(X::list(name),d)
expand((add(X)+1)^d);[op(%)]/~coeffs(%)
end:

Or, better (as you did):

mondegx:=proc(X::list(name),d)
local t;
coeffs(expand((add(X)+1)^d,X),X,t);[t]
end:

@John May 

Here is a simple example (no parameters):

SolveTools:-SemiAlgebraic([x^2+y^2<z, z^2=2, z>0], [x,y,z]);

It is correct, but I was expecting:

Note that solve (instead of SemiAlgebraic) returns the same RootOf result (because SemiAlgebraic is called) !

@John May 

Thank you.

This form of the RootOf appeared in the output of the command SolveTools:-SemiAlgebraic.
But actually it could be useful to be able to select symbolically the i-th real root.

BTW, the option real for fsolve seems to be also not documented (not needed anyway, real being the default).

It would be very useful to have a procedure to compute the area of the region given by inequalities (implicitely) using Green's formula. But I think that this would be very difficult.

It is possible to use the solve method for

rel:={(x-4)^2+y^2<=25, x^2+(y-3)^2>=9, (x+sqrt(7))^2+y^2>=16}:

and obtain

-(3/4)*55^(1/2)+(3/4)*13^(1/2)*7^(1/2)+(25/2)*arcsin(4/5)+3*7^(1/2)+(25/2)*arcsin(82/125-(9/250)*13^(1/2)*7^(1/2))+(9/2)*arcsin(6/25+(3/50)*13^(1/2)*7^(1/2))+8*arcsin((23/128)*7^(1/2)+(9/128)*55^(1/2))+(9/2)*arcsin(-(3/32)*7^(1/2)+(3/32)*55^(1/2))+8*arcsin((1/4)*7^(1/2))+12

I am going to post a method in the Application Center.

@Mariusz Iwaniuk 

For your R the result given by Area is correct because R is empty.

But in Mathematica it is +3/4 instead of -3/4. In this case Maple's solve needs a very LONG time and I had not the patience to wait for a result.

@Markiyan Hirnyk 

OK, so you are sure that I have just modified Maxim's code.
Anyway, until now you have two pieces of bad code.

@Markiyan Hirnyk 
Probably you prefer:

Area1:=proc(rel::set(relation))  # Maxim's
local x,y, xy:=indets(rel,name),cad,f,g, xmax, xmin, ymax, ymin;
if nops(xy)<>2 then error "Two vars only!" else x,y:=xy[] fi;
cad := solve(rel, [x, y]):
f := (a,b) -> if a = y then ymax = b elif b = y then ymin = a elif a = x then xmax = b else xmin = a end if:
g := s -> if subs(s, {xmax, xmin, ymax, ymin})::freeof({xmax, xmin, ymax, ymin}) then
    int(subs(s, ymax-ymin), x = subs(s, xmin) .. subs(s, xmax)) else 0 end if:
(simplify@add@curry(map, g))(applyrule(
  ['`<`(a::anything, b::anything) = f(a, b)', '`<=`(a::anything, b::anything) = f(a, b)'],
  cad)):
end:

BTW, I really like it!