Rust, 10667 bytes
Significantly faster than my last answer, and comfortably handles 4000000 on my machine. The idea is still the same as my previous answer, by first reducing everything to under 255 and then operating on bytes. I dropped nalgebra and platform intrinsics this time and added a few nightly features, mainly for portable SIMD types. Additionally, only a quarter of the grid is calculated now. Building this one with just rustc instead of cargo is fine but ensure -Ctarget-cpu=native is passed to rustc.
#![feature(portable_simd)]
#![feature(allocator_api)]
#![feature(array_chunks)]
use std::mem;
use std::ops::{Index, IndexMut};
use std::ptr::NonNull;
use std::{env::args, fmt::Debug, simd::*};
fn main() {
let arg = args()
.nth(1)
.expect("No input provided")
.trim()
.parse::<u32>()
.expect("Not parsable as an unsigned integer");
let mut sandpile = Mat::sandpile(arg);
sandpile.topple::<255>();
let mut sandpile = sandpile.map(|n| n as u8);
sandpile.topple::<3>();
sandpile.trim();
sandpile.unquarter();
sandpile.print();
}
fn max_size(n: u32) -> usize {
((n as f64).sqrt() / 2.7).ceil() as usize + 3
}
#[derive(Clone)]
struct Mat<T> {
data: Vec<T, AlignedAlloc<SIMD_ALIGN>>,
row_size: usize,
}
impl<T> Index<(usize, usize)> for Mat<T> {
type Output = T;
fn index(&self, (x, y): (usize, usize)) -> &T {
&self.data[x + y * self.row_size]
}
}
impl<T> IndexMut<(usize, usize)> for Mat<T> {
fn index_mut(&mut self, (x, y): (usize, usize)) -> &mut T {
&mut self.data[x + y * self.row_size]
}
}
impl<T> Mat<T> {
fn row(&self, y: usize) -> &[T] {
&self.data[y * self.row_size..(y + 1) * self.row_size]
}
fn row_mut(&mut self, y: usize) -> &mut [T] {
&mut self.data[y * self.row_size..(y + 1) * self.row_size]
}
fn row_count(&self) -> usize {
self.data.len() / self.row_size
}
fn row_windows(&self) -> RowWindows<'_, T> {
RowWindows {
data: &self.data,
row_size: self.row_size,
}
}
fn map<U: Default + Copy, F: Fn(T) -> U>(self, f: F) -> Mat<U> {
let mut buf = Vec::with_capacity_in(self.data.len(), AlignedAlloc);
buf.extend(self.data.into_iter().map(f));
let mut mat = Mat {
data: buf,
row_size: self.row_size,
};
mat.pad();
mat
}
}
macro_rules! make_topple {
($int:ty, $simd:ty, $lanes:literal) => {
fn topple<const MAX: $int>(&mut self) {
if !self.data.iter().any(|&x| x > MAX) {
return;
}
let mut buf = self.clone();
let four = <$simd>::splat(4);
loop {
let mut max = <$simd>::splat(0);
for ((up, center, down), buf) in self
.row_windows()
.zip(buf.data.chunks_exact_mut(self.row_size).skip(1))
{
for (((((cur, left), right), up), down), buf) in center[1..]
.array_chunks::<$lanes>()
.zip(center.array_chunks::<$lanes>())
.zip(center[2..].array_chunks::<$lanes>())
.zip(up[1..].array_chunks::<$lanes>())
.zip(down[1..].array_chunks::<$lanes>())
.zip(buf[1..].array_chunks_mut::<$lanes>())
{
let cur = <$simd>::from_array(*cur);
let left = <$simd>::from_array(*left);
let right = <$simd>::from_array(*right);
let up = <$simd>::from_array(*up);
let down = <$simd>::from_array(*down);
let sum = left / four + right / four + up / four + down / four + cur % four;
max = max.simd_max(sum);
*buf = sum.to_array();
}
}
self.row(2)[1..]
.array_chunks::<$lanes>()
.zip(buf.row_mut(1)[1..].array_chunks_mut::<$lanes>())
.for_each(|(cur, buf)| {
let cur = <$simd>::from_array(*cur);
let tmp = <$simd>::from_array(*buf);
let sum = tmp + cur / four;
max = max.simd_max(sum);
*buf = sum.to_array();
});
let col_max = (1..self.row_count() - 1)
.map(|y| {
buf[(1, y)] += self[(2, y)] / 4;
buf[(1, y)]
})
.max()
.unwrap();
mem::swap(self, &mut buf);
let max = max.reduce_max().max(col_max);
if max <= MAX {
break;
}
}
}
#[allow(dead_code)] // u32 version is unused
fn trim(&mut self) {
let mut num_rows = 0;
for row in self.data.chunks_exact(self.row_size).skip(1) {
if row.iter().all(|&x| x == 0) {
break;
}
num_rows += 1;
}
let mut buf = Vec::with_capacity_in(num_rows * num_rows, AlignedAlloc);
buf.extend(std::iter::repeat(0).take(num_rows * num_rows));
let mut new = Mat {
data: buf,
row_size: num_rows,
};
for x in 0..num_rows {
for y in 0..num_rows {
new[(x, y)] = self[(x + 1, y + 1)];
}
}
*self = new;
}
};
}
impl<T> Mat<T>
where
T: Debug + Default,
{
fn print(&mut self) {
self.data.truncate(self.row_size * self.row_size);
for n in self.data.iter() {
print!("{n:?}");
}
}
}
impl<T: Default + Copy> Mat<T> {
fn pad(&mut self) {
let round_mask = SIMD_ALIGN / mem::align_of::<T>() - 1;
let padding = SIMD_ALIGN / mem::align_of::<T>();
let size = self.row_size;
let aligned_size = if size & round_mask == 0 {
size
} else {
(size & !round_mask) + padding
};
let mut data = Vec::with_capacity_in(
aligned_size * (aligned_size + padding),
AlignedAlloc::<SIMD_ALIGN>,
);
data.extend(
std::iter::repeat_with::<T, _>(Default::default)
.take(aligned_size * (aligned_size + padding)),
);
let mut new = Mat {
data,
row_size: aligned_size + padding,
};
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[(x, y)] = self[(x, y)];
}
}
*self = new;
}
fn unquarter(&mut self) {
let mut new = Mat {
data: Vec::with_capacity_in(
self.row_size * self.row_count() * 4,
AlignedAlloc::<SIMD_ALIGN>,
),
row_size: self.row_size * 2 - 1,
};
new.data.extend(
std::iter::repeat_with::<T, _>(Default::default)
.take((self.row_size * 2 - 1).pow(2) - 1 + self.row_size * 2 - 1),
);
let size = new.row_size;
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[((size - 1) / 2 + x, (size - 1) / 2 + y)] = self[(x, y)];
}
}
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[(x, y)] = self[(self.row_size - 1 - x, self.row_count() - 1 - y)];
}
}
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[((size - 1) / 2 + x, y)] = self[(x, self.row_count() - 1 - y)];
}
}
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[(x, (size - 1) / 2 + y)] = self[(self.row_size - 1 - x, y)];
}
}
*self = new;
}
}
impl Mat<u32> {
fn sandpile(n: u32) -> Self {
let size = max_size(n);
let mut data = Vec::with_capacity_in(size * (size), AlignedAlloc::<SIMD_ALIGN>);
data.extend(std::iter::repeat(0).take(size * size));
let mut mat = Mat {
data,
row_size: size,
};
mat[(1, 1)] = n;
mat.pad();
mat
}
make_topple!(u32, u32x8, 8);
}
impl Mat<u8> {
make_topple!(u8, u8x32, 32);
}
pub struct RowWindows<'a, T> {
data: &'a [T],
row_size: usize,
}
impl<'a, T> Iterator for RowWindows<'a, T> {
type Item = (&'a [T], &'a [T], &'a [T]);
fn next(&mut self) -> Option<Self::Item> {
if self.data.len() < self.row_size * 3 {
return None;
}
let (first, rest) = self.data.split_at(self.row_size);
self.data = rest;
let (second, rest) = rest.split_at(self.row_size);
let (third, _rest) = rest.split_at(self.row_size);
Some((first, second, third))
}
}
impl<T: Debug> Debug for Mat<T> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
for i in 0..self.data.len() / self.row_size {
for j in 0..self.row_size {
write!(f, "{:?} ", self.data[i * self.row_size + j])?;
}
writeln!(f)?;
}
Ok(())
}
}
const SIMD_ALIGN: usize = mem::align_of::<u32x8>();
#[derive(Clone, Copy, Debug)]
struct AlignedAlloc<const ALIGN: usize>;
use std::alloc::{AllocError, Allocator, Global, Layout};
unsafe impl<const ALIGN: usize> Allocator for AlignedAlloc<ALIGN> {
fn allocate(&self, layout: Layout) -> Result<NonNull<[u8]>, AllocError> {
Global.allocate(layout.align_to(ALIGN).unwrap())
}
unsafe fn deallocate(&self, ptr: NonNull<u8>, layout: Layout) {
Global.deallocate(ptr, layout.align_to(ALIGN).unwrap())
}
fn allocate_zeroed(&self, layout: Layout) -> Result<NonNull<[u8]>, AllocError> {
Global.allocate_zeroed(layout.align_to(ALIGN).unwrap())
}
unsafe fn grow(
&self,
ptr: NonNull<u8>,
old_layout: Layout,
new_layout: Layout,
) -> Result<NonNull<[u8]>, AllocError> {
Global.grow(
ptr,
old_layout.align_to(ALIGN).unwrap(),
new_layout.align_to(ALIGN).unwrap(),
)
}
unsafe fn grow_zeroed(
&self,
ptr: NonNull<u8>,
old_layout: Layout,
new_layout: Layout,
) -> Result<NonNull<[u8]>, AllocError> {
Global.grow_zeroed(
ptr,
old_layout.align_to(ALIGN).unwrap(),
new_layout.align_to(ALIGN).unwrap(),
)
}
unsafe fn shrink(
&self,
ptr: NonNull<u8>,
old_layout: Layout,
new_layout: Layout,
) -> Result<NonNull<[u8]>, AllocError> {
Global.shrink(
ptr,
old_layout.align_to(ALIGN).unwrap(),
new_layout.align_to(ALIGN).unwrap(),
)
}
}
Playground
I tried parallelizing this with rayon but the code got lost somewhere in synchronization primitives and took roughly ten times longer. I might return to this later and see if I can do better.
This code is fully cross-platform, and in fact, I cross-compiled it for 32-bit arm at one point. Because this uses several nightly features it likely won't work on the latest nightly forever. If it doesn't, here's my current rustc --version:
rustc 1.66.0-nightly (57f097ea2 2022-10-01)
