From 8bb076c7fb36ae347fe2561ab71fffc69f0bb941 Mon Sep 17 00:00:00 2001 From: Juergen Stuber Date: Fri, 16 Jan 2026 21:30:35 +0100 Subject: [PATCH] Add thiele app --- README.txt | 4 + src/bin/thiele/gaussian_integers.rs | 440 ++++++++++++++++++++++++++++ src/bin/thiele/gaussian_modulo.rs | 84 ++++++ src/bin/thiele/main.rs | 145 +++++++++ 4 files changed, 673 insertions(+) create mode 100644 src/bin/thiele/gaussian_integers.rs create mode 100644 src/bin/thiele/gaussian_modulo.rs create mode 100644 src/bin/thiele/main.rs diff --git a/README.txt b/README.txt index f333a3f..83d905d 100644 --- a/README.txt +++ b/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. diff --git a/src/bin/thiele/gaussian_integers.rs b/src/bin/thiele/gaussian_integers.rs new file mode 100644 index 0000000..b6c82ec --- /dev/null +++ b/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 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 for Gaussian { + type Output = Gaussian; + + fn add(self, rhs: Gaussian) -> Gaussian { + Gaussian(self.0 + rhs.0, self.1 + rhs.1) + } +} +impl ops::Add 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 for Gaussian { + type Output = Gaussian; + + fn sub(self, rhs: Gaussian) -> Gaussian { + Gaussian(self.0 - rhs.0, self.1 - rhs.1) + } +} +impl ops::Sub 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 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 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 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 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 for Gaussian { + type Output = Gaussian; + + fn rem(self, rhs: Gaussian) -> Gaussian { + let q = self / rhs; + self - q * rhs + } +} +impl ops::Rem 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)); + } +} diff --git a/src/bin/thiele/gaussian_modulo.rs b/src/bin/thiele/gaussian_modulo.rs new file mode 100644 index 0000000..fcb4764 --- /dev/null +++ b/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))) + } +} diff --git a/src/bin/thiele/main.rs b/src/bin/thiele/main.rs new file mode 100644 index 0000000..645bee8 --- /dev/null +++ b/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, +} +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(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::>(); + eprintln!("executing {}", args[0]); + + let x_size = args[1].parse::().unwrap(); + let y_size = args[2].parse::().unwrap(); + let arg = if let Some(s) = args.get(3) { + s.parse::().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::>(); + + 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; + } +}