Your notations are not clear (at least for me).

Given the function  g = f(r,theta)*r
you can integrate it in the domain

A = { (r,theta) :  a<= r <= b,  c(r) <= theta <= d(r)}

or

 B = { (r,theta) :  a<= theta <= b,  c(theta) <= r <= d(theta)}.

In the A case, the integral is
int(int(g, theta=c(r)..d(r)), r = a..b);

and similarly for B.

Please rephrase your question in this context.

numtheory[cfrac](3/11+3*sqrt(3),'periodic', 'quotients');

The only possibility I see is to visualize.

E.g. the z-sections using

ee:=subs( [a[1]=x,a[2]=y,a[3]=z], %);
f:=unapply(rhs(ee)-lhs(ee),x,y,z);
Explore(  plots[implicitplot](f(x,y,z)>=0, x=-20..20,y=-20..20,filledregions,numpoints=15000),
parameters = [z=-20..20] );

Cf:=proc(x,n:=floor(Digits/log10(2)+2))
if not(type(x,realcons)) then
  if length(1-x)< length(x) then return 1-'Cf'(1-x) else return 'Cf'(x) fi;
fi;
if x=0 or x=1 then return x fi; if n<0 then return 0.0 fi;
piecewise(x<1/3, Cf(3*x,n-1)/2, x<=2/3, 1/2, 1/2+Cf(3*x-2,n-1)/2)
#  piecewise(x>1/2, 1-Cf(1-x,n-1), x<1/3, Cf(3*x,n-1)/2, 1/2)
end;

Cf(1-3*x)+Cf(3*x);
    1

Cf(2015/2222);
    13/16

plot(Cf,0..1);

The condition number has a standard definition, and cannot be changed. Probably there is a bug also in LU or Cholesky factorization which should be investigated.

Norm(B,infinity),Norm(B^(-1),infinity); ==> 3,2  hence OK.

Edit: for other matrices such as A:=HilbertMatrix(3,3) the results are correct.

asympt((x^2+3)/(x+2),x,8);

(there are mistakes in your computations)

convert(sin(x),FormalPowerSeries);

For two sets only, it's easier to use seq.

A:={a,b,c}; B:={x,y}:

{ seq(seq( u*v, u=A ), v=B) };

Proj:=proc(A::Matrix,b::Vector[column],v::Vector[column])
uses LinearAlgebra;
local u0,P,_xxx;
u0:=LinearSolve(A,b,free=_xxx);  #error if {x:Ax=b} is empty
u0:=eval(u0,_xxx=0);
P:=ProjectionMatrix(NullSpace(A));
#   u0+P.(v-u0)
u0+P.(Vector(v)-u0) # To work with Gradient without ConvertVector(...,free)
end;

Example

Proj(<<1|2|3>,<2|4|6>>, <1,2>, <0,0,0>);

G:=proc(k::posint, n::posint)
local f,g,x,y,u,v;
f:=(1/2)*(1-cos(k*Pi*x)*cos(n*Pi*y))/(sin((1/2)*Pi*x)^2+sin((1/2)*Pi*y)^2);
g:=simplify(subs([x=2/Pi*arctan(u),y=2/Pi*arctan(v)],f)*4/(Pi^2*(u^2+1)*(v^2+1)));
expand(g);
int(int(%,u=0..infinity),v=0..infinity)
end:

t:=time():
G(9,9);
time()-t;

   1.279

After some "reverse engineering", it is a group but for mod 70.

Simply construct the Cayley table:

G:=[5, 15, 25, 45, 55, 65];
Matrix(6, (i,j) -> G[i]*G[j] mod 70);

Alternatively, you may use Magma (see ?Magma) but the problem is too simple to need it..

phasor:=(r,a)->r*exp(I*a*Pi/180):   #simpler version
ec:=phasor(U__C, -90) = phasor(400, X)-phasor(150, 0)+phasor(1744/10, 90):
s:=solve(evalc([Re(ec),Im(ec)]),allsolutions);


uc:=eval(U__C,s):
[eval(uc,[_Z1=0,_B1=0]),eval(uc,[_Z1=0,_B1=1])];


evalf(%);


So, you will take the second (positive) solution.

Note that
fsolve(evalc([Re(ec),Im(ec)]), {U__C=0..1000, X=0..360});
does not provide a solution (it should!), but
fsolve(evalc([Re(ec),Im(ec)]), {U__C=200, X=270});

Here pi is actually numtheory[pi]; the package was loaded previously.

The function

y:=numtheory[pi](floor(x))/(x/ln(x))-1;

has indeed the graph:

Note that this function has an infinity of discontinuities and

limit(y, x=infinity) = 0.

(this is the celebrated  Prime number theorem; unfortunately Maple cannot compute the limit).

restart;
with(plottools): with(plots):
dis3d:=proc(p1,p2)
  local A,B;
  A:=getdata(p1)[3];
  B:=getdata(p2)[3];
  min( seq(seq(seq(seq(add( (A[i,j]-B[ii,jj])^2 ), i=op(2,A)[1]),j=op(2,A)[2]),ii=op(2,B)[1]),jj=op(2,B)[2]));
  sqrt(%)
end:
num:=50: np:=9:tp:=9: eps:=3:
#  increasing  np,tp  too much will slow down the program

cen:=rand(0.0..100.0):
dir:=2*rand(0..1)-1:
h:=rand(1.0..60.0):
rr:=rand(0.25..5.0):
t[0]:= tubeplot([dir()*h()*t, dir()*h()*t, dir()*h()*t], t = -1 .. 1,
        radius = rr(),style=surface,numpoints=np,tubepoints=tp):
for n to num do
  print("n=",n);
  ok:=false;
  while not(ok) do
     t[n]:= tubeplot([cen()+dir()*h()*t, cen()+dir()*h()*t, cen()+dir()*h()*t],
            t = -1 .. 1, radius = rr(), style=surface,numpoints=np,tubepoints=tp);
     for k from 0 to n-1 do
       ok := dis3d(t[k],t[n]) >eps;
       if not(ok) then break fi
     od;
  od;
od:
display( seq( t[n],n=0..num), axes=none,scaling=constrained);

 Edit. dir() was not initially used; it is intended for "more random" directions of the cylinders.
dis3d(p1,p2)  computes the distance between two 3d plots.
A new cylinder is added only if the distance to the previous ones is > eps.

Use a procedure.

f:=unapply(Determinant(M2),omega):
fsolve(f, complex);

It seems that there is no solution.

N:=Matrix( 400,400, (i,j)-> abs(f(-10+i/20 (-10+j/20)*I)) );
min(N),max(N);

So, your function is almost constant and very close to 0. You can't hope to solve it numerically; even if you find a program providing a solution, I don't think you should trust it. You should re-think the problem i.e, try another approach to solve it if you are sure that a solution (of the problem, not of this equation) exists.

P.S. I don't understand the reason for evalf[4].