Add thiele app

README.txt

@@ -31,6 +31,10 @@ primes - display primes as rod numerals (hommage to Rune Mields)
 colorcode - Show a color coded resistor
     parameter is not used yet
 
+thiele - Thieles talmønstre, quadratic residues of Gaussian integers modulo
+    parameter is the number of seconds between pattern changes
+
+
 Compile with "cargo build --release".
 
 The binary executables are in the "target/release" subdirectory.

src/bin/thiele/gaussian_integers.rs

@@ -0,0 +1,440 @@
+#![allow(clippy::suspicious_arithmetic_impl)]
+
+use std::ops;
+
+use primal::is_prime;
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
+pub struct Gaussian(pub i64, pub i64);
+impl Gaussian {
+    pub const ZERO: Gaussian = Gaussian(0, 0);
+    pub const ONE: Gaussian = Gaussian(1, 0);
+    pub const I: Gaussian = Gaussian(0, 1);
+
+    pub fn conj(&self) -> Gaussian {
+        Gaussian(self.0, -self.1)
+    }
+    pub fn norm(&self) -> i64 {
+        (self * self.conj()).0
+    }
+    pub fn re(&self) -> i64 {
+        self.0
+    }
+    pub fn im(&self) -> i64 {
+        self.1
+    }
+    pub fn is_zero(&self) -> bool {
+        *self == Self::ZERO
+    }
+    pub fn is_unit(&self) -> bool {
+        self.norm() == 1
+    }
+    /// Returns the factorization of a Gaussian integer
+    /// into a unit and its associate in the first quadrant.
+    pub fn factor_unit(&self) -> (Gaussian, Gaussian) {
+        let unit = {
+            if self.re() > 0 && self.im() >= 0 {
+                Gaussian::ONE
+            } else if self.im() > 0 && self.re() <= 0 {
+                Gaussian::I
+            } else if self.re() < 0 && self.im() <= 0 {
+                -Gaussian::ONE
+            } else if self.im() < 0 && self.re() >= 0 {
+                -Gaussian::I
+            } else {
+                // Must be zero.
+                Gaussian::ONE
+            }
+        };
+        (unit, unit.conj() * self)
+    }
+    pub fn is_prime(&self) -> bool {
+        let (_, a) = self.factor_unit();
+        if a.im() == 0 {
+            is_prime(a.re() as u64) && a.re() % 4 == 3
+        } else {
+            is_prime(a.norm() as u64)
+        }
+    }
+}
+
+impl From<i64> for Gaussian {
+    fn from(value: i64) -> Gaussian {
+        Gaussian(value, 0)
+    }
+}
+
+impl ops::Neg for Gaussian {
+    type Output = Gaussian;
+
+    fn neg(self) -> Gaussian {
+        Gaussian(-self.0, -self.1)
+    }
+}
+impl ops::Neg for &Gaussian {
+    type Output = Gaussian;
+
+    fn neg(self) -> Gaussian {
+        Gaussian(-self.0, -self.1)
+    }
+}
+
+impl ops::Add<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn add(self, rhs: Gaussian) -> Gaussian {
+        Gaussian(self.0 + rhs.0, self.1 + rhs.1)
+    }
+}
+impl ops::Add<Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn add(self, rhs: Gaussian) -> Gaussian {
+        Gaussian(self.0 + rhs.0, self.1 + rhs.1)
+    }
+}
+impl ops::Add<&Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn add(self, rhs: &Gaussian) -> Gaussian {
+        Gaussian(self.0 + rhs.0, self.1 + rhs.1)
+    }
+}
+impl ops::Add<&Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn add(self, rhs: &Gaussian) -> Gaussian {
+        Gaussian(self.0 + rhs.0, self.1 + rhs.1)
+    }
+}
+
+impl ops::Sub<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn sub(self, rhs: Gaussian) -> Gaussian {
+        Gaussian(self.0 - rhs.0, self.1 - rhs.1)
+    }
+}
+impl ops::Sub<Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn sub(self, rhs: Gaussian) -> Gaussian {
+        Gaussian(self.0 - rhs.0, self.1 - rhs.1)
+    }
+}
+impl ops::Sub<&Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn sub(self, rhs: &Gaussian) -> Gaussian {
+        Gaussian(self.0 - rhs.0, self.1 - rhs.1)
+    }
+}
+impl ops::Sub<&Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn sub(self, rhs: &Gaussian) -> Gaussian {
+        Gaussian(self.0 - rhs.0, self.1 - rhs.1)
+    }
+}
+
+impl ops::Mul<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn mul(self, rhs: Gaussian) -> Gaussian {
+        Gaussian(
+            self.0 * rhs.0 - self.1 * rhs.1,
+            self.0 * rhs.1 + self.1 * rhs.0,
+        )
+    }
+}
+impl ops::Mul<Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn mul(self, rhs: Gaussian) -> Gaussian {
+        Gaussian(
+            self.0 * rhs.0 - self.1 * rhs.1,
+            self.0 * rhs.1 + self.1 * rhs.0,
+        )
+    }
+}
+impl ops::Mul<&Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn mul(self, rhs: &Gaussian) -> Gaussian {
+        Gaussian(
+            self.0 * rhs.0 - self.1 * rhs.1,
+            self.0 * rhs.1 + self.1 * rhs.0,
+        )
+    }
+}
+impl ops::Mul<&Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn mul(self, rhs: &Gaussian) -> Gaussian {
+        Gaussian(
+            self.0 * rhs.0 - self.1 * rhs.1,
+            self.0 * rhs.1 + self.1 * rhs.0,
+        )
+    }
+}
+
+fn div_round(a: i64, b: i64) -> i64 {
+    (a + b / 2).div_euclid(b)
+}
+
+impl ops::Div<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn div(self, rhs: Gaussian) -> Gaussian {
+        let a = self * rhs.conj();
+        let n = rhs.norm();
+        Gaussian(div_round(a.re(), n), div_round(a.im(), n))
+    }
+}
+impl ops::Div<Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn div(self, rhs: Gaussian) -> Gaussian {
+        let a = self * rhs.conj();
+        let n = rhs.norm();
+        Gaussian(div_round(a.re(), n), div_round(a.im(), n))
+    }
+}
+impl ops::Div<&Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn div(self, rhs: &Gaussian) -> Gaussian {
+        let a = self * rhs.conj();
+        let n = rhs.norm();
+        Gaussian(div_round(a.re(), n), div_round(a.im(), n))
+    }
+}
+impl ops::Div<&Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn div(self, rhs: &Gaussian) -> Gaussian {
+        let a = self * rhs.conj();
+        let n = rhs.norm();
+        Gaussian(div_round(a.re(), n), div_round(a.im(), n))
+    }
+}
+
+impl ops::Rem<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn rem(self, rhs: Gaussian) -> Gaussian {
+        let q = self / rhs;
+        self - q * rhs
+    }
+}
+impl ops::Rem<Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn rem(self, rhs: Gaussian) -> Gaussian {
+        let q = self / rhs;
+        self - q * rhs
+    }
+}
+impl ops::Rem<&Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn rem(self, rhs: &Gaussian) -> Gaussian {
+        let q = self / rhs;
+        self - q * rhs
+    }
+}
+impl ops::Rem<&Gaussian> for &Gaussian {
+    type Output = Gaussian;
+
+    fn rem(self, rhs: &Gaussian) -> Gaussian {
+        let q = self / rhs;
+        self - q * rhs
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_conj() {
+        let a = Gaussian(1, 2);
+        assert_eq!(Gaussian(1, -2), a.conj());
+    }
+
+    #[test]
+    fn test_norm() {
+        let a = Gaussian(1, 2);
+        assert_eq!(5, a.norm());
+    }
+
+    #[test]
+    fn test_re() {
+        let a = Gaussian(1, 2);
+        assert_eq!(1, a.re());
+    }
+
+    #[test]
+    fn test_im() {
+        let a = Gaussian(1, 2);
+        assert_eq!(2, a.im());
+    }
+
+    #[test]
+    fn test_is_zero() {
+        assert!(Gaussian(0, 0).is_zero());
+        assert!(!Gaussian(1, 0).is_zero());
+        assert!(!Gaussian(0, 1).is_zero());
+        assert!(!Gaussian(1, 1).is_zero());
+        assert!(!Gaussian(1, -1).is_zero());
+        assert!(!Gaussian(-1, 1).is_zero());
+        assert!(!Gaussian(-1, -1).is_zero());
+    }
+
+    #[test]
+    fn test_is_unit() {
+        assert!(Gaussian(1, 0).is_unit());
+        assert!(Gaussian(-1, 0).is_unit());
+        assert!(Gaussian(0, 1).is_unit());
+        assert!(Gaussian(0, -1).is_unit());
+        assert!(!Gaussian(0, 0).is_unit());
+        assert!(!Gaussian(1, 1).is_unit());
+    }
+
+    #[test]
+    fn test_factor_unit() {
+        assert_eq!(
+            (Gaussian::ONE, Gaussian(3, 0)),
+            Gaussian(3, 0).factor_unit()
+        );
+        assert_eq!((Gaussian::I, Gaussian(3, 0)), Gaussian(0, 3).factor_unit());
+        assert_eq!(
+            (-Gaussian::ONE, Gaussian(3, 0)),
+            Gaussian(-3, 0).factor_unit()
+        );
+        assert_eq!(
+            (-Gaussian::I, Gaussian(3, 0)),
+            Gaussian(0, -3).factor_unit()
+        );
+    }
+
+    #[test]
+    fn test_is_prime() {
+        assert!(!Gaussian(0, 0).is_prime());
+
+        assert!(!Gaussian(1, 0).is_prime());
+        assert!(!Gaussian(0, 1).is_prime());
+        assert!(!Gaussian(-1, 0).is_prime());
+        assert!(!Gaussian(0, -1).is_prime());
+
+        assert!(Gaussian(1, 1).is_prime());
+        assert!(Gaussian(1, -1).is_prime());
+        assert!(Gaussian(-1, 1).is_prime());
+        assert!(Gaussian(-1, -1).is_prime());
+
+        assert!(!Gaussian(2, 0).is_prime());
+        assert!(!Gaussian(0, 2).is_prime());
+        assert!(!Gaussian(-2, 0).is_prime());
+        assert!(!Gaussian(0, -2).is_prime());
+
+        assert!(Gaussian(3, 0).is_prime());
+        assert!(Gaussian(0, 3).is_prime());
+        assert!(Gaussian(-3, 0).is_prime());
+        assert!(Gaussian(0, -3).is_prime());
+
+        assert!(!Gaussian(5, 0).is_prime());
+        assert!(Gaussian(7, 0).is_prime());
+        assert!(Gaussian(11, 0).is_prime());
+        assert!(!Gaussian(13, 0).is_prime());
+        assert!(!Gaussian(17, 0).is_prime());
+        assert!(Gaussian(19, 0).is_prime());
+        assert!(Gaussian(23, 0).is_prime());
+
+        assert!(Gaussian(2, 1).is_prime());
+        assert!(Gaussian(2, -1).is_prime());
+        assert!(Gaussian(-2, 1).is_prime());
+        assert!(Gaussian(-2, -1).is_prime());
+        assert!(Gaussian(1, 2).is_prime());
+        assert!(Gaussian(1, -2).is_prime());
+        assert!(Gaussian(-1, 2).is_prime());
+        assert!(Gaussian(-1, -2).is_prime());
+    }
+
+    #[test]
+    fn test_add() {
+        let a = Gaussian(1, 2);
+        let b = Gaussian(3, 5);
+        assert_eq!(Gaussian(4, 7), a + b);
+        assert_eq!(Gaussian(4, 7), a + &b);
+        assert_eq!(Gaussian(4, 7), &a + b);
+        assert_eq!(Gaussian(4, 7), &a + &b);
+    }
+
+    #[test]
+    fn test_sub() {
+        let a = Gaussian(1, 2);
+        let b = Gaussian(3, 5);
+        assert_eq!(Gaussian(-2, -3), a - b);
+        assert_eq!(Gaussian(-2, -3), a - &b);
+        assert_eq!(Gaussian(-2, -3), &a - b);
+        assert_eq!(Gaussian(-2, -3), &a - &b);
+    }
+
+    #[test]
+    fn test_mul() {
+        let a = Gaussian(1, 2);
+        let b = Gaussian(3, 5);
+        assert_eq!(Gaussian(-7, 11), a * b);
+        assert_eq!(Gaussian(-7, 11), a * &b);
+        assert_eq!(Gaussian(-7, 11), &a * b);
+        assert_eq!(Gaussian(-7, 11), &a * &b);
+    }
+
+    #[test]
+    fn test_div_rem() {
+        for a_re in -20..=20 {
+            for a_im in -20..=20 {
+                let a = Gaussian(a_re, a_im);
+                for b_re in -20..=20 {
+                    for b_im in -20..=20 {
+                        let b = Gaussian(b_re, b_im);
+                        if !b.is_zero() {
+                            let q = a / b;
+                            let r = a % b;
+                            assert_eq!(a, q * b + r);
+                            assert!(r.norm() <= b.norm() / 2);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn test_div() {
+        let a = Gaussian(1, 2);
+        let b = Gaussian(3, 5);
+        assert_eq!(Gaussian(0, 0), a / b);
+        assert_eq!(Gaussian(0, 0), a / &b);
+        assert_eq!(Gaussian(0, 0), &a / b);
+        assert_eq!(Gaussian(0, 0), &a / &b);
+    }
+
+    #[test]
+    fn test_rem() {
+        let a = Gaussian(1, 2);
+        let b = Gaussian(3, 5);
+        assert_eq!(Gaussian(1, 2), a % b);
+        assert_eq!(Gaussian(1, 2), a % &b);
+        assert_eq!(Gaussian(1, 2), &a % b);
+        assert_eq!(Gaussian(1, 2), &a % &b);
+    }
+
+    #[test]
+    fn test_div_round() {
+        assert_eq!(0, div_round(0, 3));
+        assert_eq!(0, div_round(1, 3));
+        assert_eq!(1, div_round(2, 3));
+        assert_eq!(1, div_round(3, 3));
+    }
+}

src/bin/thiele/gaussian_modulo.rs

@@ -0,0 +1,84 @@
+#![allow(unused)]
+
+use crate::gaussian_integers::Gaussian;
+
+pub fn mod_from(m: Gaussian, a: Gaussian) -> Gaussian {
+    a % m
+}
+
+pub fn mod_add(m: Gaussian, a: Gaussian, b: Gaussian) -> Gaussian {
+    mod_from(m, a + b)
+}
+
+pub fn mod_sub(m: Gaussian, a: Gaussian, b: Gaussian) -> Gaussian {
+    mod_from(m, a - b)
+}
+
+pub fn mod_mul(m: Gaussian, a: Gaussian, b: Gaussian) -> Gaussian {
+    mod_from(m, a * b)
+}
+
+pub fn mod_square(m: Gaussian, a: Gaussian) -> Gaussian {
+    mod_from(m, a * a)
+}
+
+#[cfg(test)]
+mod tests {
+    use std::collections::HashSet;
+
+    use crate::gaussian_integers::Gaussian;
+
+    use super::*;
+
+    #[test]
+    fn test_mod_from() {
+        for m in [
+            Gaussian(1, 1),
+            Gaussian(3, 2),
+            Gaussian(3, -2),
+            Gaussian(-3, 2),
+            Gaussian(-3, -2),
+        ] {
+            let mut representatives = HashSet::new();
+            for re in -20..=20 {
+                for im in -20..=20 {
+                    representatives.insert(mod_from(m, Gaussian(re, im)));
+                }
+            }
+            assert_eq!(m.norm() as usize, representatives.len());
+        }
+    }
+
+    #[test]
+    fn test_mod_add() {
+        let m = Gaussian(3, 2);
+        assert_eq!(
+            mod_from(m, Gaussian(4, 7)),
+            mod_add(m, Gaussian(1, 2), Gaussian(3, 5))
+        )
+    }
+
+    #[test]
+    fn test_mod_sub() {
+        let m = Gaussian(3, 2);
+        assert_eq!(
+            mod_from(m, Gaussian(-2, -3)),
+            mod_sub(m, Gaussian(1, 2), Gaussian(3, 5))
+        )
+    }
+
+    #[test]
+    fn test_mod_mul() {
+        let m = Gaussian(3, 2);
+        assert_eq!(
+            mod_from(m, Gaussian(-7, 11)),
+            mod_mul(m, Gaussian(1, 2), Gaussian(3, 5))
+        )
+    }
+
+    #[test]
+    fn test_mod_square() {
+        let m = Gaussian(3, 2);
+        assert_eq!(mod_from(m, Gaussian(-3, 4)), mod_square(m, Gaussian(1, 2)))
+    }
+}

src/bin/thiele/main.rs

@@ -0,0 +1,145 @@
+use std::collections::HashSet;
+use std::env::args;
+use std::io::stdout;
+use std::io::Write;
+use std::thread::sleep;
+use std::time::Duration;
+
+use rand::rng;
+use rand::seq::IndexedRandom;
+use rand::Rng;
+
+use lowdim::bb2d;
+use lowdim::p2d;
+use lowdim::Array2d;
+use lowdim::BBox2d;
+use lowdim::Point2d;
+
+use pixelfoo_apps::color::Color;
+
+mod gaussian_integers;
+use gaussian_integers::Gaussian;
+
+mod gaussian_modulo;
+use gaussian_modulo::mod_from;
+use gaussian_modulo::mod_square;
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+enum Square {
+    Zero,
+    Residue,
+    NonResidue,
+}
+
+#[derive(Clone, Debug)]
+struct Board {
+    map: Array2d<i64, Square>,
+}
+impl Board {
+    pub fn with(bbox: BBox2d, f: impl FnMut(Point2d) -> Square) -> Board {
+        let map = Array2d::with(bbox, f);
+        Board { map }
+    }
+    pub fn bbox(&self) -> BBox2d {
+        self.map.bbox()
+    }
+}
+
+fn send<T: Write>(w: &mut T, board: &Board) -> std::io::Result<()> {
+    for y in board.bbox().y_range() {
+        for x in board.bbox().x_range() {
+            let square = board.map[p2d(x, y)];
+            let c = match square {
+                Square::Zero => Color::blue(),
+                Square::Residue => Color::new(0xd2, 0xd4, 0xbc),
+                Square::NonResidue => Color::black(),
+            };
+            w.write_all(&c.rgb())?;
+        }
+    }
+    w.flush()
+}
+
+const DEFAULT_ARG: i64 = 60;
+
+fn main() -> std::io::Result<()> {
+    let args = args().collect::<Vec<_>>();
+    eprintln!("executing {}", args[0]);
+
+    let x_size = args[1].parse::<i64>().unwrap();
+    let y_size = args[2].parse::<i64>().unwrap();
+    let arg = if let Some(s) = args.get(3) {
+        s.parse::<i64>().unwrap_or(DEFAULT_ARG)
+    } else {
+        DEFAULT_ARG
+    };
+    eprintln!("screen size {}x{}, arg {}", x_size, y_size, arg);
+
+    let mut rng = rng();
+    let bbox = bb2d(0..x_size, 0..y_size);
+
+    let delay = Duration::from_millis(200);
+    let frame_seconds = if arg > 0 { arg } else { DEFAULT_ARG };
+    let max_time_count = frame_seconds * 5;
+
+    let mut board = Board::with(bbox, |_p| Square::NonResidue);
+
+    let mut time_count = 0;
+    loop {
+        if time_count <= 0 {
+            time_count = max_time_count;
+
+            // Pick a random gaussian prime somewhat smaller than the screen size.
+            let mut m = Gaussian::ZERO;
+            while !m.is_prime() || m.norm() < 16 {
+                let min_size = x_size.min(y_size);
+                let range = 0..=(min_size / 3).max(8);
+                let re = rng.random_range(range.clone());
+                let im = rng.random_range(range);
+                m = Gaussian(re, im);
+            }
+
+            let mut reps = HashSet::new();
+            let limit = m.re().max(m.im());
+            for re in -limit..=limit {
+                for im in -limit..=limit {
+                    let a = mod_from(m, Gaussian(re, im));
+                    reps.insert(a);
+                }
+            }
+            let reps = reps.into_iter().collect::<Vec<_>>();
+
+            let mut residues = HashSet::new();
+            for &r in &reps {
+                let a = mod_square(m, r);
+                residues.insert(a);
+            }
+
+            // Pick a random offset of the board smaller than the prime.
+            let origin = reps.choose(&mut rng).unwrap();
+
+            eprintln!("chose prime modulus {m:?}, offset {origin:?}");
+
+            board = Board::with(bbox, |p| {
+                let re = p.x();
+                let im = p.y();
+                let a = mod_from(m, origin + Gaussian(re, im));
+                if a.is_zero() {
+                    Square::Zero
+                } else if residues.contains(&a) {
+                    Square::Residue
+                } else {
+                    Square::NonResidue
+                }
+            });
+        }
+
+        let mut buf = Vec::with_capacity((x_size * y_size * 3) as usize);
+        send(&mut buf, &board)?;
+        stdout().write_all(&buf)?;
+        stdout().flush()?;
+
+        sleep(delay);
+        time_count -= 1;
+    }
+}