The incenter of a quadrilateral is a rational function of the vertex coordinates

When the system of two quadratic equations in $x$ has exactly one common root, the 1D kernel of the Sylvester matrix provides a rational expression for the root. \begin{cases}x^{2}+a_{1}x+a_{0}=0\\ x^{2}+b_{1}x+b_{0}=0\end{cases} Subtracting,$$(a_1-b_1)x+a_0-b_0=0\implies x=\frac{a_0-b_0}{b_1-a_1}$$

For two quadratic polynomials, $P(x)=x^{2}+a_{1}x+a_{0}$ and $Q(x)=x^{2}+b_{1}x+b_{0}$, the Sylvester matrix is a 4x4 matrix: $S = \begin{pmatrix} 1 & a_1 & a_0 & 0 \\ 0 & 1 & a_1 & a_0 \\ 1 & b_1 & b_0 & 0 \\ 0 & 1 & b_1 & b_0 \end{pmatrix}$

If $P(x)$ and $Q(x)$ have exactly one common root (which is a simple root of at least one polynomial), the resultant is zero, and the Sylvester matrix is singular. Its kernel will be one-dimensional. The vectors in this 1D kernel are related to the powers of the common root. Every column of the adjugate matrix $\text{Adj }S$ is proportional to the same vector $[x^3,x^2,x,1]^T$. Consecutive ratios evaluate to the same scalar$$\frac{(\text{Adj } S)_{1,i}}{(\text{Adj } S)_{2,i}}=\frac{(\text{Adj } S)_{2,i}}{(\text{Adj } S)_{3,i}}=\frac{(\text{Adj } S)_{3,i}}{(\text{Adj } S)_{4,i}}=x$$for every row $i=1,2,3,4$.

Statement being verified:

For the Sylvester matrix $$S= \begin{pmatrix} 1 & a_1 & a_0 & 0\\ 0 & 1 & a_1 & a_0\\ 1 & b_1 & b_0 & 0\\ 0 & 1 & b_1 & b_0 \end{pmatrix},$$ let $$x=\frac{a_0-b_0}{\,b_1-a_1\,}.$$ Claim. For each column $i=1,2,3,4$, $$\frac{(\operatorname{Adj}S)_{1,i}}{(\operatorname{Adj}S)_{2,i}} = \frac{(\operatorname{Adj}S)_{2,i}}{(\operatorname{Adj}S)_{3,i}} = \frac{(\operatorname{Adj}S)_{3,i}}{(\operatorname{Adj}S)_{4,i}} = x$$ in the quotient ring $$\mathbb{Q}[a_0,a_1,b_0,b_1]\big/\bigl(\det S\bigr).$$ What was computed: For each column and each consecutive ratio, we checked the cross-multiplied form $$(\operatorname{Adj}S)_{k,i}(b_1-a_1) - (\operatorname{Adj}S)_{k+1,i}(a_0-b_0),$$ and reduced this polynomial modulo the ideal $(\det S)$ using a Gröbner basis.

# Verify symbolically that (Adj S) column ratios equal x modulo det(S)
# by reducing the differences with a Groebner basis containing det(S)

import sympy as sp

# symbols
a0, a1, b0, b1 = sp.symbols('a0 a1 b0 b1')

# matrix
S = sp.Matrix([
    [1, a1, a0, 0],
    [0, 1,  a1, a0],
    [1, b1, b0, 0],
    [0, 1,  b1, b0]
])

Adj = S.adjugate()

# determinant
detS = sp.factor(S.det())

# target x
x = (a0 - b0)/(b1 - a1)

# polynomial ring generators
G = sp.groebner([detS], a0, a1, b0, b1, order='lex')

results = {}

for j in range(4):
    c = Adj.col(j)

    # three consecutive ratios
    diffs = []
    for (num, den) in [(c[0], c[1]), (c[1], c[2]), (c[2], c[3])]:
        # cross-multiplied difference: num/den - x = 0  <=>  num*(b1-a1) - den*(a0-b0)
        poly = sp.expand(num*(b1 - a1) - den*(a0 - b0))
        # reduce modulo det(S)
        reduced = G.reduce(poly)[1]
        diffs.append(sp.factor(reduced))

    results[j+1] = diffs

detS, results

SymPy returns:

{1: [0, 0, 0], 
 2: [0, 0, 0], 
 3: [0, 0, 0], 
 4: [0, 0, 0]}

Similarly, when the system of three quadratic equations in $x,y$ has exactly one common root, the 1D kernel of the Macaulay Matrix (the multivariate generalization of the Sylvester matrix) provides a rational expression for the root.

Let $L_i(x, y) = A_i x + B_i y + C_i = 0$ be the side lines and $N_i = A_i^2 + B_i^2$ be the squared side lengths. Define three quadratics with rational coefficients: $$f_1 = N_2 L_1^2 - N_1 L_2^2, \quad f_2 = N_3 L_2^2 - N_2 L_3^2, \quad f_3 = N_4 L_3^2 - N_3 L_4^2$$ Each $f_k$ has the form $a_{k1}x^2 + a_{k2}xy + a_{k3}y^2 + a_{k4}x + a_{k5}y + a_{k6} = 0$.

The Macaulay Matrix

We construct a matrix $M$ of size $9 \times 10$ by shifting each polynomial $f_k$ by $x$ and $y$. The 10 columns correspond to the monomials $[x^3, x^2y, xy^2, y^3, x^2, xy, y^2, x, y, 1]$. For each $f_k$, we generate 3 rows:

1. Row 1 ($x \cdot f_k$): $[a_{k1}, a_{k2}, a_{k3}, 0, a_{k4}, a_{k5}, 0, a_{k6}, 0, 0]$
2. Row 2 ($y \cdot f_k$): $[0, a_{k1}, a_{k2}, a_{k3}, 0, a_{k4}, a_{k5}, 0, a_{k6}, 0]$
3. Row 3 ($f_k$): $[0, 0, 0, 0, a_{k1}, a_{k2}, a_{k3}, a_{k4}, a_{k5}, a_{k6}]$

Rational Formula

For a tangential quadrilateral, the matrix $M$ has rank 9, meaning its null space (kernel) is one-dimensional. Let $\mathbf{k} = [k_9, k_8, \dots, k_1, k_0]^T$ be a non-zero vector in the kernel. By Cramer’s rule, the components $k_i$ are the $9 \times 9$ minors of $M$. The incenter $(x_I, y_I)$ is: $$x_I = \frac{k_2}{k_0} = \frac{\text{Minor}_{x}}{\text{Minor}_{1}}, \quad y_I = \frac{k_1}{k_0} = \frac{\text{Minor}_{y}}{\text{Minor}_{1}}$$

Numerical Verification

import numpy as np
from scipy.linalg import null_space

rng = np.random.default_rng(0)
Ix, Iy = rng.uniform(-3,3), rng.uniform(-3,3)
r = rng.uniform(0.5,2.0)
thetas = np.sort(rng.uniform(0, 2*np.pi, 4))

lines = []
for t in thetas:
    A, B = np.cos(t), np.sin(t)
    C = -(A*Ix + B*Iy + r)
    lines.append((A,B,C))

def get_coeffs(l1, l2):
    A1, B1, C1 = l1
    A2, B2, C2 = l2
    return [
        A1**2 - A2**2,
        2*(A1*B1 - A2*B2),
        B1**2 - B2**2,
        2*(A1*C1 - A2*C2),
        2*(B1*C1 - B2*C2),
        C1**2 - C2**2
    ]

f1 = get_coeffs(lines[0], lines[1])
f2 = get_coeffs(lines[1], lines[2])
f3 = get_coeffs(lines[2], lines[3])

def get_macaulay_rows(f):
    r1 = [f[0], f[1], f[2], 0, f[3], f[4], 0, f[5], 0, 0]
    r2 = [0, f[0], f[1], f[2], 0, f[3], f[4], 0, f[5], 0]
    r3 = [0, 0, 0, 0, f[0], f[1], f[2], f[3], f[4], f[5]]
    return [r1, r2, r3]

M = np.array(get_macaulay_rows(f1) + get_macaulay_rows(f2) + get_macaulay_rows(f3))
print("Matrix Rank:", np.linalg.matrix_rank(M))

K = null_space(M)
print("Kernel dimension (nullity):", nullity)

if K.shape[1] == 1:
    k = K[:, 0]
    xi = k[7]/k[9]
    yi = k[8]/k[9]
    print(f"Recovered Incenter: ({xi}, {yi})")
    print(f"Actual Incenter:   ({Ix}, {Iy})")

3D Generalization

For 5 planes (4 quadratic equations) in 3D, the system becomes uniquely determined (nullity 1) when expanded to degree 4. This requires 35 columns (all monomials from $x^4$ down to 1) and 40 rows (each of the 4 equations shifted by the 10 possible monomials of degree $\leq 2$).

The following numerical verification confirms in 3D, the point equidistant to 5 planes can be expressed as a rational function of the coefficients of the planes.

import numpy as np
from scipy.linalg import null_space

def get_monomials(max_deg):
    """Generates all 3D monomials (i,j,k) s.t. i+j+k <= max_deg"""
    mons = []
    for d in range(max_deg, -1, -1):
        for i in range(d, -1, -1):
            for j in range(d - i, -1, -1):
                mons.append((i, j, d - i - j))
    return mons

# 1. Generate 5 random planes tangent to a sphere
rng = np.random.default_rng(1)
Ix, Iy, Iz = rng.uniform(-2, 2, 3)
r = rng.uniform(0.5, 1.5)
normals = [n/np.linalg.norm(n) for n in rng.normal(size=(5, 3))]
planes = [np.append(n, -(n[0]*Ix + n[1]*Iy + n[2]*Iz + r)) for n in normals]

# 2. Define the 4 base quadratic equations f_i = L_0^2 - L_i^2
eqs_base = []
L0 = planes[0]
for i in range(1, 5):
    Li = planes[i]
    f = {
        (2,0,0): L0[0]**2 - Li[0]**2, (0,2,0): L0[1]**2 - Li[1]**2, (0,0,2): L0[2]**2 - Li[2]**2,
        (1,1,0): 2*(L0[0]*L0[1] - Li[0]*Li[1]), (1,0,1): 2*(L0[0]*L0[2] - Li[0]*Li[2]), 
        (0,1,1): 2*(L0[1]*L0[2] - Li[1]*Li[2]), (1,0,0): 2*(L0[0]*L0[3] - Li[0]*Li[3]),
        (0,1,0): 2*(L0[1]*L0[3] - Li[1]*Li[3]), (0,0,1): 2*(L0[2]*L0[3] - Li[2]*Li[3]),
        (0,0,0): L0[3]**2 - Li[3]**2
    }
    eqs_base.append(f)

# 3. Build the Macaulay Matrix (Degree 4 expansion)
# Basis: 35 monomials. Rows: 4 eqs * 10 shifts
mons_deg4 = get_monomials(4)
mon_to_idx = {m: i for i, m in enumerate(mons_deg4)}
shifts = get_monomials(2) # 10 shifts

M = []
for f in eqs_base:
    for s in shifts:
        row = np.zeros(35)
        for m, val in f.items():
            target = (m[0]+s[0], m[1]+s[1], m[2]+s[2])
            row[mon_to_idx[target]] = val
        M.append(row)
M = np.array(M)

# 4. Extract Incenter
K = null_space(M)
print(f"Matrix: {M.shape}, Rank: {np.linalg.matrix_rank(M)}, Nullity: {K.shape[1]}")

if K.shape[1] == 1:
    k = K[:, 0]
    idx_x, idx_y, idx_z, idx_1 = mon_to_idx[(1,0,0)], mon_to_idx[(0,1,0)], mon_to_idx[(0,0,1)], mon_to_idx[(0,0,0)]
    res = (k[idx_x]/k[idx_1], k[idx_y]/k[idx_1], k[idx_z]/k[idx_1])
    print(f"Recovered: {res}\nActual:    {(Ix, Iy, Iz)}")

This raises the question:

Can you simplify the rational expression for the point equidistant to 5 planes, expressed in a suitable coordinate system, to reveal a geometric construction analogous to the 2D case?