Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.
37 36 35 34 33 32 31 38 17 16 15 14 13 30 39 18 5 4 3 12 29 40 19 6 1 2 11 28 41 20 7 8 9 10 27 42 21 22 23 24 25 26 43 44 45 46 47 48 49
It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.
If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?
問 58 「螺旋素数」
面白いことに, 奇平方数が右下の対角線上に出現する. もっと面白いことには, 対角線上の13個の数字のうち, 8個が素数である. ここで割合は8/13 ≈ 62%である.
渦巻きに新しい層を付け加えよう. すると辺の長さが9の渦巻きが出来る. 以下, この操作を繰り返していく. 対角線上の素数の割合が10%未満に落ちる最初の辺の長さを求めよ.
fn is_prime(n: u32) -> bool { if n < 2 { return false; } if n == 2 || n == 3 || n == 5 { return true; } for d in &[2, 3, 5] { if n % *d == 0 { return false; } } let side = (n as f32).sqrt() as u32; let mut d = 5; for &i in [2, 4].iter().cycle() { d += i; if d > side { break; } if n % d == 0 { return false; } } true } struct Spiral { n: u32, diagonal_area: u32, primes: u32, } impl Spiral { fn increment(&mut self) { self.n += 1; self.diagonal_area += 4; self.primes += self.prime_corners(); } fn prime_corners(&self) -> u32 { let square = (2 * self.n + 1) * (2 * self.n + 1); (1u32..=3) .map(|i| square - 2 * self.n * i) .filter(|&n| is_prime(n)) .count() as u32 } } fn main() { let mut spiral = Spiral { n: 0, diagonal_area: 1, primes: 0, }; while { spiral.increment(); spiral.primes * 10 >= spiral.diagonal_area } {} let side = 2 * spiral.n + 1; println!("{}", side); assert_eq!(side, 26241) }
The LCG is still good enough for simple tasks like Miller-Rabin primality test, or FreeCell deals. (Linear congruential generator, rosettacode.org)
use std::time::{SystemTime, UNIX_EPOCH}; mod index { pub struct Index { pub i: usize, f: u8, } impl Index { pub fn increment(&mut self) { self.i += 2 << self.f; self.f ^= 1; } pub fn new() -> Self { Self { i: 5, f: 0 } } } } mod sieve { pub struct Sieve { sieve: [bool; 100_001], index: super::index::Index, } impl Sieve { pub fn new() -> Self { let mut s = Self { sieve: [true; 100_001], index: super::index::Index::new(), }; s.init(); s } pub fn sieve_len(&self) -> u32 { self.sieve.len() as u32 } fn init(&mut self) { let side = ((self.sieve.len() - 1) as f32).sqrt() as usize; while self.index.i <= side { if self.sieve[self.index.i] { let p = self.index.i; (p * p..self.sieve.len()) .step_by(p) .for_each(|i| self.sieve[i] = false); } self.index.increment(); } } pub fn is_prime(&self, n: u32) -> bool { assert!(n < self.sieve.len() as u32); if n < 2 { return false; } if n % 2 == 0 { return n == 2; } if n % 3 == 0 { return n == 3; } return self.sieve[n as usize]; } } } mod rand { pub struct MinStdRand { state: u64, } impl MinStdRand { const M: u64 = 2147483647; const A: u64 = 48271; const MAX: u64 = 2147483646; pub fn new(seed: u32) -> Self { Self { state: seed as u64 } } pub fn next(&mut self, partition: u32) -> u32 { let p = partition as u64; assert!(p > 0 && p <= Self::MAX); self.state = Self::A * self.state % Self::M; loop { let n = self.state * p / Self::MAX; if n < p { return n as u32; } } } } } fn mod_pow(a: u32, exp: u32, m: u32) -> u32 { let (mut a, mut exp, m) = (a as u64, exp as u64, m as u64); if m == 1 { return 0; } if exp == 0 { return 1; } let mut result = 1; a %= m; loop { if exp % 2 == 1 { result = result * a % m; } exp >>= 1; if exp == 0 { break; } a = a * a % m; } result as u32 } /// finds the k*2^e form of given n fn coefficient_and_exponent_of_two(mut n: u32) -> (u32, u32) { let mut exp = 0u32; while n % 2 == 0 { n /= 2; exp += 1; } (n, exp) } fn is_probable_prime(n: u32, rand: &mut rand::MinStdRand) -> bool { if n == 1 { return false; } if n % 2 == 0 { return n == 2; } let (d, s) = coefficient_and_exponent_of_two(n - 1); 'next_trial: for _ in 0..3 { // 2 <= a < n let a = 2 + rand.next(n - 2); let mut x = mod_pow(a, d, n); if x == 1 || x == n - 1 { continue 'next_trial; } for _ in 1..s { x = (x as u64 * x as u64 % n as u64) as u32; if x == n - 1 { continue 'next_trial; } } return false; } true } struct Spiral { n: u32, diagonal_area: u32, primes: u32, sieve: sieve::Sieve, rand: rand::MinStdRand, } impl Spiral { fn increment(&mut self) { self.n += 1; self.diagonal_area += 4; self.primes += self.prime_corners(); } fn prime_corners(&mut self) -> u32 { let square = (2 * self.n + 1) * (2 * self.n + 1); let n = self.n.clone(); (1u32..=3) .map(|i| square - 2 * n * i) .filter(|&n| { if n < self.sieve.sieve_len() { self.sieve.is_prime(n) } else { is_probable_prime(n, &mut self.rand) } }) .count() as u32 } } fn main() { let seed = SystemTime::now() .duration_since(UNIX_EPOCH) .unwrap() .subsec_nanos(); let mut spiral = Spiral { n: 0, diagonal_area: 1, primes: 0, sieve: sieve::Sieve::new(), rand: rand::MinStdRand::new(seed), }; while { spiral.increment(); spiral.primes * 10 >= spiral.diagonal_area } {} let side = 2 * spiral.n + 1; println!("{}", side); assert_eq!(side, 26241) }