The inscribed square problem, also known as the Toeplitz conjecture, is an unsolved quastion in geometry: Does every plane simple closed curve (Jordan curve) contain all four vertices of some square? This is true if the curve is convex or piecewise smooth and in other special cases. The problem was proposed by Otto Toeplitz in 1911. For detailes see  https://en.wikipedia.org/wiki/Inscribed_square_problem

The Inscribed_Square procedure finds numerically one or more solutions for a curve defined by parametric equations of its boundary or by the equation F(x,y)=0. The required parameter of procedure  L  is the list of equations of the boundary links or the equation  F(x,y)=0 . Optional parameters:  N  and  R . By default  N='onesolution' (the procedure finds one solution), if  N  is any symbol (for example  N='s'), then more solutions.  R  is the range for the length of the side of the square (by defalt  R=0.1..100 ).

The second procedure  Pic  visualizes the results obtained.

The codes of the procedures:

restart;
Inscribed_Square:=proc(L::{list(list),`=`},N::symbol:='onesolution',R::range:=0.1..100)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, M, eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, sol, Sol;
uses LinearAlgebra;
if L::list then
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
M:=<0,-1;1,0>;
eq1:=eval(C[1])=eval((B+M.(B-A))[1]);
eq2:=eval(C[2])=eval((B+M.(B-A))[2]);
eq3:=eval(D[1])=eval((C+M.(C-B))[1]);
eq4:=eval(D[2])=eval((C+M.(C-B))[2]);
eq5:=eval(DotProduct(B-A,B-A, conjugate=false))=d^2;
sol:=fsolve([eq1,eq2,eq3,eq4,eq5],{t1=op([2,2,1],L0[i])..op([2,2,2],L0[i]),t2=op([2,2,1],L0[j])..op([2,2,2],L0[j]),t3=op([2,2,1],L0[k])..op([2,2,2],L0[k]),t4=op([2,2,1],L0[m])..op([2,2,2],L0[m]),d=R});
if type(sol,set(`=`)) then if N='onesolution' then return convert~(eval([A,B,C,D],sol),list) else c:=c+1; Sol[c]:=convert~(eval([A,B,C,D],sol),list) fi;
 fi; 
od: od: od: od:
Sol:=fnormal(convert(Sol,list),7);
print(Sol);
ListTools:-Categorize((X,Y)->`and`(seq(is(convert(X,set)[i]=convert(Y,set)[i]),i=1..4)) , Sol);
return map(t->t[1],[%]);
else
A,B,C,D:=<x1,y1>,<x2,y2>,<x3,y3>,<x4,y4>:
M:=<0,-1;1,0>:
eq1:=eval(C[1])=eval((B+M.(B-A))[1]):
eq2:=eval(C[2])=eval((B+M.(B-A))[2]):
eq3:=eval(D[1])=eval((C+M.(C-B))[1]):
eq4:=eval(D[2])=eval((C+M.(C-B))[2]):
eq5:=eval(LinearAlgebra:-DotProduct((B-A,B-A), conjugate=false))=d^2:
eq6:=eval(L,[x=x1,y=y1]):
eq7:=eval(L,[x=x2,y=y2]):
eq8:=eval(L,[x=x3,y=y3]):
eq9:=eval(L,[x=x4,y=y4]):
sol:=fsolve({eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9},{seq([x||i=-2..2,y||i=-2..2][],i=1..4),d=R});
eval([[x1,y1],[x2,y2],[x3,y3],[x4,y4]], sol):
fi;
end proc:

Pic:=proc(L,Sol,R::range:=-20..20)
local P1, P2, P3, T;
uses plots, plottools;
P1:=`if`(L::list,seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L), implicitplot(L, x=R,y=R, color=blue, thickness=2, gridrefine=3));
P2:=polygon(Sol,color=yellow,thickness=0);
P3:=curve([Sol[],Sol[1]],color=red,thickness=3):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
display(P1,P2,P3,T, scaling=constrained, size=[800,500], axes=none);
end proc:

Examples of use:

The curve consists of a semicircle, a segment and a semi-ellipse (find 1 solution):

L:=[[[cos(t),sin(t)],t=0..Pi],[[t,0],t=-1..0],[[0.5+0.5*cos(t),0.8*sin(t)],t=Pi..2*Pi]]:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

The procedure finds 6 solutions for a non-convex pentagon:

 L:=[[[0,0],[9,0]],[[9,0],[8,5]],[[8,5],[5,3]],[[5,3],[0,4]],[[0,4],[0,0]]]:
Sol:=Inscribed_Square(L,'s');
plots:-display(Matrix(3,2,[seq(Pic(L,Sol[i]),i=1..6)]),size=[300,200]);

For an implicitly defined curve, only one solution can be found:

L:=abs(x)+2*abs(y)-sin((2*x-y))-cos(x+y)^2=3:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

               
See more examples in the attached file.

Inscribed_Square.mw

Hi Maplesoft Support / Community,

I've encountered a critical and bizarre bug involving Bits:-And correctness on large integers (~30 digits) derived from repeated integerdivq2exp operations.

• Maple 2023 (Linux x86_64)
• Maple 2025 (Linux x86_64)
• Maple 2025 (Windows x86_64)

The correctness of Bits:-And depends on the order of execution. 

(See attached common.mpl, bug_test2.mpl, bug_test3.mpl logic).

Case "Fail" (bug_test2.mpl):

1. Run operation (loops `integerdivq2exp`).
2. Print result num1 (semicolon).
3. Define num1_clean (hardcoded same value).
4. Bits:-And(num1) -> INCORRECT.
5. Bits:-And(num1_clean) -> INCORRECT.

Case "Pass" (bug_test3.mpl):

1. Define num1_clean.
2. Run operation (loops integerdivq2exp).
3. Bits:-And(num1) -> CORRECT.
4. Bits:-And(num1_clean) -> CORRECT.

The same behaviour can be observed in Worksheet mode using read.  (See worksheet_driver.mw)

But the result cannot be reproduced if not using read. (See worksheet_version.mw and worksheet_version2.mw)

Code below:

N := 2100:
n := 1000:
num := rand(0 .. 2^N)():
operation := proc(num, n)
    local q, k;
    q := num;
    for k from 1 to 2 do
        q := integerdivq2exp(q, n); 
    end do;
    q;
end proc:

read "common.mpl";

num1 := operation(num, n);
num1_clean := 1083029963437854242395921050992;

num1_clean_And_result := Bits:-And(num1_clean, integermul2exp(1, n) - 1);
num1_And_result := Bits:-And(num1, integermul2exp(1, n) - 1);

##################################

expected_result := irem(num1_clean, integermul2exp(1, n));

if num1 <> num1_clean then
    error "num1 does not match num1_clean";
end if;
print("num1 matches num1_clean");

if num1_And_result <> num1_clean_And_result then
    error "num1_And_result does not match num1_clean_And_result";
end if;
print("num1_And_result matches num1_clean_And_result");

if num1_And_result <> expected_result then
    error "num1_And_result does not match expected_result";
end if;
print("num1_And_result matches expected_result");

read "common.mpl";

num1_clean := 1083029963437854242395921050992:
num1 := operation(num, n):

num1_clean_And_result := Bits:-And(num1_clean, integermul2exp(1, n) - 1):
num1_And_result := Bits:-And(num1, integermul2exp(1, n) - 1);

##################################

expected_result := irem(num1_clean, integermul2exp(1, n));

if num1 <> num1_clean then
    error "num1 does not match num1_clean";
end if;
print("num1 matches num1_clean");

if num1_And_result <> num1_clean_And_result then
    error "num1_And_result does not match num1_clean_And_result";
end if;
print("num1_And_result matches num1_clean_And_result");

if num1_And_result <> expected_result then
    error "num1_And_result does not match expected_result";
end if;
print("num1_And_result matches expected_result");