Merge branch 'algebraic-stereography' into rewritten

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.*

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

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

posts/stereo/1/approx-results/.gitignore

@@ -0,0 +1,2 @@
+compare_complex
+compare_math

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

posts/stereo/1/approx-results/main.c

@@ -0,0 +1,196 @@
+#include <complex.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <unistd.h>
+
+#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;
+}

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

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

posts/stereo/1/approx-results/stereo_complex.c

@@ -0,0 +1,27 @@
+#include <complex.h>
+#include <math.h>
+
+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);
+}

posts/stereo/1/approx-results/stereo_math.c

@@ -0,0 +1,27 @@
+#include <math.h>
+
+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);
+}

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
+---
+
+<style>
+.figure-img {
+  max-width: 512px;
+  object-fit: contain;
+  height: 100%;
+}
+</style>
+
+```{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.
+
+<!-- TODO: this is gross -->
+$$
+\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.
+
+<!-- TODO: better repository link -->
+Diagrams created with GeoGebra and Matplotlib.
+Repository with approximations (as well as GeoGebra files) available [here](https://github.com/queue-miscreant/approx-trig).

posts/stereo/2/.gitignore

@@ -0,0 +1 @@
+*.mp4

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

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
+---
+
+<style>
+.figure-img {
+  max-width: 512px;
+  object-fit: contain;
+  height: 100%;
+}
+</style>
+
+```{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.

posts/stereo/_metadata.yml

@@ -0,0 +1,5 @@
+# freeze computational output
+freeze: auto
+
+# Enable banner style title blocks
+title-block-banner: true

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.