When you forget a theorem of geometry

... you should try Gröbner bases.

Suppose four points \(A,B,C,D\) lie on a common circle. How do the four quadrances to the point \(M:\equiv(AC)(BD)\) relate?



Let me Gröbner that for you.

% sage
┌──────────────────────┐
│ SageMath version 8.8 │
└──────────────────────┘
sage: R.<x0,y0,x1,y1,x2,y2,x3,y3,s,t,q0,q1,q2,q3> = PolynomialRing(QQ, order='lex')
sage: xm, ym, xm1, ym1 = x0*s+x2*(1-s), y0*s+y2*(1-s), x1*t+x3*(1-t), y1*t+y3*(1-t)

sage: ideal(x0^2+y0^2-1, x1^2+y1^2-1, x2^2+y2^2-1, x3^2+y3^2-1,

....:       xm-xm1, ym-ym1,

....:       -q0+(xm-x0)^2+(ym-y0)^2, -q1+(xm-x1)^2+(ym-y1)^2,

....:       -q2+(xm-x2)^2+(ym-y2)^2, -q3+(xm-x3)^2+(ym-y3)^2).groebner_basis()[-1]

q0*q2 - q1*q3

\[Q(A,M)Q(C,M)=Q(B,M)Q(D,M)\]

Sums of powers of roots

It is known since Newton's time that sums of the form \(\sigma_k:=\sum x_i^k\) are polynomial expressions in the coefficients of \(f=\prod(z-x_i)\). What follows is my best understanding of this. We use a formal power series approach, and therefore define \(N_f:=\sum\sigma_kz^k\). I want to prove that \[z^nf(\frac1z)N_f(z) = z^{n-1}f'(\frac1z).\label{th1}\tag{1}\] (It is written this way since \(z^nf(\frac1z)\) and \(z^{n-1}f'(\frac1z)\) are polynomials.)

The proof uses logarithmic differentiation and geometric series expansion only. \begin{align*} \frac{z^{n-1}f'(\frac1z)}{z^nf(\frac1z)} &= \frac1z(\log f)'\left(\frac1z\right) &&= \frac1z\sum\frac1{\frac1z-x_i} \\ &= \sum\frac1{1-x_iz} &&= \sum x_i^kz^k \\ &= \sum\sigma_kz^k &&= N_f(z) \end{align*} For example suppose \(f=(z-x_0)(z-x_1)=(z-2)(z-3)=z^2-5z+6\). The first few sums are \[\sigma_0 = 2,\quad\sigma_1=5,\quad\sigma_2=13,\quad\sigma_3=35,\quad\sigma_4=97\] so \(N_f\approx2+5z+13z^2+35z^3+97z^4\). Furthermore \(z^nf(\frac1z)=1-5z+6z^2\) and \(z^{n-1}f'(\frac1z)=2-5z\). We may test \((\ref{th1})\) by computing its left- and right-hand sides up to fourth order: \[(2+5z+13z^2+35z^3+97z^4)(1-5z+6z^2) = (2-5z) \mod z^5.\] Conversely, we can compute more terms \(\sigma_k\) by polynomial division: \[N_f = \frac{2-5z}{1-5z+6z^2} = \sum_{0\leq k\lt m}\sigma_kz^k + z^m\frac{\textrm{remainder polynomial}}{1-5z+6z^2}\] In python:

>>> z = sympy.var('z')
>>> Nf = sympy.fps((2-5*z)/(1-5*z+6*z*z))
>>> Nf.truncate(10)
2 + 5*z + 13*z**2 + 35*z**3 + 97*z**4 + 275*z**5
+ 793*z**6 + 2315*z**7 + 6817*z**8 + 20195*z**9 + O(z**10)

Demuxing a random variable

Suppose \(f = ag + (1-a)h\), a convex combination of pdfs. Then \(f\) is a pdf.

Suppose \(X\) has \(f\) as pdf. For any evaluation \(x\) of \(X\), produce a random variable with conditional distribution, \[Y\sim\text{Bernoulli}\left(\frac{(1-a)h(x)}{(ag(x)+(1-a)h(x)}\right).\] Then \(X|(Y=0)\sim g\) and \(X|(Y=1)\sim h\) since \begin{align} p(x|0) &= \frac{p(0|x)p(x)}{\sum_{x'}p(0|x')p(x')} = \frac{\frac{ag(x)}{f(x)}f(x)}{\sum_{x'}\frac{ag(x)}{f(x)}f(x)}\\ &= \frac{ag(x)}{\sum_{x'}ag(x')} = \frac{a\cdot g(x)}{a\cdot1} = g(x)\\ p(x|1) &= \frac{p(1|x)p(x)}{\sum_{x'}p(1|x')p(x')} = \frac{\frac{(1-a)h(x)}{f(x)}f(x)}{\sum_{x'}\frac{(1-a)h(x)}{f(x)}f(x)}\\ &= \frac{(1-a)h(x)}{\sum_{x'}(1-a)h(x')} = \frac{(1-a)\cdot h(x)}{(1-a)\cdot1} = h(x). \end{align}

The utility of this construction is as follows: It lets you turn a probabilistic experiment with definite outcomes (say the probabilities are \([0.2,\,0.6,\,0.2]\)) into an experiment with outcomes which are again probability distributions, such as \([0.4,\,0.6,\,0]\) and \([0,\,0.6,\,0.4]\), each with probability \(\frac12\). This blurs the distinction between deterministic variables (or often equivalently zero variance, zero entropy) and nondeterministic variables (positive variance, positive entropy).

In fact, it suggests that entropy can be calculated from the coefficients of an arbitrary convex combination. In the example just given, the entropy of \([0.2,\,0.6,\,0.2]\)) with respect to \([0.4,\,0.6,\,0]\) and \([0,\,0.6,\,0.4]\) is the entropy of \([0.5,\,0.5]\) with respect to \([1,\,0]\) and \([0,\,1]\), i.e. \(\frac12\log2+\frac12\log2=1\).

And why not try affine combinations! \[[0.2,\,0.6,\,0.2] = 2\cdot[\frac13,\,\frac13,\,\frac13]-1\cdot[\frac7{15},\,\frac1{15},\,\frac7{15}].\] The associated entropy in base 2 is \[-2\log2+1\log(-1) = -2+\pi i\log e = -2 + 4.53236i.\]

A proposition of trigonometry, or complex numbers

Prop. If \(0,1,z,w\) all lie on a circle centered at \(1/2\) then \(zw\) is the orthogonal projection of \(0\) onto the line through \(z\) and \(w\).

Proof. If \(z,w\) to lie on that circle, then \[(z-\frac12)(\overline{z}-\frac12)=\frac14,\text{ which simplifies to }z\overline{z}=\frac{z+\overline{z}}2\] \[(w-\frac12)(\overline{w}-\frac12)=\frac14,\text{ which simplifies to }w\overline{w}=\frac{w+\overline{w}}2.\] Introducing names \(z=a+bi,w=c+di\) these equations turn into \[a^2+b^2=a,\quad c^2+d^2=c.\] Two things need to be shown: \((1)\) \(zw\) lies on the line through \(z\) and \(w\), and \((2)\) \(zw\) is perpendicular to \(z-w\).

\((1)\) Verify that the following determinant is zero. \begin{align*} \begin{vmatrix}1&a&b\\1&c&d\\1&ac-bd&ad+bc\end{vmatrix}&=acd+bc^2-acd+bd^2-a^2d-abc+abc-b^2d+ad-bc\\ &=b(c^2+d^2)-(a^2+b^2)d+ad-bc=bc-ad+ad-bc=0. \end{align*} \((2)\) Verify that \(zw/(z-w)\) is imaginary, or equivalently that \(zw(\overline z-\overline w)\) is imaginary. \begin{align*} zw(\overline z-\overline w)+\overline{zw(\overline z-\overline w)}&=z\overline zw-zw\overline w+z\overline z\overline w-\overline zw\overline w=z\overline z(w+\overline w)-(z+\overline z)w\overline w\\ &=\frac{z+\overline z}2(w+\overline w)-(z+\overline z)\frac{w+\overline w}2=\frac12(z+\overline z)(w+\overline w)(1-1)=0. \end{align*}

Extremely rational cubic polynomials

By an extremely rational cubic polynomial, I mean a cubic polynomial \(p\) such that both \(p\) and \(p'\) are products of rational linear factors. A parametrization will now be derived.

Seek to satisfy \(p=(a+x)(b+x)(c+x),\ p'=3(d+x)(e+x)\) as polynomials.
\begin{align*}
0=p'(-e)&=ab+ac+bc-2(a+b+c)e+3e^2 \\
(3e-a-b-c)^2&=(a+b+c)^2-3(ab+ac+bc) \\
&=a^2+b^2+c^2-ab-ac-bc
\end{align*}We thus need to parametrize triples \((a,b,c)\) such that the right-hand side is a square.

By the cubic number system, I mean triples \((a,b,c)\) with the following operations.
\begin{align*}
(a,b,c)+(d,e,f)&\equiv(a+d,b+e,c+f) \\
(a,b,c)(d,e,f)&\equiv(ad+bf+ce,ae+bd+cf,af+be+cd) \\
Q(a,b,c)&\equiv a^2+b^2+c^2-ab-ac-bc
\end{align*}This system is a close sibling of the usual complex number system. Its relevance to us is in the last definition. It is the expression we want to make into a square. Verify the identity \(Q(zw)=Q(z)Q(w)\) to realize that our task can be fulfilled by setting \((a,b,c)\equiv(u,v,w)^2\). Normally, one would use fractions \(a,b,c,u,v,w\) in this number system, but we can just as well use complex parameters \(u,v,w\).

Here is the parametrization.
\begin{align*}
p&=(a+x)(b+x)(c+x),\quad p'=3(d+x)(e+x),\quad p''=6(f+x) \\
a&=u^2+2vw,\quad b=w^2+2uv,\quad c=v^2+2uw\\
d,e&=\frac13((u+v+w)^2\pm Q(u,v,w)),\quad f=\frac13(u+v+w)^2\end{align*}
Here is an example.
\begin{align*}
(u,v,w)&=(1,2,7) \\
(a,b,c)&=(1^2+2\times 2\times 7,\ 7^2+2\times 1\times 2,\ 2^2+2\times 1\times 7)=(29,53,18) \\
p&=(29+x)(53+x)(18+x)=27666+3013x+100x^2+x^3 \\
Q(u,v,w)&=1^2+2^2+7^2-1\times 2-1\times 7-2\times 7=31 \\
Q(a,b,c)&=29\times 29+53\times 53+18\times 18-29\times 53-29\times 18-53\times 18=961 \\
d,e&=\frac13({10}^2\pm 31)=\frac{131}3,23\\
p'&=3(\frac{131}3+x)(23+x)=(131+3x)(69+x)=3013+200x+3x^2\\
f&=\frac{100}3,\quad p''=6(\frac{100}3+x)=200+6x

\end{align*}
An example involving complex parameters \(u,v,w\) is best given visually and interactively.

Estimating a quarter turn in radians

Compute two sequences $\begin{matrix}a_0&\dots&a_{n-1}\\b_0&\dots&b_{n-1}\end{matrix}$ of unsigned integers with
$$\begin{align}\begin{matrix}a_0&a_{n-1}\\b_0&b_{n-1}\end{matrix} &= \begin{matrix}1&0\\0&1,\end{matrix}\\\left|\begin{matrix}a_k&a_{k+1}\\b_k&b_{k+1}\end{matrix} \right| &= 1\;\text{for all $k.$}\end{align}$$Then $$\sum\frac1{a_ka_{k+1}+b_kb_{k+1}} \approx \frac\tau4.$$

The idea is as follows: Integrate 1 from 0 to $\frac\tau4.$ A small term $\Delta\theta$ can be approximated as $\sin \Delta\theta$. Take $\sin^2 \Delta\theta$ to be the spread between two adjacent vectors $\begin{matrix}a_k\\b_k\end{matrix}$, $\begin{matrix}a_{k+1}\\b_{k+1}\end{matrix}.$ This spread is $$\frac{(a_kb_{k+1}-a_{k+1}b_k)^2}{(a_k^2+b_k^2)(a_{k+1}^2+b_{k+1}^2)} = \frac1{(a_ka_{k+1}+b_kb_{k+1})^2+1} \approx \frac1{(a_ka_{k+1}+b_kb_{k+1})^2}.$$Thus $\Delta\theta\approx\frac1{a_ka_{k+1}+b_kb_{k+1}}.$

Here's an algorithm:

#include <iostream>
#include <iomanip>
using namespace std;

long double integrate(unsigned long a0, unsigned long b0,
                      unsigned long a1, unsigned long b1) {   
    if (a0*a1+b0*b1 > 100000)
        return 1.0L/(a0*a1+b0*b1);
    else
        return integrate(a0+a1, b0+b1, a1, b1) + integrate(a0, b0, a0+a1, b0+b1);
}

int main() {
    cout << "approx " << setprecision(20) << integrate(1, 0, 0, 1) << endl;
}

The output is 1.5707963268188070548, which equals $\frac\tau4$ to 9 decimal places.

Continuous matrices

This will be a naive discussion about functions roughly of the form $A: \mathbb{R}^2\to\mathbb{R}$ but allowing delta behaviour which I only know of informally so far. $A$ is a continuous matrix. Addition $A+B$ is done componentwise and multiplication $AB$ is carried out by integrating
$$(AB)(x,y)=\int_{-\infty}^{+\infty}A(x,z)B(z,y)dz.$$

The plot operator $P$, defined through $P_f(x,y)=\delta(x-f(y))$, creates a continuous matrix. One interesting property, when combined with the usual associativity and distributivity of matrix operations, is that
$$ P_fP_g=\int_{-\infty}^{+\infty}\delta(x-f(z))\delta(z-g(y))dz=\delta(x-f(g(y)))=P_{f\circ g}.$$ This property makes it possible to interpret an arbitrary continuous matrix as a superposition of graphs $A = \sum \alpha_iP_{f_i}$ of single variable functions $f_i$, and to interpret multiplication of two matrices $A,B$ as a weaving together of all possible compositions: $AB = \sum\alpha_i\beta_jP_{f_i\circ g_j}.$ But all this is just conceptual speculation.


Here is a new light on fourier transforms: $F(x,y)=e^{-ixy}$. The inversion formula $FFf(x)=\tau f(-x)$ can be written as $F^2=\tau\delta(x+y)$ and verified using the definition of matrix multiplication.


In simple cases of signal analysis a linear operator $T$ transforming the signal is characterized through its impulse response $h=T\delta$ in the fashion $Tx(t)=\int_{-\infty}^{+\infty}h(t-u)x(u)du$, i.e. through the convolution $Tx=h*x$. The $h*$ part can be written as a matrix. Indeed, for any functions $x,y$,
$$ x*y(t) = \int_{-\infty}^{+\infty}x(t-u)y(u)du = C_xy(t) $$ where $C_f(x,y)=f(x-y)$ in general. In this case we have $Tx=h*x=C_hx$ so that $T=C_h$ in the same sense that a finite-dimensional linear transformation is equal to its transformation matrix.


This reminds me of the familiar way to define complex multiplication:
$$\left(\begin{matrix}a\\b\end{matrix}\right)\left(\begin{matrix}c\\d\end{matrix}\right)=\left(\begin{matrix}a&-b\\b&a\end{matrix}\right)\left(\begin{matrix}c\\d\end{matrix}\right).$$ Maybe it's more common to use matrices in both steads but we see that with this formulation $zw=M_zw$ where $M_z$ is chosen in the obvious way. The similarity between $C$ and $M$ is that both are obtained by repeating data across diagonals.

Another similarity: In the same way that matrix models verify the associativity $M_uM_vM_w$ of complex number multiplication $uvw$, they verify the associativity $C_fC_gC_h$ of convolutions $f*g*h$. And the same holds for the distributive property.