The shape of a quadrilateral inscribed in the unit circle and having an inscribed circle.

f := (a,b) -> 1/2*(b+arccos(sin(2*a-b)/tan(b))):
dis:=(A,B)-> sqrt(inner(A-B,A-B)):
bisA:=proc(A,B,C)
  local b,c,M;
  b,c:=dis(A,C),dis(A,B);
  M:=(b*B+c*C)/(b+c);
  x*(A[2]-M[2])+y*(M[1]-A[1]) + A[1]*M[2]-A[2]*M[1]
end:
P:=proc(a0,b0)
  local a:=evalf(a0), b:=evalf(b0), c:=f(a,b), p1,p2,p3,p4, r, II;
  if a0>b0-0.001 then return plot() fi;
  p1,p2,p3,p4:=[1,0],[cos(4*a),sin(4*a)],[cos(4*b),sin(4*b)],[cos(4*c),sin(4*c)];
  II:= eval([x,y], solve( [bisA(p2,p1,p3), bisA(p3,p2,p4)]) );
  r:=evalf(eval(dis(II,p2)*cos(b))); 
  plots:-display(plot([p1,p2,p3,p4,p1]), plottools:-circle(II,r, color=blue), plottools:-circle([0,0],1));
end:
Explore(P(a,b), parameters=[a=0.01..Pi/2-0.02, b=0.02..Pi/2-0.01], initialvalues=[a=Pi/8,b=Pi/4]);

###    0 < a < b < Pi/2

restart;
PP:=proc(n::posint, kmin::posint, L::list(prime), nSolMax::posint)
# Prime partitions with elements >= kmin, starting with L
local Sol:=table(), nSol:=0, E;

E:=proc(L)
local p,sL:=add(L);
if L<>[] then p:=L[-1] else p:=nextprime(kmin-1) fi;
while p+sL<=n and nSol<=nSolMax do
  if p+sL=n and nSol<nSolMax then Sol[++nSol]:=[L[],p] fi;
  if nSol>=nSolMax then return fi;  
  E([L[],p]);
  p:=nextprime(p);
od;
end:

if L<>[] and add(L)=n then  Sol[++nSol]:=L fi;   
E(L);
entries(Sol,nolist);
end:

SingPrime:=proc(p::prime,n::posint)  
# there exist only 1 PP starting with p
evalb(nops([PP(n,p,[p],2)])=1)
end:

Q := proc(n)    # list of SingPrime's
local p, P:=select(isprime, [seq(1..n)]);
select(SingPrime, P, n)
end:

SingPrime(3,8), SingPrime(5,8);
                          true, false
for n from 2 to 500 do
    if Q(n)=[] then printf("%d ",n) fi
od;
63 161 195 235 253 425

(The procedures are not optimized.)

An implementation of Lehmer's algorithm.

Lpi := proc(x)  # Lehmer pi
  local phi, Lpi1, M:=10^7 ;
    
  phi := proc(x, a)
  local t;
  if a <= 1 then  return iquo(x + 1, 2) fi;
  t := thisproc(x, a-1) - thisproc(iquo(x, ithprime(a)), a-1);
  procname(x, a) := t;
  t;
  end;

  Lpi1 := proc(x)
  local a,b,c,s,w,i,j,L, xx:=trunc(x);
  if (trunc(x) < M) then return NumberTheory:-pi(x) fi;
  a := thisproc(iroot(xx, 4));;
  b := thisproc(isqrt(xx));
  c := thisproc(iroot(xx, 3));
  s := phi(xx,a) + (b+a-2) * (b-a+1) / 2;
  for i from a+1 to b do
    w := iquo(xx , ithprime(i));
    L := thisproc(isqrt(w));
    s := s - thisproc(w);
    if (i <= c) then
      for j from i to L do
        s := s - thisproc(iquo(w , ithprime(j))) + j - 1
      od
    fi
  od;
  s
  end;

Lpi1(x)
end:

CodeTools:-Usage(Lpi(6469693230));
memory used=6.11GiB, alloc change=80.00MiB, cpu time=55.19s, real time=50.11s, gc time=10.83s
                           300369796

Download IntEq-vv.mw

Use add instead of sum i.e.

c := k -> piecewise(k = 0,b(0)^N,1 <= k,1/k/b(0)*add((q*(N+1)-k)*b(q)*c(k-q),q = 1.. k));

restart: 
with(GraphTheory):
L:=ImportGraph("yourpath/graph5c.g6", graph6, output=list):
DrawGraph(L);

The integrals are highly oscillating. t1 is convergent but not absolutely (so it does not exist as a Lebesgue integral).
Maple needs help from the user for such integrals.
I doubt Maple V was able to calculate them directly (and correctly).

Here is a special method to compute t1 (t2 goes similarly).

restart;
f1:= unapply( eval((1-exp(-(d-y)*I*x))/I/x*exp(y*I*x/(1-I*x)),  [d=3, y=9/10]), x):
F1:= u -> f1(1/u)/u^2:
h := t -> t+I*t*(1-t): # special change of variables
G1:=unapply(simplify(F1(h(t))*D(h)(t)), t):
t1:=evalf[25](2*Int(Re@f1, 0..1) + 2*Int(Re@G1, 0 .. 1));

       t1 := 3.418878800322797147227936

Edit: added a forgotten term.

Seems to be an old bug.

with(GraphTheory):
G:=Graph( Trail(1,2,3,4,5) ):
DrawGraph(G): # ok
DrawPlanar(G); # error

Probably the graph G being "collinear", it's "too planar"  :-)

Actually you do not have integral equations, but expressions containing integrals depending on the unknowns.
You should proceed as in the next example:

restart;
Digits:=15:
eq1 := proc(x,y) local t; int(sin(x*t)/(t^2+7), t=0..(x+y), numeric)  - 1/10. end:
eq2 := proc(x,y) local t; int(exp(x*t)/(t^2+3), t=0..(x-y), numeric) + 0.13 end:
solxy:=fsolve([eq1,eq2], [1 .. 2.5,  1.5 .. 4]);
          solxy := [1.474881121279997081543, 2.084760793500574599373]
eq1(solxy[]), eq2(solxy[]); # check
         -.4e-15, .1e-14

Download diff_conj.mw

plots:-display(seq(seq(plottools:-polygon([[m,n],[m+1/2,n+1/2],[m+1,n],[m+1/2,n-1/2]],
               color=`if`(m::odd,red,blue)), m=0..3),n=0..3), axes=none);

Note that ignoring the color, another possibility is

plots:-inequal(abs(x-floor(x)-1/2) + abs(y-floor(y)-1/2)<=1/2, x=0..4, y=0..4);

Yes, ChainTools  and also SolveTools:-SemiAlgebraic  are too slow for this problem.
But after inspection and some numeric experiments it seems that the system has solution iff:

( a=1  and  1 <= b <= sqrt(4/3) )    or    (a=b>1).

In general, such a brute-force solution is used only if a direct one is not available. This is not the case.
For t=0 the number of the roots of the equation is either <= 4 or infinity (so 6 is out of the question).

It remains t=1 which was solved in your previous post [or, just form the 4 quadratic equations (obtained from abs(u) = +- u) and find the conditions to have integer roots]. The case t=2 is similar.

I see two problems in your worksheet.

1. The distribution corresponding to the RV MIXTURE is not correctly defined. It should be:

pdf := (p,x) -> p*PDF(U1, x)*Heaviside(x+3) + (1-p)*PDF(U2, x)*Heaviside(x-1);
cdf := (p,x) -> int(pdf(p,y), y=-infinity..x);

MIXTURE := p -> RandomVariable( Distribution( PDF = (x -> pdf(p,x)), CDF = (x -> cdf(p,x)), Conditions=[p >=0, p <=1] )):

2. It seems that Sample works only for distributions for which the support is an interval and the PDF is differentiable here.
The support of  MIXTURE is (-3,-1) union (1,3).
[I don't know how to declare such a support, or even if it's possible].
So,  Sample(MIXTURE(1/3), 10);  will fail, but

Sample(MIXTURE(1/3), 10, method=[envelope, range=1..3]);  # works.

Probably it's possible to use somehow method=custom.