dxdy  :=   y=-sqrt(1-x^2) .. sqrt(1-x^2), x=-1..1 ;

int(sin(x^2) *cos(x^2) , dxdy) =

1/2*int( sin(x^2+y^2) + sin(x^2 - y^2), dxdy ) =

1/2*int( sin(x^2+y^2) , dxdy)  + 0 =

1/2* int( sin(r^2)*r, r=0..1, t=0..2*Pi ) =

1/2*(1 - cos(1))*Pi =

Pi * sin(1/2)^2

It's a very difficult task to compute u3 even at a single point.
For example, u3(0) can be simplified to the triple integral of
f := -(2.533317577*10^7*I)*z*exp((-.64-12.56637062*I)*x^2-(7.853982892*10^7*I)*z^2-(25.13274124*I)*y^2)*y*x*BesselJ(0, 25.13274123*y*x)*BesselJ(0, 62831.85308*z*y)

But f is very oscillatory and to compute its triple integral over [0,1]^3 is almost impossible even using MonteCarlo method.
Probably some special methods will be needed.
Just look at this plot:
plot( eval(Re(f),[x=1/2,y=1/2]),z=0..1 );

P.S. Was this integral invented, or it appears in a concrete problem?

eq:=proc(L::listlist,x::name,y::name,z::name)
  primpart(LinearAlgebra:-Determinant(<Matrix([[x,y,z],L[]])|<1,1,1,1>>))=0
end:

map(eq,L,x,y,z);

TriOnUnitSphere:=proc(A::Vector[column](3),B::Vector[column](3),C::Vector[column](3)) #A,B,C = directions
uses LinearAlgebra,plots,plottools;
local a:=Normalize(A,2),b:=Normalize(B,2),c:=Normalize(C,2),u,v;
display(sphere(),  plot3d(1.01*Normalize(c+v*(a+u*(b-a)-c),2), u=0..1,v=0..1,color=red,style=surface,_rest),
pointplot3d([1.01*a,1.01*b,1.01*c],symbolsize=20,color=blue,symbol=solidsphere)  )
end:

TriOnUnitSphere(<1,2,3>, <-1,1,-2>, <-1,-1,3>);

TriOnUnitSphere(<1,0,0>, <-10,1,0>, <0,1,5>, grid=[200,200]);

with(LinearAlgebra):with(plots):with(plottools):
A1:=<0,0,1>: A2:=<0,3/5,4/5>: A3:=<4/5,3/5,0>: A4:=<5/13,0,12/13>:
seg:= (a,b)->spacecurve(Normalize(a+t*(b-a),2),t=0..1, thickness=3, color=red):
display(sphere([0,0,0],1),seg(A1,A2),seg(A2,A3),seg(A3,A4),seg(A1,A4));

You can use e.g.

binaryvariables = indets(L, specindex(integer,x));

where L is a list containing the constraints (and the function).
Or, you may generate their sequence using seq.
 

f:=(x,y) ->exp(-x^2-y^2):
L:=4: n:=64:
r1:=rand(0..1.0): r:=rand(-1..1.0):
a:=Vector(n,L*r): b:=Vector(n,L*r):
c:=Vector(n,0.3*r):
da:=Vector(n, i->0.5*''r1()''):
P:=proc(aa)
  global n,a,b,c,da,h,L;
  a:=a+eval~(da);
  map[inplace](frem,a,L);
  h:=add(c[i]*f(x+a[i], y+b[i]),i=1..n):
  plot3d(h, x=-L..L,y=-L..L, color="Aquamarine", scaling=constrained, view=-2..2, axes=none, style=surface, lightmodel=light1);
end:
plots:-display(seq(P(aa),aa=1..50),insequence=true);
## or ...
Explore(P(aa), aa=0..1, animate,loop,size=[800,800]);

Download req-int.mw

restart;
f:=r^2 *cos(theta)+r*sin(theta):
r1:= 0.3+0.1*cos(theta):
r2:= 0.5+0.1*cos(theta):
plots:-densityplot(f, theta=0..2*Pi, r=r1 .. r2,  
coords=polar, colorstyle = HUE, style = patchnogrid);

If you need labels:

plots:-display(%, labels=["x","y"]);

c:=Matrix([    # result = 22.613, [4, 2, 5, 3, 6, 1, 4]
[0., 3.60555127546399, 2.00000000000000, 8.06225774829855, 5.09901951359278, 1.41421356237310],
[3.60555127546399, 0., 2.23606797749979, 6.32455532033676, 2.23606797749979, 3.60555127546399],
[2.00000000000000, 2.23606797749979, 0., 8.06225774829855, 3.16227766016838, 1.41421356237310],
[8.06225774829855, 6.32455532033676, 8.06225774829855, 0., 8.06225774829855, 9.00000000000000],
[5.09901951359278, 2.23606797749979, 3.16227766016838, 8.06225774829855, 0., 4.47213595499958],
[1.41421356237310, 3.60555127546399, 1.41421356237310, 9.00000000000000, 4.47213595499958, 0.]]):

n := sqrt(numelems(c)):  N:={seq(1..n)}:
z:=add(add( c[i,j]*x[i,j],j=N minus {i}),i=N): 
conx:=seq( add(x[i,j], i=N minus {j})=1,j=N),  seq( add(x[i,j], j=N minus {i})=1,i=N):
conu:= seq(seq(u[i]-u[j]+n*x[i,j] <= n-1, i=N minus {1,j}), j=N minus {1}):

r := Optimization[LPSolve](z, {conx, conu}, #integervariables={seq(u[i],i=N minus {1})},
                           binaryvariables={ seq(seq(x[i,j],i=N minus {j}),j=N)} );

X:=eval(Matrix(n, symbol=x),{r[2][], seq(x[i,i]=0,i=1..n)}):
ord:=1: k:=1: to n do k:=max[index](X[k]); ord:=ord,k od:'order'=ord;

r := [22.6135858310497, [u[2] = -.999999998967403, u[3] = -2.99999997900383, u[4] = -8.88178419700125*10^(-16), u[5] = -2.00000000722817, u[6] = -3.99999998106902, x[1, 2] = 0, x[1, 3] = 0, x[1, 4] = 0, x[1, 5] = 0, x[1, 6] = 1, x[2, 1] = 0, x[2, 3] = 0, x[2, 4] = 1, x[2, 5] = 0, x[2, 6] = 0, x[3, 1] = 0, x[3, 2] = 0, x[3, 4] = 0, x[3, 5] = 1, x[3, 6] = 0, x[4, 1] = 1, x[4, 2] = 0, x[4, 3] = 0, x[4, 5] = 0, x[4, 6] = 0, x[5, 1] = 0, x[5, 2] = 1, x[5, 3] = 0, x[5, 4] = 0, x[5, 6] = 0, x[6, 1] = 0, x[6, 2] = 0, x[6, 3] = 1, x[6, 4] = 0, x[6, 5] = 0]]

order = (1, 6, 3, 5, 2, 4, 1)

If you need to read the same data several times you could do this:

1. Import once using ImportMatrix into A
2.    save A, "hugeA.m";
3.  Whenever you need A, use
     read "hugeA.m";
     It will be executed much faster (being in internal format .m).

The result is wrong, the implicitplot algorithm does not work well here.




Probably some "clever" options could be added, but always a better idea is to try a parametric representation of the curve:

restart;
f:=(1/2)*x*tan(log((1/5)*(2*(x^2+y^2)))) = y:
g:=simplify(eval(f, [x=r*cos(t), y=r*sin(t)])) assuming r>0:
u:=solve(g,t);

plot( [[r*cos(u),r*sin(u), r=0..5],[-r*cos(u),-r*sin(u), r=0..5]]);