Repository: SpectralSequences/sseq (language: python)

Issue:
Revisit the saving mechanism
As we push into higher stems, I think it might be useful to rethink the way we store resolutions. Right now the `Resolution` struct is written to a file as-is. While that has its advantages, I can think of a few downsides.

- Very sensitive to the data layout. Any minor change is susceptible to make all prior archives unreadable, unless we're careful.
- Impossible to inspect manually.
- No support for incremental saving. This is the big one. When we start testing multi-day/week runs it will be inevitable that some crashes will happen, and it would be great to not lose the quasi-inverses that we will have spent many CPU-hours calculating. I think it would help to fragment the archive into smaller files. Maybe quasi-inverses could be stored in their own files?
- No compression. As the matrices get bigger the storage space becomes an issue. I've calculated that computing into the 200+ stems will require on the order of 10TB of disk space at the very least. It would be helpful to have the ability of writing/reading compressed data, especially since the archives are empirically very compressible. (I just checked and the resolution of S_2 up to t = 100 can be compressed by 60% using gzip.)

I don't have a concrete suggestion for a replacement right now, but I think it's something we should consider before pushing computations further (see #30).


Discussion / hints:
Here are some miscellaneous ideas. I think the main question is what the
performance costs are.

  - Instead of saving the full quasi-inverses, how costly is it to just
    save the quasi-inverse on the generators.

  - Can we use something like sqlite to store the data? It would be
    pretty useful to have atomic writes in case we have to kill the
    process halfway through the computation.

  - Simply piping the output through gzip should be easy using flate2.

I'll add that the resolution up to the 140th stem can be compressed by 50% using gzip with default settings.

I think the IO costs are going to be minimal compared to the rest of the computation, and there's also always the option of dispatching the save to a dedicated thread.

I think sqlite might still be too cumbersome.  I spoke with a friend today and he said that geneticists, who use more than 1 PB of storage on the WSU grid, have tools that exclusively read in / write out `.tar.gz` tarballs. That might be worth considering. It would really prefer having dedicated files for the quasi-inverses, so that we can mmap them and use them as-is without deserializing first. In large stems a single quasi-inverse can occupy _several_ GB.

I hadn't thought about storing the quasi-inverses only on the generators. That would certainly cut down the storage requirements by a lot, but at the cost of making it harder to compute lifts down the line. I would say that more experimentation is needed.
It depends on exactly what you are doing. For example, if you want to compute a single Massey product in a high stem, you would most likely spend most of the time loading the resolution.
I see your point. Even if the quasi-inverses are individually compressed, a lot of them would have to be decompressed to compute the lifts. However I still think that to compute them in the first place compression is essential. Maybe have the output be a tarball, and enable compressed IO through a feature flag, something like "compressed"? Obviously we should also use the fastest/most parallelized compression algorithm available, since that should already save a great deal of storage.

Another option would be to have the in-memory representation of the matrices be already compressed, so that we still can mmap them without deserializing. That's pretty similar to your idea of only storing the images on module generators.
This is my plan of attack:

The resolution is stored in a directory, which we shall call `/` for the purposes of this discussion. This has three subdirectories: `/differentials/`, `/quasi-inverses/` and `/kernels/`. In each directory, the data of each bidegree is stored in a file called `s_t.ext` for some extension to be determined, and each file is compressed.

For the kernels directory, we delete the file once we no longer need it.

This allows us to save data incrementally, and also it is now possible to load a smaller resolution given a save data.

Saving the quasi-inverses separately from the differentials makes it easy to implement #29.

Each save file will be compressed. I've experimented with several compression algorithms, trying to compress the data for `s = 3 to 7` and `t = 140` with the Milnor basis. The results are as follows (the numbers following an algorithm represent the compression level):
```
algorithm: write_time read_time total size
raw: 292.87629ms 245.78584ms 795806013
snap: 697.168531ms 854.239285ms 415882110
zstd 1: 2.275886699s 1.026356308s 365239503
zstd 3: 4.799971583s 881.428894ms 359660237
zstd 21: 168.854727843s 1.253918986s 339565081
gzip 1: 5.901635364s 2.558595974s 387680655
gzip 9: 137.415523131s 2.902248096s 378109490
brotli 1: 4.165027028s 4.57096398s 368688292
brotli 8: 48.285620746s 7.624755351s 353941345
```

From this data I think zstd is a clear winner. Only snap is faster but we get a ~20% size improvement. During write time, we can use compression level 1, and then after computations are done we can re-compress the data to a higher compression level to save space. My main concern with zstd is that it uses C++ bindings so it cannot be used with wasm.

An alternative approach is to always write the data uncompressed but be able to read both uncompressed and compressed data. We can then compress separately when we are idle.

I'll probably start working on this tomorrow. Comments are welcome.

NB. Compressing each bidegree individually is necessarily less efficient for compression purposes, but the ability to work incrementally is too useful to forgo imo.

I agree with mostly everything. Right now the main bottleneck is the space usage, both in storage and in RAM, so it makes sense to sacrifice some loading time. It's also a good idea to be more modular so that we avoid problems when changing the saving scheme for a given struct.

However I think it is essential to be able to work with uncompressed as well as compressed data. In the near future it will be possible to load vectors and matrices directly from disk using `mmap` and zero-copy deserialization, so it is crucial that the data be stored uncompressed in that case.

For the compression algorithm, have you tried lz4 using something like [lz4_flex](https://crates.io/crates/lz4_flex)? They claim to be the fastest compression scheme in pure Rust.
lz4 data is as follows:
```
lz4: 1.277215386s 417.09343ms 398264945
```

This gives slightly smaller sizes than snap and significantly faster decompression, but compression is slower. I think compression/decompression speed are less important, and compression time is specially unimportant. So it may make sense to use this over snap.

mmap is very tricky to do safely, since users would be able to modify the underlying files while the program is running, which gives instant UB. If we are happy to ignore this possibility then this should be easy to implement with uncompressed data.

It should not be too difficult to support multiple data formats.

Are there particular computations you want to do where storing even a single quasi-inverse in memory is problematic? In these cases we might even be able to get away with reading the quasi-inverse line by line. We can save the quasi-inverse in the format

```
<1st basis vector of image>
<lift of 1st basis vector>
<2nd basis vector of image>
<lift of 2nd basis vector>
...
```
and compute the quasi-inverse on a set of vector in a single pass through the file. Then at any point in time we only have to deal with two rows (plus whatever we are working with).

This seems to be a pretty niche use case, though; we need to have had enough memory to generate the quasi-inverse in the first place, so this is most likely only relevant if we want to use the pre-computed values on a different machine that is significantly less powerful.

> Are there particular computations you want to do where storing even a single quasi-inverse in memory is problematic?

Not right now, but if we minimize the memory usage we'll be able to spread the work over several machines in a cluster. I'm imagining that should be feasible in the not-too-distant future.

I think that storing the quasi-inverses in a zero-copy-deserializable format will make it much easier to compute extra data in the future. Once we have a large resolution, we will want to use the quasi-inverses to compute the products and Massey products in Ext, and the d2 differentials in the E_2 page. If we implement efficient matrix-matrix multiplication that should not take too much time, but that's only if they are accessible as matrices. I feel like that's the point of saving the quasi-inverses anyway, to make generating lifts easier in the next step.

Edit: Also I think it's a fair assumption that we have exclusive access to the save files. I'm thinking of using an immutable mmap so we can even access them via multiple threads. New data will always be written to a new save anyway.

Fix patch:
diff --git a/ext/Cargo.toml b/ext/Cargo.toml
index e66d23019..b7521b9de 100644
--- a/ext/Cargo.toml
+++ b/ext/Cargo.toml
@@ -16,6 +16,7 @@ once = { path = "crates/once"}
 query = { path = "crates/query" }
 saveload = { path = "crates/saveload" }
 thread-token = { path = "crates/thread-token", optional = true }
+adler = "1"
 
 anyhow = "1.0.0"
 crossbeam-channel = { version = "0.5.0", optional = true }
@@ -25,11 +26,16 @@ serde_json = "1.0.0"
 itertools = { version = "0.10.0", default-features = false, features = ["use_alloc"] }
 chart = { path = "crates/chart" }
 dashmap = "4.0.0"
+byteorder = "1.4.3"
+
+[target.'cfg(not(target_arch = "wasm32"))'.dependencies]
+zstd = "0.9.0"
 
 [dev-dependencies]
 expect-test = "1.1.0"
 rstest = "0.11.0"
 iai = "0.1.1"
+tempfile = "3.0.0"
 
 [lib]
 crate-type = ["cdylib", "rlib"]
@@ -67,3 +73,7 @@ harness = false
 [[bench]]
 name = "secondary"
 harness = false
+
+[[bench]]
+name = "load_resolution"
+harness = false
diff --git a/ext/Makefile b/ext/Makefile
index b3368147d..a4592c10d 100644
--- a/ext/Makefile
+++ b/ext/Makefile
@@ -4,6 +4,7 @@ all:
 	cargo build
 
 test:
+	cargo test
 	cargo test --features concurrent --lib --tests --workspace
 	cargo test --features concurrent --doc --workspace
 
diff --git a/ext/SAVE-FORMAT.md b/ext/SAVE-FORMAT.md
new file mode 100644
index 000000000..e60cb15ff
--- /dev/null
+++ b/ext/SAVE-FORMAT.md
@@ -0,0 +1,179 @@
+# Save file format
+
+This file documents how we save data pertaining to a resolution. We let `/`
+denote the base directory containing the data of a single resolution.
+
+All integers are stored in little-endian.
+
+## Design philosophy
+We have the following goals:
+
+ - Modularity: Each piece of data pertinent to each bidegree is contained in
+   its own file. This gives the possibility of sharding large workloads and
+   only loading the relevant data for each computation.
+
+ - Zero-copy deserialization: For large data structures (e.g. quasi-inverses),
+   it should be possible to `mmap` the save file and use the pointer directly.
+   In particular, the data should be 64 bit-aligned. This use case ties in with
+   modularity, since each such data structure should be contained in exactly
+   one file in order to do RAII properly.
+
+   Of course, we only get true zero-copy deserializaion on little-endian 64-bit
+   machines.
+
+ - Incrementality: We should be able to save data incrementally in case the
+   program crashes, and resume old computations to push them further without
+   too much additional cost.
+
+ - Robustness: It should be easy to detect and correct for data corruption,
+   possibly due to the program crashing halfway through a write. Here "correct"
+   would largely mean deleting the corrupted file and recomputing the data, but
+   modularity helps with minimizing the amount we need to recompute.
+
+ - Space and time efficiency: We should be able to compress the saved data to
+   limit disk usage. If we run this over a cluster, we would be reading files
+   over the network, and space efficiency would be even more important. On the
+   other hand, decompression should not be too slow either. This is necessarily
+   at odds with zero-copy deserialization, and the user should be able to
+   choose whether or not to compress.
+
+## Overview
+In general, we have many "kinds" of data associated to a resolution, such as
+the differential itself and the quasi-inverses. Each kind of data has a name,
+and the data for the bidegree `(s, t)` is contained in
+`/{name}s/{s}_{t}_{name}`.
+
+We support reading compressed files. At the moment, we only support zstd, and
+the compressed file should be saved in `/{name}s/{s}_{t}_{name}.zst`. When
+seeking saved data, we always try the uncompressed version first, as it tends
+to be faster.
+
+(when dealing with compressed files, the `_{name}` suffix makes it easier to
+wildcard match uncompressed files via `/{name}s/*_{name}`)
+
+In general, we always write uncompressed data files, and the compression can be
+performed with other programs separately. This should not pose storage
+problems, since the raw data has to be small enough to be held in memory prior
+to being written.
+
+### File headers
+In addition to a name, each kind of data has a 4-byte magic number to ensure we
+do not mix them up. Further, since it is easy to forget to specify which
+algebra to use, we have a 2-byte magic number for each algebra. These are
+defined as follows:
+```
+adem algebra: 0x0000
+milnor algebra without profile: 0x8000
+milnor algebra with profile: 0x8001
+```
+
+
+Each data file starts with a 16-byte header of the form
+```
+struct {
+    magic: u32,
+    algebra_magic: u16,
+    prime: u16,
+    s: u32,
+    t: i32
+}
+```
+
+Apart from the prime, all the other data is also present in the file name, and
+the header serves as a sanity check that we did not mess up our files. Note
+that it is important that the length of the header is a multiple of 64 bits, so
+that the remaining data is 64 bit-aligned.
+
+### File footer
+At the end of the file, an adler32 checksum of the contents (including the
+header) is appended as a u32.
+
+## Data types
+
+### Differentials
+This has magic `0xD1FF0000`.
+
+This contains the differential in the resolution as well as the augmentation
+map. There is not much reason to put these two together apart from the fact
+that they are both quite small and are useful for most use cases.
+
+For each bidegree, we store the outputs of the maps on the free module
+generators. In particular, this data file records the number of generators in
+each bidegree.
+
+Let `x_1, ..., x_n` be the generators and `d, ε` be the differential and
+augmentation map. The format of the save file is then as follows:
+
+```
+struct {
+    num_gens: u64,
+    len_of_d(x_i): u64,
+    len_of_ε(x_i): u64,
+    d(x_i): [[u64; num_limbs(len_of_d(x_i))]; num_gens]
+    ε(x_i): [[u64; num_limbs(len_of_d(x_i))]; num_gens]
+}
+```
+
+The length is the vector space dimension of the target, while the vectors
+themselves are stored as the limbs of the FpVectors. This will be the case for
+all other data types.
+
+Note that the length of the value of the differential may be shorter than the
+dimension of the target if we use `compute_through_stem`, since we do not yet
+know how many new generators will be added.
+
+### Subspaces
+After computing the differential in bidegree `(s, t)`, we save the kernel of
+the differrential, which will be used when we compute `(s + 1, t)`. We then
+delete the saved file after we are done with computing `(s + 1, t)`. This has
+magic `0x0000D1FF`.
+
+In general, given any subspace, we can choose a canonical basis by taking a
+matrix of spanning (row) vectors and putting it in reduced row echelon form. We
+will always use this basis when dealing with subspaces of a vector space.
+
+Having found such a basis, we have a vector of pivots, whose length is the
+dimension of the ambient space. The kth entry of this vector is index of the
+basis element whose pivot column is `k`, and `-1` if there is no such basis
+element. For example, if the basis is
+```
+[1, 1, 0, 1, 0, 0]
+[0, 0, 1, 0, 0, 0]
+[0, 0, 0, 0, 1, 0]
+[0, 0, 0, 0, 0, 1]
+```
+then the pivot is
+```
+[0, -1, 1, -1, 2, 3]
+```
+
+The format of a subspace is as follows:
+```
+struct {
+    subspace_dimension: u64,
+    total_space_dimension: u64,
+    basis_vectors: [[u64; num_limbs(total_space_dimension)]; subspace_dimension]
+    pivots: [i64; total_space_dimension]
+}
+```
+As usual, the basis elements of a subspace are chosen to be in row reduced
+echelon form.
+
+### Quasi-inverses
+We also store the quasi-inverses to the differential and the augmentation map.
+The associated magics are `0x0100D1FF` and `0x0100A000` respectively.
+
+The format of a quasi-inverse is as follows:
+```
+struct {
+    source_dimension: u64,
+    target_dimension: u64,
+    image_dimension: u64,
+    pivots_of_image: [i64; target_dimension],
+    lift_of_basis_elements: [[u64; num_limbs(source_dimension)]; image_dimension]
+}
+```
+
+Note that it is not necessary to know what the image is, as long as the user of
+the quasi-inverse guarantees their vector is in the image. This is due to the
+fact that our basis is in *reduced* row ecehlon form.
diff --git a/ext/benches/load_resolution.rs b/ext/benches/load_resolution.rs
new file mode 100644
index 000000000..df142643d
--- /dev/null
+++ b/ext/benches/load_resolution.rs
@@ -0,0 +1,11 @@
+use algebra::Algebra;
+use ext::chain_complex::ChainComplex;
+use ext::utils::construct;
+
+fn main() {
+    let resolution = construct("S_2@milnor", Some("/tmp/.ext_bench/S_2_milnor".into())).unwrap();
+    resolution.algebra().compute_basis(100);
+    let start = std::time::Instant::now();
+    resolution.compute_through_bidegree(50, 100);
+    println!("Time: {:?}", start.elapsed());
+}
diff --git a/ext/benches/secondary.rs b/ext/benches/secondary.rs
index 7da04216e..e27a77289 100644
--- a/ext/benches/secondary.rs
+++ b/ext/benches/secondary.rs
@@ -9,7 +9,12 @@ use std::time::Instant;
 fn main() {
     // Attempt to load a resolution of S_2 from resolution_milnor.save, and generates one from
     // scratch if it isn't available.
-    let save_file = std::fs::File::open("resolution_milnor.save").ok();
+    let save_file = std::path::PathBuf::from("S_2_milnor");
+    let save_file = if save_file.exists() {
+        Some(save_file)
+    } else {
+        None
+    };
     let resolution = construct("S_2@milnor", save_file).unwrap();
 
     // Compute the minimal resolution R_{s, t}
diff --git a/ext/crates/algebra/Cargo.toml b/ext/crates/algebra/Cargo.toml
index 11301fb81..baa2205f4 100644
--- a/ext/crates/algebra/Cargo.toml
+++ b/ext/crates/algebra/Cargo.toml
@@ -19,7 +19,6 @@ rustc-hash = "1.1.0"
 bivec = { path = "../bivec"}
 fp = { path = "../fp", default-features = false }
 once = { path = "../once" }
-saveload = { path = "../saveload" }
 
 [dev-dependencies]
 rstest = "0.10.0"
diff --git a/ext/crates/algebra/src/algebra/adem_algebra.rs b/ext/crates/algebra/src/algebra/adem_algebra.rs
index 7b5a8e38a..88c8e25a0 100644
--- a/ext/crates/algebra/src/algebra/adem_algebra.rs
+++ b/ext/crates/algebra/src/algebra/adem_algebra.rs
@@ -204,6 +204,10 @@ impl fmt::Display for AdemAlgebra {
 }
 
 impl Algebra for AdemAlgebra {
+    fn magic(&self) -> u16 {
+        0
+    }
+
     fn prime(&self) -> ValidPrime {
         self.p
     }
diff --git a/ext/crates/algebra/src/algebra/algebra_trait.rs b/ext/crates/algebra/src/algebra/algebra_trait.rs
index 455b65230..25b36f63d 100644
--- a/ext/crates/algebra/src/algebra/algebra_trait.rs
+++ b/ext/crates/algebra/src/algebra/algebra_trait.rs
@@ -18,6 +18,16 @@ use fp::vector::FpVector;
 ///
 /// Algebras may have a distinguished set of generators; see [`GeneratedAlgebra`].
 pub trait Algebra: std::fmt::Display + Send + Sync {
+    /// A magic constant used to identify the algebra in save files. When working with the
+    /// Milnor algebra, it is easy to forget to specify the algebra and load Milnor save files
+    /// with the Adem basis. If we somehow manage to resume computation, this can have
+    /// disasterous consequences. So we store the magic in the save files.
+    ///
+    /// This defaults to 0 for other kinds of algebra that don't care about this problem.
+    fn magic(&self) -> u16 {
+        0
+    }
+
     /// Returns the prime the algebra is over.
     fn prime(&self) -> ValidPrime;
 
@@ -252,6 +262,7 @@ macro_rules! dispatch_algebra {
     ($struct:ty, $dispatch_macro: ident) => {
         impl Algebra for $struct {
             $dispatch_macro! {
+                fn magic(&self) -> u16;
                 fn prime(&self) -> ValidPrime;
                 fn compute_basis(&self, degree: i32);
                 fn dimension(&self, degree: i32, excess: i32) -> usize;
diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs
index 6350e3afa..e78267cc7 100644
--- a/ext/crates/algebra/src/algebra/milnor_algebra.rs
+++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs
@@ -207,6 +207,14 @@ impl MilnorAlgebra {
 }
 
 impl Algebra for MilnorAlgebra {
+    fn magic(&self) -> u16 {
+        if self.profile.is_trivial() {
+            0x8000
+        } else {
+            0x8001
+        }
+    }
+
     fn prime(&self) -> ValidPrime {
         self.p
     }
diff --git a/ext/crates/algebra/src/module/free_module.rs b/ext/crates/algebra/src/module/free_module.rs
index c7e124cc6..2915ec107 100644
--- a/ext/crates/algebra/src/module/free_module.rs
+++ b/ext/crates/algebra/src/module/free_module.rs
@@ -2,7 +2,6 @@ use std::sync::Arc;
 
 use crate::algebra::Algebra;
 use crate::module::Module;
-use bivec::BiVec;
 use fp::vector::{Slice, SliceMut};
 use once::{OnceBiVec, OnceVec};
 
@@ -414,39 +413,6 @@ impl<A: JsonAlgebra> FreeModule<A> {
         Value::from(result)
     }
 }
-
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
-
-impl<A: Algebra> Save for FreeModule<A> {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.num_gens.save(buffer)
-    }
-}
-
-impl<A: Algebra> Load for FreeModule<A> {
-    type AuxData = (Arc<A>, i32);
-
-    fn load(buffer: &mut impl Read, data: &Self::AuxData) -> io::Result<Self> {
-        let algebra = Arc::clone(&data.0);
-        let min_degree = data.1;
-
-        let result = FreeModule::new(algebra, "".to_string(), min_degree);
-
-        let num_gens: BiVec<usize> = Load::load(buffer, &(min_degree, ()))?;
-        result.algebra().compute_basis(num_gens.len() - min_degree);
-
-        for (degree, num) in num_gens.iter_enum() {
-            result.add_generators(degree, *num, None);
-        }
-        // We extend to one degree beyond the number of generators added, which is needed for
-        // resolving to stem. It is always safe to extend more than we "need".
-        result.extend_table_entries(num_gens.len());
-        Ok(result)
-    }
-}
-
 /*
 #[cfg(not(feature = "cache-multiplication"))]
 impl<A: Algebra> FreeModule<A> {
diff --git a/ext/crates/algebra/src/module/free_unstable_module.rs b/ext/crates/algebra/src/module/free_unstable_module.rs
index 5a16fe590..fb046ce99 100644
--- a/ext/crates/algebra/src/module/free_unstable_module.rs
+++ b/ext/crates/algebra/src/module/free_unstable_module.rs
@@ -3,7 +3,6 @@ use std::sync::Arc;
 use crate::algebra::{AdemAlgebra, AdemAlgebraT};
 use crate::module::Module;
 use crate::module::OperationGeneratorPair;
-use bivec::BiVec;
 use fp::vector::SliceMut;
 use once::{OnceBiVec, OnceVec};
 
@@ -324,34 +323,6 @@ impl<A: JsonAlgebra + AdemAlgebraT> FreeUnstableModule<A> {
     }
 }
 
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
-
-impl<A: AdemAlgebraT> Save for FreeUnstableModule<A> {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.num_gens.save(buffer)
-    }
-}
-
-impl<A: AdemAlgebraT> Load for FreeUnstableModule<A> {
-    type AuxData = (Arc<A>, i32);
-
-    fn load(buffer: &mut impl Read, data: &Self::AuxData) -> io::Result<Self> {
-        let algebra = Arc::clone(&data.0);
-        let min_degree = data.1;
-
-        let result = FreeUnstableModule::new(algebra, "".to_string(), min_degree);
-
-        let num_gens: BiVec<usize> = Load::load(buffer, &(min_degree, ()))?;
-        for (degree, num) in num_gens.iter_enum() {
-            result.add_generators(degree, *num, None);
-        }
-        result.extend_table_entries(num_gens.max_degree());
-        Ok(result)
-    }
-}
-
 /*
 #[cfg(not(feature = "cache-multiplication"))]
 impl<A: AdemAlgebraT> FreeUnstableModule<A> {
diff --git a/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs b/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs
index 5c2611244..774507a9e 100644
--- a/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs
+++ b/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs
@@ -246,48 +246,3 @@ impl<A: Algebra> FreeModuleHomomorphism<FreeModule<A>> {
         result
     }
 }
-
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
-
-impl<M: Module> Save for FreeModuleHomomorphism<M> {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.outputs.save(buffer)?;
-        self.images.save(buffer)?;
-        self.kernels.save(buffer)?;
-        self.quasi_inverses.save(buffer)?;
-
-        Ok(())
-    }
-}
-
-impl<M: Module> Load for FreeModuleHomomorphism<M> {
-    type AuxData = (Arc<FreeModule<M::Algebra>>, Arc<M>, i32);
-
-    fn load(buffer: &mut impl Read, data: &Self::AuxData) -> io::Result<Self> {
-        let source: Arc<FreeModule<M::Algebra>> = Arc::clone(&data.0);
-        let target: Arc<M> = Arc::clone(&data.1);
-        let degree_shift = data.2;
-        let min_degree = std::cmp::max(source.min_degree(), target.min_degree() + degree_shift);
-        let p = source.prime();
-
-        let outputs: OnceBiVec<Vec<FpVector>> = Load::load(buffer, &(min_degree, p))?;
-
-        let images = OnceBiVec::<Option<Subspace>>::load(buffer, &(min_degree, Some(p)))?;
-        let kernels = OnceBiVec::<Option<Subspace>>::load(buffer, &(min_degree, Some(p)))?;
-        let quasi_inverses =
-            OnceBiVec::<Option<QuasiInverse>>::load(buffer, &(min_degree, Some(p)))?;
-
-        Ok(Self {
-            source,
-            target,
-            outputs,
-            images,
-            kernels,
-            quasi_inverses,
-            min_degree,
-            degree_shift,
-        })
-    }
-}
diff --git a/ext/crates/bivec/Cargo.toml b/ext/crates/bivec/Cargo.toml
index adb59463b..1fd9e72ac 100644
--- a/ext/crates/bivec/Cargo.toml
+++ b/ext/crates/bivec/Cargo.toml
@@ -6,7 +6,6 @@ edition = "2021"
 
 [dependencies]
 serde = { version = "1.0.0", optional = true }
-saveload = { path = "../saveload" }
 
 [features]
 json = ["serde"]
diff --git a/ext/crates/bivec/src/lib.rs b/ext/crates/bivec/src/lib.rs
index 5756d4f5a..c6485d8b3 100644
--- a/ext/crates/bivec/src/lib.rs
+++ b/ext/crates/bivec/src/lib.rs
@@ -254,23 +254,3 @@ impl<T> IndexMut<i32> for BiVec<T> {
         &mut (self.data[(i - self.min_degree) as usize])
     }
 }
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
-
-impl<T: Save> Save for BiVec<T> {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.data.save(buffer)
-    }
-}
-
-impl<T: Load> Load for BiVec<T> {
-    type AuxData = (i32, T::AuxData);
-
-    fn load(buffer: &mut impl Read, data: &Self::AuxData) -> io::Result<Self> {
-        let min_degree = data.0;
-        let data = Load::load(buffer, &data.1)?;
-
-        Ok(Self { data, min_degree })
-    }
-}
diff --git a/ext/crates/fp/Cargo.toml b/ext/crates/fp/Cargo.toml
index 3334cafd3..7b719e306 100644
--- a/ext/crates/fp/Cargo.toml
+++ b/ext/crates/fp/Cargo.toml
@@ -13,6 +13,7 @@ saveload = { path = "../saveload" }
 itertools = { version = "0.10.0", default-features = false }
 cfg-if = "1.0.0"
 build_const = "0.2.2"
+byteorder = "1.4.3"
 
 [dev-dependencies]
 rand = "0.8.4"
diff --git a/ext/crates/fp/src/matrix/matrix_inner.rs b/ext/crates/fp/src/matrix/matrix_inner.rs
index 51a99dd2f..7311dfb2a 100644
--- a/ext/crates/fp/src/matrix/matrix_inner.rs
+++ b/ext/crates/fp/src/matrix/matrix_inner.rs
@@ -4,6 +4,7 @@ use crate::prime::{self, ValidPrime};
 use crate::vector::{FpVector, Slice, SliceMut};
 
 use std::fmt;
+use std::io::{Read, Write};
 use std::ops::{Index, IndexMut};
 
 /// A matrix! In particular, a matrix with values in F_p. The way we store matrices means it is
@@ -17,7 +18,7 @@ pub struct Matrix {
     /// The pivot columns of the matrix. `pivots[n]` is `k` if column `n` is the `k`th pivot
     /// column, and a negative number otherwise. Said negative number is often -1 but this is not
     /// guaranteed.
-    pivots: Vec<isize>,
+    pub(crate) pivots: Vec<isize>,
 }
 
 impl PartialEq for Matrix {
@@ -62,6 +63,70 @@ impl Matrix {
             pivots: Vec::new(),
         }
     }
+
+    pub fn from_bytes(
+        p: ValidPrime,
+        rows: usize,
+        columns: usize,
+        data: &mut impl Read,
+    ) -> std::io::Result<Matrix> {
+        let mut vectors: Vec<FpVector> = Vec::with_capacity(rows);
+        for _ in 0..rows {
+            vectors.push(FpVector::from_bytes(p, columns, data)?);
+        }
+        Ok(Matrix {
+            p,
+            columns,
+            vectors,
+            pivots: Vec::new(),
+        })
+    }
+
+    pub fn to_bytes(&self, data: &mut impl Write) -> std::io::Result<()> {
+        for v in &self.vectors {
+            v.to_bytes(data)?;
+        }
+        Ok(())
+    }
+
+    /// Read a vector of `isize`
+    pub(crate) fn write_pivot(v: &[isize], buffer: &mut impl Write) -> std::io::Result<()> {
+        cfg_if::cfg_if! {
+            if #[cfg(all(target_endian = "little", target_pointer_width="64"))] {
+                unsafe {
+                    let buf: &[u8] = std::slice::from_raw_parts(v.as_ptr() as *const u8, v.len() * 8);
+                    buffer.write_all(buf).unwrap();
+                }
+            } else {
+                use byteorder::{LittleEndian, WriteBytesExt};
+                for &i in v {
+                    buffer.write_i64::<LittleEndian>(i as i64)?;
+                }
+            }
+        }
+        Ok(())
+    }
+
+    /// Read a vector of `isize` of length `dim`.
+    pub(crate) fn read_pivot(dim: usize, data: &mut impl Read) -> std::io::Result<Vec<isize>> {
+        cfg_if::cfg_if! {
+            if #[cfg(all(target_endian = "little", target_pointer_width="64"))] {
+                let mut image = vec![0; dim];
+                unsafe {
+                    let buf: &mut [u8] = std::slice::from_raw_parts_mut(image.as_mut_ptr() as *mut u8, dim * 8);
+                    data.read_exact(buf).unwrap();
+                }
+                Ok(image)
+            } else {
+                use byteorder::{LittleEndian, ReadBytesExt};
+                let mut image = Vec::with_capacity(dim);
+                for _ in 0..dim {
+                    image.push(data.read_i64::<LittleEndian>()? as isize);
+                }
+                Ok(image)
+            }
+        }
+    }
 }
 
 impl Matrix {
@@ -1071,47 +1136,6 @@ impl<'a> MatrixSliceMut<'a> {
     }
 }
 
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
-
-impl Save for Matrix {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.columns.save(buffer)?;
-        self.vectors.save(buffer)?;
-        self.pivots.save(buffer)?;
-        Ok(())
-    }
-}
-
-impl Load for Matrix {
-    type AuxData = ValidPrime;
-
-    fn load(buffer: &mut impl Read, p: &ValidPrime) -> io::Result<Self> {
-        let columns = usize::load(buffer, &())?;
-        let vectors: Vec<FpVector> = Load::load(buffer, p)?;
-        let mut result = Matrix::from_rows(*p, vectors, columns);
-        result.pivots = Load::load(buffer, &())?;
-        Ok(result)
-    }
-}
-
-impl Save for Subspace {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.matrix.save(buffer)?;
-        Ok(())
-    }
-}
-
-impl Load for Subspace {
-    type AuxData = ValidPrime;
-
-    fn load(buffer: &mut impl Read, p: &ValidPrime) -> io::Result<Self> {
-        let matrix: Matrix = Matrix::load(buffer, p)?;
-        Ok(Subspace { matrix })
-    }
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;
diff --git a/ext/crates/fp/src/matrix/quasi_inverse.rs b/ext/crates/fp/src/matrix/quasi_inverse.rs
index 1d4b6ec82..a6e7ed73a 100644
--- a/ext/crates/fp/src/matrix/quasi_inverse.rs
+++ b/ext/crates/fp/src/matrix/quasi_inverse.rs
@@ -1,6 +1,8 @@
 use super::Matrix;
 use crate::prime::ValidPrime;
-use crate::vector::{Slice, SliceMut};
+use crate::vector::{FpVector, Slice, SliceMut};
+use byteorder::{LittleEndian, ReadBytesExt, WriteBytesExt};
+use std::io::{Read, Write};
 
 /// Given a matrix M, a quasi-inverse Q is a map from the co-domain to the domain such that xQM = x
 /// for all x in the image (recall our matrices act on the right).
@@ -21,6 +23,112 @@ impl QuasiInverse {
         Self { image, preimage }
     }
 
+    pub fn image_dimension(&self) -> usize {
+        self.preimage.rows()
+    }
+
+    pub fn source_dimension(&self) -> usize {
+        self.preimage.columns()
+    }
+
+    pub fn target_dimension(&self) -> usize {
+        match self.image.as_ref() {
+            Some(v) => v.len(),
+            None => self.image_dimension(),
+        }
+    }
+
+    pub fn to_bytes(&self, buffer: &mut impl Write) -> std::io::Result<()> {
+        buffer.write_u64::<LittleEndian>(self.source_dimension() as u64)?;
+        buffer.write_u64::<LittleEndian>(self.target_dimension() as u64)?;
+        buffer.write_u64::<LittleEndian>(self.image_dimension() as u64)?;
+
+        match self.image.as_ref() {
+            None => {
+                for i in 0..self.preimage.rows() {
+                    buffer.write_i64::<LittleEndian>(i as i64)?;
+                }
+            }
+            Some(v) => {
+                Matrix::write_pivot(v, buffer)?;
+            }
+        }
+        self.preimage.to_bytes(buffer)
+    }
+
+    pub fn from_bytes(p: ValidPrime, data: &mut impl Read) -> std::io::Result<Self> {
+        let source_dim = data.read_u64::<LittleEndian>()? as usize;
+        let target_dim = data.read_u64::<LittleEndian>()? as usize;
+        let image_dim = data.read_u64::<LittleEndian>()? as usize;
+
+        let image = Matrix::read_pivot(target_dim, data)?;
+        let preimage = Matrix::from_bytes(p, image_dim, source_dim, data)?;
+        Ok(Self {
+            image: Some(image),
+            preimage,
+        })
+    }
+
+    /// Given a data file containing a quasi-inverse, apply it to all the vectors in `input`
+    /// and write the results to the corresponding vectors in `results`. This reads in the
+    /// quasi-inverse row by row to minimize memory usage.
+    pub fn stream_quasi_inverse(
+        p: ValidPrime,
+        data: &mut impl Read,
+        results: &mut [SliceMut],
+        inputs: &[Slice],
+    ) -> std::io::Result<()> {
+        use crate::limb::Limb;
+
+        let source_dim = data.read_u64::<LittleEndian>()? as usize;
+        let target_dim = data.read_u64::<LittleEndian>()? as usize;
+        let _image_dim = data.read_u64::<LittleEndian>()? as usize;
+
+        let image = Matrix::read_pivot(target_dim, data)?;
+        let mut v = FpVector::new(p, source_dim);
+
+        assert_eq!(results.len(), inputs.len());
+        for result in &*results {
+            assert_eq!(result.as_slice().len(), source_dim);
+        }
+        for input in inputs {
+            // The quasi-inverse might be computed with a smaller target dimension when computing
+            // through stem.
+            assert!(input.len() >= target_dim);
+        }
+
+        let num_limbs = FpVector::num_limbs(p, source_dim);
+
+        for (i, r) in image.into_iter().enumerate() {
+            if r < 0 {
+                continue;
+            }
+
+            let limbs = v.limbs_mut();
+            // Read in data
+            cfg_if::cfg_if! {
+                if #[cfg(target_endian = "little")] {
+                    let num_bytes = num_limbs * std::mem::size_of::<Limb>();
+                    unsafe {
+                        let buf: &mut [u8] = std::slice::from_raw_parts_mut(limbs.as_mut_ptr() as *mut u8, num_bytes);
+                        data.read_exact(buf).unwrap();
+                    }
+                } else {
+                    for limb in limbs.iter_mut() {
+                        let mut bytes: [u8; size_of::<Limb>()] = [0; size_of::<Limb>()];
+                        data.read_exact(&mut bytes)?;
+                        *limb = Limb::from_le_bytes(bytes);
+                    }
+                }
+            };
+            println!("{v}");
+            for (input, result) in inputs.iter().zip(&mut *results) {
+                result.add(v.as_slice(), input.entry(i));
+            }
+        }
+        Ok(())
+    }
+
     pub fn preimage(&self) -> &Matrix {
         &self.preimage
     }
@@ -57,25 +165,52 @@ impl QuasiInverse {
     }
 }
 
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
+#[cfg(test)]
+mod test {
+    use super::*;
 
-impl Save for QuasiInverse {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.image.save(buffer)?;
-        self.preimage.save(buffer)?;
-        Ok(())
-    }
-}
+    #[test]
+    fn test_stream_qi() {
+        let p = ValidPrime::new(2);
+        let qi = QuasiInverse {
+            image: Some(vec![0, -1, 1, -1, 2, 3]),
+            preimage: Matrix::from_vec(
+                p,
+                &[
+                    vec![1, 0, 1, 1],
+                    vec![1, 1, 0, 0],
+                    vec![0, 1, 0, 1],
+                    vec![1, 1, 1, 0],
+                ],
+            ),
+        };
+        let v0 = FpVector::from_slice(p, &[1, 1, 0, 0, 1, 0]);
+        let v1 = FpVector::from_slice(p, &[0, 0, 1, 0, 1, 1]);
 
-impl Load for QuasiInverse {
-    type AuxData = ValidPrime;
+        let mut out0 = FpVector::new(p, 4);
+        let mut out1 = FpVector::new(p, 4);
 
-    fn load(buffer: &mut impl Read, p: &ValidPrime) -> io::Result<Self> {
-        Ok(Self {
-            image: Option::<Vec<isize>>::load(buffer, &Some(()))?,
-            preimage: Matrix::load(buffer, p)?,
-        })
+        let mut cursor = std::io::Cursor::new(Vec::<u8>::new());
+        qi.to_bytes(&mut cursor).unwrap();
+        cursor.set_position(0);
+
+        QuasiInverse::stream_quasi_inverse(
+            p,
+            &mut cursor,
+            &mut [out0.as_slice_mut(), out1.as_slice_mut()],
+            &[v0.as_slice(), v1.as_slice()],
+        )
+        .unwrap();
+
+        let mut bench0 = FpVector::new(p, 4);
+        let mut bench1 = FpVector::new(p, 4);
+
+        qi.apply(bench0.as_slice_mut(), 1, v0.as_slice());
+        qi.apply(bench1.as_slice_mut(), 1, v1.as_slice());
+
+        assert_eq!(out0, bench0, "{out0} != {bench0}");
+        assert_eq!(out1, bench1, "{out1} != {bench1}");
+
+        assert_eq!(cursor.position() as usize, cursor.get_ref().len());
     }
 }
diff --git a/ext/crates/fp/src/matrix/subspace.rs b/ext/crates/fp/src/matrix/subspace.rs
index a3603b5b2..51dd2b7d4 100644
--- a/ext/crates/fp/src/matrix/subspace.rs
+++ b/ext/crates/fp/src/matrix/subspace.rs
@@ -2,6 +2,9 @@ use super::Matrix;
 use crate::prime::ValidPrime;
 use crate::vector::{FpVector, Slice, SliceMut};
 
+use byteorder::{LittleEndian, ReadBytesExt, WriteBytesExt};
+use std::io::{Read, Write};
+
 /// A subspace of a vector space.
 /// # Fields
 ///  * `matrix` - A matrix in reduced row echelon, whose number of columns is the dimension of the
@@ -36,6 +39,25 @@ impl Subspace {
         Self { matrix }
     }
 
+    pub fn from_bytes(p: ValidPrime, data: &mut impl Read) -> std::io::Result<Self> {
+        let rows = data.read_u64::<LittleEndian>()? as usize;
+        let ambient_dimension = data.read_u64::<LittleEndian>()? as usize;
+
+        let mut matrix = Matrix::from_bytes(p, rows, ambient_dimension, data)?;
+
+        matrix.pivots = Matrix::read_pivot(matrix.columns(), data)?;
+
+        Ok(Self { matrix })
+    }
+
+    pub fn to_bytes(&self, buffer: &mut impl Write) -> std::io::Result<()> {
+        buffer.write_u64::<LittleEndian>(self.rows() as u64)?;
+        buffer.write_u64::<LittleEndian>(self.ambient_dimension() as u64)?;
+
+        self.matrix.to_bytes(buffer)?;
+        Matrix::write_pivot(self.matrix.pivots(), buffer)
+    }
+
     pub fn empty_space(p: ValidPrime, dim: usize) -> Self {
         Self::new(p, 0, dim)
     }
diff --git a/ext/crates/fp/src/vector.rs b/ext/crates/fp/src/vector.rs
index 186ca841b..3af81d93c 100644
--- a/ext/crates/fp/src/vector.rs
+++ b/ext/crates/fp/src/vector.rs
@@ -16,6 +16,8 @@ use itertools::Itertools;
 #[cfg(feature = "json")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 use std::convert::TryInto;
+use std::io::{Read, Write};
+use std::mem::size_of;
 
 macro_rules! dispatch_vector_inner {
     // other is a type, but marking it as a :ty instead of :tt means we cannot use it to access its
@@ -191,9 +193,49 @@ impl FpVector {
         match_p!(p, FpVectorP::from(&slice))
     }
 
-    pub fn padded_len(p: ValidPrime, len: usize) -> usize {
+    pub fn num_limbs(p: ValidPrime, len: usize) -> usize {
         let entries_per_limb = entries_per_limb(p);
-        ((len + entries_per_limb - 1) / entries_per_limb) * entries_per_limb
+        (len + entries_per_limb - 1) / entries_per_limb
+    }
+    pub fn padded_len(p: ValidPrime, len: usize) -> usize {
+        Self::num_limbs(p, len) * entries_per_limb(p)
+    }
+
+    pub fn from_bytes(p: ValidPrime, len: usize, data: &mut impl Read) -> std::io::Result<Self> {
+        let num_limbs = Self::num_limbs(p, len);
+
+        let mut limbs: Vec<Limb>;
+        cfg_if::cfg_if! {
+            if #[cfg(target_endian = "little")] {
+                limbs = vec![0; num_limbs];
+                let num_bytes = num_limbs * size_of::<Limb>();
+                unsafe {
+                    let buf: &mut [u8] = std::slice::from_raw_parts_mut(limbs.as_mut_ptr() as *mut u8, num_bytes);
+                    data.read_exact(buf).unwrap();
+                }
+            } else {
+                limbs = Vec::with_capacity(num_limbs);
+
+                for _ in 0..num_limbs {
+                    let mut bytes: [u8; size_of::<Limb>()] = [0; size_of::<Limb>()];
+                    data.read_exact(&mut bytes)?;
+                    limbs.push(Limb::from_le_bytes(bytes));
+                }
+            }
+        };
+        Ok(match_p!(p, FpVectorP::from_raw_parts(len, limbs)))
+    }
+
+    pub fn to_bytes(&self, buffer: &mut impl Write) -> std::io::Result<()> {
+        let num_limbs = Self::num_limbs(self.prime(), self.len());
+        // self.limbs is allowed to have more limbs than necessary, but we only save the
+        // necessary ones.
+
+        for limb in &self.limbs()[0..num_limbs] {
+            let bytes = limb.to_le_bytes();
+            buffer.write_all(&bytes)?;
+        }
+        Ok(())
     }
 
     dispatch_vector! {
@@ -372,7 +414,6 @@ impl<'de> Deserialize<'de> for FpVector {
 
 use saveload::{Load, Save};
 use std::io;
-use std::io::{Read, Write};
 
 impl Save for FpVector {
     fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
@@ -547,6 +588,20 @@ mod test {
     }
 
     test_dim! {
+        fn test_serialize(p: ValidPrime, dim: usize) {
+            use std::io::{Seek, Cursor, SeekFrom};
+
+            let v_arr = random_vector(p, dim);
+            let v = FpVector::from_slice(p, &v_arr);
+
+            let mut cursor = Cursor::new(Vec::<u8>::new());
+            v.to_bytes(&mut cursor).unwrap();
+            cursor.seek(SeekFrom::Start(0)).unwrap();
+
+            let w = FpVector::from_bytes(p, dim, &mut cursor).unwrap();
+            v.assert_vec_eq(&w);
+        }
+
         fn test_add(p: ValidPrime, dim: usize) {
             let mut v_arr = random_vector(p, dim);
             let w_arr = random_vector(p, dim);
@@ -558,7 +613,6 @@ mod test {
                 v_arr[i] = (v_arr[i] + w_arr[i]) % *p;
             }
             v.assert_list_eq(&v_arr);
-
         }
 
         fn test_scale(p: ValidPrime, dim: usize) {
diff --git a/ext/crates/fp/src/vector_2.rs b/ext/crates/fp/src/vector_2.rs
index 4d3adf949..e991217ec 100644
--- a/ext/crates/fp/src/vector_2.rs
+++ b/ext/crates/fp/src/vector_2.rs
@@ -8,6 +8,9 @@ use itertools::Itertools;
 #[cfg(feature = "json")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
+use std::io::{Read, Write};
+use std::mem::size_of;
+
 pub type FpVector = FpVectorP<2>;
 pub type Slice<'a> = SliceP<'a, 2>;
 pub type SliceMut<'a> = SliceMutP<'a, 2>;
@@ -26,9 +29,37 @@ impl FpVector {
         Self::from(&slice)
     }
 
-    pub fn padded_len(_p: ValidPrime, len: usize) -> usize {
+    pub fn num_limbs(_p: ValidPrime, len: usize) -> usize {
         let entries_per_limb = entries_per_limb_const::<2>();
-        ((len + entries_per_limb - 1) / entries_per_limb) * entries_per_limb
+        (len + entries_per_limb - 1) / entries_per_limb
+    }
+
+    pub fn padded_len(p: ValidPrime, len: usize) -> usize {
+        Self::num_limbs(p, len) * entries_per_limb_const::<2>()
+    }
+
+    pub fn from_bytes(p: ValidPrime, len: usize, data: &mut impl Read) -> std::io::Result<Self> {
+        let num_limbs = Self::num_limbs(p, len);
+        let mut limbs = Vec::with_capacity(num_limbs);
+
+        for _ in 0..num_limbs {
+            let mut bytes: [u8; size_of::<Limb>()] = [0; size_of::<Limb>()];
+            data.read_exact(&mut bytes)?;
+            limbs.push(Limb::from_le_bytes(bytes));
+        }
+        Ok(Self::from_raw_parts(len, limbs))
+    }
+
+    pub fn to_bytes(&self, buffer: &mut impl Write) -> std::io::Result<()> {
+        let num_limbs = Self::num_limbs(self.prime(), self.len());
+        // self.limbs is allowed to have more limbs than necessary, but we only save the
+        // necessary ones.
+
+        for limb in &self.limbs()[0..num_limbs] {
+            let bytes = limb.to_le_bytes();
+            buffer.write_all(&bytes)?;
+        }
+        Ok(())
     }
 }
 
@@ -83,7 +114,6 @@ impl<'de> Deserialize<'de> for FpVector {
 
 use saveload::{Load, Save};
 use std::io;
-use std::io::{Read, Write};
 
 impl Save for FpVector {
     fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
diff --git a/ext/crates/fp/src/vector_inner.rs b/ext/crates/fp/src/vector_inner.rs
index 55ad193b7..626040652 100644
--- a/ext/crates/fp/src/vector_inner.rs
+++ b/ext/crates/fp/src/vector_inner.rs
@@ -51,6 +51,11 @@ impl<const P: u32> FpVectorP<P> {
         }
     }
 
+    pub fn from_raw_parts(len: usize, limbs: Vec<Limb>) -> Self {
+        debug_assert_eq!(limbs.len(), limb::number::<P>(len));
+        Self { len, limbs }
+    }
+
     pub fn new_with_capacity_(len: usize, capacity: usize) -> Self {
         let mut limbs = Vec::with_capacity(limb::number::<P>(capacity));
         limbs.resize(limb::number::<P>(len), 0);
diff --git a/ext/crates/once/Cargo.toml b/ext/crates/once/Cargo.toml
index a407b4117..2919a35f0 100644
--- a/ext/crates/once/Cargo.toml
+++ b/ext/crates/once/Cargo.toml
@@ -8,4 +8,3 @@ edition = "2021"
 
 [dependencies]
 bivec = { path = "../bivec"}
-saveload = { path = "../saveload"}
diff --git a/ext/crates/once/src/lib.rs b/ext/crates/once/src/lib.rs
index b89bd58d7..28ecbbb3c 100644
--- a/ext/crates/once/src/lib.rs
+++ b/ext/crates/once/src/lib.rs
@@ -315,33 +315,6 @@ impl<T> IndexMut<u32> for OnceVec<T> {
 unsafe impl<T: Send> Send for OnceVec<T> {}
 unsafe impl<T: Sync> Sync for OnceVec<T> {}
 
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
-
-impl<T: Save> Save for OnceVec<T> {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.len().save(buffer)?;
-        for x in self.iter() {
-            x.save(buffer)?;
-        }
-        Ok(())
-    }
-}
-
-impl<T: Load> Load for OnceVec<T> {
-    type AuxData = T::AuxData;
-
-    fn load(buffer: &mut impl Read, data: &Self::AuxData) -> io::Result<Self> {
-        let len = usize::load(buffer, &())?;
-        let result: OnceVec<T> = OnceVec::new();
-        for _ in 0..len {
-            result.push(T::load(buffer, data)?);
-        }
-        Ok(result)
-    }
-}
-
 #[derive(Clone, PartialEq, Eq)]
 pub struct OnceBiVec<T> {
     pub data: OnceVec<T>,
@@ -493,22 +466,6 @@ impl<T> IndexMut<i32> for OnceBiVec<T> {
     }
 }
 
-impl<T: Save> Save for OnceBiVec<T> {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        self.data.save(buffer)
-    }
-}
-
-impl<T: Load> Load for OnceBiVec<T> {
-    type AuxData = (i32, T::AuxData);
-
-    fn load(buffer: &mut impl Read, data: &Self::AuxData) -> io::Result<Self> {
-        let min_degree = data.0;
-        let data = Load::load(buffer, &data.1)?;
-        Ok(Self { data, min_degree })
-    }
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -531,36 +488,4 @@ mod tests {
             assert_eq!(v[i], i);
         }
     }
-
-    #[test]
-    fn test_saveload() {
-        use std::io::{Cursor, Seek, SeekFrom};
-
-        let v: OnceVec<u32> = OnceVec::new();
-        v.push(6);
-        v.push(3);
-        v.push(4);
-        v.push(2);
-
-        let mut cursor: Cursor<Vec<u8>> = Cursor::new(Vec::new());
-        v.save(&mut cursor).unwrap();
-
-        cursor.seek(SeekFrom::Start(0)).unwrap();
-        let v_saved_then_loaded: OnceVec<u32> = Load::load(&mut cursor, &()).unwrap();
-        assert_eq!(v, v_saved_then_loaded);
-        assert_eq!(0, cursor.bytes().count());
-
-        // let mut w = BiVec::new(-3);
-        // w.push(2);
-        // w.push(6);
-        // w.push(2);
-        // w.push(7);
-
-        // let mut cursor2 : Cursor<Vec<u8>> = Cursor::new(Vec::new());
-        // w.save(&mut cursor2).unwrap();
-        // cursor2.seek(SeekFrom::Start(0)).unwrap();
-        // let w_saved_then_loaded : BiVec<u32> = Load::load(&mut cursor, &(-3, ())).unwrap();
-
-        // assert_eq!(w, w_saved_then_loaded);
-    }
 }
diff --git a/ext/examples/benchmarks/lift_hom-Cnu-S_2 b/ext/examples/benchmarks/lift_hom-Cnu-S_2
index a7ce206db..ba7039f0f 100644
--- a/ext/examples/benchmarks/lift_hom-Cnu-S_2
+++ b/ext/examples/benchmarks/lift_hom-Cnu-S_2
@@ -1,4 +1,4 @@
-lift_hom -- Cnu S_2 0 -4 "" 15 15 [1]
+lift_hom -- Cnu "" no S_2 "" 0 -4 15 15 [1]
 F(x_(0, 0, 0)) = []
 F(x_(0, 1, 0)) = []
 F(x_(0, 2, 0)) = []
diff --git a/ext/examples/benchmarks/lift_hom-S_2-Cnu b/ext/examples/benchmarks/lift_hom-S_2-Cnu
index 23e553fcc..f99648d20 100644
--- a/ext/examples/benchmarks/lift_hom-S_2-Cnu
+++ b/ext/examples/benchmarks/lift_hom-S_2-Cnu
@@ -1,4 +1,4 @@
-lift_hom -- S_2 Cnu 0 0 "" 15 15 [1]
+lift_hom -- S_2 "" no Cnu "" 0 0 15 15 [1]
 F(x_(0, 0, 0)) = [1]
 F(x_(0, 1, 0)) = [1]
 F(x_(0, 2, 0)) = [1]
diff --git a/ext/examples/benchmarks/lift_hom-g-S_2 b/ext/examples/benchmarks/lift_hom-g-S_2
index b23763335..9a0517b4c 100644
--- a/ext/examples/benchmarks/lift_hom-g-S_2
+++ b/ext/examples/benchmarks/lift_hom-g-S_2
@@ -1,4 +1,4 @@
-lift_hom -- S_2 S_2 4 20 "" 5 5 [1]
+lift_hom -- S_2 "" y 4 20 5 5 [1]
 F(x_(0, 0, 0)) = [1]
 F(x_(0, 1, 0)) = [1]
 F(x_(0, 2, 0)) = [1]
diff --git a/ext/examples/lift_hom.rs b/ext/examples/lift_hom.rs
index 464598e34..21ed6414d 100644
--- a/ext/examples/lift_hom.rs
+++ b/ext/examples/lift_hom.rs
@@ -42,23 +42,25 @@
 use algebra::module::{BoundedModule, Module};
 use ext::chain_complex::ChainComplex;
 use ext::resolution_homomorphism::ResolutionHomomorphism;
-use ext::utils::{construct, Config};
+use ext::utils;
 use fp::matrix::Matrix;
 
-use anyhow::anyhow;
-use std::convert::TryInto;
-use std::{fs::File, sync::Arc};
+use std::sync::Arc;
 
 fn main() -> anyhow::Result<()> {
-    let target: Config = query::with_default("Target module", "S_2", |name| {
-        let target: Config = name.try_into()?;
-        match target.module["type"].as_str() {
-            Some("finite dimensional module") => Ok(target),
-            _ => Err(anyhow!("Target must be finite dimensional")),
-        }
-    });
+    let target = Arc::new(utils::query_module_only("Target module", None)?);
+
+    let source_equal_target = query::yes_no("Source equal to target?");
+    let source = if source_equal_target {
+        Arc::clone(&target)
+    } else {
+        let mut s = utils::query_module_only("Source module", None)?;
+        s.load_quasi_inverse = false;
+        Arc::new(s)
+    };
 
-    let source: Config = query::with_default("Source module", "Cnu", |name| name.try_into());
+    assert_eq!(source.prime(), target.prime());
+    let p = source.prime();
 
     let shift_s: u32 = query::with_default("s of Ext class", "0", str::parse);
     let shift_n: i32 = query::with_default("n of Ext class", "0", str::parse);
@@ -67,66 +69,28 @@ fn main() -> anyhow::Result<()> {
     #[cfg(feature = "concurrent")]
     let bucket = ext::utils::query_bucket();
 
-    #[allow(clippy::redundant_closure)]
-    let target_save_file = query::optional("Target save file", |s| File::open(s));
+    let s: u32 = query::with_default("Max target s", "10", str::parse);
+    let n: i32 = query::with_default("Max target n", "10", str::parse);
 
-    let (target, source) = if target == source {
-        let target = match target_save_file {
-            Some(f) => construct(target, Some(f))?,
-            None => {
-                let s: u32 = query::with_default("Max target s", "10", str::parse);
-                let n: i32 = query::with_default("Max target n", "10", str::parse);
-                let target = construct(target, None)?;
+    if source_equal_target {
+        #[cfg(feature = "concurrent")]
+        target.compute_through_stem_concurrent(s + shift_s, n + std::cmp::max(0, shift_n), &bucket);
 
-                #[cfg(feature = "concurrent")]
-                target.compute_through_stem_concurrent(
-                    s + shift_s,
-                    n + std::cmp::max(0, shift_n),
-                    &bucket,
-                );
-
-                #[cfg(not(feature = "concurrent"))]
-                target.compute_through_stem(s + shift_s, n + std::cmp::max(0, shift_n));
-
-                target
-            }
-        };
-        let target = Arc::new(target);
-        (Arc::clone(&target), target)
+        #[cfg(not(feature = "concurrent"))]
+        target.compute_through_stem(s + shift_s, n + std::cmp::max(0, shift_n));
     } else {
-        match target_save_file {
-            Some(f) => (
-                Arc::new(construct(target, Some(f))?),
-                Arc::new(construct(
-                    source,
-                    #[allow(clippy::redundant_closure)]
-                    Some(query::raw("Source save file", |s| File::open(s))),
-                )?),
-            ),
-            None => {
-                let s = query::with_default("Max target s", "10", str::parse);
-                let n: i32 = query::with_default("Max target n", "10", str::parse);
-                let (target, source) = (construct(target, None)?, construct(source, None)?);
-
-                #[cfg(feature = "concurrent")]
-                {
-                    source.compute_through_stem_concurrent(s + shift_s, n + shift_n, &bucket);
-                    target.compute_through_stem_concurrent(s, n, &bucket);
-                }
-
-                #[cfg(not(feature = "concurrent"))]
-                {
-                    source.compute_through_stem(s + shift_s, n + shift_n);
-                    target.compute_through_stem(s, n);
-                }
-
-                (Arc::new(target), Arc::new(source))
-            }
+        #[cfg(feature = "concurrent")]
+        {
+            source.compute_through_stem_concurrent(s + shift_s, n + shift_n, &bucket);
+            target.compute_through_stem_concurrent(s, n, &bucket);
         }
-    };
 
-    assert_eq!(source.prime(), target.prime());
-    let p = source.prime();
+        #[cfg(not(feature = "concurrent"))]
+        {
+            source.compute_through_stem(s + shift_s, n + shift_n);
+            target.compute_through_stem(s, n);
+        }
+    }
 
     let target_module = target.complex().module(0);
     let hom = ResolutionHomomorphism::new(String::new(), source, target, shift_s, shift_t);
diff --git a/ext/examples/resolve.rs b/ext/examples/resolve.rs
index b16d69ce1..e0343ebfa 100644
--- a/ext/examples/resolve.rs
+++ b/ext/examples/resolve.rs
@@ -16,22 +16,23 @@
 //!
 use ext::chain_complex::{ChainComplex, FreeChainComplex};
 use ext::utils::construct;
-use saveload::Save;
-use std::fs::File;
+use std::path::PathBuf;
 
 fn main() -> anyhow::Result<()> {
-    let res = query::with_default("Module", "S_2", |name| {
+    let mut res = query::with_default("Module", "S_2", |name| {
         construct(
             name,
-            query::optional("Load from save?", |filename| File::open(filename)),
+            query::optional("Save directory", |filename| {
+                core::result::Result::<PathBuf, std::convert::Infallible>::Ok(PathBuf::from(
+                    filename,
+                ))
+            }),
         )
     });
+    res.load_quasi_inverse = false;
 
     let max_s = query::with_default("Max s", "15", str::parse);
     let max_t = query::with_default("Max t", "30", str::parse);
-    // Clippy false positive
-    #[allow(clippy::redundant_closure)]
-    let save_file: Option<File> = query::optional("Save file", |s| File::create(s));
 
     #[cfg(not(feature = "concurrent"))]
     res.compute_through_bidegree(max_s, max_t);
@@ -43,10 +44,5 @@ fn main() -> anyhow::Result<()> {
     }
 
     println!("{}", res.graded_dimension_string());
-
-    if let Some(file) = save_file {
-        let mut file = std::io::BufWriter::new(file);
-        res.save(&mut file)?;
-    }
     Ok(())
 }
diff --git a/ext/examples/resolve_through_stem.rs b/ext/examples/resolve_through_stem.rs
index 013fcde23..dccf6d305 100644
--- a/ext/examples/resolve_through_stem.rs
+++ b/ext/examples/resolve_through_stem.rs
@@ -21,22 +21,23 @@
 
 use ext::chain_complex::FreeChainComplex;
 use ext::utils::construct;
-use saveload::Save;
-use std::fs::File;
+use std::path::PathBuf;
 
 fn main() -> anyhow::Result<()> {
-    let res = query::with_default("Module", "S_2", |name| {
+    let mut res = query::with_default("Module", "S_2", |name| {
         construct(
             name,
-            query::optional("Load from save?", |filename| File::open(filename)),
+            query::optional("Save directory", |filename| {
+                core::result::Result::<PathBuf, std::convert::Infallible>::Ok(PathBuf::from(
+                    filename,
+                ))
+            }),
         )
     });
+    res.load_quasi_inverse = false;
 
     let max_s = query::with_default("Max s", "15", str::parse);
     let max_n = query::with_default("Max n", "30", str::parse);
-    // Clippy false positive
-    #[allow(clippy::redundant_closure)]
-    let save_file: Option<File> = query::optional("Save file", |s| File::create(s));
 
     #[cfg(not(feature = "concurrent"))]
     res.compute_through_stem(max_s, max_n);
@@ -49,10 +50,5 @@ fn main() -> anyhow::Result<()> {
 
     println!("{}", res.graded_dimension_string());
 
-    if let Some(file) = save_file {
-        let mut file = std::io::BufWriter::new(file);
-        res.save(&mut file)?;
-    }
-
     Ok(())
 }
diff --git a/ext/examples/steenrod.rs b/ext/examples/steenrod.rs
index 1a17d2334..84c9fe153 100644
--- a/ext/examples/steenrod.rs
+++ b/ext/examples/steenrod.rs
@@ -4,7 +4,7 @@ use algebra::module::homomorphism::{
 use algebra::module::{BoundedModule, Module};
 use ext::chain_complex::{ChainComplex, TensorChainComplex};
 use ext::resolution_homomorphism::ResolutionHomomorphism;
-use ext::utils::construct;
+use ext::utils;
 use ext::yoneda::yoneda_representative_element;
 use fp::matrix::Matrix;
 use fp::prime::ValidPrime;
@@ -21,8 +21,10 @@ use std::{thread, thread::JoinHandle};
 use crossbeam_channel::{unbounded, Receiver};
 
 fn main() -> anyhow::Result<()> {
-    let resolution =
-        Arc::new(construct("S_2", std::fs::File::open("resolution_adem.save").ok()).unwrap());
+    let resolution = Arc::new(utils::query_module_only(
+        "Module",
+        Some(algebra::AlgebraType::Adem),
+    )?);
 
     let complex = resolution.complex();
     let module = complex.module(0);
diff --git a/ext/src/chain_complex/mod.rs b/ext/src/chain_complex/mod.rs
index 4faa7f29b..7e4d5f1b3 100644
--- a/ext/src/chain_complex/mod.rs
+++ b/ext/src/chain_complex/mod.rs
@@ -8,7 +8,7 @@ use algebra::module::{FreeModule, Module};
 use algebra::Algebra;
 use fp::matrix::Subquotient;
 use fp::prime::ValidPrime;
-use fp::vector::FpVector;
+use fp::vector::{FpVector, Slice, SliceMut};
 use std::sync::Arc;
 
 // pub use hom_complex::HomComplex;
@@ -162,6 +162,31 @@ pub trait ChainComplex: Send + Sync {
             max_s: self.next_homological_degree(),
         }
     }
+
+    /// Apply the quasi-inverse of the (s, t)th differential to the list of inputs and results.
+    /// This defaults to applying `self.differentials(s).quasi_inverse(t)`, but in some cases
+    /// the quasi-inverse might be stored separately on disk.
+    ///
+    /// This returns whether the application was successful
+    #[must_use]
+    fn apply_quasi_inverse(
+        &self,
+        results: &mut [SliceMut],
+        s: u32,
+        t: i32,
+        inputs: &[Slice],
+    ) -> bool {
+        assert_eq!(results.len(), inputs.len());
+
+        if let Some(qi) = self.differential(s).quasi_inverse(t) {
+            for (input, result) in inputs.iter().zip(results) {
+                qi.apply(result.copy(), 1, *input);
+            }
+            true
+        } else {
+            false
+        }
+    }
 }
 
 /// An iterator returned by [`ChainComplex::iter_stem`]
diff --git a/ext/src/resolution.rs b/ext/src/resolution.rs
index 0ba65d0f8..1493101c5 100644
--- a/ext/src/resolution.rs
+++ b/ext/src/resolution.rs
@@ -1,16 +1,24 @@
 use std::sync::{Arc, Mutex};
 
 use crate::chain_complex::{AugmentedChainComplex, ChainComplex};
+use crate::utils;
 use algebra::module::homomorphism::{FreeModuleHomomorphism, ModuleHomomorphism};
 use algebra::module::{FreeModule, Module};
 use algebra::Algebra;
-use fp::matrix::{AugmentedMatrix, Subspace};
+use fp::matrix::{AugmentedMatrix, QuasiInverse, Subspace};
 use fp::prime::ValidPrime;
-use fp::vector::FpVector;
+use fp::vector::{FpVector, Slice, SliceMut};
 use once::OnceVec;
 
+use std::fs::File;
+use std::io::{BufReader, BufWriter, Read, Write};
+use std::path::PathBuf;
+
+use anyhow::Context;
 use dashmap::DashMap;
 
+use byteorder::{LittleEndian, ReadBytesExt, WriteBytesExt};
+
 #[cfg(feature = "concurrent")]
 use crossbeam_channel::{unbounded, Receiver};
 
@@ -23,6 +31,40 @@ use thread_token::TokenBucket;
 /// number if needs be, but up to the 140th stem we only see at most 8 new generators.
 const MAX_NEW_GENS: usize = 10;
 
+#[derive(Debug, Copy, Clone, Eq, PartialEq)]
+pub enum SaveData {
+    Kernel,
+    Differential,
+    ResQi,
+    AugmentationQi,
+}
+
+impl SaveData {
+    pub fn magic(self) -> u32 {
+        match self {
+            Self::Kernel => 0x0000D1FF,
+            Self::Differential => 0xD1FF0000,
+            Self::ResQi => 0x0100D1FF,
+            Self::AugmentationQi => 0x0100A000,
+        }
+    }
+
+    pub fn name(self) -> &'static str {
+        match self {
+            Self::Kernel => "kernel",
+            Self::Differential => "differential",
+            Self::ResQi => "res_qi",
+            Self::AugmentationQi => "augmentation_qi",
+        }
+    }
+
+    pub fn resolution_data() -> impl Iterator<Item = SaveData> {
+        use SaveData::*;
+        static KINDS: [SaveData; 4] = [Kernel, Differential, ResQi, AugmentationQi];
+        KINDS.iter().copied()
+    }
+}
+
 /// A resolution of a chain complex.
 pub struct Resolution<CC: ChainComplex> {
     lock: Mutex<()>,
@@ -37,24 +79,126 @@ pub struct Resolution<CC: ChainComplex> {
     ///  returned by `generate_old_kernel_and_compute_new_kernel`, to be used if we run
     ///  compute_through_degree again.
     kernels: DashMap<(u32, i32), Subspace>,
+    save_dir: Option<PathBuf>,
+    pub should_save: bool,
+    pub load_quasi_inverse: bool,
 }
 
 impl<CC: ChainComplex> Resolution<CC> {
     pub fn new(complex: Arc<CC>) -> Self {
+        // It doesn't error if the save file is None
+        Self::new_with_save(complex, None).unwrap()
+    }
+
+    pub fn new_with_save(complex: Arc<CC>, mut save_dir: Option<PathBuf>) -> anyhow::Result<Self> {
         let algebra = complex.algebra();
         let min_degree = complex.min_degree();
         let zero_module = Arc::new(FreeModule::new(algebra, "F_{-1}".to_string(), min_degree));
 
-        Self {
+        if let Some(p) = save_dir.as_mut() {
+            use std::io::{Error, ErrorKind};
+
+            if !p.exists() {
+                std::fs::create_dir_all(&p)
+                    .context("Failed to create save directory for resolution")?;
+            } else if !p.is_dir() {
+                return Err(Error::new(
+                    ErrorKind::InvalidInput,
+                    format!("{p:?} is not a directory"),
+                )
+                .into());
+            }
+            for subdir in SaveData::resolution_data() {
+                p.push(format!("{}s", subdir.name()));
+                if !p.exists() {
+                    std::fs::create_dir_all(&p)
+                        .with_context(|| format!("Failed to create directory {p:?}"))?;
+                } else if !p.is_dir() {
+                    return Err(Error::new(
+                        ErrorKind::InvalidInput,
+                        format!("{p:?} is not a directory"),
+                    )
+                    .into());
+                }
+                p.pop();
+            }
+        }
+
+        Ok(Self {
             complex,
             zero_module,
+            should_save: save_dir.is_some(),
+            save_dir,
             lock: Mutex::new(()),
 
             chain_maps: OnceVec::new(),
             modules: OnceVec::new(),
             differentials: OnceVec::new(),
             kernels: DashMap::new(),
+            load_quasi_inverse: true,
+        })
+    }
+
+    /// This panics if there is no save dir
+    fn get_save_path(&self, kind: SaveData, s: u32, t: i32) -> PathBuf {
+        let name = kind.name();
+        let mut p = self.save_dir.clone().unwrap();
+        p.push(format!("{name}s/{s}_{t}_{name}"));
+        p
+    }
+
+    // When compiling to wasm we don't mutate path
+    #[allow(unused_mut)]
+    fn search_file(mut path: PathBuf) -> Option<Box<dyn Read>> {
+        // We should try in decreasing order of access speed.
+        match File::open(&path) {
+            Ok(f) => return Some(Box::new(utils::ChecksumReader::new(BufReader::new(f)))),
+            Err(e) => {
+                if e.kind() != std::io::ErrorKind::NotFound {
+                    panic!("Error when opening {path:?}");
+                }
+            }
         }
+
+        #[cfg(not(target_arch = "wasm32"))]
+        {
+            path.set_extension("zst");
+            match File::open(&path) {
+                Ok(f) => {
+                    return Some(Box::new(utils::ChecksumReader::new(
+                        zstd::stream::Decoder::new(f).unwrap(),
+                    )))
+                }
+                Err(e) => {
+                    if e.kind() != std::io::ErrorKind::NotFound {
+                        panic!("Error when opening {path:?}");
+                    }
+                }
+            }
+        }
+
+        None
+    }
+
+    /// This panics if there is no save dir
+    fn open_save_file(&self, kind: SaveData, s: u32, t: i32) -> Option<Box<dyn Read>> {
+        let mut f = Self::search_file(self.get_save_path(kind, s, t))?;
+        utils::validate_header(kind.magic(), &*self.algebra(), self.prime(), s, t, &mut f).unwrap();
+        Some(f)
+    }
+
+    fn create_save_file(&self, kind: SaveData, s: u32, t: i32) -> impl Write {
+        let p = self.get_save_path(kind, s, t);
+
+        let f = std::fs::OpenOptions::new()
+            .write(true)
+            .create_new(true)
+            .open(&p)
+            .with_context(|| format!("Failed to create save file {p:?}"))
+            .unwrap();
+        let mut f = utils::ChecksumWriter::new(BufWriter::new(f));
+        utils::write_header(kind.magic(), &*self.algebra(), self.prime(), s, t, &mut f).unwrap();
+        f
     }
 
     /// This function prepares the Resolution object to perform computations up to the
@@ -110,6 +254,14 @@ impl<CC: ChainComplex> Resolution<CC> {
 
         let p = self.prime();
 
+        if self.save_dir.is_some() {
+            if let Some(mut f) = self.open_save_file(SaveData::Kernel, s, t) {
+                return Subspace::from_bytes(p, &mut f)
+                    .with_context(|| format!("Failed to read kernel at ({s}, {t})"))
+                    .unwrap();
+            }
+        }
+
         let complex = self.complex();
         complex.compute_through_bidegree(s, t);
 
@@ -138,7 +290,16 @@ impl<CC: ChainComplex> Resolution<CC> {
         matrix.segment(2, 2).add_identity();
         matrix.row_reduce();
 
-        matrix.compute_kernel()
+        let kernel = matrix.compute_kernel();
+
+        if self.should_save && self.save_dir.is_some() {
+            let mut f = self.create_save_file(SaveData::Kernel, s, t);
+            kernel
+                .to_bytes(&mut f)
+                .with_context(|| format!("Failed to write kernel at ({s}, {t})"))
+                .unwrap();
+        }
+        kernel
     }
 
     /// Call our resolution $X$, and the chain complex to resolve $C$. This is a legitimate
@@ -259,6 +420,88 @@ impl<CC: ChainComplex> Resolution<CC> {
         target_res.extend_table_entries(t);
         let target_res_dimension = target_res.dimension(t);
 
+        if self.save_dir.is_some() {
+            if let Some(mut f) = self.open_save_file(SaveData::Differential, s, t) {
+                let num_new_gens = f.read_u64::<LittleEndian>().unwrap() as usize;
+                // It can be smaller than target_res_dimension if we resolved through stem
+                let saved_target_res_dimension = f.read_u64::<LittleEndian>().unwrap() as usize;
+                assert!(
+                    target_res_dimension >= saved_target_res_dimension,
+                    "Malformed data: mismatched resolution target dimension"
+                );
+                assert_eq!(
+                    target_cc_dimension,
+                    f.read_u64::<LittleEndian>().unwrap() as usize,
+                    "Malformed data: mismatched augmentation target dimension"
+                );
+
+                source.add_generators(t, num_new_gens, None);
+
+                let mut d_targets = Vec::with_capacity(num_new_gens);
+                let mut a_targets = Vec::with_capacity(num_new_gens);
+
+                for _ in 0..num_new_gens {
+                    d_targets
+                        .push(FpVector::from_bytes(p, saved_target_res_dimension, &mut f).unwrap());
+                }
+                for _ in 0..num_new_gens {
+                    a_targets.push(FpVector::from_bytes(p, target_cc_dimension, &mut f).unwrap());
+                }
+                drop(f);
+                current_differential.add_generators_from_rows(t, d_targets);
+                current_chain_map.add_generators_from_rows(t, a_targets);
+
+                // res qi
+                if self.load_quasi_inverse {
+                    if let Some(mut f) = self.open_save_file(SaveData::ResQi, s, t) {
+                        let res_qi = QuasiInverse::from_bytes(p, &mut f).unwrap();
+
+                        assert!(
+                            target_res_dimension >= res_qi.target_dimension(),
+                            "Malformed data: mismatched resolution target dimension in qi at ({s}, {t})"
+                            );
+                        assert_eq!(
+                            res_qi.source_dimension(),
+                            source_dimension + num_new_gens,
+                            "Malformed data: mismatched source dimension in resolution qi at ({s}, {t})"
+                            );
+
+                        current_differential.set_quasi_inverse(t, Some(res_qi));
+                    } else {
+                        current_differential.set_quasi_inverse(t, None);
+                    }
+
+                    if let Some(mut f) = self.open_save_file(SaveData::AugmentationQi, s, t) {
+                        let cm_qi = QuasiInverse::from_bytes(p, &mut f).unwrap();
+
+                        assert_eq!(
+                        cm_qi.target_dimension(),
+                        target_cc_dimension,
+                        "Malformed data: mismatched augmentation target dimension in qi at ({s}, {t})"
+                    );
+                        assert_eq!(
+                        cm_qi.source_dimension(),
+                        source_dimension + num_new_gens,
+                        "Malformed data: mismatched source dimension in augmentation qi at ({s}, {t})"
+                    );
+
+                        current_chain_map.set_quasi_inverse(t, Some(cm_qi));
+                    } else {
+                        current_chain_map.set_quasi_inverse(t, None);
+                    }
+                } else {
+                    current_differential.set_quasi_inverse(t, None);
+                    current_chain_map.set_quasi_inverse(t, None);
+                }
+                current_differential.set_kernel(t, None);
+                current_differential.set_image(t, None);
+
+                current_chain_map.set_kernel(t, None);
+                current_chain_map.set_image(t, None);
+                return;
+            }
+        }
+
         let mut matrix = AugmentedMatrix::<3>::new_with_capacity(
             p,
             source_dimension,
@@ -275,7 +518,17 @@ impl<CC: ChainComplex> Resolution<CC> {
         matrix.row_reduce();
 
         if !self.has_computed_bidegree(s + 1, t) {
-            self.kernels.insert((s, t), matrix.compute_kernel());
+            let kernel = matrix.compute_kernel();
+            if self.should_save && self.save_dir.is_some() {
+                let mut f = self.create_save_file(SaveData::Kernel, s, t);
+
+                kernel
+                    .to_bytes(&mut f)
+                    .with_context(|| format!("Failed to write kernel at ({s}, {t})"))
+                    .unwrap();
+            }
+
+            self.kernels.insert((s, t), kernel);
         }
 
         // Now add generators to surject onto C_{s, t}.
@@ -413,13 +666,60 @@ impl<CC: ChainComplex> Resolution<CC> {
                 }
             }
         }
+        let (cm_qi, res_qi) =
+            if (self.should_save && self.save_dir.is_some()) || self.load_quasi_inverse {
+                let (c, r) = matrix.compute_quasi_inverses();
+                (Some(c), Some(r))
+            } else {
+                (None, None)
+            };
+
+        if self.should_save && self.save_dir.is_some() {
+            // Write differentials
+            let mut f = self.create_save_file(SaveData::Differential, s, t);
+
+            f.write_u64::<LittleEndian>(num_new_gens as u64).unwrap();
+            f.write_u64::<LittleEndian>(target_res_dimension as u64)
+                .unwrap();
+            f.write_u64::<LittleEndian>(target_cc_dimension as u64)
+                .unwrap();
+
+            for n in 0..num_new_gens {
+                current_differential.output(t, n).to_bytes(&mut f).unwrap();
+            }
+            for n in 0..num_new_gens {
+                current_chain_map.output(t, n).to_bytes(&mut f).unwrap();
+            }
+            drop(f);
+
+            // Write resolution qi
+            let mut f = self.create_save_file(SaveData::ResQi, s, t);
+            res_qi.as_ref().unwrap().to_bytes(&mut f).unwrap();
+            drop(f);
+
+            // Write augmentation qi
+            let mut f = self.create_save_file(SaveData::AugmentationQi, s, t);
+            cm_qi.as_ref().unwrap().to_bytes(&mut f).unwrap();
+            drop(f);
+
+            // Delete kernel
+            if self.should_save && s > 0 {
+                let ker_path = self.get_save_path(SaveData::Kernel, s - 1, t);
+                if ker_path.exists() {
+                    std::fs::remove_file(ker_path).unwrap();
+                }
+            }
+        }
 
-        let (cm_qi, res_qi) = matrix.compute_quasi_inverses();
-
-        current_chain_map.set_quasi_inverse(t, Some(cm_qi));
+        if self.load_quasi_inverse {
+            current_chain_map.set_quasi_inverse(t, cm_qi);
+            current_differential.set_quasi_inverse(t, res_qi);
+        } else {
+            current_chain_map.set_quasi_inverse(t, None);
+            current_differential.set_quasi_inverse(t, None);
+        }
         current_chain_map.set_kernel(t, None);
         current_chain_map.set_image(t, None);
-        current_differential.set_quasi_inverse(t, Some(res_qi));
         current_differential.set_kernel(t, None);
         current_differential.set_image(t, None);
     }
@@ -606,10 +906,21 @@ impl<CC: ChainComplex> Resolution<CC> {
                             token = bucket.recv_or_release(token, &last_receiver);
                             if !self.has_computed_bidegree(s, t) {
                                 if s as i32 + max_n + 1 == t {
-                                    // We compute the kernel now and save it. If we don't call
-                                    // this, then the step_resolution stage will compute it, but we
-                                    // may be unnecessarily idling in that case.
-                                    self.kernels.insert((s, t), self.get_kernel(s, t));
+                                    // This is the bidegree just beyond max_n. We are not computing
+                                    // this bidegree, but if we have to compute the next s, we have
+                                    // to compute the kernel of this bidegree.
+                                    //
+                                    // We can wait until the next step to compute the kernel, but
+                                    // we are already ready to do so, so let's do it now while we
+                                    // have time.
+                                    if !self.has_computed_bidegree(s + 1, t)
+                                        && (self.save_dir.is_none()
+                                            || !self
+                                                .get_save_path(SaveData::Differential, s + 1, t)
+                                                .exists())
+                                    {
+                                        self.kernels.insert((s, t), self.get_kernel(s, t));
+                                    }
                                     // The next t cannot be computed yet
                                     continue;
                                 } else {
@@ -685,6 +996,32 @@ impl<CC: ChainComplex> ChainComplex for Resolution<CC> {
     fn next_homological_degree(&self) -> u32 {
         self.modules.len() as u32
     }
+
+    fn apply_quasi_inverse(
+        &self,
+        results: &mut [SliceMut],
+        s: u32,
+        t: i32,
+        inputs: &[Slice],
+    ) -> bool {
+        assert_eq!(results.len(), inputs.len());
+
+        if let Some(qi) = self.differential(s).quasi_inverse(t) {
+            for (input, result) in inputs.iter().zip(results) {
+                qi.apply(result.copy(), 1, *input);
+            }
+            true
+        } else if self.save_dir.is_some() {
+            if let Some(mut f) = self.open_save_file(SaveData::ResQi, s, t) {
+                QuasiInverse::stream_quasi_inverse(self.prime(), &mut f, results, inputs).unwrap();
+                true
+            } else {
+                false
+            }
+        } else {
+            false
+        }
+    }
 }
 
 impl<CC: ChainComplex> AugmentedChainComplex for Resolution<CC> {
@@ -700,73 +1037,6 @@ impl<CC: ChainComplex> AugmentedChainComplex for Resolution<CC> {
     }
 }
 
-use saveload::{Load, Save};
-use std::io;
-use std::io::{Read, Write};
-
-impl<CC: ChainComplex> Save for Resolution<CC> {
-    fn save(&self, buffer: &mut impl Write) -> io::Result<()> {
-        // This is no longer used but we keep it around for backwards compatibility
-        let max_algebra_dim = self.module(0).max_computed_degree() - self.min_degree();
-        max_algebra_dim.save(buffer)?;
-
-        self.modules.save(buffer)?;
-        self.kernels.save(buffer)?;
-        self.differentials.save(buffer)?;
-        self.chain_maps.save(buffer)?;
-        Ok(())
-    }
-}
-
-impl<CC: ChainComplex> Load for Resolution<CC> {
-    type AuxData = Arc<CC>;
-
-    fn load(buffer: &mut impl Read, cc: &Self::AuxData) -> io::Result<Self> {
-        // This is no longer used but we keep it around for backwards compatibility
-        let _max_algebra_dim = i32::load(buffer, &())?;
-
-        let mut result = Resolution::new(Arc::clone(cc));
-
-        let algebra = result.algebra();
-        let p = result.prime();
-        let min_degree = result.min_degree();
-
-        result.modules = Load::load(buffer, &(Arc::clone(&algebra), min_degree))?;
-        result.kernels = Load::load(buffer, &(((), ()), p))?;
-
-        let max_s = result.modules.len();
-        assert!(max_s > 0, "cannot load uninitialized resolution");
-
-        let len = usize::load(buffer, &())?;
-        assert_eq!(len, max_s);
-
-        result.differentials.push(Load::load(
-            buffer,
-            &(result.module(0), result.zero_module(), 0),
-        )?);
-        for s in 1..max_s as u32 {
-            let d: Arc<FreeModuleHomomorphism<FreeModule<CC::Algebra>>> =
-                Load::load(buffer, &(result.module(s), result.module(s - 1), 0))?;
-            result.differentials.push(d);
-        }
-
-        let len = usize::load(buffer, &())?;
-        assert_eq!(len, max_s);
-
-        for s in 0..max_s as u32 {
-            let c: Arc<FreeModuleHomomorphism<CC::Module>> =
-                Load::load(buffer, &(result.module(s), result.complex().module(s), 0))?;
-            result.chain_maps.push(c);
-        }
-
-        result
-            .zero_module
-            .extend_by_zero(result.module(0).max_computed_degree());
-
-        Ok(result)
-    }
-}
-
 #[cfg(test)]
 mod test {
     use super::*;
@@ -802,4 +1072,22 @@ mod test {
         "#]]
         .assert_eq(&res.graded_dimension_string());
     }
+
+    #[test]
+    fn test_apply_quasi_inverse() {
+        let tempdir = tempfile::TempDir::new().unwrap();
+
+        let mut res = construct("S_2", Some(tempdir.path().into())).unwrap();
+        res.load_quasi_inverse = false;
+
+        res.compute_through_bidegree(8, 8);
+
+        assert!(res.differential(8).quasi_inverse(8).is_none());
+
+        let v = FpVector::new(res.prime(), res.module(7).dimension(8));
+        let mut w = FpVector::new(res.prime(), res.module(8).dimension(8));
+
+        assert!(res.apply_quasi_inverse(&mut [w.as_slice_mut()], 8, 8, &[v.as_slice()]));
+        assert!(w.is_zero());
+    }
 }
diff --git a/ext/src/utils.rs b/ext/src/utils.rs
index b2cebf083..821516f31 100644
--- a/ext/src/utils.rs
+++ b/ext/src/utils.rs
@@ -2,14 +2,15 @@ use crate::chain_complex::{ChainComplex, FiniteChainComplex, FreeChainComplex};
 use crate::resolution::Resolution;
 use crate::CCC;
 use algebra::module::{FiniteModule, Module};
-use algebra::{AlgebraType, SteenrodAlgebra};
-use saveload::Load;
+use algebra::{Algebra, AlgebraType, SteenrodAlgebra};
+use fp::prime::ValidPrime;
 
 use anyhow::{anyhow, Context};
+use byteorder::{LittleEndian, ReadBytesExt, WriteBytesExt};
 use serde_json::Value;
 
 use std::convert::{TryFrom, TryInto};
-use std::fs::File;
+use std::io::{Read, Write};
 use std::path::PathBuf;
 use std::sync::Arc;
 
@@ -116,10 +117,7 @@ impl<T: TryInto<AlgebraType>> TryFrom<(Value, T)> for Config {
 ///       `algebra` are as above.
 ///  - `save_file`: The save file for the module. If it points to an invalid save file, an error is
 ///    returned.
-pub fn construct<T, E>(
-    module_spec: T,
-    mut save_file: Option<File>,
-) -> anyhow::Result<Resolution<CCC>>
+pub fn construct<T, E>(module_spec: T, save_dir: Option<PathBuf>) -> anyhow::Result<Resolution<CCC>>
 where
     anyhow::Error: From<E>,
     T: TryInto<Config, Error = E>,
@@ -169,17 +167,7 @@ where
         chain_complex = Arc::new(yoneda);
     }
 
-    if let Some(f) = json["save_file"].as_str() {
-        save_file = Some(File::open(f)?);
-    }
-
-    Ok(match save_file {
-        Some(f) => {
-            let mut f = std::io::BufReader::new(f);
-            Resolution::load(&mut f, &chain_complex)?
-        }
-        None => Resolution::new(Arc::clone(&chain_complex)),
-    })
+    Resolution::new_with_save(chain_complex, save_dir)
 }
 
 pub fn load_module_json(name: &str) -> anyhow::Result<Value> {
@@ -226,7 +214,6 @@ pub fn print_resolution_color<C: FreeChainComplex, S: std::hash::BuildHasher>(
     max_s: u32,
     highlight: &std::collections::HashMap<(u32, i32), u32, S>,
 ) {
-    use std::io::Write;
     let stderr = std::io::stderr();
     let mut stderr = stderr.lock();
     for s in (0..max_s).rev() {
@@ -259,37 +246,37 @@ pub struct QueryModuleResult {
     pub bucket: thread_token::TokenBucket,
 }
 
-pub fn query_module(algebra: Option<AlgebraType>) -> anyhow::Result<QueryModuleResult> {
-    let module: Config = query::with_default("Module", "S_2", |s| match algebra {
+pub fn query_module_only(
+    prompt: &str,
+    algebra: Option<AlgebraType>,
+) -> anyhow::Result<Resolution<CCC>> {
+    let module: Config = query::with_default(prompt, "S_2", |s| match algebra {
         Some(algebra) => (s, algebra).try_into(),
         None => s.try_into(),
     });
 
-    // Clippy false positive
-    #[allow(clippy::redundant_closure)]
-    let save_file = query::optional("Resolution save file", |x| File::open(x));
+    let save_dir = query::optional(&format!("{prompt} save directory"), |x| {
+        core::result::Result::<PathBuf, std::convert::Infallible>::Ok(PathBuf::from(x))
+    });
+
+    construct(module, save_dir).context("Failed to load module from save file")
+}
+
+pub fn query_module(algebra: Option<AlgebraType>) -> anyhow::Result<QueryModuleResult> {
+    let resolution = query_module_only("Module", algebra)?;
 
     #[cfg(feature = "concurrent")]
     let bucket = query_bucket();
 
-    let resolution: Resolution<CCC> = match save_file {
-        Some(save_file) => {
-            construct(module, Some(save_file)).context("Failed to load module from save file")?
-        }
-        None => {
-            let max_s: u32 = query::with_default("Max s", "7", str::parse);
-            let max_n: i32 = query::with_default("Max n", "30", str::parse);
+    let max_s: u32 = query::with_default("Max s", "7", str::parse);
+    let max_n: i32 = query::with_default("Max n", "30", str::parse);
 
-            let resolution = construct(module, None).context("Failed to construct module")?;
-            #[cfg(not(feature = "concurrent"))]
-            resolution.compute_through_stem(max_s, max_n);
+    #[cfg(not(feature = "concurrent"))]
+    resolution.compute_through_stem(max_s, max_n);
 
-            #[cfg(feature = "concurrent")]
-            resolution.compute_through_stem_concurrent(max_s, max_n, &bucket);
+    #[cfg(feature = "concurrent")]
+    resolution.compute_through_stem_concurrent(max_s, max_n, &bucket);
 
-            resolution
-        }
-    };
     Ok(QueryModuleResult {
         resolution,
         #[cfg(feature = "concurrent")]
@@ -356,6 +343,168 @@ pub fn print_element(v: fp::vector::Slice, n: i32, s: u32) {
     }
 }
 
+pub fn write_header<A: Algebra>(
+    magic: u32,
+    algebra: &A,
+    p: ValidPrime,
+    s: u32,
+    t: i32,
+    buffer: &mut impl Write,
+) -> std::io::Result<()> {
+    buffer.write_u32::<LittleEndian>(magic)?;
+    buffer.write_u16::<LittleEndian>(algebra.magic())?;
+    buffer.write_u16::<LittleEndian>(*p as u16)?;
+    buffer.write_u32::<LittleEndian>(s)?;
+    buffer.write_i32::<LittleEndian>(t)
+}
+
+pub fn validate_header<A: Algebra>(
+    magic: u32,
+    algebra: &A,
+    p: ValidPrime,
+    s: u32,
+    t: i32,
+    buffer: &mut impl Read,
+) -> std::io::Result<()> {
+    use std::io::{Error, ErrorKind};
+
+    let data_magic = buffer.read_u32::<LittleEndian>()?;
+    if data_magic != magic {
+        return Err(Error::new(
+            ErrorKind::InvalidData,
+            format!("Invalid magic {data_magic:#010x}; expected {magic:#010x}"),
+        ));
+    }
+
+    let algebra_magic = buffer.read_u16::<LittleEndian>()?;
+    if algebra_magic != algebra.magic() {
+        return Err(Error::new(
+            ErrorKind::InvalidData,
+            format!(
+                "Invalid algebra magic {algebra_magic:#06x}; expected {:#06x}",
+                algebra.magic()
+            ),
+        ));
+    }
+
+    let data_p = buffer.read_u16::<LittleEndian>()? as u32;
+    if data_p != *p {
+        return Err(Error::new(
+            ErrorKind::InvalidData,
+            format!("Invalid prime {data_p}; expected {p}"),
+        ));
+    }
+
+    let data_s = buffer.read_u32::<LittleEndian>()?;
+    if data_s != s {
+        return Err(Error::new(
+            ErrorKind::InvalidData,
+            format!("Invalid s {data_s}; expected {s}"),
+        ));
+    }
+
+    let data_t = buffer.read_i32::<LittleEndian>()?;
+    if data_t != t {
+        return Err(Error::new(
+            ErrorKind::InvalidData,
+            format!("Invalid s {data_t}; expected {t}"),
+        ));
+    }
+
+    Ok(())
+}
+
+pub struct ChecksumWriter<T: Write> {
+    writer: T,
+    adler: adler::Adler32,
+}
+
+impl<T: Write> ChecksumWriter<T> {
+    pub fn new(writer: T) -> Self {
+        Self {
+            writer,
+            adler: adler::Adler32::new(),
+        }
+    }
+}
+
+/// We only implement the functions required and the ones we actually use.
+impl<T: Write> Write for ChecksumWriter<T> {
+    fn write(&mut self, buf: &[u8]) -> std::io::Result<usize> {
+        let bytes_written = self.writer.write(buf)?;
+        self.adler.write_slice(&buf[0..bytes_written]);
+        Ok(bytes_written)
+    }
+
+    fn flush(&mut self) -> std::io::Result<()> {
+        self.writer.flush()
+    }
+
+    fn write_all(&mut self, buf: &[u8]) -> std::io::Result<()> {
+        self.writer.write_all(buf)?;
+        self.adler.write_slice(buf);
+        Ok(())
+    }
+}
+
+impl<T: Write> std::ops::Drop for ChecksumWriter<T> {
+    fn drop(&mut self) {
+        if !std::thread::panicking() {
+            // We may not have finished writing, so the data is wrong. It should not be given a
+            // valid checksum
+            self.writer
+                .write_u32::<LittleEndian>(self.adler.checksum())
+                .unwrap();
+        }
+    }
+}
+
+pub struct ChecksumReader<T: Read> {
+    reader: T,
+    adler: adler::Adler32,
+}
+
+impl<T: Read> ChecksumReader<T> {
+    pub fn new(reader: T) -> Self {
+        Self {
+            reader,
+            adler: adler::Adler32::new(),
+        }
+    }
+}
+
+/// We only implement the functions required and the ones we actually use.
+impl<T: Read> Read for ChecksumReader<T> {
+    fn read(&mut self, buf: &mut [u8]) -> std::io::Result<usize> {
+        let bytes_read = self.reader.read(buf)?;
+        self.adler.write_slice(&buf[0..bytes_read]);
+        Ok(bytes_read)
+    }
+
+    fn read_exact(&mut self, buf: &mut [u8]) -> std::io::Result<()> {
+        self.reader.read_exact(buf)?;
+        self.adler.write_slice(buf);
+        Ok(())
+    }
+}
+
+impl<T: Read> std::ops::Drop for ChecksumReader<T> {
+    fn drop(&mut self) {
+        if !std::thread::panicking() {
+            // If we are panicking, we may not have read everything, and panic in panic
+            // is bad.
+            assert_eq!(
+                self.adler.checksum(),
+                self.reader.read_u32::<LittleEndian>().unwrap(),
+                "Invalid file checksum"
+            );
+            let mut buf = [0];
+            // Check EOF
+            assert_eq!(self.reader.read(&mut buf).unwrap(), 0, "EOF not reached");
+        }
+    }
+}
+
 #[cfg(test)]
 mod test {
     use super::*;
diff --git a/web_ext/sseq_gui/examples/charts.rs b/web_ext/sseq_gui/examples/charts.rs
index 70fabc352..f456aa7f7 100644
--- a/web_ext/sseq_gui/examples/charts.rs
+++ b/web_ext/sseq_gui/examples/charts.rs
@@ -6,7 +6,7 @@ use sseq_gui::actions::SseqChoice;
 use sseq_gui::sseq::SseqWrapper;
 use std::fs::File;
 use std::io::{BufRead, BufReader};
-use std::path::Path;
+use std::path::{Path, PathBuf};
 use std::time::Instant;
 
 const TWO: ValidPrime = ValidPrime::new(2);
@@ -19,7 +19,9 @@ fn main() -> anyhow::Result<()> {
 
     // Clippy false positive
     #[allow(clippy::redundant_closure)]
-    let save_file: Option<File> = query::optional("Resolution save file", |s| File::open(s));
+    let save_file: Option<PathBuf> = query::optional("Resolution save file", |s| {
+        std::result::Result::<PathBuf, std::convert::Infallible>::Ok(PathBuf::from(s))
+    });
 
     #[cfg(feature = "concurrent")]
     let bucket = ext::utils::query_bucket();


Test patch:
diff --git a/ext/tests/milnor_vs_adem.rs b/ext/tests/milnor_vs_adem.rs
index 4d927b51d..d8ffd08bc 100644
--- a/ext/tests/milnor_vs_adem.rs
+++ b/ext/tests/milnor_vs_adem.rs
@@ -17,8 +17,18 @@ fn compare(#[case] module_name: &str, #[case] max_degree: i32) {
     let a = construct((module_name, "adem"), None).unwrap();
     let b = construct((module_name, "milnor"), None).unwrap();
 
-    a.compute_through_bidegree(max_degree as u32, max_degree);
-    b.compute_through_bidegree(max_degree as u32, max_degree);
+    #[cfg(feature = "concurrent")]
+    {
+        let bucket = thread_token::TokenBucket::new(2.try_into().unwrap());
+        a.compute_through_bidegree_concurrent(max_degree as u32, max_degree, &bucket);
+        b.compute_through_bidegree_concurrent(max_degree as u32, max_degree, &bucket);
+    }
+
+    #[cfg(not(feature = "concurrent"))]
+    {
+        a.compute_through_bidegree(max_degree as u32, max_degree);
+        b.compute_through_bidegree(max_degree as u32, max_degree);
+    }
 
     assert_eq!(a.graded_dimension_string(), b.graded_dimension_string());
 }
diff --git a/ext/tests/save_load_resolution.rs b/ext/tests/save_load_resolution.rs
index a91edb622..c96569da7 100644
--- a/ext/tests/save_load_resolution.rs
+++ b/ext/tests/save_load_resolution.rs
@@ -1,30 +1,91 @@
 use algebra::module::homomorphism::ModuleHomomorphism;
 use ext::chain_complex::{ChainComplex, FreeChainComplex};
-use ext::{resolution::Resolution, utils::construct};
-use saveload::{Load, Save};
-use std::io::{Cursor, Read, Seek, SeekFrom};
+use ext::resolution::{Resolution, SaveData};
+use ext::utils::construct;
+use ext::CCC;
+
+use std::path::{Path, PathBuf};
+
+fn set_readonly(p: &Path, readonly: bool) {
+    let mut perm = p.metadata().unwrap().permissions();
+    perm.set_readonly(readonly);
+    std::fs::set_permissions(p, perm).unwrap();
+}
+
+fn lock_tempdir(dir: &Path) {
+    let mut dir: PathBuf = dir.into();
+    for kind in SaveData::resolution_data() {
+        dir.push(format!("{}s", kind.name()));
+        set_readonly(&dir, true);
+        dir.pop();
+    }
+    set_readonly(&dir, true);
+}
+
+/// Should unlock after the test so that cleanup can be performed
+fn unlock_tempdir(dir: &Path) {
+    set_readonly(dir, false);
+
+    let mut dir: PathBuf = dir.into();
+    for kind in SaveData::resolution_data() {
+        dir.push(format!("{}s", kind.name()));
+        set_readonly(&dir, false);
+        dir.pop();
+    }
+}
+
+#[test]
+#[should_panic]
+fn test_tempdir_lock() {
+    let tempdir = tempfile::TempDir::new().unwrap();
+    let resolution1 = construct("S_2", Some(tempdir.path().into())).unwrap();
+    resolution1.compute_through_bidegree(5, 5);
+
+    lock_tempdir(tempdir.path());
+    resolution1.compute_through_bidegree(6, 6);
+}
+
+#[test]
+fn test_tempdir_unlock() {
+    let tempdir = tempfile::TempDir::new().unwrap();
+    let resolution1 = construct("S_2", Some(tempdir.path().into())).unwrap();
+    resolution1.compute_through_bidegree(5, 5);
+
+    lock_tempdir(tempdir.path());
+    unlock_tempdir(tempdir.path());
+    resolution1.compute_through_bidegree(6, 6);
+}
 
 #[test]
 fn test_save_load() {
-    let resolution1 = construct("S_2", None).unwrap();
-    resolution1.compute_through_bidegree(10, 6);
-    resolution1.compute_through_bidegree(6, 10);
+    let tempdir = tempfile::TempDir::new().unwrap();
+    let mut resolution1 = construct("S_2", Some(tempdir.path().into())).unwrap();
 
-    let mut cursor: Cursor<Vec<u8>> = Cursor::new(Vec::new());
-    resolution1.save(&mut cursor).unwrap();
+    #[cfg(feature = "concurrent")]
+    let bucket = thread_token::TokenBucket::new(2.try_into().unwrap());
 
-    cursor.seek(SeekFrom::Start(0)).unwrap();
+    let resolve = |resolution: &Resolution<CCC>, s, t| {
+        #[cfg(feature = "concurrent")]
+        resolution.compute_through_bidegree_concurrent(s, t, &bucket);
+        #[cfg(not(feature = "concurrent"))]
+        resolution.compute_through_bidegree(s, t);
+    };
 
-    let resolution2 = Resolution::load(&mut cursor, &resolution1.complex()).unwrap();
-    assert_eq!(0, cursor.bytes().count());
+    resolve(&resolution1, 10, 6);
+    resolve(&resolution1, 6, 10);
+    resolution1.should_save = false;
 
-    assert_eq!(
-        resolution1.graded_dimension_string(),
-        resolution2.graded_dimension_string()
-    );
+    let mut resolution2 = construct("S_2", Some(tempdir.path().into())).unwrap();
 
-    resolution1.compute_through_bidegree(20, 20);
-    resolution2.compute_through_bidegree(20, 20);
+    // Check that we are not writing anything new.
+    lock_tempdir(tempdir.path());
+    resolve(&resolution2, 10, 6);
+    resolve(&resolution2, 6, 10);
+
+    resolution2.should_save = false;
+
+    resolve(&resolution1, 20, 20);
+    resolve(&resolution2, 20, 20);
 
     assert_eq!(
         resolution1.graded_dimension_string(),
@@ -35,20 +96,41 @@ fn test_save_load() {
         resolution1.differential(5).quasi_inverse(7),
         resolution2.differential(5).quasi_inverse(7)
     );
+    unlock_tempdir(tempdir.path());
+}
+
+#[test]
+#[should_panic]
+fn wrong_algebra() {
+    let tempdir = tempfile::TempDir::new().unwrap();
+    let resolution1 = construct("S_2", Some(tempdir.path().into())).unwrap();
+    resolution1.compute_through_bidegree(2, 2);
+
+    let resolution2 = construct("S_2@milnor", Some(tempdir.path().into())).unwrap();
+    resolution2.compute_through_bidegree(2, 2);
 }
 
 #[test]
 fn test_save_load_stem() {
-    let resolution1 = construct("S_2", None).unwrap();
-    resolution1.compute_through_stem(10, 10);
+    let tempdir = tempfile::TempDir::new().unwrap();
+
+    #[cfg(feature = "concurrent")]
+    let bucket = thread_token::TokenBucket::new(2.try_into().unwrap());
 
-    let mut cursor: Cursor<Vec<u8>> = Cursor::new(Vec::new());
-    resolution1.save(&mut cursor).unwrap();
+    let resolution1 = construct("S_2", Some(tempdir.path().into())).unwrap();
+
+    #[cfg(feature = "concurrent")]
+    resolution1.compute_through_stem_concurrent(10, 10, &bucket);
+    #[cfg(not(feature = "concurrent"))]
+    resolution1.compute_through_stem(10, 10);
 
-    cursor.seek(SeekFrom::Start(0)).unwrap();
+    let resolution2 = construct("S_2", Some(tempdir.path().into())).unwrap();
+    lock_tempdir(tempdir.path());
 
-    let resolution2 = Resolution::load(&mut cursor, &resolution1.complex()).unwrap();
-    assert_eq!(0, cursor.bytes().count());
+    #[cfg(feature = "concurrent")]
+    resolution2.compute_through_stem_concurrent(10, 10, &bucket);
+    #[cfg(not(feature = "concurrent"))]
+    resolution2.compute_through_stem(10, 10);
 
     assert_eq!(
         resolution1.graded_dimension_string(),
@@ -59,26 +141,66 @@ fn test_save_load_stem() {
         resolution1.differential(5).quasi_inverse(7),
         resolution2.differential(5).quasi_inverse(7)
     );
+    unlock_tempdir(tempdir.path());
 }
 
 #[test]
 fn test_save_load_resume() {
-    let resolution1 = construct("S_2", None).unwrap();
-    resolution1.compute_through_stem(8, 14);
+    let tempdir = tempfile::TempDir::new().unwrap();
 
-    let mut cursor: Cursor<Vec<u8>> = Cursor::new(Vec::new());
-    resolution1.save(&mut cursor).unwrap();
+    #[cfg(feature = "concurrent")]
+    let bucket = thread_token::TokenBucket::new(2.try_into().unwrap());
+    let resolve = |resolution: &Resolution<CCC>, s, n| {
+        #[cfg(feature = "concurrent")]
+        resolution.compute_through_stem_concurrent(s, n, &bucket);
+        #[cfg(not(feature = "concurrent"))]
+        resolution.compute_through_stem(s, n);
+    };
 
-    cursor.seek(SeekFrom::Start(0)).unwrap();
+    let resolution1 = construct("S_2", Some(tempdir.path().into())).unwrap();
+    resolve(&resolution1, 8, 14);
 
-    let resolution2 = Resolution::load(&mut cursor, &resolution1.complex()).unwrap();
-    assert_eq!(0, cursor.bytes().count());
+    let resolution2 = construct("S_2", Some(tempdir.path().into())).unwrap();
+    lock_tempdir(tempdir.path());
+    resolve(&resolution2, 8, 14);
+    unlock_tempdir(tempdir.path());
 
-    resolution1.compute_through_stem(5, 19);
-    resolution2.compute_through_stem(5, 19);
+    resolve(&resolution1, 5, 19);
+    lock_tempdir(tempdir.path());
+    resolve(&resolution2, 5, 19);
 
     assert_eq!(
         resolution1.graded_dimension_string(),
         resolution2.graded_dimension_string()
     );
+    unlock_tempdir(tempdir.path());
+}
+
+#[test]
+#[should_panic]
+fn test_checksum() {
+    use std::fs::OpenOptions;
+    use std::io::{Seek, SeekFrom, Write};
+
+    let tempdir = tempfile::TempDir::new().unwrap();
+
+    construct("S_2", Some(tempdir.path().into()))
+        .unwrap()
+        .compute_through_bidegree(2, 2);
+
+    let mut path = tempdir.path().to_owned();
+    path.push("differentials/2_2_differential");
+
+    let mut file = OpenOptions::new()
+        .read(true)
+        .write(true)
+        .open(path)
+        .unwrap();
+
+    file.seek(SeekFrom::Start(41)).unwrap();
+    file.write_all(&[1]).unwrap();
+
+    construct("S_2", Some(tempdir.path().into()))
+        .unwrap()
+        .compute_through_bidegree(2, 2);
 }