It is obvious that b should not be allowed to be too large, otherwise the graph may become almost "dense" and will fit any data.

For example, DirectSearch ==>

DataFit(a*sin(b*x+c)+d, [b<=2,b>=0], Vector(X), Vector(Y), x, fitmethod = lms, strategy = globalsearch, method = quadratic);

gave

which seems resonable.

The arctan function seems to be not very "smart".

For example,

arctan( cos(t/2),sin(t/2)) assuming t>0,t<Pi/2;

simplifies to Pi/2 - t/2

but equivalent expressions

arctan( 1+cos(t),sin(t)) assuming t>0,t<Pi/2;

simplify(  arctan( 1+cos(2*t),sin(2*t))  ) assuming t>0,t<Pi/4;

do not.

In such circumstances I think that it would be better to use an inert variant; at least Maple will not try to simplify.

roundcoeffs:=proc(p,x)  local t;
 add(round~([coeffs(p,x,t)]) *~ [t])
end:

p:=3.5*x^2+3.2*x-6.5+88.3*x*y-y^3+a*y:
roundcoeffs(p,[x,y]);

Use

f := unapply(evalf(x^2/(sin(x)+x+1)), x);

restart;
f:=(x,y)->evalf(2*x*Int(sqrt(1+y^2*(t*x-1)^2/(1-(t*x-1)^2)), t = 0 .. 1)):
m:=19;n:=9;
G:=Matrix(m*n,3,datatype=float);
k:=0;
for i to m  do   for j to n do
k:=k+1;
G[k,1]:=2*i/(m+1); G[k,2]:=j/(n+1); G[k,3]:=f(2*i/(m+1),j/(n+1))
od od:

f_fit:=a1*x*sqrt(1+y^2*((a2*x-a3)^2/(1-(a4*x-1)^2)));
Statistics[NonlinearFit](f_fit, G, [x, y], initialvalues = [a1 = 1, a2 = 0, a3 = 0, a4=1/2],
output = [leastsquaresfunction, residuals]);
f1:=unapply(%[1],x,y);

printf("x    y    f(x,y)    f_fit(x,y)  \n");    
for x_ from 0.1 to 1.9 by 0.1 do
for y_ from 0.1 to 0.9 by 0.1 do
printf("%2.1f  %2.1f  %6f  %6f\n",x_,y_,f(x_,y_), f1(x_,y_)      );
od; od;

You may replace arcsin by wrongarcsin where:

wrongarcsin := x -> piecewise(x<0, arcsin(x)+2*Pi, arcsin(x));

You may simply "cut" the polynomial, or use:

hor:= (f,x,i) -> numapprox[hornerform](select(u -> degree(u,x)<=i, f) ,x);

f:=add( (10*i+1)*x^i, i=0..6);
hor(f,x,3);

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