@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)));

But as I can see, it is not possible to obtain the "singular" value c=3  just using infolevel.
You will have to use the suggested method in the answer.

@codell 

The variables having xk=true or xk=false  are fixed. They were determined in the first step because BooleanSimplify cannot simplify the whole expression (being too large).
The rest of the variables appear in `simp`; this is the simplified version of your expression (after fixing the previous variables). From here you can imediately see for what values of the remaining variables the expression (ee) is true.

@roubeur 

@Preben Alsholm 

What I meant is that it seems that Basis does not call primpart any more, so the infolevel method does not work.