diff --git a/_freeze/posts/stereo/2/index/execute-results/html.json b/_freeze/posts/stereo/2/index/execute-results/html.json new file mode 100644 index 0000000..c657e0f --- /dev/null +++ b/_freeze/posts/stereo/2/index/execute-results/html.json @@ -0,0 +1,12 @@ +{ + "hash": "e723e799442b3e52f9543809a686e1f9", + "result": { + "engine": "jupyter", + "markdown": "---\ntitle: \"Further Notes on Algebraic Stereography\"\ndescription: |\n How do you rotate in 2D and 3D without standard trigonometry?\nformat:\n html:\n html-math-method: katex\njupyter: python3\ndate: \"2021-10-10\"\ndate-modified: \"2025-06-30\"\ncategories:\n - algebra\n - complex analysis\n - polar roses\n - generating functions\n---\n\n\n\n\n\nIn my previous post, I discussed the stereographic projection of a circle as it pertains\n to complex numbers, as well as its applications in 2D and 3D rotation.\nIn an effort to document more interesting facts about this mathematical object\n (of which scarce information is immediately available online),\n I will now elaborate on more of its properties.\n\n\nChebyshev Polynomials\n---------------------\n\n[Previously](/posts/chebyshev/1), I derived the\n [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials)\n with the archetypal complex exponential.\nThese polynomials express the sines and cosines of a multiple of an angle from\n the sine and cosine of the base angle.\nWhere $T_n(t)$ are Chebyshev polynomials of the first kind and $U_n(t)$ are those of the second kind,\n\n$$\n\\begin{gather*}\n \\cos(n \\theta) = T_n(\\cos(\\theta))\n \\\\\n \\sin(n \\theta) = U_{n - 1}(\\cos(\\theta)) \\sin(\\theta)\n\\end{gather*}\n$$\n\nThe complex exponential derivation begins by squaring and developing a second-order recurrence.\n\n$$\n\\begin{align*}\n (e^{i\\theta})^2 &= (\\cos + i\\sin)^2\n \\\\\n &= \\cos^2 + 2i\\cos \\cdot \\sin - \\sin^2 + (0 = \\cos^2 + \\sin^2 - 1)\n \\\\\n &= 2\\cos^2 + 2i\\cos \\cdot \\sin - 1\n \\\\\n &= 2\\cos \\cdot (\\cos + i\\sin) - 1\n \\\\\n &= 2\\cos(\\theta)e^{i\\theta} - 1\n \\\\\n (e^{i\\theta})^{n+2} &= 2\\cos(\\theta)(e^{i\\theta})^{n+1} - (e^{i\\theta})^n\n\\end{align*}\n$$\n\nThis recurrence relation can then be used to obtain the Chebyshev polynomials, and hence,\n the expressions using sine and cosine above.\nPresented this way with such a simple derivation, it appears as though these relationships\n are inherently trigonometric.\n\nHowever, these polynomials actually have *nothing* to do with sine and cosine on their own.\nFor one, [they appear in graph theory](/posts/chebyshev/2), and for two,\n the importance of the complex exponential is overstated.\n$e^{i\\theta}$ really just specifies a point on the complex unit circle.\nThis property is used on the second line to coax the equation into a quadratic in $e^{i\\theta}$.\nThis is also the *only* property upon which the recurrence depends; all else is algebraic manipulation.\n\n\n### Back to the Stereograph\n\nKnowing this, let's start over with the stereographic projection of the circle:\n\n$$\no_1(t) = {1 + it \\over 1 - it}\n = {1 - t^2 \\over 1 + t^2} + i {2t \\over 1 + t^2}\n = \\text{c}_1 + i\\text{s}_1\n$$\n\nThe subscript \"1\" is because as *t* ranges over $(-\\infty, \\infty)$, the function loops once\n around the unit circle.\nTaking this to higher powers keeps points on the circle since all points on the circle\n have a norm of 1.\nIt also makes more loops around the circle, which we can denote by larger subscripts:\n\n$$\n\\begin{align*}\n o_n &= (o_1)^n\n = \\left( {1 + it \\over 1 - it} \\right)^n\n \\\\\n \\text{c}_n + i\\text{s}_n\n &= (\\text{c}_1 + i\\text{s}_1)^n\n\\end{align*}\n$$\n\nThis mirrors raising the complex exponential to a power\n (which loops over the range $(-\\pi, \\pi)$ instead).\nThe final line is analogous to de Moivre's formula, but in a form where everything is\n a ratio of polynomials in *t*.\nThis means that the Chebyshev polynomials can be obtained directly from these rational expressions:\n\n$$\n\\begin{align*}\n o_2 = (o_1)^2 &= (\\text{c}_1 + i\\text{s}_1)^2\n \\\\\n &= \\text{c}_1^2 + 2i\\text{c}_1\\text{s}_1 - \\text{s}_1^2\n + (0 = \\text{c}_1^2 + \\text{s}_1^2 - 1)\n \\\\\n &= 2\\text{c}_1^2 + 2i\\text{c}_1\\text{s}_1 - 1\n \\\\\n &= 2\\text{c}_1(\\text{c}_1 + i\\text{s}_1) - 1\n \\\\\n &= 2\\text{c}_1 o_1 - 1\n \\\\\n o_2 \\cdot (o_1)^n &= 2\\text{c}_1 o_1 \\cdot (o_1)^n - (o_1)^n\n \\\\\n o_{n+2} &= 2\\text{c}_1 o_{n+1} - o_n\n\\end{align*}\n$$\n\nThis matches the earlier recurrence relation with the complex exponential and therefore\n the recurrence relation of the Chebyshev polynomials.\nIt also means that the the rational functions obey the same relationship as sine and cosine:\n\n$$\n\\begin{matrix}\n \\begin{gather*}\n \\text{c}_n = T_n(\\text{c}_1)\n \\\\\n \\text{s}_n = U_{n-1}(\\text{c}_1) \\text{s}_1\n \\end{gather*}\n & \\text{where }\n \\text{c}_1 = {1 - t^2 \\over 1 + t^2}, &\n \\text{s}_1 = {2t \\over 1 + t^2}\n\\end{matrix}\n$$\n\nThus, the Chebyshev polynomials are tied to (coordinates on) circles,\n rather than explicitly to the trigonometric functions.\nIt is a bit strange that these polynomials are in terms of rational functions, but no stranger\n than them being in terms of *ir*rational functions like sine and cosine.\n\n\nCalculus\n--------\n\nSince these functions behave similarly to sine and cosine, one might wonder about\n the nature of these expressions in the context of calculus.\n\nFor comparison, the complex exponential (as it is a parallel construction) has a simple derivative[^1].\nSince the exponential function is its own derivative, the expression acquires\n an imaginary coefficient through the chain rule.\n\n[^1]: This is forgoing the fact that complex derivatives require more care than their real counterparts.\n It matters slightly less in this case since this function is complex-valued, but has a real parameter.\n\n$$\n\\begin{align*}\n e^{it} &= \\cos(t) + i\\sin(t)\n \\\\\n {d \\over dt} e^{it}\n &= {d \\over dt} \\cos(t) + {d \\over dt} i\\sin(t)\n \\\\\n i e^{it} &= -\\sin(t) + i\\cos(t)\n \\\\\n i[\\cos(t) + i\\sin(t)]\n &\\stackrel{\\checkmark}{=} -\\sin(t) + i\\cos(t)\n\\end{align*}\n$$\n\nMeanwhile, the complex stereograph has derivative\n\n$$\n\\begin{align*}\n {d \\over dt} o_1(t) &= {d \\over dt} {1 + it \\over 1 - it}\n = {i(1 - it) + i(1 + it) \\over (1 - it)^2}\n \\\\\n &= {2i \\over (1 - it)^2}\n = {2i(1 + it)^2 \\over (1 + t^2)^2}\n = {2i(1 - t^2 + 2it) \\over (1 + t^2)^2}\n \\\\\n &= {-4t \\over (1 + t^2)^2} + i {2(1 - t^2) \\over (1 + t^2)^2}\n \\\\\n &= {-2 \\over 1 + t^2}s_1 + i {2 \\over 1 + t^2}c_1\n \\\\\n &= -(1 + c_1)s_1 + i(1 + c_1)c_1\n \\\\\n &= i(1 + c_1)o_1\n\\end{align*}\n$$\n\nJust like the complex exponential, an imaginary coefficient falls out.\nHowever, the expression also accrues a $1 + c_1$ term, almost like an adjustment factor\n for its failure to be the complex exponential.\nSine and cosine obey a simpler relationship with respect to the derivative,\n and thus need no adjustment.\n\n\n### Complex Analysis\n\nSince $o_n$ is a curve which loops around the unit circle *n* times, that possibly suits it\n to showing a simple result from complex analysis.\nIntegrating along a contour which wraps around a sufficiently nice function's pole\n (i.e., where its magnitude grows without bound) yields a familiar value.\nThis is easiest to see with $f(z) = 1 / z$:\n\n$$\n\\oint_\\Gamma {1 \\over z} dz\n = \\int_a^b {\\gamma'(t) \\over \\gamma(t)} dt\n = 2\\pi i\n$$\n\nIn this example, Γ is a counterclockwise curve parametrized by γ which loops once around\n the pole at *z* = 0.\nMore loops will scale this by a factor according to the number of loops.\n\nNormally this equality is demonstrated with the complex exponential, but will $o_1$ work just as well?\nIf Γ is the unit circle, the integral is:\n\n$$\n\\oint_\\Gamma {1 \\over z} dz\n = \\int_{-\\infty}^\\infty {o_1'(t) \\over o_1(t)} dt\n = \\int_{-\\infty}^\\infty i(1 + c_1(t)) dt\n = 2i\\int_{-\\infty}^\\infty {1 \\over 1 + t^2} dt\n$$\n\nIf one has studied their integral identities, the indefinite version of the final integral\n will be obvious as $\\arctan(t)$, which has horizontal asymptotes of $\\pi / 2$ and $-\\pi / 2$.\nTherefore, the value of the integral is indeed $2\\pi i$.\n\nIf there are *n* loops, then naturally there are *n* of these $2\\pi i$s.\nSince powers of *o* are more loops around the circle, the chain and power rules show:\n\n$$\n\\begin{gather*}\n {d \\over dt} (o_1)^n = n(o_1)^{n-1} {d \\over dt} o_1\n \\\\[14pt]\n \\oint_\\Gamma {1 \\over z} dz\n = \\int_{-\\infty}^\\infty {n o_1(t)^{n-1} o_1'(t) \\over o_1(t)^n} dt\n = n \\int_{-\\infty}^\\infty {o_1'(t) \\over o_1(t)} dt\n = 2 \\pi i n\n\\end{gather*}\n$$\n\nIt is certainly possible to perform these contour integrals along straight lines;\n in fact, integrating along lines from 1 to *i* to -1 to -*i* deals with a\n similar integral involving arctangent.\nHowever, the best one can do to construct more loops with lines is to count each line\n multiple times, which isn't extraordinarily convincing.\n\nPerhaps the use of $\\infty$ in the integral bounds is also unconvincing.\nThe integral can be shifted back into the realm of plausibility by considering simpler bounds on $o_2$:\n\n$$\n\\begin{align*}\n \\oint_\\Gamma {1 \\over z} dz\n &= \\int_{-1}^1 {2 o_1(t) o_1'(t) \\over o_1(t)^2} dt\n \\\\\n &= 2 \\int_{-1}^1 {o_1'(t) \\over o_1(t)} dt\n \\\\\n &= 2(2i\\arctan(1) - 2i\\arctan(-1))\n \\\\\n &= 2\\pi i\n\\end{align*}\n$$\n\nThis has an additional benefit: using the series form of $1 / (1 + t^2)$ and integrating,\n one obtains the series form of the arctangent.\nThis series converges for $-1 \\le t \\le 1$, which happens to match the bounds of integration.\nThe convergence of this series is fairly important, since it is tied to formulas for π,\n in particular [Leibniz's formula](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80).\n\nWere one to integrate with the complex exponential, we would instead use the bounds $(0, 2\\pi)$,\n since at this point a full loop has been made.\nBut think to yourself -- how do you know the period of the complex exponential?\nHow do you know that 2π radians is equivalent to 0 radians?\nThe result using stereography relies on neither of these prior results and is directly pinned\n to a formula for π instead an apparent detour through the number *e*.\n\n\nPolar Curves\n------------\n\nPolar coordinates are useful for expressing for which the distance from the origin is\n a function of the angle with respect to the positive *x*-axis.\nThey can also be readily converted to parametric forms:\n\n$$\n\\begin{gather*}\n r(\\theta) &\\Longleftrightarrow&\n \\begin{matrix}\n x(\\theta) = r \\cos(\\theta) \\\\\n y(\\theta) = r \\sin(\\theta)\n \\end{matrix}\n\\end{gather*}\n$$\n\nPolar curves frequently loop in on themselves, and so it is necessary to choose appropriate bounds\n for θ (usually as multiples of π) when plotting.\nEvidently, this is due to the use of sine and cosine in the above parametrization.\nFortunately, $s_n$ and $c_n$ (as shown by the calculus above) have much simpler bounds.\nSo what happens when one substitutes the rational functions in place of the trig ones?\n\n\n### Polar Roses\n\n[Polar roses](https://en.wikipedia.org/wiki/Rose_(mathematics)) are beautiful shapes which have\n a simple form when expressed in polar coordinates.\n\n$$\nr(\\theta) = \\cos \\left( {p \\over q} \\cdot \\theta \\right)\n$$\n\nThe ratio $p/q$ in least terms uniquely determines the shape of the curve.\n\nIf you weren't reading this post, you might assume this curve is transcendental since it uses cosine,\n but you probably know better at this point.\nThe Chebyshev examples above demonstrate the resemblance between $c_n$ and $\\cos(n\\theta)$.\nThe subscript of $c$ is easiest to work with as an integer, so let $q = 1$.\n\n$$\nx(t) = c_p(t) c_1(t) \\qquad y(t) = c_p(t) s_1(t)\n$$\n\nwill plot a $p/1$ polar rose as t ranges over $(-\\infty, \\infty)$.\n\n\n\n::: {#fig-polar-roses-1}\n{{< video \"./polar_roses_1.mp4\" >}}\n\np/1 polar roses as rational curves.\nSince *t* never reaches infinity, a bite appears to be taken out of the graphs near (-1, 0).\"\n:::\n\n$q = 1$ happens to match the subscript *c* term of *x* and *s* term of *y*, so one might wonder\n whether the other polar curves can be obtained by allowing it to vary as well.\nAnd you'd be right.\n\n$$\nx(t) = c_p(t) c_q(t) \\qquad y(t) = c_p(t) s_q(t)\n$$\n\nwill plot a $p/q$ polar rose as t ranges over $(-\\infty, \\infty)$.\n\n\n\n::: {#fig-polar-roses-2}\n{{< video \"./polar_roses_2.mp4\" >}}\n\np/q polar roses as rational curves\n:::\n\nJust as with the prior calculus examples, doubling all subscripts of *c* and *s* will\n only require *t* to range over $(-1, 1)$, which removes the ugly bite mark.\nPerhaps it is also slightly less satisfying, since the fraction $p/q$ directly appears in the\n typical polar incarnation with cosine.\nOn the other hand, it exposes an important property of these curves: they are all rational.\n\nThis approach lends additional precision to a prospective pseudo-polar coordinate system.\nIn the next few examples, I will be using the following notation for compactness:\n\n$$\n \\begin{gather*}\n R_n(t) = f(t) &\\Longleftrightarrow&\n \\begin{matrix}\n x(t) = f(t) c_n(t) \\\\\n y(t) = f(t) s_n(t)\n \\end{matrix}\n\\end{gather*}\n$$\n\n\n### Conic Sections\n\nThe polar equation for a conic section (with a particular unit length occurring somewhere)\n in terms of its eccentricity $\\varepsilon$ is:\n\n$$\nr(\\theta) = {1 \\over 1 - \\varepsilon \\cos(\\theta)}\n$$\n\nCorrespondingly, the rational polar form can be expressed as\n\n$$\nR_1(t) = {1 \\over 1 - \\varepsilon c_1}\n$$\n\nSince polynomial arithmetic is easier to work with than trigonometric identities,\n it is a matter of pencil-and-paper algebra to recover the implicit form from a parametric one.\n\n\n#### Parabola ($|\\varepsilon| = 1$)\n\nThe conic section with the simplest implicit equation is the parabola.\nSince $c_n$ is a simple ratio of polynomials in *t*, it is much simpler to recover the implicit equation.\nFor $\\varepsilon = 1$,\n\n:::: {layout-ncol=\"2\"}\n\n::: {#39eea2a2 .cell execution_count=5}\n\n::: {.cell-output .cell-output-display}\n![$x = {c_1 \\over 1 - c_1} \\quad y = {s_1 \\over 1 - c_1}$](index_files/figure-html/cell-5-output-1.png){}\n:::\n:::\n\n\n::: {}\n$$\n\\begin{align*}\n 1 - c_1 &= 1 - {1 - t^2 \\over 1 + t^2}\n = {2 t^2 \\over 1 + t^2}\n \\\\\n y &= {s_1 \\over 1 - c_1}\n = {2t \\over 1 + t^2} {1 + t^2 \\over 2 t^2}\n = {1 \\over t}\n \\\\\n x &= {c_1 \\over 1 - c_1}\n = {1 - t^2 \\over 1 + t^2} \\cdot {1 + t^2 \\over 2 t^2}\n = {1 - t^2 \\over 2t^2}\n \\\\\n &= {1 \\over 2t^2} - {1 \\over 2}\n = {y^2 \\over 2} - {1 \\over 2}\n\\end{align*}\n$$\n:::\n::::\n\n*x* is a quadratic polynomial in *y*, so trivially the figure formed is a parabola.\nTechnically it is missing the point where $y = 0 ~ (t = \\infty)$, and this is not a circumstance\n where using a higher $c_n$ would help.\nIt is however, similar to the situation where we allow $o_1(\\infty) = -1$, and an argument\n can be made to waive away any concerns one might have.\n\n\n#### Ellipse ($|\\varepsilon| < 1$)\n\nEllipses are next.\nThe simplest fraction between zero and one is 1/2, so for $\\varepsilon = 1/2$,\n\n:::: {layout-ncol = \"2\"}\n\n::: {#f9484e49 .cell execution_count=6}\n\n::: {.cell-output .cell-output-display}\n![$x = {c_1 \\over 1 - c_1 / 2} \\quad y = {s_1 \\over 1 - c_1 / 2}$](index_files/figure-html/cell-6-output-1.png){}\n:::\n:::\n\n\n::: {}\n$$\n\\begin{align*}\n 1 - {1 \\over 2}c_1 &= 1 - {1 \\over 2} \\cdot {1 - t^2 \\over 1 + t^2}\n = {3 t^2 + 1 \\over 2 + 2t^2}\n \\\\\n y &= {s_1 \\over 1 - {1 \\over 2}c_1}\n = {4t \\over 3t^2 + 1}\n \\\\\n x &= {c_1 \\over 1 - {1 \\over 2}c_1}\n = {2 - 2t^2 \\over 3t^2 + 1}\n\\end{align*}\n$$\n:::\n::::\n\nThere isn't an obvious way to combine products of *x* and *y* into a single equation.\nThe general form of a conic section is $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$, so\n we know that the implicit equation for the curve almost certainly involves $x^2$ and $y^2$.\n\n$$\nx^2 = {4 - 8t^2 + 4t^4 \\over (3t^2 + 1)^2} \\qquad\n y^2 = {16t^2 \\over (3t^2 + 1)^2}\n$$\n\nSquaring produces some $t^4$ terms which cannot exist outside of these terms and *xy*.\nA linear combination of $x^2$ and $y^2$ never includes any cubic terms in the numerator\n which would appear in *xy*, so $B = 0$.\nSince all remaining terms are linear in *x* and *y*, their denominator must appear as a factor\n in the numerator of $Ax^2 + Cy^2$, whatever *A* and *C* are.\n\nSince the coefficient of $t^4$ in $x^2$ is 4, *A* must be multiple of 3.\nThrough trial and error, $A = 3, C = 4$ gives:\n\n$$\n\\begin{align*}\n 3x^2 + 4y^2\n &= {12 - 24t^2 + 12t^4 + 64t^2 \\over (3t^2 + 1)^2}\n \\\\\n &= {12 - 40t^2 + 12t^4 \\over (3t^2 + 1)^2}\n \\\\\n &= {(4t^2 + 12) (3t^2 + 1) \\over (3t^2 + 1)^2}\n \\\\\n &= {4t^2 + 12 \\over 3t^2 + 1}\n\\end{align*}\n$$\n\nSince the numerator of *y* has a *t*, this is clearly some combination of *x* and a constant.\nBy the previous line of thought, the constant term must be a multiple of 4, and picking the smallest\n option finally results in the implicit form:\n\n$$\n\\begin{align*}\n 4 &= {4(3t^2 + 1) \\over 3t^2 + 1}\n = {12t^2 + 4 \\over 3t^2 + 1}\n \\\\\n {4t^2 + 12 \\over 3t^2 + 1} - 4\n &= {8 - 8t^2 \\over 3t^2 + 1} = 4x\n \\\\[14pt]\n 3x^2 + 4y^2 &= 4x + 4\n \\\\\n 3x^2 + 4y^2 - 4x - 4 &= 0\n\\end{align*}\n$$\n\nNotably, the coefficients of *x* and *y* are 3 and 4.\nSimultaneously, $o_1(\\varepsilon) = o_1(1/2) = {3 \\over 5} + i{4 \\over 5}$.\nThis binds together three concepts: the simplest case of the Pythagorean theorem,\n the 3-4-5 right triangle; the coefficients of the implicit form; and the role of eccentricity\n with respect to stereography.\n\n\n#### Hyperbola ($|\\varepsilon| > 1$)\n\nAs evidenced by the bound on the eccentricity above, hyperbolae are in some way the inverses of ellipses.\nSince $o_1(2)$ is a reflection of $o_1(1/2)$, you might think the implicit equation for\n $\\varepsilon = 2$ to be the same, but with a flipped sign or two.\nUnfortunately, you'd be wrong.\n\n:::: {layout-ncol=\"2\"}\n\n::: {#226c9935 .cell execution_count=7}\n\n::: {.cell-output .cell-output-display}\n![$x = {c_1 \\over 1 - 2c_1} \\quad y = {s_1 \\over 1 - 2c_1}$](index_files/figure-html/cell-7-output-1.png){}\n:::\n:::\n\n\n::: {}\n$$\n\\begin{gather*}\n \\begin{align*}\n x &= {c_1 \\over 1 - 2c_1} = {1 - t^2 \\over 3t^2 - 1}\n \\\\\n y &= {s_1 \\over 1 - 2c_1} = {2t \\over 3t^2 - 1}\n \\\\[14pt]\n 3x^2 - y^2 &= {3 - 6t^2 + 3t^4 - 4t^2 \\over 3t^2 - 1}\n \\\\\n &= {(t^2 - 3)(3t^2 - 1) \\over (3t^2 - 1)^2 }\n \\\\\n &= {t^2 - 3 \\over 3t^2 - 1 }\n = ... = -4x - 1\n \\end{align*}\n \\\\[14pt]\n 3x^2 - y^2 + 4x + 1 = 0\n\\end{gather*}\n$$\n:::\n::::\n\nAt the very least, the occurrences of 1 in the place of 4 have a simple explanation: 1 = 4 - 3.\n\n\n### Archimedean Spiral\n\nArguably the simplest (non-circular) polar curve is $r(\\theta) = \\theta$, the unit\n [Archimedean spiral](https://en.wikipedia.org/wiki/Archimedean_spiral).\nSince the curve is defined by a constant turning, this is a natural application of the properties\n of sine and cosine.\nThe closest equivalent in rational polar coordinates is $R_1(t) = t$.\nBut this can be converted to an implicit form:\n\n$$\n\\begin{gather*}\n x = tc_1 \\qquad y = ts_1\n \\\\[14pt]\n x^2 + y^2 = t^2(c_1^2 + s_1^2) = t^2\n \\\\\n y = {2t^2 \\over 1 + t^2} = {2(x^2 + y^2) \\over 1 + (x^2 + y^2)}\n \\\\[14pt]\n (1 + x^2 + y^2)y = 2(x^2 + y^2)\n\\end{gather*}\n$$\n\nThe curve produced by this equation is a\n [right strophoid](https://mathworld.wolfram.com/RightStrophoid.html)\n with a node at (0, 1) and asymptote $y = 2$.\nThis form suggests something interesting about this curve: it approximates the Archimedean spiral\n (specifically the one with polar equation $r(\\theta) = \\theta/2$).\nIndeed, the sequence of curves with parametrization $R_n(t) = 2nt$ approximate the (unit) spiral\n for larger *n*, as can be seen in the following video.\n\n\n\n::: {#fig-approx-archimedes}\n{{< video ./approximate_archimedes.mp4 >}}\n\nApproximations to the Archimedean spiral\n:::\n\n\nSince R necessarily defines a rational curve, the curves will never be equal,\n just as any stretching of $c_n$ will never exactly become cosine.\n\n\nClosing\n-------\n\nSine, cosine, and the exponential function, are useful in a calculus setting precisely\n because of their constant \"velocity\" around the circle.\nAlso, nearly every modern scientific calculator in the world features buttons\n for trigonometric functions, so there seems to be no reason *not* to use them.\n\nWe can however be misled by their apparent omnipresence.\nStereographic projection has been around for *millennia*, and not every formula needs to be rewritten\n in its language.\nFor example (and as previously mentioned), defining the Chebyshev polynomials really only requires\n understanding the multiplication of two complex numbers whose norm cannot grow,\n not trigonometry and dividing angles.\nMany other instances of sine and cosine merely rely on a number (or ratio) of loops around a circle.\nWhen velocity does not factor, it will obviously do to \"stay rational\".\n\nOne of my favorite things to plot as a kid were polar roses, so I was somewhat intrigued\n to see that they are, in fact, rational curves.\nOn the other hand, their rationality follows immediately from the rationality of the circle\n (which itself follows from the existence of Pythagorean triples).\nIf I were more experienced with manipulating Chebyshev polynomials or willing to set up a\n linear system in (way too) many terms, I might have considered attempting to find\n an implicit form for them as well.\n\nDiagrams created with Sympy and Matplotlib.\n\n", + "supporting": [ + "index_files" + ], + "filters": [], + "includes": {} + } +} \ No newline at end of file diff --git a/_freeze/posts/stereo/2/index/figure-html/cell-5-output-1.png b/_freeze/posts/stereo/2/index/figure-html/cell-5-output-1.png new file mode 100644 index 0000000..35d0d33 Binary files /dev/null and b/_freeze/posts/stereo/2/index/figure-html/cell-5-output-1.png differ diff --git a/_freeze/posts/stereo/2/index/figure-html/cell-6-output-1.png b/_freeze/posts/stereo/2/index/figure-html/cell-6-output-1.png new file mode 100644 index 0000000..bb71449 Binary files /dev/null and b/_freeze/posts/stereo/2/index/figure-html/cell-6-output-1.png differ diff --git a/_freeze/posts/stereo/2/index/figure-html/cell-7-output-1.png b/_freeze/posts/stereo/2/index/figure-html/cell-7-output-1.png new file mode 100644 index 0000000..9aa671b Binary files /dev/null and b/_freeze/posts/stereo/2/index/figure-html/cell-7-output-1.png differ diff --git a/index.qmd b/index.qmd index 3a63b58..d0e3b2b 100644 --- a/index.qmd +++ b/index.qmd @@ -3,6 +3,7 @@ title: "Posts by topic" listing: contents: - posts/permutations/index.* + - posts/stereo/index.* - posts/chebyshev/index.* - posts/pentagons/index.* - posts/polycount/index.* diff --git a/posts/polycount/sand-1/todo.txt b/posts/polycount/sand-1/todo.txt index b9bd8d8..a9b5086 100644 --- a/posts/polycount/sand-1/todo.txt +++ b/posts/polycount/sand-1/todo.txt @@ -3,3 +3,5 @@ mp4s saved do not get copied or linked over to the rendering folder center figures make sure to object-fit: scale + +use SympyAnimationWrapper from stereo/2 diff --git a/posts/stereo/1/approx-results/.clang-format b/posts/stereo/1/approx-results/.clang-format new file mode 100644 index 0000000..120882a --- /dev/null +++ b/posts/stereo/1/approx-results/.clang-format @@ -0,0 +1,9 @@ +BasedOnStyle: llvm +IndentWidth: 4 + +AlignAfterOpenBracket: BlockIndent +BinPackParameters: OnePerLine +BreakBeforeBraces: Custom +BraceWrapping: + AfterFunction: true +PointerAlignment: Left diff --git a/posts/stereo/1/approx-results/.gitignore b/posts/stereo/1/approx-results/.gitignore new file mode 100644 index 0000000..07337da --- /dev/null +++ b/posts/stereo/1/approx-results/.gitignore @@ -0,0 +1,2 @@ +compare_complex +compare_math diff --git a/posts/stereo/1/approx-results/Makefile b/posts/stereo/1/approx-results/Makefile new file mode 100644 index 0000000..757487b --- /dev/null +++ b/posts/stereo/1/approx-results/Makefile @@ -0,0 +1,14 @@ +default: report + +clean: + rm -rf compare_complex compare_math results_complex.txt results_math.txt + +report: compare_complex compare_math + ./compare_complex > results_complex.txt + ./compare_math > results_math.txt + +compare_complex: + gcc main.c stereo_complex.c -lm -DUSE_COMPLEX -o compare_complex + +compare_math: + gcc main.c stereo_math.c -lm -o compare_math diff --git a/posts/stereo/1/approx-results/main.c b/posts/stereo/1/approx-results/main.c new file mode 100644 index 0000000..db1e705 --- /dev/null +++ b/posts/stereo/1/approx-results/main.c @@ -0,0 +1,196 @@ +#include +#include +#include +#include +#include +#include + +#define STR_RED "\x1b[31m" +#define STR_GREEN "\x1b[32m" +#define STR_NORM "\x1b[m" + +#define SECONDS_PER_NANOSECOND 1000000000 +#define NUM_STATS 100000 +#define SQRT_NUM_STATS 100 +#define NUM_LOOPS 10000000 + +struct circle { + double c; + double s; +}; + +#ifdef USE_COMPLEX +#define EXTRACT_COSINE(a) creal(a) +#define EXTRACT_SINE(a) cimag(a) +#define CIRCLE_TYPE double complex +#else +#define EXTRACT_COSINE(a) (a.c) +#define EXTRACT_SINE(a) (a.s) +#define CIRCLE_TYPE struct circle +#endif + +void turn_update(double turn, CIRCLE_TYPE* result); +void approx_turn_update(double turn, CIRCLE_TYPE* result); + +void print_errors( + const double* inputs, + const CIRCLE_TYPE* ideals, + const CIRCLE_TYPE* approxs, + int n +) +{ + double c_error, s_error; + double largest_c_error = 0, largest_s_error = 0; + size_t largest_c_index, largest_s_index; + + double total_c_error = 0, total_s_error = 0; + + size_t i; + CIRCLE_TYPE ideal, approx; + + for (i = 0; i < n; i++) { + ideal = ideals[i]; + approx = approxs[i]; + + // squared error in c components + c_error = EXTRACT_COSINE(ideal) - EXTRACT_COSINE(approx); + c_error *= c_error; + // squared error in s components + s_error = EXTRACT_SINE(ideal) - EXTRACT_SINE(approx); + s_error *= s_error; + + if (largest_c_error < c_error) { + largest_c_error = c_error; + largest_c_index = i; + } + + if (largest_s_error < s_error) { + largest_s_error = s_error; + largest_s_index = i; + } + + total_c_error += c_error; + total_s_error += s_error; + } + // these now contain the *average* squared error + total_c_error /= (double)n; + total_s_error /= (double)n; + + printf( + "Squared error in cosines (%d runs): \n" + "\tAverage: %f (%f%% error)\n" + "\tLargest: %f (%f%% error)\n" + "\t\tInput:\t\t%f\n" + "\t\tValue:\t\t%f\n" + "\t\tApproximation:\t%f\n", + NUM_STATS, total_c_error, sqrt(total_c_error) * SQRT_NUM_STATS, + largest_c_error, sqrt(largest_c_error) * SQRT_NUM_STATS, + inputs[largest_c_index], EXTRACT_COSINE(ideals[largest_c_index]), + EXTRACT_COSINE(approxs[largest_c_index]) + ); + printf( + "Squared error in sines (%d runs): \n" + "\tAverage: %f (%f%% error)\n" + "\tLargest: %f (%f%% error)\n" + "\t\tInput:\t\t%f\n" + "\t\tValue:\t\t%f\n" + "\t\tApproximation:\t%f\n", + NUM_STATS, total_s_error, sqrt(total_s_error) * SQRT_NUM_STATS, + largest_s_error, sqrt(largest_s_error) * SQRT_NUM_STATS, + inputs[largest_s_index], EXTRACT_SINE(ideals[largest_s_index]), + EXTRACT_SINE(approxs[largest_s_index]) + ); +} + +// time the length of the computation `f` in nanoseconds +// using a macro rather than a function taking a function pointer for +// more accurate time +#define TIME_COMPUTATION(f, inputs, results, n_stats, n_loop, time) \ + do { \ + size_t i; \ + struct timespec tp1; \ + struct timespec tp2; \ + \ + for (i = 0; i < n_stats; i++) { \ + f(inputs[i], results + i); \ + } \ + CIRCLE_TYPE temp; \ + clock_gettime(CLOCK_MONOTONIC, &tp1); \ + for (i = 0; i < n_stats; i++) { \ + f(rand() / (double)RAND_MAX, &temp); \ + } \ + clock_gettime(CLOCK_MONOTONIC, &tp2); \ + \ + time = SECONDS_PER_NANOSECOND * (tp2.tv_sec - tp1.tv_sec) + \ + (tp2.tv_nsec - tp1.tv_nsec); \ + } while (0) + +int main(int argn, char** args) +{ + long trig_time, rat_time; + + double rands[NUM_STATS]; + CIRCLE_TYPE trigs[NUM_STATS]; + CIRCLE_TYPE rats[NUM_STATS]; + + size_t i; + for (i = 0; i < NUM_STATS; i++) { + rands[i] = rand() / (double)RAND_MAX; + } + + TIME_COMPUTATION( + turn_update, rands, trigs, NUM_STATS, NUM_LOOPS, trig_time + ); + printf( +#ifdef USE_COMPLEX + "Timing for %d complex.h cexp:\t%ldns\n", +#else + "Timing for %d math.h sin and cos:\t%ldns\n", +#endif + NUM_LOOPS, trig_time + ); + + TIME_COMPUTATION( + approx_turn_update, rands, rats, NUM_STATS, NUM_LOOPS, rat_time + ); + printf("Timing for %d approximations:\t%ldns\n", NUM_LOOPS, rat_time); + + // Report results + long diff = rat_time - trig_time; + double frac_speed; + if (diff > 0) { + frac_speed = rat_time / (double)trig_time; + + // Disable colors for non-terminal output + if (isatty(STDOUT_FILENO)) { + printf( + STR_RED "stdlib" STR_NORM " faster, speedup: %ldns (%2.2fx)\n", + diff, frac_speed + ); + } else { + printf( + "stdlib faster, speedup: %ldns (%2.2fx)\n", diff, frac_speed + ); + } + } else { + frac_speed = trig_time / (double)rat_time; + + // Disable colors for non-terminal output + if (isatty(STDOUT_FILENO)) { + printf( + STR_GREEN "Approximation" STR_NORM + " faster, speedup: %ldns (%2.2fx)\n", + -diff, frac_speed + ); + } else { + printf( + "Approximation faster, speedup: %ldns (%2.2fx)\n", -diff, + frac_speed + ); + } + + print_errors(rands, trigs, rats, NUM_STATS); + } + + return 0; +} diff --git a/posts/stereo/1/approx-results/results_complex.txt b/posts/stereo/1/approx-results/results_complex.txt new file mode 100644 index 0000000..27b4634 --- /dev/null +++ b/posts/stereo/1/approx-results/results_complex.txt @@ -0,0 +1,15 @@ +Timing for 10000000 complex.h cexp: 2864133ns +Timing for 10000000 approximations: 1321049ns +Approximation faster, speedup: 1543084ns (2.17x) +Squared error in cosines (100000 runs): + Average: 0.000051 (0.713743% error) + Largest: 0.000174 (1.320551% error) + Input: 0.729202 + Value: -0.659428 + Approximation: -0.672634 +Squared error in sines (100000 runs): + Average: 0.000070 (0.835334% error) + Largest: 0.000288 (1.698413% error) + Input: 0.842206 + Value: 0.475669 + Approximation: 0.458685 diff --git a/posts/stereo/1/approx-results/results_math.txt b/posts/stereo/1/approx-results/results_math.txt new file mode 100644 index 0000000..11dbb60 --- /dev/null +++ b/posts/stereo/1/approx-results/results_math.txt @@ -0,0 +1,15 @@ +Timing for 10000000 math.h sin and cos: 2822266ns +Timing for 10000000 approximations: 1223832ns +Approximation faster, speedup: 1598434ns (2.31x) +Squared error in cosines (100000 runs): + Average: 0.000051 (0.713743% error) + Largest: 0.000174 (1.320551% error) + Input: 0.729202 + Value: -0.659428 + Approximation: -0.672634 +Squared error in sines (100000 runs): + Average: 0.000070 (0.835334% error) + Largest: 0.000288 (1.698413% error) + Input: 0.842206 + Value: 0.475669 + Approximation: 0.458685 diff --git a/posts/stereo/1/approx-results/stereo_complex.c b/posts/stereo/1/approx-results/stereo_complex.c new file mode 100644 index 0000000..c63e4b3 --- /dev/null +++ b/posts/stereo/1/approx-results/stereo_complex.c @@ -0,0 +1,27 @@ +#include +#include + +double complex complex_approx_turn(double turn) +{ + const static double a = 8 * M_SQRT2 / 3 - 3; + const static double b = 4 - 8 * M_SQRT2 / 3; + + double p = turn * (b * turn * turn + a); + double q = p * p; + double r = 1 + q; + double c = (1 - q) / r, s = 2 * p / r; + + return (c * c - s * s) + I * (2 * c * s); +} + +double complex complex_turn(double turn) { return cexp(I * M_PI * turn); } + +void turn_update(double turn, double complex* result) +{ + *result = complex_turn(turn); +} + +void approx_turn_update(double turn, double complex* result) +{ + *result = complex_approx_turn(turn); +} diff --git a/posts/stereo/1/approx-results/stereo_math.c b/posts/stereo/1/approx-results/stereo_math.c new file mode 100644 index 0000000..a7025f6 --- /dev/null +++ b/posts/stereo/1/approx-results/stereo_math.c @@ -0,0 +1,27 @@ +#include + +struct circle { + double c; + double s; +}; + +void approx_turn_update(double turn, struct circle* ret) +{ + static const double a = 8 * M_SQRT2 / 3 - 3; + static const double b = 4 - 8 * M_SQRT2 / 3; + + double p = turn * (b * turn * turn + a); + double q = p * p; + double r = 1 + q; + double c = (1 - q) / r, s = 2 * p / r; + + ret->c = c * c - s * s; + ret->s = 2 * c * s; +} + +void turn_update(double turn, struct circle* ret) +{ + double arg = M_PI * turn; + ret->c = cos(arg); + ret->s = sin(arg); +} diff --git a/posts/stereo/1/index.qmd b/posts/stereo/1/index.qmd new file mode 100644 index 0000000..7e00740 --- /dev/null +++ b/posts/stereo/1/index.qmd @@ -0,0 +1,707 @@ +--- +title: "Algebraic Rotations and Stereography" +description: | + How do you rotate in 2D and 3D without standard trigonometry? +format: + html: + html-math-method: katex +jupyter: python3 +date: "2021-09-26" +date-modified: "2025-06-28" +categories: + - geometry + - symmetry +--- + + + +```{python} +#| echo: false + +from IPython.display import display, Code +import sympy +from sympy.abc import t +``` + + +Expressing rotation mathematically can be rather challenging, even in two dimensional space. +Typically, the subject is approached from complicated-looking rotation matrices, + with the occasional accompanying comment about complex numbers. +The challenge only increases when examining three dimensional space, + where everyone has different ideas about the best implementations. +In what follows, I attempt to provide a derivation of 2D and 3D rotations based in intuition + and without the use of trigonometry. + + +Complex Numbers: 2D Rotations +----------------------------- + +As points in the plane, multiplication between complex numbers is frequently analyzed + as a combination of scaling (a ratio of two scales) and rotation (a point on the complex unit circle). +However, the machinery used in these analyses usually expresses complex numbers in a polar form + involving the complex exponential. +This leverages prior knowledge of trigonometry and some calculus, and in fact obfuscates + the algebraic properties of complex numbers. +Neither concept is truly necessary to understand points on the complex unit circle, + or complex number multiplication in general. + + +### A review of complex multiplication + +A complex number $z = a + bi$ multiplies with another complex number + $w = c + di$ in the obvious way: by distributing over addition. + +$$ +zw = (a + bi)(c + di) = ac + bci + adi + bdi^2 = (ac - bd) + i(ad + bc) +$$ + +*z* also has associated to it a conjugate $z^* = a - bi$. +The product $zz^*$ is the *norm* $a^2 + b^2$, a positive real quantity whose square root + is typically called the *magnitude* of the complex number. +Taking the square root is frequently unnecessary, as many of the properties of + this quantity do not depend on this normalization. + +One such property is that for two complex numbers *z* and *w*, + the norm of the product is the product of the norm. + +$$ +\begin{align*} + (zw)(zw)^* &= (ac - bd)^2 + (ad + bc)^2 + \\ + &= a^2c^2 - 2abcd + b^2d^2 + a^2d^2 + 2abcd + b^2c^2 + \\ + &= a^2c^2 + b^2d^2 + a^2d^2 + b^2c^2 + = a^2(c^2 + d^2) + b^2(c^2 + d^2) + \\ + &= (a^2 + b^2)(c^2 + d^2) +\end{align*} +$$ + +Conjugation also distributes over multiplication, and complex numbers possess multiplicative inverses. +Also, the norm of the multiplicative inverse is the inverse of the norm. + +$$ +\begin{gather*} + z^*w^* = (a - bi)(c - di) = ac - bci - adi + bdi^2 + \\ + = (ac - bd) - i(ad + bc) = (zw)^{*} + \\ + {1 \over z} = {z^{*} \over z z^{*}} = {a - bi \over a^2 + b^2} + \\[14pt] + \left({1 \over z}\right) + \left({1 \over z}\right)^{*} + = {1 \over zz^{*}} = {1 \over a^2 + b^2} +\end{gather*} +$$ + +A final interesting note about complex multiplication is the product $z^{*} w$ + +$$ +z^{*} w = (a - bi)(c + di) = ac - bci + adi - bdi^2 = (ac + bd) + i(ad - bc) +$$ + +The real part is the same as the dot product $(a, c) \cdot (b, d)$ and the imaginary part + looks like the determinant of $\begin{pmatrix}a & b \\ c & d\end{pmatrix}$. +This shows the inherent value of complex arithmetic, as both of these quantities + have important interpretations in vector algebra. + + +### Onto the Circle + +If $zz^{*} = a^2 + b^2 = 1$, then *z* lies on the complex unit circle. +Since the norm is 1, a general complex number *w* has the same norm as its product with *z*, *zw*. +Specifically, if *w* is 1, multiplication by *z* can be seen as the rotation that maps 1 to *z*. +But how do you describe a point on the unit circle? + +Simple. +For any complex *w*, there are three other complex numbers that are guaranteed + to share its norm: $-w, w^{*}$, and $-w^{*}$. +Dividing any one of these numbers by another results in a complex number on the unit circle. +Dividing *w* by -*w* is boring, since that just results in -1. +However, division by $w^*$ turns out to be important: + +$$ +\begin{align*} + z &= {w \over w^*} = + \left( {c + di \over c - di} \right) + \left( {c + di \over c + di} \right) = 1 + \\[8pt] + &= {(c^2 - d^2) + (2cd)i \over c^2 + d^2} +\end{align*} +$$ + +Dividing the numerator and denominator of this expression through by $c^2$ yields an expression + in $t = d/c$, where t can be interpreted as the slope of the line joining 0 and *w*. +Slopes range from $-\infty$ to $\infty$, which means that the variable *t* does as well. + +$$ +\begin{align*} + {(c^2 - d^2) + (2cd)i \over c^2 + d^2} + &= {1/c^2 \over 1/c^2} \cdot {(c^2 - d^2) + (2cd)i \over c^2 + d^2} + \\ + &= {(1 - d^2/c^2) + (2d/c)i \over 1 + d^2/c^2} + = {(1 - t^2) + (2t)i \over 1 + t^2} + \\ + &= {1 - t^2 \over 1 + t^2} + i{2t \over 1 + t^2} = z(t) +\end{align*} +$$ + +This gives us a point *z* on the unit circle. +Its reciprocal is: + +$$ +{1 \over z} = {z^{*} \over zz^{*}} = z^{*} +$$ + +Applying the same operation to *z* as we did to *w* gives us $z / z^{*} = z^2$, + or as an action on points, the rotation from 1 to $z^2$. + +This demonstrates the decomposition of the rotation into two parts. +For a general complex number *w* which may not lie on the unit circle, multiplication by + this number dilates the complex plane as well rotating 1 toward the direction of *w*. +Then, division by $w^{*}$ undoes the dilation and performs the rotation again. +General (i.e., non-unit) rotations *must* occur in pairs. + + +### Two Halves, Twice Over + +Plugging in -1, 0, and 1 for *t* reveals some interesting behavior: + $z(0) = 1 + 0i$, $z(1) = 0 + 1i$, and $z(-1) = 0 - 1i$. +Assuming continuity, this shows that $t \in (-1, 1)$ plots the semicircle in the right half-plane. +The other half of the circle comes from *t* with magnitude greater than 1, + and the point where $z(t) = -1$ exists if we admit the extra point $t = \pm \infty$. + +Is it possible to get the other half into a nicer interval? +No problem, just double the rotation. + +$$ +\begin{align*} + z(t)^2 &= (\Re[z]^2 - \Im[z]^2) + i(2\Re[z]\Im[z]) + \\[8pt] + &= {(1 - t^2)^2 - (2t)^2 \over (1 + t^2)^2} + i{2(1 - t^2)(2t) \over (1 + t^2)^2} + \\[8pt] + &= {1 - 6t^2 + t^4 \over 1 + 2t^2 + t^4} + + i {4t - 4t^3 \over 1 + 2t^2 + t^4} +\end{align*} +$$ + +Now $z(-1)^2 = z(1)^2 = -1$ and the entire circle fits in the interval \[-1, 1\]. +This action of squaring *z* means that as *t* ranges from $-\infty$ to $\infty$, the point + $z^2$ loops around the circle twice, since the unit rotation *z* is applied twice. + + +### Plotting Transcendence + +Let's go on a brief tangent and compare these expressions to their + transcendental trigonometric counterparts. +Graphing the real and imaginary parts separately shows how much they resemble + $\cos(\pi t)$ and $\sin(\pi t)$ around zero[^1]. + +[^1]: An approximation of a function as the ratio of two polynomials is called a + [Padé approximant](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant). + Specifically, the real part (which approximates cosine) has order \[4/4\] and the imaginary part + (which approximates sine) has order \[3/4\]. + +```{python} +#| echo: false +#| layout-ncol: 2 + +z = (1 + sympy.I * t) / (1 - sympy.I * t) + +sympy.plot( + sympy.re(z**2), + sympy.im(z**2), + (t, -2, 2), + ylim=(-2, 2), + label=[ + r"$\Re \left [ \left ( \frac{1 + it}{1 - it} \right )^2 \right ]$", + r"$\Im \left [ \left ( \frac{1 + it}{1 - it} \right )^2 \right ]$", + ], + legend=True, + xlabel=False, + ylabel=False, +) + +sympy.plot( + sympy.re(z**2) - sympy.cos(sympy.pi * t), + sympy.im(z**2) - sympy.sin(sympy.pi * t), + (t, -2, 2), + ylim=(-2, 2), + label=[ + r"Cosine error", + r"Sine error", + ], + legend=True, + xlabel=False, + ylabel=False, +) +``` + +Ideally, the zeroes of the real part should occur at $t = \pm 1/2$, since $\cos(\pm \pi/2) = 0$. +Instead, they occur at + +$$ +\begin{align*} + 0 &= 1 - 6t^2 + t^4 + = (t^2 - 2t - 1)(t^2 + 2t - 1) + \\ + &\implies t = \pm {1 \over \delta_s} = \pm(\sqrt 2 - 1) \approx \pm 0.414 +\end{align*} +$$ + +With a bit of polynomial interpolation, this can be rectified. +A cubic interpolation between the expected and actual values is given by: + +$$ +\begin{gather*} + P(x) = ax^3 + bx^2 + cx + d + \\ + \begin{pmatrix} + 1 & 0 & 0 & 0 \\ + 1 & 1/2 & 1/4 & 1/8 \\ + 1 & -1/2 & 1/4 & -1/8 \\ + 1 & 1 & 1 & 1 + \end{pmatrix} \begin{pmatrix} + d \\ + c \\ + b \\ + a + \end{pmatrix} = \begin{pmatrix} + 0 \\ + \sqrt 2 - 1 \\ + -\sqrt 2 + 1 \\ + 1 + \end{pmatrix} + \\ + \implies \begin{pmatrix} + d \\ + c \\ + b \\ + a + \end{pmatrix} = \begin{pmatrix} + 0 \\ + -3 + 8\sqrt 2/3 \\ + 0 \\ + 4 - 8\sqrt 2/3 + \end{pmatrix} + \implies P(x) \approx 0.229 x^3 + 0.771x +\end{gather*} +$$ + +::: {} +```{python} +#| echo: false +#| layout-ncol: 2 + +interpolate = (-3 + 8 * sympy.sqrt(2)/3) * t + (4 - 8 * sympy.sqrt(2)/3) * t**3 + +sympy.plot( + sympy.re(z**2).subs(t, interpolate), + sympy.im(z**2).subs(t, interpolate), + (t, -2, 2), + ylim=(-2, 2), + label=[ + r"$\Re \left [ z(P(t))^2 \right ]$", + r"$\Im \left [ z(P(t))^2 \right ]$", + ], + legend=True, + xlabel=False, + ylabel=False, +) + +sympy.plot( + sympy.re(z**2).subs(t, interpolate) - sympy.cos(sympy.pi * t), + sympy.im(z**2).subs(t, interpolate) - sympy.sin(sympy.pi * t), + (t, -2, 2), + ylim=(-2, 2), + label=[ + r"Cosine error", + r"Sine error", + ], + legend=True, + xlabel=False, + ylabel=False, +) +``` + +The error over the interval \[-1, 1\] in these approximations is around 1% (RMS). +::: + +Notably, this interpolating polynomial is entirely odd, which gives it some symmetry about 0. +Along with being a very good approximation, it has the feature that a rotation by an + *n*^th^ of a turn is about $z(P(1/n))^2$. +The approximation be improved further by taking *z* to higher powers and deriving a higher-order + interpolating polynomial. + +It is impossible to shrink this error to 0 because the derivative of the imaginary part at + $t = 1$ would be π. +But the imaginary part is a rational expression, so its derivative is also a rational expression. +Since π is transcendental, there must always be some error[^2]. + +[^2]: Even with the interpolating polynomial, the quantity is a ratio of two algebraic numbers, + which is also algebraic and not transcendental. + +Even though it is arguably easier to calculate, there probably isn't a definite benefit to this approximation. +For example, + +- Other accurate approximants for sine and cosine can be calculated directly from their power series. +- There are no inverse functions like `acos` or `asin` to complement these approximations. + - This isn't that bad, since I have refrained from describing rotations with an angle + (i.e., the quantity returned by these functions). + Even so, the inverse functions have their use cases. +- Trigonometric functions can be hardware-accelerated (at least in some FPUs). +- If no such hardware exists, approximations of `sin` and `cos` can be calculated by software beforehand + and stored in a lookup table (which is also probably involved at some stage of the FPU). + +Elaborating on the latter two points, using a lookup table means evaluation happens in roughly constant time. +A best-case analysis of the above approximation, given some value *t* is... + + +| Expression | Additions | Multiplications | Divisions | +|-----------------------------------------------------------|-----------|-----------------|-----------| +| $p = t \cdot (0.228763834 \cdot t \cdot t + 0.771236166)$ | 1 | 3 | | +| $q = p \cdot p$ | | 1 | | +| $r = 1 + q$ | 1 | | | +| $c = { 1 - q \over r}$ | 1 | | 1 | +| $s = { p + p \over r}$ | 1 | | 1 | +| real = $c \cdot c - s \cdot s$ | 1 | 2 | | +| imag = $c \cdot s + c \cdot s$ | 1 | 1 | | +| Total | 6 | 7 | 2 | + +...or 15 FLOPs. + +On a more optimistic note, a Monte Carlo test of the above approximation on my computer yields + promising results when compared with GCC's `libm` implementations of `sin`, `cos`, and `cexp`. + +```{python} +#| echo: false + +for file in [ + "approx-results/results_math.txt", + "approx-results/results_complex.txt", +]: + with open(file) as a: + display(Code(a.read())) +``` + +For the source which I used to generate this output, see the repository linked at the bottom + of this article. + +Even though results can be inaccurate, this exercise in the algebraic manipulation + of complex numbers is fairly interesting since it requires no calculus to define, + unlike sine, cosine, and their Padé approximants (and to a degree, π). + + +From Circles to Spheres +----------------------- + +The expression for complex points on the unit circle coincides with the stereographic projection of a circle. +This method is achieved by selecting a point on the circle, fixing *t* along a line + (such as the *y*-axis), and letting *t* range over all possible values. + +![ + 2D stereographic projection +](./stereo_circle.png) + +The equations of the line and circle appear in the diagram above. +Through much algebra, expressions for *x* and *y* as functions of *t* can be formed. + + +$$ +\begin{align*} + y &= t(x + 1) + \\ + 1 &= x^2 + y^2 = x^2 + (tx + t)^2 + \\ + &= x^2 + t^2x^2 + 2t^2 x + t^2 + \\[14pt] + {1 \over t^2 + 1} -\ {t^2 \over t^2 + 1} + &= x^2 + {t^2 \over t^2 + 1} 2x + \\ + {1 -\ t^2 \over t^2 + 1} + \left( {t^2 \over t^2 + 1} \right)^2 + &= \left(x + {t^2 \over t^2 + 1} \right)^2 + {t^2 \over t^2 + 1} + \\ + {1 -\ t^4 \over (t^2 + 1)^2} + {t^4 \over (t^2 + 1)^2} + &= \left(x + {t^2 \over t^2 + 1} \right)^2 + \\ + {1 \over (t^2 + 1)^2} + &= \left(x + {t^2 \over t^2 + 1} \right)^2 + \\ + {1 \over t^2 + 1} + &= x + {t^2 \over t^2 + 1} + \\ + x &= {1 -\ t^2 \over 1 + t^2} + \\ + y &= t(x + 1) + = t\left( {1 -\ t^2 \over 1 + t^2} + {1 + t^2 \over 1 + t^2} \right) + \\ + &= {2t \over 1 + t^2} +\end{align*} +$$ + +These are exactly the same expressions which appear in the real and imaginary components + of the ratio between $1 + ti$ and its conjugate. + +Compared to the algebra using complex numbers, this is quite a bit more work. +One might ask whether, given a proper arithmetic setting, this can be extended to three dimensions. +In other words, we want to find the projection of a sphere using some new number system + analogously to the circle with respect to complex numbers. + +![ + What we are *not* doing: intersecting the sphere with a line through a pole and a point in the xy plane +](./stereo_sphere.png) + + +### Algebraic Projection + +The simplest thing to do is try it out. +Let's say an extended number $h = 1 + si + tj$ has a conjugate of $h^{*} = 1 - si - tj$. +Then, following our noses[^3]: + +[^3]: We don't know whether *i* and *j* form a + [field](https://en.wikipedia.org/wiki/Field_%28mathematics%29) + or [division ring](https://en.wikipedia.org/wiki/Division_ring)), + so it's a bit of an assumption for division to be possible. + We'll ignore that. + +$$ +\begin{gather*} + {h \over h^{*}} = {1 + si + tj \over 1 - si - tj} + = \left({1 + si + tj \over 1 - si - tj}\right) + \left({1 + si + tj \over 1 + si + tj}\right) + \\[8pt] + = { + 1 + si + tj + si + s^2i^2 + stij + tj + stji + t^2j^2 + \over 1 - si - tj + si - s^2 i^2 - stij + tj - stji - t^2j^2 + } +\end{gather*} +$$ + +While there is some cancellation in the denominator, both products are rather messy. +To get more cancellation, we can add some nice algebraic properties between *i* and *j*. +If we let i and j be *anticommutative* (meaning that $ij = -ji$) then $stij$ cancels with $stji$. +So that the denominator is totally real (and therefore can be guaranteed to divide), + we can also assert that $i^2$ and $j^2$ are both real. +Then the expression becomes + +$$ +{1 + si + tj \over 1 - si - tj} + = {1 + s^2i^2 + t^2j^2 \over 1 - s^2 i^2 - t^2 j^2} + + i{2s \over 1 - s^2 i^2 - t^2 j^2} + + j{2t \over 1 - s^2 i^2 - t^2 j^2} +$$ + +If it is chosen that $i^2 = j^2 = -1$, this produces the correct equations parametrizing the unit sphere + (see [Wikipedia](https://en.wikipedia.org/wiki/Stereographic_projection#First_formulation)). + +$$ +{h \over h^{*}} = {1 - s^2 - t^2 \over 1 + s^2 + t^2} + + i{2s \over 1 + s^2 + t^2} + + j{2t \over 1 + s^2 + t^2} +$$ + +One can also check that this makes the squares of the real, *i*, and *j* components sum to 1. + +If both *i* and *j* anticommute, then their product also anticommutes with both *i* and *j*. + +$$ +\textcolor{red}{ij}j = -j\textcolor{red}{ij} + ~,~ i\textcolor{red}{ij} = -\textcolor{red}{ij}i +$$ + +Calling this product *k* and noticing that $k^2 = (ij)(ij) = -ijji = i^2 = -1$ + completely characterizes the [quaternions](https://en.wikipedia.org/wiki/Quaternion). + + +### Other projections + +Choosing different values for $i^2$ and $j^2$ yield different shapes than a sphere. +If you know a little group theory, you might know there are only two nonabelian + (noncommutative) groups of order 8: + +- the [quaternion group](https://en.wikipedia.org/wiki/Quaternion_group) +- and the [ + dihedral group of degree 4 + ](https://en.wikipedia.org/wiki/Examples_of_groups#dihedral_group_of_order_8). + +In the latter group, *j* and *k* are both imaginary, but square to 1 (*i* still squares to -1)[^4]. + +[^4]: I'm being a bit careless with the meanings of "1" and "-1" here. + Properly, these are the group identity and another group element which commutes with all others. + +Changing the sign of one (or both) of the imaginary squares in the expression $h / h^{*}$ above + switches the multiplicative structure from quaternions to the dihedral group. +In this group, picking $i^2 = -j^2 = -1$ parametrizes a hyperboloid of one sheet, + and picking $i^2 = j^2 = 1$ parametrizes a hyperboloid of two sheets. + + +Quaternions and Rotation +------------------------ + +We've already established that complex numbers are useful for describing 2D rotations. +They also elegantly describe the stereographic projection of a circle. +Consequently, since quaternions elegantly describe the stereographic projection of a sphere, + they are useful for 3D rotations. + +Since the imaginary units *i*, *j*, and *k* are anticommutative, a general quaternion does not + commute with other quaternions like complex numbers do. +This embodies a difficulty with 3D rotations in general: unlike 2D rotations, they do not commute. + +Only three of the four components in a quaternion are necessary to describe a point in 2D space. +The *ijk* subspace is ideal since the three imaginary axes are symmetric (i.e., they all square to -1). +Quaternions in this subspace are called *vectors*. +Another reason for using the *ijk* subspace comes from considering a point + $u = ai + bj + ck$ on the unit sphere ($a^2 + b^2 + c^2 = 1$). +Then the square of this point is: + +$$ +\begin{align*} + u^2 &= (ai + bj + ck)(ai + bj + ck) + \\ + &= a^2i^2 + abij + acik + abji + b^2j^2 + bcjk + acki + bckj + c^2k^2 + \\ + &= (a^2 + b^2 + c^2)(-1) + ab(ij + ji) + ac(ik + ki) + bc(jk + kj) + \\ + &= -1 + 0 + 0 + 0 +\end{align*} +$$ + +This property means that *u* behaves similarly to a typical imaginary unit. +If we form pseudo-complex numbers of the form $a + b u$, then ${1 + tu \over 1 - tu} = \alpha + \beta u$ + specifies *some kind* of rotation in terms of the parameter *t*. +As with complex numbers, the inverse rotation is the conjugate $1 - tu \over 1 + tu$. + +Another useful feature of quaternion algebra involves symmetry transformations of the imaginary units. +If an imaginary unit (one of *i*, *j*, *k*) is left-multiplied by a quaternion $q$ + and right-multiplied by its conjugate $q^{*}$, then the result is still imaginary. +In other words, if *p* is a vector, then $qpq^{*}$ is also a vector since + *i*, *j*, and *k* form a basis and multiplication distributes over addition. +I will not demonstrate this fact, as it requires a large amount of algebra. + +These two features combine to characterize a transformation for a vector quaternion. +For a vector *p*, if $q_u(t)$ is a "rotation" as above for a point *u* on the unit sphere, + then another vector (the image of *p*) can be described by + +$$ +\begin{gather*} + qpq^* = qpq^{-1} + = \left({1 + tu \over 1 - tu}\right) p \left({1 - tu \over 1 + tu}\right) + = (\alpha + \beta u)p(\alpha - \beta u) + \\[4pt] + = (\alpha + \beta u)(\alpha p - \beta pu) + = \alpha^2 p - \alpha \beta pu + \alpha \beta up - \beta^2 upu +\end{gather*} +$$ + + +### Testing a Transformation + +3D space contains 2D subspaces. +If this transformation is a rotation, then a 2D subspace will be affected in a way + which can be described more simply by complex numbers. +For example, if $u = k$ and $p = xi + yj$, then this can be expanded as: + +$$ +\begin{align*} + qpq^{-1} &= (\alpha + \beta k)(xi + yj)(\alpha - \beta k) + \\[10pt] + &= \alpha^2(xi + yj) - \alpha \beta (xi + yj)k + \alpha \beta k(xi + yj) - \beta^2 k(xi + yj)k + \vphantom{\over} + \\[10pt] + &= \alpha^2(xi + yj) - \alpha \beta(-xj + yi) + \alpha \beta(xj - yi) - \beta^2(xi + yj) + \\[10pt] + &= (\alpha^2 - \beta^2)(xi + yj) + 2\alpha \beta(xj - yi) + \\[10pt] + &= [(\alpha^2 - \beta^2)x - 2\alpha \beta y]i + [(\alpha^2 - \beta^2)y + 2\alpha \beta x]j +\end{align*} +$$ + +If *p* is rewritten as a vector (in the linear algebra sense), then this final expression + can be rewritten as a linear transformation of *p*: + +$$ +\begin{align*} + \begin{pmatrix} + (\alpha^2 - \beta^2)x - 2 \alpha \beta y \\ + (\alpha^2 - \beta^2)y + 2 \alpha \beta x + \end{pmatrix} + &= \begin{pmatrix} + \alpha^2 - \beta^2 & -2\alpha \beta \\ + 2\alpha \beta & \alpha^2 - \beta^2 + \end{pmatrix} + \begin{pmatrix} + x \\ + y + \end{pmatrix} + \\ + &= \begin{pmatrix} + \alpha & -\beta \\ + \beta & \alpha + \end{pmatrix}^2 + \begin{pmatrix} + x \\ + y + \end{pmatrix} +\end{align*} +$$ + +The complex number $\alpha + \beta i$ is isomorphic to the matrix on the second line[^5]. +Since $\alpha$ and $\beta$ were chosen so to lie on "a" complex unit circle, + this means that the vector (*x*, *y*) is rotated by $\alpha + \beta i$ twice. +This also demonstrates that *u* specifies the axis of rotation, or equivalently, + the plane of a great circle. +This means that the "some kind of rotation" specified by *q* is a rotation around + the great circle normal to *u*. + +[^5]: Observe this by comparing the entries of the square of the matrix and the square of the complex number. + +The double rotation resembles the half-steps in the complex numbers, + where $1 + ti$ performs a dilation, ${1 \over 1 - ti}$ reverses it, and together they describe a rotation. +The form $qpq^{-1}$ is similar to this -- the $q$ on the left performs some transformation which is (partially) + undone by the $q^{-1}$ on the right. +In this case, since *q*'s conjugate is its inverse and quaternions do not commute, the only thing to do without + the transformation being trivial is to sandwich *p* between them. + +As shown above, the product of a vector with itself should be totally real + and equal to the negative of the norm of *p*, the sum of the squares of its components. +The rotated *p* shares the same norm as *p*, so it should equal $p^2$. +Indeed, this is the case, as + +$$ +(qpq^{-1})^2 = (qpq^{-1})(qpq^{-1}) = q p p q^{-1} = p^2 q q^{-1} = p^2 +$$ + +As a final remark, due to rotation through quaternions being doubled, the interpolating polynomial + used above to smooth out double rotations can also be used here. +That is, $q_u(P(t))pq_u^{-1}(P(t))$ rotates *p* by (approximately) the fraction of a turn specified + by *t* through the great circle which intersects the plane normal to *u*. + + +Closing +------- + +It is difficult to find a treatment of rotations which does not hesitate + to use the complex exponential $e^{ix} = \cos(x) + i \sin(x)$. +The true criterion for rotations is simply that a point lies on the unit circle. +Perhaps this contributes to why understanding the 3D situation with quaternions is challenging for many. +Past the barrier to entry, I believe them to be rather intuitive. +I have outlined some futher benefits of this approach in this post. + +As a disclaimer, this post was (at least subconsciously) inspired by + [this video](https://www.youtube.com/watch?v=d4EgbgTm0Bg) by 3blue1brown on YouTube + (though I had not seen the video in years before checking that I wasn't just plagiarizing it). +It *also* uses stereographic projections of the circle and sphere to describe rotations + in the complex plane and quaternion vector space. +However, I feel like it fails to provide an algebraic motivation for quaternions, + or even stereography in the first place. +Hopefully, my remarks on the algebraic approach can be used to augment the information in the video. + + +Diagrams created with GeoGebra and Matplotlib. +Repository with approximations (as well as GeoGebra files) available [here](https://github.com/queue-miscreant/approx-trig). diff --git a/posts/stereo/1/stereo_circle.png b/posts/stereo/1/stereo_circle.png new file mode 100644 index 0000000..fc967b2 Binary files /dev/null and b/posts/stereo/1/stereo_circle.png differ diff --git a/posts/stereo/1/stereo_sphere.png b/posts/stereo/1/stereo_sphere.png new file mode 100644 index 0000000..52ec591 Binary files /dev/null and b/posts/stereo/1/stereo_sphere.png differ diff --git a/posts/stereo/2/.gitignore b/posts/stereo/2/.gitignore new file mode 100644 index 0000000..2641667 --- /dev/null +++ b/posts/stereo/2/.gitignore @@ -0,0 +1 @@ +*.mp4 diff --git a/posts/stereo/2/anim.py b/posts/stereo/2/anim.py new file mode 100644 index 0000000..357714f --- /dev/null +++ b/posts/stereo/2/anim.py @@ -0,0 +1,67 @@ +import logging +import time +from pathlib import Path +from threading import Thread +from typing import Callable + +from matplotlib import animation, pyplot as plt + + +def animate( + fig, frames: list[int] | None, **kwargs +) -> Callable[..., animation.FuncAnimation]: + """Decorator-style wrapper for FuncAnimation.""" + # Why isn't this how the function is exposed by default? + return lambda func: animation.FuncAnimation(fig, func, frames, **kwargs) + +# writer = animation.writers["ffmpeg"](fps=15, metadata={"artist": "Me"}) + +def bindfig(fig): + """ + Create a function which is compatible with the signature of `pyplot.figure`, + but alters the already-existing figure `fig` instead. + """ + def ret(**kwargs): + for i, j in kwargs.items(): + if i == "figsize": + if j is not None: + fig.set_figwidth(j[0]) + fig.set_figheight(j[1]) + continue + fig.__dict__["set_" + i](j) + return fig + return ret + +class SympyAnimationWrapper: + """ + Context manager for binding a figure for animation. + """ + def __init__(self, filename: Path | str): + self._filename = filename if isinstance(filename, Path) else Path(filename) + self._current_figure = plt.gcf() + + self._figure_temp = plt.figure + + def __enter__(self): + plt.figure = bindfig(self._current_figure) + return self.animate + + def __exit__(self, *_): + plt.figure = self._figure_temp + + def animate(self, *args, **kwargs): + def wrapper(func): + ret = animate(self._current_figure, *args, **kwargs)(func) + current_save = ret.save + def save(*args, **kwargs): + if self._filename.exists(): + logging.error("Ignoring saving animation '%s': file already exists", self._filename) + return + thread = Thread(target=lambda: time.sleep(1) or plt.close()) + thread.run() + current_save(self._filename, *args, **kwargs) + ret.save = save + + return ret + + return wrapper diff --git a/posts/stereo/2/index.qmd b/posts/stereo/2/index.qmd new file mode 100644 index 0000000..7631557 --- /dev/null +++ b/posts/stereo/2/index.qmd @@ -0,0 +1,783 @@ +--- +title: "Further Notes on Algebraic Stereography" +description: | + How do you rotate in 2D and 3D without standard trigonometry? +format: + html: + html-math-method: katex +jupyter: python3 +date: "2021-10-10" +date-modified: "2025-06-30" +categories: + - algebra + - complex analysis + - polar roses + - generating functions +--- + + + +```{python} +#| echo: false + +from dataclasses import dataclass + +import matplotlib.pyplot as plt +import sympy +from sympy.abc import t + +from anim import SympyAnimationWrapper + +# rational circle function +o = (1 + sympy.I*t) / (1 - sympy.I*t) + +# cache some values for plotting later +def generate_os(amt): + ret = sympy.sympify(1) + for _ in range(amt): + yield ret + ret = (ret * o).expand().cancel().simplify() + +os = list(generate_os(11)) +cs = [sympy.re(i).simplify() for i in os] +ss = [sympy.im(i).simplify() for i in os] +``` + + +In my previous post, I discussed the stereographic projection of a circle as it pertains + to complex numbers, as well as its applications in 2D and 3D rotation. +In an effort to document more interesting facts about this mathematical object + (of which scarce information is immediately available online), + I will now elaborate on more of its properties. + + +Chebyshev Polynomials +--------------------- + +[Previously](/posts/chebyshev/1), I derived the + [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) + with the archetypal complex exponential. +These polynomials express the sines and cosines of a multiple of an angle from + the sine and cosine of the base angle. +Where $T_n(t)$ are Chebyshev polynomials of the first kind and $U_n(t)$ are those of the second kind, + +$$ +\begin{gather*} + \cos(n \theta) = T_n(\cos(\theta)) + \\ + \sin(n \theta) = U_{n - 1}(\cos(\theta)) \sin(\theta) +\end{gather*} +$$ + +The complex exponential derivation begins by squaring and developing a second-order recurrence. + +$$ +\begin{align*} + (e^{i\theta})^2 &= (\cos + i\sin)^2 + \\ + &= \cos^2 + 2i\cos \cdot \sin - \sin^2 + (0 = \cos^2 + \sin^2 - 1) + \\ + &= 2\cos^2 + 2i\cos \cdot \sin - 1 + \\ + &= 2\cos \cdot (\cos + i\sin) - 1 + \\ + &= 2\cos(\theta)e^{i\theta} - 1 + \\ + (e^{i\theta})^{n+2} &= 2\cos(\theta)(e^{i\theta})^{n+1} - (e^{i\theta})^n +\end{align*} +$$ + +This recurrence relation can then be used to obtain the Chebyshev polynomials, and hence, + the expressions using sine and cosine above. +Presented this way with such a simple derivation, it appears as though these relationships + are inherently trigonometric. + +However, these polynomials actually have *nothing* to do with sine and cosine on their own. +For one, [they appear in graph theory](/posts/chebyshev/2), and for two, + the importance of the complex exponential is overstated. +$e^{i\theta}$ really just specifies a point on the complex unit circle. +This property is used on the second line to coax the equation into a quadratic in $e^{i\theta}$. +This is also the *only* property upon which the recurrence depends; all else is algebraic manipulation. + + +### Back to the Stereograph + +Knowing this, let's start over with the stereographic projection of the circle: + +$$ +o_1(t) = {1 + it \over 1 - it} + = {1 - t^2 \over 1 + t^2} + i {2t \over 1 + t^2} + = \text{c}_1 + i\text{s}_1 +$$ + +The subscript "1" is because as *t* ranges over $(-\infty, \infty)$, the function loops once + around the unit circle. +Taking this to higher powers keeps points on the circle since all points on the circle + have a norm of 1. +It also makes more loops around the circle, which we can denote by larger subscripts: + +$$ +\begin{align*} + o_n &= (o_1)^n + = \left( {1 + it \over 1 - it} \right)^n + \\ + \text{c}_n + i\text{s}_n + &= (\text{c}_1 + i\text{s}_1)^n +\end{align*} +$$ + +This mirrors raising the complex exponential to a power + (which loops over the range $(-\pi, \pi)$ instead). +The final line is analogous to de Moivre's formula, but in a form where everything is + a ratio of polynomials in *t*. +This means that the Chebyshev polynomials can be obtained directly from these rational expressions: + +$$ +\begin{align*} + o_2 = (o_1)^2 &= (\text{c}_1 + i\text{s}_1)^2 + \\ + &= \text{c}_1^2 + 2i\text{c}_1\text{s}_1 - \text{s}_1^2 + + (0 = \text{c}_1^2 + \text{s}_1^2 - 1) + \\ + &= 2\text{c}_1^2 + 2i\text{c}_1\text{s}_1 - 1 + \\ + &= 2\text{c}_1(\text{c}_1 + i\text{s}_1) - 1 + \\ + &= 2\text{c}_1 o_1 - 1 + \\ + o_2 \cdot (o_1)^n &= 2\text{c}_1 o_1 \cdot (o_1)^n - (o_1)^n + \\ + o_{n+2} &= 2\text{c}_1 o_{n+1} - o_n +\end{align*} +$$ + +This matches the earlier recurrence relation with the complex exponential and therefore + the recurrence relation of the Chebyshev polynomials. +It also means that the the rational functions obey the same relationship as sine and cosine: + +$$ +\begin{matrix} + \begin{gather*} + \text{c}_n = T_n(\text{c}_1) + \\ + \text{s}_n = U_{n-1}(\text{c}_1) \text{s}_1 + \end{gather*} + & \text{where } + \text{c}_1 = {1 - t^2 \over 1 + t^2}, & + \text{s}_1 = {2t \over 1 + t^2} +\end{matrix} +$$ + +Thus, the Chebyshev polynomials are tied to (coordinates on) circles, + rather than explicitly to the trigonometric functions. +It is a bit strange that these polynomials are in terms of rational functions, but no stranger + than them being in terms of *ir*rational functions like sine and cosine. + + +Calculus +-------- + +Since these functions behave similarly to sine and cosine, one might wonder about + the nature of these expressions in the context of calculus. + +For comparison, the complex exponential (as it is a parallel construction) has a simple derivative[^1]. +Since the exponential function is its own derivative, the expression acquires + an imaginary coefficient through the chain rule. + +[^1]: This is forgoing the fact that complex derivatives require more care than their real counterparts. + It matters slightly less in this case since this function is complex-valued, but has a real parameter. + +$$ +\begin{align*} + e^{it} &= \cos(t) + i\sin(t) + \\ + {d \over dt} e^{it} + &= {d \over dt} \cos(t) + {d \over dt} i\sin(t) + \\ + i e^{it} &= -\sin(t) + i\cos(t) + \\ + i[\cos(t) + i\sin(t)] + &\stackrel{\checkmark}{=} -\sin(t) + i\cos(t) +\end{align*} +$$ + +Meanwhile, the complex stereograph has derivative + +$$ +\begin{align*} + {d \over dt} o_1(t) &= {d \over dt} {1 + it \over 1 - it} + = {i(1 - it) + i(1 + it) \over (1 - it)^2} + \\ + &= {2i \over (1 - it)^2} + = {2i(1 + it)^2 \over (1 + t^2)^2} + = {2i(1 - t^2 + 2it) \over (1 + t^2)^2} + \\ + &= {-4t \over (1 + t^2)^2} + i {2(1 - t^2) \over (1 + t^2)^2} + \\ + &= {-2 \over 1 + t^2}s_1 + i {2 \over 1 + t^2}c_1 + \\ + &= -(1 + c_1)s_1 + i(1 + c_1)c_1 + \\ + &= i(1 + c_1)o_1 +\end{align*} +$$ + +Just like the complex exponential, an imaginary coefficient falls out. +However, the expression also accrues a $1 + c_1$ term, almost like an adjustment factor + for its failure to be the complex exponential. +Sine and cosine obey a simpler relationship with respect to the derivative, + and thus need no adjustment. + + +### Complex Analysis + +Since $o_n$ is a curve which loops around the unit circle *n* times, that possibly suits it + to showing a simple result from complex analysis. +Integrating along a contour which wraps around a sufficiently nice function's pole + (i.e., where its magnitude grows without bound) yields a familiar value. +This is easiest to see with $f(z) = 1 / z$: + +$$ +\oint_\Gamma {1 \over z} dz + = \int_a^b {\gamma'(t) \over \gamma(t)} dt + = 2\pi i +$$ + +In this example, Γ is a counterclockwise curve parametrized by γ which loops once around + the pole at *z* = 0. +More loops will scale this by a factor according to the number of loops. + +Normally this equality is demonstrated with the complex exponential, but will $o_1$ work just as well? +If Γ is the unit circle, the integral is: + +$$ +\oint_\Gamma {1 \over z} dz + = \int_{-\infty}^\infty {o_1'(t) \over o_1(t)} dt + = \int_{-\infty}^\infty i(1 + c_1(t)) dt + = 2i\int_{-\infty}^\infty {1 \over 1 + t^2} dt +$$ + +If one has studied their integral identities, the indefinite version of the final integral + will be obvious as $\arctan(t)$, which has horizontal asymptotes of $\pi / 2$ and $-\pi / 2$. +Therefore, the value of the integral is indeed $2\pi i$. + +If there are *n* loops, then naturally there are *n* of these $2\pi i$s. +Since powers of *o* are more loops around the circle, the chain and power rules show: + +$$ +\begin{gather*} + {d \over dt} (o_1)^n = n(o_1)^{n-1} {d \over dt} o_1 + \\[14pt] + \oint_\Gamma {1 \over z} dz + = \int_{-\infty}^\infty {n o_1(t)^{n-1} o_1'(t) \over o_1(t)^n} dt + = n \int_{-\infty}^\infty {o_1'(t) \over o_1(t)} dt + = 2 \pi i n +\end{gather*} +$$ + +It is certainly possible to perform these contour integrals along straight lines; + in fact, integrating along lines from 1 to *i* to -1 to -*i* deals with a + similar integral involving arctangent. +However, the best one can do to construct more loops with lines is to count each line + multiple times, which isn't extraordinarily convincing. + +Perhaps the use of $\infty$ in the integral bounds is also unconvincing. +The integral can be shifted back into the realm of plausibility by considering simpler bounds on $o_2$: + +$$ +\begin{align*} + \oint_\Gamma {1 \over z} dz + &= \int_{-1}^1 {2 o_1(t) o_1'(t) \over o_1(t)^2} dt + \\ + &= 2 \int_{-1}^1 {o_1'(t) \over o_1(t)} dt + \\ + &= 2(2i\arctan(1) - 2i\arctan(-1)) + \\ + &= 2\pi i +\end{align*} +$$ + +This has an additional benefit: using the series form of $1 / (1 + t^2)$ and integrating, + one obtains the series form of the arctangent. +This series converges for $-1 \le t \le 1$, which happens to match the bounds of integration. +The convergence of this series is fairly important, since it is tied to formulas for π, + in particular [Leibniz's formula](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80). + +Were one to integrate with the complex exponential, we would instead use the bounds $(0, 2\pi)$, + since at this point a full loop has been made. +But think to yourself -- how do you know the period of the complex exponential? +How do you know that 2π radians is equivalent to 0 radians? +The result using stereography relies on neither of these prior results and is directly pinned + to a formula for π instead an apparent detour through the number *e*. + + +Polar Curves +------------ + +Polar coordinates are useful for expressing for which the distance from the origin is + a function of the angle with respect to the positive *x*-axis. +They can also be readily converted to parametric forms: + +$$ +\begin{gather*} + r(\theta) &\Longleftrightarrow& + \begin{matrix} + x(\theta) = r \cos(\theta) \\ + y(\theta) = r \sin(\theta) + \end{matrix} +\end{gather*} +$$ + +Polar curves frequently loop in on themselves, and so it is necessary to choose appropriate bounds + for θ (usually as multiples of π) when plotting. +Evidently, this is due to the use of sine and cosine in the above parametrization. +Fortunately, $s_n$ and $c_n$ (as shown by the calculus above) have much simpler bounds. +So what happens when one substitutes the rational functions in place of the trig ones? + + +### Polar Roses + +[Polar roses](https://en.wikipedia.org/wiki/Rose_(mathematics)) are beautiful shapes which have + a simple form when expressed in polar coordinates. + +$$ +r(\theta) = \cos \left( {p \over q} \cdot \theta \right) +$$ + +The ratio $p/q$ in least terms uniquely determines the shape of the curve. + +If you weren't reading this post, you might assume this curve is transcendental since it uses cosine, + but you probably know better at this point. +The Chebyshev examples above demonstrate the resemblance between $c_n$ and $\cos(n\theta)$. +The subscript of $c$ is easiest to work with as an integer, so let $q = 1$. + +$$ +x(t) = c_p(t) c_1(t) \qquad y(t) = c_p(t) s_1(t) +$$ + +will plot a $p/1$ polar rose as t ranges over $(-\infty, \infty)$. + +```{python} +#| echo: false +#| output: false + +@dataclass +class Rose: + x: sympy.Expr + y: sympy.Expr + title: str | None = None + + @classmethod + def from_rational(cls, p: int, q: int): + return cls( + cs[p]*cs[q], + cs[p]*ss[q], + title=f"$x = c_{{{p}}}c_{{{q}}}; y = c_{{{p}}}s_{{{q}}}$", + ) + + +def animate_roses(filename: str, arguments: list[Rose], interval=500): + with SympyAnimationWrapper(filename) as animate: + @animate(list(range(len(arguments))), interval=interval) + def ret(fr): + plt.clf() + argument = arguments[fr] + + sympy.plot_parametric( + argument.x, + argument.y, + xlim=(-2,2), + ylim=(-2,2), + title=argument.title, + backend="matplotlib", + ) + + ret.save() # type: ignore + +animate_roses( + "polar_roses_1.mp4", + [ Rose.from_rational(i, 1) for i in range(1, 11) ], +) +``` + +::: {#fig-polar-roses-1} +{{< video "./polar_roses_1.mp4" >}} + +p/1 polar roses as rational curves. +Since *t* never reaches infinity, a bite appears to be taken out of the graphs near (-1, 0)." +::: + +$q = 1$ happens to match the subscript *c* term of *x* and *s* term of *y*, so one might wonder + whether the other polar curves can be obtained by allowing it to vary as well. +And you'd be right. + +$$ +x(t) = c_p(t) c_q(t) \qquad y(t) = c_p(t) s_q(t) +$$ + +will plot a $p/q$ polar rose as t ranges over $(-\infty, \infty)$. + +```{python} +#| echo: false +#| output: false + +animate_roses( + "polar_roses_2.mp4", + [ + Rose.from_rational(i,j) + for i,j in [(1,1),(1,2),(2,1),(3,1),(3,2),(2,3),(1,3),(1,4),(3,4),(4,3),(4,1)] + ], +) +``` + +::: {#fig-polar-roses-2} +{{< video "./polar_roses_2.mp4" >}} + +p/q polar roses as rational curves +::: + +Just as with the prior calculus examples, doubling all subscripts of *c* and *s* will + only require *t* to range over $(-1, 1)$, which removes the ugly bite mark. +Perhaps it is also slightly less satisfying, since the fraction $p/q$ directly appears in the + typical polar incarnation with cosine. +On the other hand, it exposes an important property of these curves: they are all rational. + +This approach lends additional precision to a prospective pseudo-polar coordinate system. +In the next few examples, I will be using the following notation for compactness: + +$$ + \begin{gather*} + R_n(t) = f(t) &\Longleftrightarrow& + \begin{matrix} + x(t) = f(t) c_n(t) \\ + y(t) = f(t) s_n(t) + \end{matrix} +\end{gather*} +$$ + + +### Conic Sections + +The polar equation for a conic section (with a particular unit length occurring somewhere) + in terms of its eccentricity $\varepsilon$ is: + +$$ +r(\theta) = {1 \over 1 - \varepsilon \cos(\theta)} +$$ + +Correspondingly, the rational polar form can be expressed as + +$$ +R_1(t) = {1 \over 1 - \varepsilon c_1} +$$ + +Since polynomial arithmetic is easier to work with than trigonometric identities, + it is a matter of pencil-and-paper algebra to recover the implicit form from a parametric one. + + +#### Parabola ($|\varepsilon| = 1$) + +The conic section with the simplest implicit equation is the parabola. +Since $c_n$ is a simple ratio of polynomials in *t*, it is much simpler to recover the implicit equation. +For $\varepsilon = 1$, + +:::: {layout-ncol="2"} +```{python} +#| echo: false +#| fig-cap: $x = {c_1 \over 1 - c_1} \quad y = {s_1 \over 1 - c_1}$ + +sympy.plot_parametric( + 1 / (1 - sympy.re(o)) * sympy.re(o), + 1 / (1 - sympy.re(o)) * sympy.im(o), + xlim=(-5,5), + ylim=(-5,5), +) +``` + +::: {} +$$ +\begin{align*} + 1 - c_1 &= 1 - {1 - t^2 \over 1 + t^2} + = {2 t^2 \over 1 + t^2} + \\ + y &= {s_1 \over 1 - c_1} + = {2t \over 1 + t^2} {1 + t^2 \over 2 t^2} + = {1 \over t} + \\ + x &= {c_1 \over 1 - c_1} + = {1 - t^2 \over 1 + t^2} \cdot {1 + t^2 \over 2 t^2} + = {1 - t^2 \over 2t^2} + \\ + &= {1 \over 2t^2} - {1 \over 2} + = {y^2 \over 2} - {1 \over 2} +\end{align*} +$$ +::: +:::: + +*x* is a quadratic polynomial in *y*, so trivially the figure formed is a parabola. +Technically it is missing the point where $y = 0 ~ (t = \infty)$, and this is not a circumstance + where using a higher $c_n$ would help. +It is however, similar to the situation where we allow $o_1(\infty) = -1$, and an argument + can be made to waive away any concerns one might have. + + +#### Ellipse ($|\varepsilon| < 1$) + +Ellipses are next. +The simplest fraction between zero and one is 1/2, so for $\varepsilon = 1/2$, + +:::: {layout-ncol = "2"} +```{python} +#| echo: false +#| fig-cap: $x = {c_1 \over 1 - c_1 / 2} \quad y = {s_1 \over 1 - c_1 / 2}$ + +sympy.plot_parametric( + 1 / (1 - (1/2)*sympy.re(o)) * sympy.re(o), + 1 / (1 - (1/2)*sympy.re(o)) * sympy.im(o), + xlim=(-5,5), + ylim=(-5,5), +) +``` + +::: {} +$$ +\begin{align*} + 1 - {1 \over 2}c_1 &= 1 - {1 \over 2} \cdot {1 - t^2 \over 1 + t^2} + = {3 t^2 + 1 \over 2 + 2t^2} + \\ + y &= {s_1 \over 1 - {1 \over 2}c_1} + = {4t \over 3t^2 + 1} + \\ + x &= {c_1 \over 1 - {1 \over 2}c_1} + = {2 - 2t^2 \over 3t^2 + 1} +\end{align*} +$$ +::: +:::: + +There isn't an obvious way to combine products of *x* and *y* into a single equation. +The general form of a conic section is $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$, so + we know that the implicit equation for the curve almost certainly involves $x^2$ and $y^2$. + +$$ +x^2 = {4 - 8t^2 + 4t^4 \over (3t^2 + 1)^2} \qquad + y^2 = {16t^2 \over (3t^2 + 1)^2} +$$ + +Squaring produces some $t^4$ terms which cannot exist outside of these terms and *xy*. +A linear combination of $x^2$ and $y^2$ never includes any cubic terms in the numerator + which would appear in *xy*, so $B = 0$. +Since all remaining terms are linear in *x* and *y*, their denominator must appear as a factor + in the numerator of $Ax^2 + Cy^2$, whatever *A* and *C* are. + +Since the coefficient of $t^4$ in $x^2$ is 4, *A* must be multiple of 3. +Through trial and error, $A = 3, C = 4$ gives: + +$$ +\begin{align*} + 3x^2 + 4y^2 + &= {12 - 24t^2 + 12t^4 + 64t^2 \over (3t^2 + 1)^2} + \\ + &= {12 - 40t^2 + 12t^4 \over (3t^2 + 1)^2} + \\ + &= {(4t^2 + 12) (3t^2 + 1) \over (3t^2 + 1)^2} + \\ + &= {4t^2 + 12 \over 3t^2 + 1} +\end{align*} +$$ + +Since the numerator of *y* has a *t*, this is clearly some combination of *x* and a constant. +By the previous line of thought, the constant term must be a multiple of 4, and picking the smallest + option finally results in the implicit form: + +$$ +\begin{align*} + 4 &= {4(3t^2 + 1) \over 3t^2 + 1} + = {12t^2 + 4 \over 3t^2 + 1} + \\ + {4t^2 + 12 \over 3t^2 + 1} - 4 + &= {8 - 8t^2 \over 3t^2 + 1} = 4x + \\[14pt] + 3x^2 + 4y^2 &= 4x + 4 + \\ + 3x^2 + 4y^2 - 4x - 4 &= 0 +\end{align*} +$$ + +Notably, the coefficients of *x* and *y* are 3 and 4. +Simultaneously, $o_1(\varepsilon) = o_1(1/2) = {3 \over 5} + i{4 \over 5}$. +This binds together three concepts: the simplest case of the Pythagorean theorem, + the 3-4-5 right triangle; the coefficients of the implicit form; and the role of eccentricity + with respect to stereography. + + +#### Hyperbola ($|\varepsilon| > 1$) + +As evidenced by the bound on the eccentricity above, hyperbolae are in some way the inverses of ellipses. +Since $o_1(2)$ is a reflection of $o_1(1/2)$, you might think the implicit equation for + $\varepsilon = 2$ to be the same, but with a flipped sign or two. +Unfortunately, you'd be wrong. + +:::: {layout-ncol="2"} +```{python} +#| echo: false +#| fig-cap: $x = {c_1 \over 1 - 2c_1} \quad y = {s_1 \over 1 - 2c_1}$ + +# Ignore discontinuity at sqrt(1/3) +disc = (1/3)**(1/2) +sympy.plot_parametric( + ( + 1 / (1 - 2*sympy.re(o)) * sympy.re(o), + 1 / (1 - 2*sympy.re(o)) * sympy.im(o), + (t, -disc + 1e-5, disc - 1e-5) + ), + ( + 1 / (1 - 2*sympy.re(o)) * sympy.re(o), + 1 / (1 - 2*sympy.re(o)) * sympy.im(o), + (t, -5, -disc - 1e-5) + ), + ( + 1 / (1 - 2*sympy.re(o)) * sympy.re(o), + 1 / (1 - 2*sympy.re(o)) * sympy.im(o), + (t, disc + 1e-5, 5) + ), + xlim=(-5,5), + ylim=(-5,5), + line_color = "blue", +) +``` + +::: {} +$$ +\begin{gather*} + \begin{align*} + x &= {c_1 \over 1 - 2c_1} = {1 - t^2 \over 3t^2 - 1} + \\ + y &= {s_1 \over 1 - 2c_1} = {2t \over 3t^2 - 1} + \\[14pt] + 3x^2 - y^2 &= {3 - 6t^2 + 3t^4 - 4t^2 \over 3t^2 - 1} + \\ + &= {(t^2 - 3)(3t^2 - 1) \over (3t^2 - 1)^2 } + \\ + &= {t^2 - 3 \over 3t^2 - 1 } + = ... = -4x - 1 + \end{align*} + \\[14pt] + 3x^2 - y^2 + 4x + 1 = 0 +\end{gather*} +$$ +::: +:::: + +At the very least, the occurrences of 1 in the place of 4 have a simple explanation: 1 = 4 - 3. + + +### Archimedean Spiral + +Arguably the simplest (non-circular) polar curve is $r(\theta) = \theta$, the unit + [Archimedean spiral](https://en.wikipedia.org/wiki/Archimedean_spiral). +Since the curve is defined by a constant turning, this is a natural application of the properties + of sine and cosine. +The closest equivalent in rational polar coordinates is $R_1(t) = t$. +But this can be converted to an implicit form: + +$$ +\begin{gather*} + x = tc_1 \qquad y = ts_1 + \\[14pt] + x^2 + y^2 = t^2(c_1^2 + s_1^2) = t^2 + \\ + y = {2t^2 \over 1 + t^2} = {2(x^2 + y^2) \over 1 + (x^2 + y^2)} + \\[14pt] + (1 + x^2 + y^2)y = 2(x^2 + y^2) +\end{gather*} +$$ + +The curve produced by this equation is a + [right strophoid](https://mathworld.wolfram.com/RightStrophoid.html) + with a node at (0, 1) and asymptote $y = 2$. +This form suggests something interesting about this curve: it approximates the Archimedean spiral + (specifically the one with polar equation $r(\theta) = \theta/2$). +Indeed, the sequence of curves with parametrization $R_n(t) = 2nt$ approximate the (unit) spiral + for larger *n*, as can be seen in the following video. + +```{python} +#| echo: false +#| output: false + +with SympyAnimationWrapper("approximate_archimedes.mp4") as animate: + @animate(list(range(10)), interval=500) + def ret(fr): + plt.clf() + + p = sympy.plot_parametric( + t*sympy.cos(t), + t*sympy.sin(t), + xlim=(-4,4), + ylim=(-4,4), + label="Archimedean Spiral", + backend="matplotlib", + show=False + ) + i = fr + 1 + p.extend( + sympy.plot_parametric( + 2*i*t*cs[i], + 2*i*t*ss[i], + (t, -5, 5), + line_color="black", + label=f"$R_{{{i}}}(t) = {2*i}t$", + backend="matplotlib", + show=False + ) + ) + p.show() + plt.legend() + + ret.save() # type: ignore +``` + +::: {#fig-approx-archimedes} +{{< video ./approximate_archimedes.mp4 >}} + +Approximations to the Archimedean spiral +::: + + +Since R necessarily defines a rational curve, the curves will never be equal, + just as any stretching of $c_n$ will never exactly become cosine. + + +Closing +------- + +Sine, cosine, and the exponential function, are useful in a calculus setting precisely + because of their constant "velocity" around the circle. +Also, nearly every modern scientific calculator in the world features buttons + for trigonometric functions, so there seems to be no reason *not* to use them. + +We can however be misled by their apparent omnipresence. +Stereographic projection has been around for *millennia*, and not every formula needs to be rewritten + in its language. +For example (and as previously mentioned), defining the Chebyshev polynomials really only requires + understanding the multiplication of two complex numbers whose norm cannot grow, + not trigonometry and dividing angles. +Many other instances of sine and cosine merely rely on a number (or ratio) of loops around a circle. +When velocity does not factor, it will obviously do to "stay rational". + +One of my favorite things to plot as a kid were polar roses, so I was somewhat intrigued + to see that they are, in fact, rational curves. +On the other hand, their rationality follows immediately from the rationality of the circle + (which itself follows from the existence of Pythagorean triples). +If I were more experienced with manipulating Chebyshev polynomials or willing to set up a + linear system in (way too) many terms, I might have considered attempting to find + an implicit form for them as well. + +Diagrams created with Sympy and Matplotlib. diff --git a/posts/stereo/_metadata.yml b/posts/stereo/_metadata.yml new file mode 100644 index 0000000..eaf1857 --- /dev/null +++ b/posts/stereo/_metadata.yml @@ -0,0 +1,5 @@ +# freeze computational output +freeze: auto + +# Enable banner style title blocks +title-block-banner: true diff --git a/posts/stereo/index.qmd b/posts/stereo/index.qmd new file mode 100644 index 0000000..82c3696 --- /dev/null +++ b/posts/stereo/index.qmd @@ -0,0 +1,8 @@ +--- +title: "Algebraic Stereography" +listing: + contents: . + sort: "date" +--- + +Articles about the use of rational circle functions and some applications.