@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!

@tomleslie 

The intention was not to write efficiently, but to exhibit a bug.

@williamdiao 

There are only 2 possibilities:
1. You did something wrong.
2. It is a bug.

This is all that can be said without seeing the worksheet.

@Carl Love 

Nice solution, vote up!

I think that the module should also contain a function to remove the "noise".
Compare:

P:=plot(sin(x)*sin(40*x),x=0..5*Pi);
and
P:=plot(sin(x)*sin(40*x)+sin(1000*x)/100,x=0..5*Pi);

1. See  https://en.wikipedia.org/wiki/Extreme_point

2. By the definition of the extreme points

@Carl Love 

Yes, I know this. What I meant is that the result would be HUGE or impossible to compute (think at the degree of the denominator).

@Carl Love 

I don't think that the extension is viable. Try:

A:= LinearAlgebra:-RandomMatrix((n$2), generator= (()-> randpoly(x, degree= 2)/randpoly(x, degree= 2)));