Multiplication via parabola

Posted on 29 June 2023 by John

Here’s a graphical way to multiply two positive numbers a and b using the parabola y = x².

Start at the origin, move a units to the left, then go up vertically to the parabola, and draw a point.
Go back to the origin, move b units to the right, go up vertically to the parabola, and draw another point.
Connect the points and see where they cross the y-axis. That point is ab.

Here’s an example multiplying 3 and 5.

Multiplying numbers geometrically via a parabola

Here’s why this works. The slope of the line is the change in y over the change in x which is

m = (b² − a²)/(b − (−a)) = b − a.

Use the equation of a line

y − y₀ = m(x – x₀)

with x₀ = b and y₀ = b² to get

y − b² = (b − a)(x − b).

Stick in x = 0 and you get y = ab.

Equation of an ellipse or hyperbola through five points

Posted on 19 June 2023 by John

Yesterday I wrote about the equation of a circle through three points. This post will be similar, looking at the equation of an ellipse or hyperbola through five points. As with the earlier post, we can write down an elegant equation right away. Making the equation practical will take more work.

The equation of a general conic section through five points

$\begin{vmatrix} x^2 & y^2 & xy & x & y & 1 \\ x_1^2 & y_1^2 & x_1 y_1 & x_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 y_2 & x_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 y_3 & x_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 y_4 & x_4 & y_4 & 1 \\ x_5^2 & y_5^2 & x_5 y_5 & x_5 & y_5 & 1 \\ \end{vmatrix} = 0$

Direct but inefficient solution

The most direct (though not the most efficient) way to turn this into an equation

is to expand the determinant by minors of the first row.

$A = \begin{vmatrix} y_1^2 & x_1 y_1 & x_1 & y_1 & 1 \\ y_2^2 & x_2 y_2 & x_2 & y_2 & 1 \\ y_3^2 & x_3 y_3 & x_3 & y_3 & 1 \\ y_4^2 & x_4 y_4 & x_4 & y_4 & 1 \\ y_5^2 & x_5 y_5 & x_5 & y_5 & 1 \\ \end{vmatrix} ,\, B = - \begin{vmatrix} x_1^2 & x_1 y_1 & x_1 & y_1 & 1 \\ x_2^2 & x_2 y_2 & x_2 & y_2 & 1 \\ x_3^2 & x_3 y_3 & x_3 & y_3 & 1 \\ x_4^2 & x_4 y_4 & x_4 & y_4 & 1 \\ x_5^2 & x_5 y_5 & x_5 & y_5 & 1 \\ \end{vmatrix}$

$C = \begin{vmatrix} x_1^2 & y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 & y_4 & 1 \\ x_5^2 & y_5^2 & x_5 & y_5 & 1 \\ \end{vmatrix} ,\, D = - \begin{vmatrix} x_1^2 & y_1^2 & x_1 y_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 y_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 y_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 y_4 & y_4 & 1 \\ x_5^2 & y_5^2 & x_5 y_5 & y_5 & 1 \\ \end{vmatrix}$

$E = \begin{vmatrix} x_1^2 & y_1^2 & x_1 y_1 & x_1 & 1 \\ x_2^2 & y_2^2 & x_2 y_2 & x_2 & 1 \\ x_3^2 & y_3^2 & x_3 y_3 & x_3 & 1 \\ x_4^2 & y_4^2 & x_4 y_4 & x_4 & 1 \\ x_5^2 & y_5^2 & x_5 y_5 & x_5 & 1 \\ \end{vmatrix} ,\, F = - \begin{vmatrix} x_1^2 & y_1^2 & x_1 y_1 & x_1 & y_1\\ x_2^2 & y_2^2 & x_2 y_2 & x_2 & y_2\\ x_3^2 & y_3^2 & x_3 y_3 & x_3 & y_3\\ x_4^2 & y_4^2 & x_4 y_4 & x_4 & y_4\\ x_5^2 & y_5^2 & x_5 y_5 & x_5 & y_5\\ \end{vmatrix}$

Python example

The following code will generate five points and set up the system of equations Ms = b whose solution s will give our coefficients.

    import numpy as np

    np.random.seed(20230619)

    M = np.zeros((5,5))
    b = np.zeros(5)
    for i in range(5):
        x = np.random.random()
        y = np.random.random()    
        M[i, 0] = y**2
        M[i, 1] = x*y
        M[i, 2] = x
        M[i, 3] = y
        M[i, 4] = 1
        b[i] = -x*x

Now we solve for the coefficients.

    s = np.linalg.solve(M, b)
    A, B, C, D, E, F = 1, s[0], s[1], s[2], s[3], s[4]

Next we verify that the points we generated indeed lie on the conic section whose equation we just found.

    for i in range(5):
        x = M[i, 2]
        y = M[i, 3]
        assert( abs(A*x*x + B*y*y + C*x*y + D*x + E*y + F) < 1e-10 )

Ellipse or hyperbola?

It turns out that the five points generated above lie on a hyperbola.

If we had set the random number generator seed to 20230620 we would have gotten an ellipse.

When we generate points at random, as we did above, we will almost certainly get an ellipse or a hyperbola. Since we’re generating pseudorandom numbers, there is some extremely small chance the generated points would lie on a line or on a parabola. In theory, this almost never happens, in the technical sense of “almost never.”

If 4AB − C² is positive we have a an ellipse. If it is negative we have a hyperbola. If it is zero, we have a parabola.

A parabola, a triangle, and a circle

Posted on 18 June 2023 by John

Let P be a parabola. Draw tangents to P at three different points. These three lines intersect to form a triangle.

Theorem: The circumcircle of this triangle passes through the focus of P. [1]

In the image above, the dashed lines are tangents and the black dot is the focus of the parabola.

(See this post for an explanation of what focus and directrix mean.)

By the way, creating the image above required finding the circle through the three intersection points. See the previous post for how to find such a circle.

[1] Ross Honsberger. Episodes in Nineteenth and Twentieth Century Euclidean Geometry. Mathematical Association of America. 1995. Page 47.

Equation of a circle through three points

Posted on 18 June 2023 by John

Three given points on a circle with unknown center

A few months ago I wrote about several equations that come up in geometric calculations that all have the form of a determinant equal to 0, where one column of the determinant is all 1s.

The equation of a circle through three points (x₁, y₁), (x₂, y₂), and (x₃, y₃) has this form:

$\begin{vmatrix} x^2 + y^2 & x & y & 1 \\ x_1^2 + y_1^2 & x_1 & y_1 & 1 \\ x_2^2 + y_2^2 & x_2 & y_2 & 1 \\ x_3^2 + y_3^2 & x_3 & y_3 & 1 \\ \end{vmatrix} = 0$

This is elegant equation, but it’s not in a form that’s ready to program. This post will start with the equation above and conclude with code.

Let M_ij be the determinant of the minor of the (i, j) element of the matrix above. That is, M_ij is the determinant of the 3 × 3 matrix obtained by deleting the ith row and jth column. Then we can expand the determinant above by minors of the last column to get

$M_{11}(x^2 + y^2) - M_{12}x + M_{13}y - M_{14} = 0$

which we can rewrite as

$x^2 + y^2 -\frac{M_{12}}{M_{11}} x + \frac{M_{13}}{M_{11}} y - \frac{M_{14}}{M_{11}} = 0.$

A circle with center (x₀, y₀) and radius r₀ has equation

By equating the coefficients x and y in the two previous equations we have

$\begin{align*} x_0 &= \phantom{-}\frac{1}{2} \frac{M_{12}}{M_{11}} \\ y_0 &= -\frac{1}{2} \frac{M_{13}}{M_{11}} \\ \end{align*}$

Next we expand the determinant minors in more explicit form.

$\begin{align*} M_{11} &= \begin{vmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \\ \end{vmatrix} \\ &= (x_2y_3 - x_3y_2) - (x_1y_3 - x_3y_1) + (x_1y_2 - x_2y_1) \\ &= (x1y_2 + x_2y_3 + x_3y_1) - (x_2y_1 + x_3y_2 + x_1y_3) \end{align*}$

Having gone to the effort of finding a nice expression for M₁₁ we can reuse it to find M₁₂ by the substitution

$x_i \to s_i \equiv x_i^2 + y_i^2$

as follows:

$\begin{align*} M_{12} &= \begin{vmatrix} x_1^2 + y_1^2 & y_1 & 1 \\ x_2^2 + y_2^2 & y_2 & 1 \\ x_3^2 + y_3^2 & y_3 & 1 \\ \end{vmatrix} = \begin{vmatrix} s_1^2 & y_1 & 1 \\ s_2^2 & y_2 & 1 \\ s_3^2 & y_3 & 1 \\ \end{vmatrix} \\ &= (s_1y_2 + s_2y_3 + s_3y_1) - (s_2y_1 + s_3y_2 + s_1y_3) \end{align*}$

And we can find M₁₃ from M₁₂ by reversing the role of the xs and ys.

$\begin{align*} M_{13} &= \begin{vmatrix} x_1^2 + y_1^2 & x_1 & 1 \\ x_2^2 + y_2^2 & x_2 & 1 \\ x_3^2 + y_3^2 & x_3 & 1 \\ \end{vmatrix} = \begin{vmatrix} s_1^2 & x_1 & 1 \\ s_2^2 & x_2 & 1 \\ s_3^2 & x_3 & 1 \\ \end{vmatrix} \\ &= (s_1x_2 + s_2x_3 + s_3x_1) - (s_2x_1 + s_3x_2 + s_1x_3) \end{align*}$

Now we have concrete expressions for M₁₁, M₁₂, and M₁₃, and these give us concrete expressions for x₀ and y₀.

To find the radius of the circle we simply find the distance from the center (x₀, y₀) to any of the points on the circle, say (x₁, y₁).

The following Python code incorporates the derivation above.

    def circle_thru_pts(x1, y1, x2, y2, x3, y3):
        s1 = x1**2 + y1**2
        s2 = x2**2 + y2**2
        s3 = x3**2 + y3**2
        M11 = x1*y2 + x2*y3 + x3*y1 - (x2*y1 + x3*y2 + x1*y3)
        M12 = s1*y2 + s2*y3 + s3*y1 - (s2*y1 + s3*y2 + s1*y3)
        M13 = s1*x2 + s2*x3 + s3*x1 - (s2*x1 + s3*x2 + s1*x3)
        x0 =  0.5*M12/M11
        y0 = -0.5*M13/M11
        r0 = ((x1 - x0)**2 + (y1 - y0)**2)**0.5
        return (x0, y0, r0)

The code doesn’t use any loops or arrays, so it ought to be easy to port to any programming language.

For some inputs the value of M₁₁ can be (nearly) zero, in which case the three points (nearly) lie on a line. In practice M₁₁ can be nonzero and yet zero for practical purposes and so that’s something to watch out for when using the code.

The determinant at the top of the post is zero when the points are colinear. There’s a way to make rigorous the notion that in this case the points line on a circle centered at infinity.

Update: Equation of an ellipse or hyperbola through five points

Similar triangles and complex numbers

Posted on 5 June 2023 by John

Suppose the vertices of two triangles are given by complex numbers a, b, c and x, y, z. The two triangles are similar if

$\left| \begin{matrix} 1 & 1 & 1 \\ a & b & c \\ x & y & z \\ \end{matrix} \right| = 0$

This can be found in G. H. Hardy’s classic A Course of Pure Mathematics. It’s on page 93 in the 10th edition.

Corollary

The theorem above generalizes a result from an earlier post giving a test for whether points determine an equilateral triangle. A triangle ABC is equilateral if it is similar to BCA. For triangles to be similar, corresponding vertex angles must be equal. If you can rotate the order of the angles in one triangle and the two triangles remain similar, the angles must all be equal.

Proof

Hardy’s proof of the equation above is direct and terse. I’ll quote it here in its entirety.

The condition required is that

$\overline{AB}\,/\,\overline{AC} = \overline{XY}\,/\,\overline{XZ}$

(large letters denoting the points whose arguments are the corresponding small letters), or

which is the same as the condition given.

Commentary

To expand on this a little bit, Hardy is saying that we need to show that if we take two corresponding vertices, A and X, the ratio of the lengths of corresponding sides is the same in both triangles. If memory serves, this is called the “SAS postulate.”

The determinant condition implies the last equation in Hardy’s proof. This is not obvious, but it’s easy to work out with pencil and paper. The reader may be left thinking “OK, I can see it’s true, but why would anyone think of that?” A more motivated derivation would require more background than the reader of the book is presumed to have.

You might object that Hardy’s last equation is a ratio of complex numbers, not a ratio of lengths. But you can take the absolute value of both sides. And the absolute value of a fraction is the absolute value of the numerator divided by the absolute value of the denominator. And now you’re taking a ratio of side lengths. In symbols,

$\frac{b-a}{c-a} = \frac{y-x}{z-x} \implies \frac{|b-a|}{|c-a|} = \frac{|y-x|}{|z-x|}$

Comparing approximations for ellipse perimeter

Posted on 28 May 2023 by John

This post will compare the accuracy of approximations for the perimeter of an ellipse.

The exact perimeter is given in terms of an elliptic integral. (That’s where the elliptic integrals gets their name.) And so an obvious way to approximate the perimeter would be to expand the elliptic integral in a power series. Unfortunately this series converges slowly for highly eccentric ellipses.

Ernst Kummer found a series that converges more quickly [1]. Let a and b be the semi-major and semi-minor axes of an ellipse and let

$x = \left(\frac{a-b}{a+b} \right )^2$

Then Kummer’s series is

$\pi(a+b)\left(1 + \frac{x}{4} + \frac{x^2}{64} + \frac{x^3}{256} + \frac{25x^4}{16384} + \cdots \right )$

The following plot compares the accuracy of truncating Kummer’s series after powers of x equal to 1 through 4.

The noise at the bottom of the chart comes from the limits of machine precision. For moderate values of eccentricity, three terms of Kummer’s series gives an approximation so accurate that all the error is coming from the limitations of floating point resolution.

I’ve written about approximations for the perimeter of an ellipse before. The approximation based on the 3/2 norm is about as accurate as Kummer’s series truncated after the x term. The error in Ramanujan’s approximation is between that of Kummer’s series with the cubic and quartic terms.

For the orbit of Pluto (e = 0.2488) the first-order Kummer approximation is good to 9 figures and the higher-order approximations are essentially good to machine precision.

For the orbit of Haley’s comet (e = 0.9671) the successive Kummer approximations are good to 2, 3, 4, and 5 figures.

[1] Carl E. Linderholm and Arthur C. Segal. An Overlooked Series for the Elliptic Perimeter. Mathematics Magazine, Vol. 68, No. 3 (Jun., 1995), pp. 216–220.

Contraharmonic mean

Posted on 20 May 2023 by John

I’ve mentioned the harmonic mean multiple times here, most recently last week. The harmonic mean pops up in many contexts.

The contraharmonic mean is a variation on the harmonic mean that comes up occasionally, though not as often as its better known sibling.

Definition

The contraharmonic mean of two positive numbers a and b is

$C = \frac{a^2 + b^2}{a + b}$

and more generally, the contraharmonic mean of a sequence of positive numbers is the sum of their squares over their sum.

Why the name?

What is “contra” about the contraharmonic mean? The harmonic mean and contraharmonic mean have a sort of reciprocal relationship. The ratio

$\frac{a - m}{m - b}$

equals a/b when m is the harmonic mean and it equals b/a when m is the contraharmonic mean.

Geometric visualization

Construct a trapezoid with two vertical sides, one of length a and another of length b. The vertical slice through the intersection of the diagonals, the dashed blue line in the diagram below, has length equal to the harmonic mean of a and b.

The vertical slice midway between the two vertical sides, the solid red line in the diagram, has length equal to the arithmetic mean of a and b.

The vertical slice on the opposite side of the red line, as from the red line on one side as the blue line is on the other, the dash-dot green line above, has length equal to the contraharmonic mean of a and b.

Connection to statistics

The contraharmonic mean of a list of numbers is the ratio of the sum of the squares to the sum. If we divide the numerator and denominator by the number of terms n we see this is the ratio of the second moment to the first moment, i.e. the mean of the squared values divided the mean of the values.

This smells like the ratio of variance to mean, known as the index of dispersion D, except variance is the second central moment rather than the second moment. The second moment equals variance plus the square of the mean, and so

$D = \frac{\sigma^2}{\mu}$

and

$C = \frac{\text{E}(X^2)}{\text{E}(X)} = \frac{\sigma^2 + \mu^2}{\mu} = D + \mu.$

How Albrecht Dürer drew an 11-sided figure

Posted on 14 April 2023 by John

You cannot exactly construct an 11-sided regular polygon (called a hendecagon or an undecagon) using only a straight edge and compass. Gauss fully classified which regular n-gons can be constructed, and this isn’t one of them [1].

However, Albrecht Dürer [2] came up with a good approximate construction for a hendecagon.

To construct an eleven-sided figure by means of a compass, I take a quarter of a circle’s diameter, extend it by one eighth of its length, and use this for the construction of the eleven-sided figure.

It took some trial and error for me to understand what Dürer meant.

In more explicit terms, draw a circle of radius R centered at the origin. Then draw a circle of radius 9R/16 centered at (R, 0). Then draw a line from the origin to the intersection of the two circles. (There are two intersections, and it doesn’t matter which one you pick.) Then this line makes an angle of approximately (360/11) degrees with the x-axis.

Without loss of generality we can assume R = 1. A point at the intersection of both circles satisfies the equation of each circle, and from these two equations in two unknowns we can determine that the x coordinate of the two intersection points is 431/512. We can then use the Pythagorean theorem to solve for y and find the inverse tangent of y/x. This shows that we’ve constructed an angle of 32.6696°. We were after an angle of 32.7272°, and so the relative error in Dürer’s approximation is less than 0.2%.

Diagram of Durer's construction

To be even more explicit, here is the Python code that created the image above.

    import matplotlib.pyplot as plt
    from numpy import linspace, sin, cos, pi

    r = 9/16
    t = linspace(0, 2*pi)
    plt.plot(cos(t), sin(t), 'r-')
    plt.plot(r*cos(t) + 1, r*sin(t), 'b-')

    x = 431/512
    y = (1 - x*x)**0.5
    plt.plot([0, 1], [0, 0], 'k')
    plt.plot([0, x], [0, y], 'k')

    plt.gca().set_aspect("equal")
    plt.axis("off")
    plt.savefig("durer11.png")

More Dürer posts

[1] The Guass-Wantzel theorem says a regular n-gon can be constructed with straight edge and compass if n factors into a power of 2 and distinct Fermat primes.

[2] I found this in A Panoply of Polygons which in turn cites “The Polygons of Albrecht Durer -1525” By Gordon Hughes, available on arXiv.

A pentagon, hexagon, and decagon walk into a bar …

Posted on 14 April 2023 by John

The new book A Panoply of Polygons cites a theorem Euclid (Proposition XIII.10) saying

If a regular pentagon, a regular hexagon, and a regular decagon are inscribed in congruent circles, then their side lengths form a right triangle.

This isn’t exactly what Euclid said, but it’s an easy deduction from what he did say.

Here’s an illustration.

Illustrating Euclid Proposition XIII.10

Update: Ralph Corderoy suggested in the comments that it would be good to add the congruent circles through the polygon vertices. Here is a graphic with the circles added.

Illustrating Euclid Proposition XIII.10

Connection with golden ratio

I learned about this theorem from the book cited above, which arrived in the mail yesterday. The figure looked familiar because Cliff Pickover had posted essentially the same illustration on Twitter recently, saying that not only is there a right triangle in the middle of the diagram, this triangle is half of a golden rectangle.

The length of the side of a regular polygon with n sides is 2R sin(π/n) where R is the radius of the circle through the vertices. So to prove that the ratio between the side of the hexagon and the side of the decagon is the golden ratio φ, it suffices to prove

sin(π/6) / sin(π/10) = φ

You could verify this in Mathematica by seeing that the following returns True.

    Sin[Pi/6]/ Sin[Pi/10] == GoldenRatio

Python code

Here’s the code I used to create the image above.

    from numpy import exp, pi, arange
    import matplotlib.pyplot as plt
    
    # Vertices of a hexagon
    hexa = exp(2j*pi*arange(6)/6)
    
    # Vertices of a decagon, shifted and rotated
    # to have a side perpendicular to the hexagon
    deca = exp(2j*pi*(0.5+arange(10))/10)
    deca += hexa[1] - deca[5]
    
    # Shift and rotate a pentagon to share one vertex 
    # with the hexagon and one vertex with the decagon
    omega = exp(2j*pi/5)
    c = (hexa[2]*omega - deca[4])/(omega - 1)
    penta = c + (hexa[2] - c)*exp(2j*pi*arange(5)/5) 
    
    def connect(p, q, c):
        plt.plot([p.real, q.real], [p.imag, q.imag], '-', color=c)
    
    plt.gca().set_aspect("equal")
    
    for n in range(-1, 5):
        connect(hexa[n], hexa[n+1], 'blue')
    for n in range(-1, 9):
        connect(deca[n], deca[n+1], 'green')
    for n in range(-1, 4):
        connect(penta[n], penta[n+1], 'red')
    
    plt.show()

The ability to use −1 to as the index to the last element of an array simplifies the code that connects the vertices. In a language like C you’d have to write an extra statement after your loop to join the last vertex to the first.

When I first saw the illustration in Panoply of Pentagons I though that the top of the pentagon was parallel to the top of the hexagon. When I wrote the Python code first assuming this orientation of the pentagon, things didn’t line up. Writing the code made me pay closer attention, as is often the case.

A new trig identity

Posted on 13 April 2023 by John

This evening I ran across a trig identity I hadn’t seen before. I doubt it’s new to the world, but it’s new to me.

Let A, B, and C be the angles of an arbitrary triangle. Then

sin² A + sin² B + sin² C = 2 + 2 cos A cos B cos C.

This looks a little like the Pythagorean theorem, but the Pythagorean theorem involves the sides of a triangle, not the angles. (I expect there’s an interesting generalization of the identity above to spherical geometry where sides are angles.)

This identity also looks a little like the law of cosines, but the law of cosines mixes sides and angles, and this identity only involves angles.

Source: Lemma 4.4.3 in A Panoply of Polygons by Claudi Alsina and Roger B. Nelsen.

Geometry

Multiplication via parabola

Equation of an ellipse or hyperbola through five points

Direct but inefficient solution

More efficient solution

Python example

Ellipse or hyperbola?

Related posts

A parabola, a triangle, and a circle

Equation of a circle through three points

Similar triangles and complex numbers

Corollary

Proof

Commentary

Related posts

Comparing approximations for ellipse perimeter

Contraharmonic mean

Definition

Why the name?

Geometric visualization

Connection to statistics

How Albrecht Dürer drew an 11-sided figure

More Dürer posts

A pentagon, hexagon, and decagon walk into a bar …

Connection with golden ratio

Python code

Related posts

A new trig identity