Merge pull request #2 from NathanaelLane/master

Refactor checksum algorithm to process 4 bytes at a time

src/lib.rs

@@ -12,12 +12,14 @@
 #![doc(test(attr(deny(unused_imports, unused_must_use))))]
 #![cfg_attr(docsrs, feature(doc_cfg))]
 #![warn(missing_debug_implementations)]
+#![forbid(unsafe_code)]
 #![cfg_attr(not(feature = "std"), no_std)]
 
 #[cfg(not(feature = "std"))]
 extern crate core as std;
 
 use std::hash::Hasher;
+use std::ops::{AddAssign, MulAssign, RemAssign};
 
 #[cfg(feature = "std")]
 use std::io::{self, BufRead};
@@ -108,30 +110,92 @@ impl Adler32 {
  // b = n×D1 + (n−1)×D2 + (n−2)×D3 + ... + Dn + n*1 (mod 65521)
  // Where n is the number of bytes and Di is the i-th Byte. We need to change this to account
  // for the previous values of a and b, as well as treat every input Byte as being 255:
- // b_inc = n×255 + (n-1)×255 + ... + 255 + n*65521
+ // b_inc = n×255 + (n-1)×255 + ... + 255 + n*65520
  // Or in other words:
- // b_inc = n*65521 + n(n+1)/2*255
+ // b_inc = n*65520 + n(n+1)/2*255
  // The max chunk size is thus the largest value of n so that b_inc <= 2^32-65521.
- // 2^32-65521 = n*65521 + n(n+1)/2*255
+ // 2^32-65521 = n*65520 + n(n+1)/2*255
  // Plugging this into an equation solver since I can't math gives n = 5552.18..., so 5552.
+ //
+ // On top of the optimization outlined above, the algorithm can also be parallelized with a
+ // bit more work:
+ //
+ // Note that b is a linear combination of a vector of input bytes (D1, ..., Dn).
+ //
+ // If we fix some value k<N and rewrite indices 1, ..., N as
+ //
+ // 1_1, 1_2, ..., 1_k, 2_1, ..., 2_k, ..., (N/k)_k,
+ //
+ // then we can express a and b in terms of sums of smaller sequences kb and ka:
+ //
+ // ka(j) := D1_j + D2_j + ... + D(N/k)_j where j <= k
+ // kb(j) := (N/k)*D1_j + (N/k-1)*D2_j + ... + D(N/k)_j where j <= k
+ //
+ // a = ka(1) + ka(2) + ... + ka(k) + 1
+ // b = k*(kb(1) + kb(2) + ... + kb(k)) - 1*ka(2) - ... - (k-1)*ka(k) + N
+ //
+ // We use this insight to unroll the main loop and process k=4 bytes at a time.
+ // The resulting code is highly amenable to SIMD acceleration, although the immediate speedups
+ // stem from increased pipeline parallelism rather than auto-vectorization.
+ //
+ // This technique is described in-depth (here:)[https://software.intel.com/content/www/us/\
+ // en/develop/articles/fast-computation-of-fletcher-checksums.html]
 
  const MOD: u32 = 65521;
- const CHUNK_SIZE: usize = 5552;
+ const CHUNK_SIZE: usize = 5552 * 4;
 
  let mut a = u32::from(self.a);
  let mut b = u32::from(self.b);
- for chunk in bytes.chunks(CHUNK_SIZE) {
- for byte in chunk {
- let val = u32::from(*byte);
- a += val;
- b += a;
- }
+ let mut a_vec = U32X4([0; 4]);
+ let mut b_vec = a_vec;
+
+ let (bytes, remainder) = bytes.split_at(bytes.len() - bytes.len() % 4);
 
- a %= MOD;
+ // iterate over 4 bytes at a time
+ let chunk_iter = bytes.chunks_exact(CHUNK_SIZE);
+ let remainder_chunk = chunk_iter.remainder();
+ for chunk in chunk_iter {
+ for byte_vec in chunk.chunks_exact(4) {
+ let val = U32X4::from(byte_vec);
+ a_vec += val;
+ b_vec += a_vec;
+ }
+ b += CHUNK_SIZE as u32 * a;
+ a_vec %= MOD;
+ b_vec %= MOD;
  b %= MOD;
  }
- self.a = a as u16;
- self.b = b as u16;
+ // special-case the final chunk because it may be shorter than the rest
+ for byte_vec in remainder_chunk.chunks_exact(4) {
+ let val = U32X4::from(byte_vec);
+ a_vec += val;
+ b_vec += a_vec;
+ }
+ b += remainder_chunk.len() as u32 * a;
+ a_vec %= MOD;
+ b_vec %= MOD;
+ b %= MOD;
+
+ // combine the sub-sum results into the main sum
+ b_vec *= 4;
+ b_vec.0[1] += MOD - a_vec.0[1];
+ b_vec.0[2] += (MOD - a_vec.0[2]) * 2;
+ b_vec.0[3] += (MOD - a_vec.0[3]) * 3;
+ for &av in a_vec.0.iter() {
+ a += av;
+ }
+ for &bv in b_vec.0.iter() {
+ b += bv;
+ }
+
+ // iterate over the remaining few bytes in serial
+ for &byte in remainder.iter() {
+ a += u32::from(byte);
+ b += a;
+ }
+
+ self.a = (a % MOD) as u16;
+ self.b = (b % MOD) as u16;
  }
 }
 
@@ -160,6 +224,44 @@ pub fn adler32_slice(data: &[u8]) -> u32 {
  h.checksum()
 }
 
+#[derive(Copy, Clone)]
+struct U32X4([u32; 4]);
+
+impl U32X4 {
+ fn from(bytes: &[u8]) -> Self {
+ U32X4([
+ u32::from(bytes[0]),
+ u32::from(bytes[1]),
+ u32::from(bytes[2]),
+ u32::from(bytes[3]),
+ ])
+ }
+}
+
+impl AddAssign<Self> for U32X4 {
+ fn add_assign(&mut self, other: Self) {
+ for (s, o) in self.0.iter_mut().zip(other.0.iter()) {
+ *s += o;
+ }
+ }
+}
+
+impl RemAssign<u32> for U32X4 {
+ fn rem_assign(&mut self, quotient: u32) {
+ for s in self.0.iter_mut() {
+ *s %= quotient;
+ }
+ }
+}
+
+impl MulAssign<u32> for U32X4 {
+ fn mul_assign(&mut self, rhs: u32) {
+ for s in self.0.iter_mut() {
+ *s *= rhs;
+ }
+ }
+}
+
 /// Calculates the Adler-32 checksum of a `BufRead`'s contents.
 ///
 /// The passed `BufRead` implementor will be read until it reaches EOF.