Just one region occludes the other.

Workaround:

with(plots):

I1 := inequal({x < 1, y < 1}, x = 0 .. 1, y = 0 .. 1, optionsfeasible = [color = "Brown", thickness = 2], view = [0 .. 2, 0 .. 2]):

I2 := inequal({y > 1, x < 1}, x = 0 .. 1, y = 1 .. 2, optionsfeasible = [color = "Yellow", thickness = 2], view = [0 .. 2, 0 .. 2]):

display(I1);

display(I2);

display({I1, I2});

restart;

A:=(1/a[x])^(-k-1+epsilon)/(1/a[I])^(-k-1+epsilon);

Rule:=(x::anything)^(z::anything)/(y::anything)^(z::anything)=(x/y)^z:

B:=applyrule(Rule, A);

subs({op(1,B)=op(1,B)^(-1), op(2,B)=-op(2,B)}, B);  # Result

For example

simplify((a+b)^2, {a^2+b^2=1, a*b=1/2});

                                 2

Addition. The above example is trivial. Here is less obvious example:

simplify(a^20+b^20, {a+b=5, a*b=1});

                       40648568638127

Idea of ​​the proof:
If  x*(x +1)  is divisible by 100!  then because the numbers  x  and  x+1  are relatively prime,  in the expansion of  100!  into prime factors any  p^r  can belong to only one of the factors  x  or  x+1.

restart;

A:={op(ifactors(100!)[2])}:

B:=[seq(op(combinat[choose](A,k)), k=1..7)]:

C:=[seq([B[i], A minus B[i]], i=1..nops(B))]:

for c in C do

a:=mul(c[1,i,1]^c[1,i,2], i=1..nops(c[1])): b:=mul(c[2,i,1]^c[2,i,2], i=1..nops(c[2])):

isolve(a*u-b*v=1): x:=tcoeff(b*rhs(%[2])): if x<4*99! then break fi:

isolve(b*v-a*u=1): x:=tcoeff(a*rhs(%[1])): if x<4*99! then break fi:

od:

x;

Verification:

is(irem(x*(x+1), 100!)=0 and x<4*99!);

                               true

This method lets you easily find the smallest  x  satisfying the given conditions, but it will take much more time.

You can use  plot3d  command:

A:=mtaylor(sin(x*y), [x=1,y=2], 6);

plot3d(A, x=0..3, y=0..3, view=-2..2, axes=normal, scaling=constrained);

Here is an independent solution of the logical problem. Code is quite complicated, because based on the phase-solving, but it works fast. Carl, do  you get the same answers in your package?

restart;

Position:=[$ 1..5]:

Names:= [Bob, Keeley, Rachael, Eilish, Amy]:

TV:=[Simpsons, Coronation, Eastenders, Desperate, Neighbours]:

Dest:= [Fra, Aus, Eng, Afr,Ita]:

Ages:= [14, 21, 46, 52, 81]:

Hair:=[afro, long, straight, curly , bald]:

Lives:= [town, city, village, farm, youth]:

P:=proc(A,B)

if (A implies B) and (B implies A) then true else false fi;

end:

L:=[]:

for i in Position do

for j in Names do

for k in TV do

for p in Dest do

for q in Ages do

for s in Hair do

for t in Lives do

if convert([P(i=3,k=Desperate), P(i=1,j=Bob), (k=Simpsons implies t<>youth) and (t=youth implies k<>Simpsons), (p=Afr implies j<>Rachael) and (j=Rachael implies p<>Afr), P(t=village,q=52), P(p=Aus,s=straight), P(p=Afr,k=Desperate), P(q=14,i=5), P(j=Amy, k=Eastenders), P(p=Ita,s=long), P(j=Keeley,t=village), P(q=46,s=bald), P(i=4, p=Eng),  (k=Coronation implies s<>afro) and (s=afro implies k<>Coronation), (j=Rachael implies s<>afro) and (s=afro implies j<>Rachael), P(q=21,t=youth), P(k=Coronation, s=long), P(q=81,t=farm), P(p=Fra,t=town), (j=Eilish implies s<>straight) and (s=straight implies j<>Eilish)], `and`) then L:=[op(L), [i,j,k,p,q,s,t]] fi:

od: od: od: od: od: od: od:

M:=[ListTools[Categorize]((x,y)->x[1]=y[1], L)]:

R:=[]:

for i in M[1] do

for j in M[2] do

for k in M[3] do

for p in M[4] do

for q in M[5] do

if ((i[3]=Simpsons and j[7]=youth) or (j[3]=Simpsons and (i[7]=youth or k[7]=youth)) or (k[3]=Simpsons and (j[7]=youth or p[7]=youth)) or (p[3]=Simpsons and (k[7]=youth or q[7]=youth)) or (q[3]=Simpsons and p[7]=youth)) and ((i[2]=Rachael and (j[4]=Afr or k[4]=Afr or p[4]=Afr or q[4]=Afr)) or (j[2]=Rachael and (k[4]=Afr or p[4]=Afr or q[4]=Afr)) or (k[2]=Rachael and (p[4]=Afr or q[4]=Afr)) or (p[2]=Rachael and q[4]=Afr)) and ((i[3]=Desperate and j[3]=Neighbours) or (j[3]=Desperate and (i[3]=Neighbours or k[3]=Neighbours)) or (k[3]=Desperate and (j[3]=Neighbours or p[3]=Neighbours)) or (p[3]=Desperate and (k[3]=Neighbours or q[3]=Neighbours)) or (q[3]=Desperate and p[3]=Neighbours)) and ((i[3]=Coronation and j[6]=afro) or (j[3]=Coronation and (i[6]=afro or k[6]=afro)) or (k[3]=Coronation and (j[6]=afro or p[6]=afro)) or (p[3]=Coronation and (k[6]=afro or q[6]=afro)) or (q[3]=Coronation and p[6]=afro)) and ((i[2]=Rachael and j[6]=afro) or (j[2]=Rachael and (i[6]=afro or k[6]=afro)) or (k[2]=Rachael and (j[6]=afro or p[6]=afro)) or (p[2]=Rachael and (k[6]=afro or q[6]=afro)) or (q[2]=Rachael and p[6]=afro)) and not ((i[2]=Eilish and j[6]=straight) or (j[2]=Eilish and (i[6]=straight or k[6]=straight)) or (k[2]=Eilish and (j[6]=straight or p[6]=straight)) or (p[2]=Eilish and (k[6]=straight or q[6]=straight)) or (q[2]=Eilish and p[6]=straight))  then

R:=[op(R), [i,j,k,p,q]] fi:

od: od: od: od: od:

Result:=select(x->convert([seq(nops({seq(x[i,j], i=1..5)})=5, j=1..7)],`and`), R):

print(`First solution`);

Matrix(Result[1])^+;

print(``);

print(`Second solution`);

Matrix(Result[2])^+; 

Puzzle.mws

with(LinearAlgebra):

a, b, c:=<22,11,5>, <13,6,3>, <-5,-2,-1>:

a1:=a/norm(a,2):

b0:=b-(a1.b)*a1:

b1:=b0/norm(b0,2):

c0:=c-((a1.c)*a1+(b1.c)*b1):

c1:=c0/norm(c0,2):

a1,  b1,  c1;

Vectors  a1, b1, c1  are unit and mutually perpendicular:

a1.a1, b1.b1, c1.c1, a1.b1, a1.c1, b1.c1;

            1, 1, 1, 0, 0, 0

convert(op(0.25)[1], symbol)/(10^(-op(0.25)[2]));

 ``(op(0.25)[1])/``(10^(-op(0.25)[2]));

The second method allows you to use this fraction in subsequent calculations:

expand(%*3/5);

              3/20

Solution without  dsolve (only for rectangular tables):

restart;

A:=plots[animate](plot,[[2/Pi*abs(arcsin(sin(Pi/2*(0.5+0.2*t)))), 2/Pi*abs(arcsin(sin(Pi/2*(0.5+sqrt(3)/5*t)))), t=0..a], thickness=5, color=red, numpoints=1000], a=0..20, frames=200):

B:=plottools[rectangle]([0,1], [1,0], thickness=2, color=green):

plots[display](A, B, scaling=constrained, axes=none);

with(RealDomain):

 plot(x^(1/3), x=-8 .. 8, thickness=5, scaling=constrained);

For Classic M 12:

A:=1/2+(1/3)*sqrt(2);

numer(A)/(``*denom(A));

For Standard:

A := 1/2+(1/3)*sqrt(2);

numer(A)/convert(denom(A), symbol);

Another variant:

A := 1/2+(1/3)*sqrt(2);

``(numer(A))/denom(A);

For your example it's more reliable to use a linear model, which is obtained if to apply  ln  function to your equation:

ln(x)=ln(a)+b*ln(y)+c*ln(z)  or  x1=a1+b*y1+c*z1 , where  x1=ln(x) ,  a1=ln(a) ,  y1=ln(y) ,  z1=ln(z) 

xlist := [0.1e-1, 0.1e-1, 0.1e-1, 0.113e-2, 0.144e-2, 0.171e-2, 0.21e-2, 0.23e-2, 0.183e-2, 0.233e-2, 0.28e-2, 0.331e-2, 0.409e-2, 0.304e-2, 0.377e-2, 0.459e-2, 0.576e-2, 0.708e-2, 0.424e-2, 0.502e-2, 0.623e-2, 0.763e-2, 0.918e-2, 0.529e-2, 0.63e-2, 0.778e-2, 0.965e-2, 0.116e-1, 0.63e-2, 0.747e-2, 0.899e-2, 0.113e-1, 0.137e-1, 0.724e-2, 0.872e-2, 0.105e-1, 0.13e-1, 0.158e-1, 0.848e-2, 0.1e-1, 0.122e-1, 0.15e-1, 0.182e-1]:

ylist := [16.529, 13.333, 10.929, 2.66, 2.75, 2.77, 2.69, 2.58, 4.4, 4.67, 4.62, 4.51, 4.48, 7.5, 7.58, 7.64, 7.75, 7.77, 10.1, 10.1, 10.1, 10.2, 10.2, 12.6, 12.7, 12.8, 12.7, 12.8, 14.9, 15, 15, 15, 15, 17.3, 17.5, 17.4, 17.3, 17.3, 20, 19.9, 19.9, 19.9, 20]:

zlist := [100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000]: 

X:=convert(map(ln,xlist), Vector):

Y:=convert(map(ln,ylist), Vector):

Z:=convert(map(ln,zlist), Vector):

YZ:=Matrix([Y,Z]):

Statistics[LinearFit]([1,y1,z1], YZ, X, [y1,z1]):

a=exp(tcoeff(%)), b=coeff(%,y1), c=coeff(%,z1);

            a = 0.0004189912543,  b =0 .996002908082386118, c = 0.0849219222885009006

L1:=x + y - 1 = 0, x - z - 1 = 0:

L2:=x + y - 7 = 0, x - y + 1 = 0:

A:=plots[intersectplot](L1, x=-5..5, y=-5..5, z=-3..5, thickness=2, color=red):

B:=plots[intersectplot](L2, x=-5..5, y=-5..5, z=-3..5, thickness=2, color=blue):

plots[display](A, B, scaling=constrained, axes=normal);

For your data, you can find exactly unknown a, b, c:

X:=[0.01, 0.01, 0.01]:

Y:=[16.526, 13.333, 10.929]:

Z:=[100, 1000, 10000]:

f:=x=a*y^b*z^c:

solve({seq(eval(f, {x=X[i], y=Y[i], z=Z[i]}), i=1..3)});

                             {a = 0.01, b = 0., c = 0.}

Let the numbers assigned to each chair from 1 to 8. 

Seating for 8 people means to establish a one-to-one correspondence between these people and these numbers. Total we get  96 solutions, including 6 unique solutions:

restart;

Names := [Alan, Amy, Brad, Beth, Charles, Debbie, Emily, Frances]:

P := proc (n1, n2)

`if`(irem(abs(n2-n1), 7) <= 1, true, false) end proc:

Permutations := combinat[permute]([`$`(1 .. 8)]):

Arrangements := []:

for p in Permutations do

if eval(`and`(`and`(`and`(`and`(`and`(`and`(`and`(`and`(P(Amy, Alan), P(Brad, Beth)), `or`(P(Charles, Debbie), P(Charles, Emily))), P(Frances, Debbie)), `and`(`and`(`and`(`not`(P(Amy, Brad)), `not`(P(Amy, Beth))), `not`(P(Alan, Brad))), `not`(P(Alan, Beth)))), `and`(`not`(P(Brad, Charles)), `not`(P(Brad, Frances)))), `not`(P(Debbie, Emily))), `and`(`not`(P(Alan, Debbie)), `not`(P(Alan, Emily)))), `not`(P(Amy, Charles))), {seq(Names[i] = p[i], i = 1 .. 8)})then Arrangements := [op(Arrangements), p] end if end do:

seq([seq([Names[i], p[i]], i = 1 .. 8)], p = Arrangements);

Arrangement_Birthda.mw