ext/
secondary.rs

1use std::{io, sync::Arc};
2
3use algebra::{
4    Algebra,
5    module::{
6        FreeModule, Module,
7        homomorphism::{FreeModuleHomomorphism, ModuleHomomorphism},
8    },
9    pair_algebra::PairAlgebra,
10};
11use bivec::BiVec;
12use byteorder::{LittleEndian, ReadBytesExt, WriteBytesExt};
13use dashmap::DashMap;
14use fp::{
15    matrix::Matrix,
16    prime::ValidPrime,
17    vector::{FpSlice, FpSliceMut, FpVector},
18};
19use itertools::Itertools;
20use maybe_rayon::prelude::*;
21use once::OnceBiVec;
22use sseq::coordinates::{Bidegree, BidegreeElement, BidegreeGenerator, BidegreeRange};
23use tracing::Level;
24
25use crate::{
26    chain_complex::{ChainComplex, ChainHomotopy, FreeChainComplex},
27    resolution_homomorphism::ResolutionHomomorphism,
28    save::{SaveDirectory, SaveFile, SaveKind},
29};
30
31pub static LAMBDA_BIDEGREE: Bidegree = Bidegree::n_s(0, 1);
32
33pub type CompositeData<A> = Vec<(
34    u32,
35    Arc<FreeModuleHomomorphism<FreeModule<A>>>,
36    Arc<FreeModuleHomomorphism<FreeModule<A>>>,
37)>;
38
39/// A homotopy of a map A -> M of pair modules. We assume this map does not hit generators.
40pub struct SecondaryComposite<A: PairAlgebra> {
41    target: Arc<FreeModule<A>>,
42    degree: i32,
43    /// The component of the map on the R_B portion.
44    /// gen_deg -> gen_idx -> coefficient
45    composite: BiVec<Vec<A::Element>>,
46}
47
48impl<A: PairAlgebra> SecondaryComposite<A> {
49    pub fn algebra(&self) -> Arc<A> {
50        self.target.algebra()
51    }
52
53    pub fn new(target: Arc<FreeModule<A>>, degree: i32, hit_generator: bool) -> Self {
54        let algebra = target.algebra();
55        let min_degree = target.min_degree();
56
57        let mut composite = BiVec::with_capacity(min_degree, degree);
58
59        let end = if hit_generator { degree + 1 } else { degree };
60        for t_ in min_degree..end {
61            let num_gens = target.number_of_gens_in_degree(t_);
62            let mut c = Vec::with_capacity(num_gens);
63            c.resize_with(num_gens, || algebra.new_pair_element(degree - t_));
64            composite.push(c);
65        }
66
67        Self {
68            target,
69            degree,
70            composite,
71        }
72    }
73
74    pub fn to_bytes(&self, buffer: &mut impl io::Write) -> io::Result<()> {
75        let algebra = self.target.algebra();
76        for composites in self.composite.iter() {
77            buffer.write_u64::<LittleEndian>(composites.len() as u64)?;
78            for composite in composites {
79                algebra.element_to_bytes(composite, buffer)?;
80            }
81        }
82        Ok(())
83    }
84
85    pub fn from_bytes(
86        target: Arc<FreeModule<A>>,
87        degree: i32,
88        hit_generator: bool,
89        buffer: &mut impl io::Read,
90    ) -> io::Result<Self> {
91        let min_degree = target.min_degree();
92        let algebra = target.algebra();
93        let mut composite = BiVec::with_capacity(min_degree, degree);
94
95        let end = if hit_generator { degree + 1 } else { degree };
96        for t in min_degree..end {
97            let num_gens = buffer.read_u64::<LittleEndian>()? as usize;
98            let mut c = Vec::with_capacity(num_gens);
99            for _ in 0..num_gens {
100                c.push(algebra.element_from_bytes(degree - t, buffer)?);
101            }
102            composite.push(c);
103        }
104
105        Ok(Self {
106            target,
107            degree,
108            composite,
109        })
110    }
111
112    pub fn finalize(&mut self) {
113        for r in self.composite.iter_mut() {
114            for r in r.iter_mut() {
115                A::finalize_element(r);
116            }
117        }
118    }
119
120    pub fn add_composite(
121        &mut self,
122        coeff: u32,
123        gen_degree: i32,
124        gen_idx: usize,
125        d1: &FreeModuleHomomorphism<FreeModule<A>>,
126        d0: &FreeModuleHomomorphism<FreeModule<A>>,
127    ) {
128        assert!(Arc::ptr_eq(&d1.target(), &d0.source()));
129        assert!(Arc::ptr_eq(&d0.target(), &self.target));
130
131        let middle = d1.target();
132        let dx = d1.output(gen_degree, gen_idx);
133        let algebra = self.algebra();
134
135        for (gen_deg1, gen_idx1, op_deg1, slice1) in
136            middle.iter_slices(gen_degree - d1.degree_shift(), dx.as_slice())
137        {
138            if slice1.is_zero() {
139                continue;
140            }
141            if gen_deg1 < d0.degree_shift() {
142                continue;
143            }
144            let dy = d0.output(gen_deg1, gen_idx1);
145
146            for (gen_deg2, gen_idx2, op_deg2, slice2) in self
147                .target
148                .iter_slices(gen_deg1 - d0.degree_shift(), dy.as_slice())
149            {
150                if slice2.is_zero() {
151                    continue;
152                }
153                algebra.sigma_multiply(
154                    &mut self.composite[gen_deg2][gen_idx2],
155                    coeff,
156                    op_deg1,
157                    slice1,
158                    op_deg2,
159                    slice2,
160                )
161            }
162        }
163    }
164
165    pub fn act(&self, mut result: FpSliceMut, coeff: u32, op_degree: i32, op: FpSlice) {
166        let algebra = self.algebra();
167        for (gen_deg, row) in self.composite.iter_enum() {
168            let module_op_deg = self.degree - gen_deg;
169            for (gen_idx, c) in row.iter().enumerate() {
170                if gen_deg > self.target.max_computed_degree() {
171                    // If we are resolving up to a stem then the target might be missing some
172                    // degrees. This is fine but we want to assert that c is zero.
173                    assert!(A::element_is_zero(c));
174                    continue;
175                }
176
177                let offset =
178                    self.target
179                        .generator_offset(self.degree + op_degree - 1, gen_deg, gen_idx);
180                let len = algebra.dimension(module_op_deg + op_degree - 1);
181
182                algebra.a_multiply(
183                    result.slice_mut(offset, offset + len),
184                    coeff,
185                    op_degree,
186                    op,
187                    module_op_deg,
188                    c,
189                );
190            }
191        }
192    }
193}
194
195pub struct SecondaryHomotopy<A: PairAlgebra> {
196    pub source: Arc<FreeModule<A>>,
197    pub target: Arc<FreeModule<A>>,
198    /// output_t = input_t - shift_t
199    pub shift_t: i32,
200
201    /// gen_deg -> gen_idx -> composite
202    pub(crate) composites: OnceBiVec<Vec<SecondaryComposite<A>>>,
203
204    /// gen_deg -> gen_idx -> homotopy
205    pub homotopies: FreeModuleHomomorphism<FreeModule<A>>,
206
207    hit_generator: bool,
208}
209
210impl<A: PairAlgebra + Send + Sync> SecondaryHomotopy<A> {
211    pub fn new(
212        source: Arc<FreeModule<A>>,
213        target: Arc<FreeModule<A>>,
214        shift_t: i32,
215        hit_generator: bool,
216    ) -> Self {
217        Self {
218            composites: OnceBiVec::new(std::cmp::max(
219                source.min_degree(),
220                target.min_degree() + shift_t,
221            )),
222            homotopies: FreeModuleHomomorphism::new(
223                Arc::clone(&source),
224                Arc::clone(&target),
225                shift_t + 1,
226            ),
227            source,
228            target,
229            shift_t,
230            hit_generator,
231        }
232    }
233
234    /// Add composites up to and including the specified degree
235    #[tracing::instrument(skip(self, maps, dir), fields(source = %self.source, target = %self.target))]
236    pub fn add_composite(&self, s: i32, degree: i32, maps: CompositeData<A>, dir: &SaveDirectory) {
237        for (_, d1, d0) in &maps {
238            assert!(Arc::ptr_eq(&d1.target(), &d0.source()));
239            assert!(Arc::ptr_eq(&d0.target(), &self.target));
240            assert_eq!(d1.degree_shift() + d0.degree_shift(), self.shift_t);
241        }
242
243        let tracing_span = tracing::Span::current();
244        let f = |t, idx| {
245            let _tracing_guard = tracing_span.enter();
246            let g = BidegreeGenerator::s_t(s, t, idx);
247            let save_file = SaveFile {
248                algebra: self.target.algebra(),
249                kind: SaveKind::SecondaryComposite,
250                b: g.degree(),
251                idx: Some(g.idx()),
252            };
253            if let Some(dir) = dir.read()
254                && let Some(mut f) = save_file.open_file(dir.to_owned())
255            {
256                return SecondaryComposite::from_bytes(
257                    Arc::clone(&self.target),
258                    g.t() - self.shift_t,
259                    self.hit_generator,
260                    &mut f,
261                )
262                .unwrap();
263            }
264
265            let mut composite = SecondaryComposite::new(
266                Arc::clone(&self.target),
267                g.t() - self.shift_t,
268                self.hit_generator,
269            );
270
271            tracing::info_span!("Computing composite", %g).in_scope(|| {
272                for (coef, d1, d0) in &maps {
273                    composite.add_composite(*coef, g.t(), g.idx(), d1, d0);
274                }
275                composite.finalize();
276            });
277
278            if let Some(dir) = dir.write() {
279                let mut f = save_file.create_file(dir.to_owned(), false);
280                composite.to_bytes(&mut f).unwrap();
281            }
282
283            composite
284        };
285
286        self.composites.maybe_par_extend(degree, |t| {
287            (0..self.source.number_of_gens_in_degree(t))
288                .into_maybe_par_iter()
289                .map(|i| f(t, i))
290                .collect()
291        });
292    }
293
294    /// Compute the image of an element in the source under the homotopy, writing the result in
295    /// `result`. It is assumed that the coefficients of generators are zero in `op`.
296    ///
297    /// # Arguments
298    ///  - full: Whether to include the action of the homotopy part as well
299    pub fn act(
300        &self,
301        mut result: FpSliceMut,
302        coeff: u32,
303        elt_degree: i32,
304        elt: FpSlice,
305        full: bool,
306    ) {
307        for (gen_deg, gen_idx, op_deg, slice) in self.source.iter_slices(elt_degree, elt) {
308            if gen_deg < self.composites.min_degree() {
309                continue;
310            }
311            // This is actually necessary. We don't have the homotopies on the
312            // generators at the edge of the resolution, but we don't need them since they never
313            // get hit.
314            if slice.is_zero() {
315                continue;
316            }
317            self.composites[gen_deg][gen_idx].act(result.copy(), coeff, op_deg, slice);
318        }
319
320        if full {
321            self.homotopies.apply(result, coeff, elt_degree, elt);
322        }
323    }
324
325    pub fn composite(&self, gen_deg: i32, gen_idx: usize) -> &SecondaryComposite<A> {
326        &self.composites[gen_deg][gen_idx]
327    }
328}
329
330/// Logic that is common to all secondary lifts.
331///
332/// When lifting a thing to its secondary version, often what we have to do is to specify an
333/// explicit homotopy to witnesses that some equation holds. For example, to lift a chain complex,
334/// we need a homotopy witnessing the fact that $d^2 \simeq 0$. This homotopy in turn is required
335/// to satisfy certain recursive relations.
336///
337/// To specify this lifting problem, one needs to supply two pieces of data. First is the equation
338/// that we are trying to witness, which is usually of the form
339///
340/// $$ \sum_i c_i f_i g_i = 0, $$
341///
342/// where $f_i$ and $g_i$ are free module homomorphisms and $c_i$ are constants. This is specified
343/// by [`SecondaryLift::composite`].
344///
345/// The next is a compatibility equation, which restricts the λ part of the null-homotopy, and is
346/// usually of the form
347///
348/// $$ dh = hd + \mathrm{stuff} $$
349///
350/// The λ part of $hd + \mathrm{stuff}$ is known as the intermediate data, and is what
351/// [`SecondaryLift::compute_intermediate`] returns.
352pub trait SecondaryLift: Sync + Sized {
353    type Algebra: PairAlgebra;
354    type Source: FreeChainComplex<Algebra = Self::Algebra>;
355    type Target: FreeChainComplex<Algebra = Self::Algebra>;
356    type Underlying;
357
358    /// Whether the composite can hit generators. This is true for `SecondaryChainHomotopy` and
359    /// false for the rest. This is important because for [`SecondaryResolution`], we don't
360    /// actually know all the generators if we resolve up to a stem. So in composites for
361    /// [`SecondaryResolution`], we need to ignore target generators of the same degree uniformly.
362    const HIT_GENERATOR: bool = false;
363
364    fn underlying(&self) -> Arc<Self::Underlying>;
365    fn algebra(&self) -> Arc<Self::Algebra>;
366    fn prime(&self) -> ValidPrime {
367        self.algebra().prime()
368    }
369
370    fn source(&self) -> Arc<Self::Source>;
371    fn target(&self) -> Arc<Self::Target>;
372    fn shift(&self) -> Bidegree;
373
374    fn max(&self) -> BidegreeRange<'_, Self>;
375
376    fn homotopies(&self) -> &OnceBiVec<SecondaryHomotopy<Self::Algebra>>;
377    fn intermediates(&self) -> &DashMap<BidegreeGenerator, FpVector>;
378
379    fn save_dir(&self) -> &SaveDirectory;
380
381    fn compute_intermediate(&self, g: BidegreeGenerator) -> FpVector;
382    fn composite(&self, s: i32) -> CompositeData<Self::Algebra>;
383
384    #[tracing::instrument(skip(self))]
385    fn initialize_homotopies(&self) {
386        let shift = self.shift();
387        let max = self.max();
388
389        self.homotopies().extend(max.s() - 1, |s| {
390            SecondaryHomotopy::new(
391                self.source().module(s),
392                self.target().module(s - shift.s()),
393                shift.t(),
394                Self::HIT_GENERATOR,
395            )
396        });
397    }
398
399    #[tracing::instrument(skip(self))]
400    fn compute_composites(&self) {
401        let tracing_span = tracing::Span::current();
402        let f = |s| {
403            let _tracing_guard = tracing_span.enter();
404            self.homotopies()[s].add_composite(
405                s,
406                self.max().t(s) - 1,
407                self.composite(s),
408                self.save_dir(),
409            );
410        };
411
412        self.homotopies().range().into_maybe_par_iter().for_each(f);
413    }
414
415    #[tracing::instrument(skip(self), ret(Display, level = Level::DEBUG), fields(%g))]
416    fn get_intermediate(&self, g: BidegreeGenerator) -> FpVector {
417        if let Some((_, v)) = self.intermediates().remove(&g) {
418            return v;
419        }
420
421        let save_file = SaveFile {
422            algebra: self.algebra(),
423            kind: SaveKind::SecondaryIntermediate,
424            b: g.degree(),
425            idx: Some(g.idx()),
426        };
427
428        if let Some(dir) = self.save_dir().read()
429            && let Some(mut f) = save_file.open_file(dir.to_owned())
430        {
431            // The target dimension can depend on whether we resolved to stem
432            let dim = f.read_u64::<LittleEndian>().unwrap() as usize;
433            return FpVector::from_bytes(self.prime(), dim, &mut f).unwrap();
434        }
435
436        let result = self.compute_intermediate(g);
437
438        if let Some(dir) = self.save_dir().write() {
439            let mut f = save_file.create_file(dir.to_owned(), false);
440            f.write_u64::<LittleEndian>(result.len() as u64).unwrap();
441            result.to_bytes(&mut f).unwrap();
442        }
443
444        result
445    }
446
447    #[tracing::instrument(skip(self))]
448    fn compute_partial(&self, s: i32) {
449        self.initialize_homotopies();
450        let homotopies = self.homotopies();
451        let tracing_span = tracing::Span::current();
452
453        if s < homotopies.min_degree() {
454            eprintln!(
455                "Computing partial for s = {s} when minimum degree is {}",
456                homotopies.min_degree()
457            );
458            return;
459        }
460
461        homotopies[s].add_composite(s, self.max().t(s) - 1, self.composite(s), self.save_dir());
462
463        if let Some(homotopy) = homotopies.get(s + 1) {
464            (0..self.max().t(s + 1))
465                .into_maybe_par_iter()
466                .for_each(|t| {
467                    (0..homotopy.source.number_of_gens_in_degree(t))
468                        .into_maybe_par_iter()
469                        .for_each(|i| {
470                            let _tracing_guard = tracing_span.enter();
471                            self.get_intermediate(BidegreeGenerator::s_t(s + 1, t, i));
472                        })
473                });
474        }
475    }
476
477    #[tracing::instrument(skip(self))]
478    fn compute_intermediates(&self) {
479        let tracing_span = tracing::Span::current();
480        let f = |g: BidegreeGenerator| {
481            let _tracing_guard = tracing_span.enter();
482
483            // If we already have homotopies, we don't need to compute intermediate
484            if self.homotopies()[g.s()].homotopies.next_degree() >= g.t() {
485                return;
486            }
487            // Check if we have a saved homotopy
488            if let Some(dir) = self.save_dir().read() {
489                let save_file = SaveFile {
490                    algebra: self.algebra(),
491                    kind: SaveKind::SecondaryHomotopy,
492                    b: g.degree(),
493                    idx: None,
494                };
495
496                if save_file.exists(dir.to_owned()) {
497                    return;
498                }
499            }
500            self.intermediates().insert(g, self.get_intermediate(g));
501        };
502
503        self.homotopies()
504            .maybe_par_iter()
505            .skip(1)
506            .for_each(|(s, homotopy)| {
507                homotopy
508                    .composites
509                    .range()
510                    .into_maybe_par_iter()
511                    .for_each(|t| {
512                        (0..homotopy.source.number_of_gens_in_degree(t))
513                            .into_maybe_par_iter()
514                            .for_each(|i| f(BidegreeGenerator::s_t(s, t, i)))
515                    })
516            })
517    }
518
519    #[tracing::instrument(skip(self), fields(%b))]
520    fn compute_homotopy_step(&self, b: Bidegree) -> std::ops::Range<i32> {
521        let homotopy = &self.homotopies()[b.s()];
522        if homotopy.homotopies.next_degree() > b.t() {
523            return b.t()..b.t() + 1;
524        }
525        let p = self.prime();
526        let shift = self.shift();
527        let target_b = b - shift - Bidegree::s_t(0, 1);
528
529        let d = self.source().differential(b.s());
530        let source = self.source().module(b.s());
531        let target = self.target();
532        let num_gens = source.number_of_gens_in_degree(b.t());
533        let target_dim = target.module(target_b.s()).dimension(target_b.t());
534
535        if let Some(dir) = self.save_dir().read() {
536            let save_file = SaveFile {
537                algebra: self.algebra(),
538                kind: SaveKind::SecondaryHomotopy,
539                b,
540                idx: None,
541            };
542
543            if let Some(mut f) = save_file.open_file(dir.to_owned()) {
544                let mut results = Vec::with_capacity(num_gens);
545                for _ in 0..num_gens {
546                    results.push(FpVector::from_bytes(p, target_dim, &mut f).unwrap());
547                }
548                return self.homotopies()[b.s()]
549                    .homotopies
550                    .add_generators_from_rows_ooo(b.t(), results);
551            }
552        }
553
554        let tracing_span = tracing::Span::current();
555        let get_intermediate = |i| {
556            let _tracing_guard = tracing_span.enter();
557
558            let g = BidegreeGenerator::new(b, i);
559            let mut v = self.get_intermediate(g);
560            if g.s() > shift.s() + 1 {
561                self.homotopies()[g.s() - 1].homotopies.apply(
562                    v.as_slice_mut(),
563                    1,
564                    g.t(),
565                    d.output(g.t(), g.idx()).as_slice(),
566                );
567            }
568            v
569        };
570
571        let mut intermediates: Vec<FpVector> = (0..num_gens)
572            .into_maybe_par_iter()
573            .map(get_intermediate)
574            .collect();
575
576        let mut results = vec![FpVector::new(p, target_dim); num_gens];
577
578        assert!(target.apply_quasi_inverse(&mut results, target_b, &intermediates,));
579
580        if b.s() == shift.s() + 1 {
581            // Check that we indeed had a lift
582            let d = target.differential(target_b.s());
583            for (src, tgt) in std::iter::zip(&results, &mut intermediates) {
584                d.apply(tgt.as_slice_mut(), p - 1, target_b.t(), src.as_slice());
585                assert!(
586                    tgt.is_zero(),
587                    "secondary: Failed to lift at {b}. This likely indicates an invalid input."
588                );
589            }
590        }
591
592        if let Some(dir) = self.save_dir().write() {
593            let save_file = SaveFile {
594                algebra: self.algebra(),
595                kind: SaveKind::SecondaryHomotopy,
596                b,
597                idx: None,
598            };
599
600            let mut f = save_file.create_file(dir.to_owned(), false);
601            for row in &results {
602                row.to_bytes(&mut f).unwrap();
603            }
604            drop(f);
605
606            let mut save_file = SaveFile {
607                algebra: self.algebra(),
608                kind: SaveKind::SecondaryIntermediate,
609                b,
610                idx: None,
611            };
612
613            for i in 0..num_gens {
614                save_file.idx = Some(i);
615                save_file.delete_file(dir.to_owned()).unwrap();
616            }
617        }
618
619        homotopy
620            .homotopies
621            .add_generators_from_rows_ooo(b.t(), results)
622    }
623
624    #[tracing::instrument(skip(self))]
625    fn compute_homotopies(&self) {
626        let shift = self.shift();
627
628        // When s = shift_s, the homotopies are just zero
629        {
630            let h = &self.homotopies()[shift.s()];
631            h.homotopies.extend_by_zero(h.composites.max_degree());
632        }
633
634        let min_t = self.homotopies()[shift.s()].homotopies.min_degree();
635        let s_range = self.homotopies().range();
636        let min = Bidegree::s_t(s_range.start + 1, min_t);
637        let max = self.max().restrict(s_range.end);
638        sseq::coordinates::iter_s_t(&|b| self.compute_homotopy_step(b), min, max);
639    }
640
641    #[tracing::instrument(skip(self))]
642    fn extend_all(&self) {
643        self.initialize_homotopies();
644        self.compute_composites();
645        self.compute_intermediates();
646        self.compute_homotopies();
647    }
648}
649
650pub struct SecondaryResolution<CC: FreeChainComplex>
651where
652    CC::Algebra: PairAlgebra,
653{
654    underlying: Arc<CC>,
655    /// s -> t -> idx -> homotopy
656    pub(crate) homotopies: OnceBiVec<SecondaryHomotopy<CC::Algebra>>,
657    intermediates: DashMap<BidegreeGenerator, FpVector>,
658}
659
660impl<CC: FreeChainComplex> SecondaryLift for SecondaryResolution<CC>
661where
662    CC::Algebra: PairAlgebra,
663{
664    type Algebra = CC::Algebra;
665    type Source = CC;
666    type Target = CC;
667    type Underlying = CC;
668
669    fn underlying(&self) -> Arc<CC> {
670        Arc::clone(&self.underlying)
671    }
672
673    fn algebra(&self) -> Arc<Self::Algebra> {
674        self.underlying.algebra()
675    }
676
677    fn source(&self) -> Arc<Self::Source> {
678        Arc::clone(&self.underlying)
679    }
680
681    fn target(&self) -> Arc<Self::Target> {
682        Arc::clone(&self.underlying)
683    }
684
685    fn shift(&self) -> Bidegree {
686        Bidegree::s_t(2, 0)
687    }
688
689    fn max(&self) -> BidegreeRange<'_, Self> {
690        BidegreeRange::new(
691            self,
692            self.underlying.next_homological_degree(),
693            &|selff, s| {
694                std::cmp::min(
695                    selff.underlying.module(s).max_computed_degree(),
696                    selff.underlying.module(s - 2).max_computed_degree() + 1,
697                ) + 1
698            },
699        )
700    }
701
702    fn homotopies(&self) -> &OnceBiVec<SecondaryHomotopy<CC::Algebra>> {
703        &self.homotopies
704    }
705
706    fn intermediates(&self) -> &DashMap<BidegreeGenerator, FpVector> {
707        &self.intermediates
708    }
709
710    fn save_dir(&self) -> &SaveDirectory {
711        self.underlying.save_dir()
712    }
713
714    fn composite(&self, s: i32) -> CompositeData<CC::Algebra> {
715        let d1 = self.underlying.differential(s);
716        let d0 = self.underlying.differential(s - 1);
717        vec![(1, d1, d0)]
718    }
719
720    fn compute_intermediate(&self, g: BidegreeGenerator) -> FpVector {
721        let p = self.prime();
722        let target = self.underlying.module(g.s() - 3);
723        let mut result = FpVector::new(p, target.dimension(g.t() - 1));
724        let d = self.underlying.differential(g.s());
725        self.homotopies[g.s() - 1].act(
726            result.as_slice_mut(),
727            1,
728            g.t(),
729            d.output(g.t(), g.idx()).as_slice(),
730            false,
731        );
732        result
733    }
734}
735
736impl<CC: FreeChainComplex> SecondaryResolution<CC>
737where
738    CC::Algebra: PairAlgebra,
739{
740    pub fn new(cc: Arc<CC>) -> Self {
741        if let Some(p) = cc.save_dir().write() {
742            for subdir in SaveKind::secondary_data() {
743                subdir.create_dir(p).unwrap();
744            }
745        }
746
747        Self {
748            underlying: cc,
749            homotopies: OnceBiVec::new(2),
750            intermediates: DashMap::new(),
751        }
752    }
753
754    pub fn homotopy(&self, s: i32) -> &SecondaryHomotopy<CC::Algebra> {
755        &self.homotopies[s]
756    }
757
758    pub fn e3_page(&self) -> sseq::Sseq<2, sseq::Adams> {
759        let p = self.prime();
760
761        let mut sseq = self.underlying.to_sseq();
762
763        let mut source_vec = FpVector::new(p, 0);
764        let mut target_vec = FpVector::new(p, 0);
765
766        for b in self.underlying.iter_stem() {
767            if b.t() > 0
768                && self
769                    .underlying
770                    .has_computed_bidegree(b + Bidegree::n_s(-1, 2))
771            {
772                let m = self.homotopy(b.s() + 2).homotopies.hom_k(b.t());
773                if m.is_empty() || m[0].is_empty() {
774                    continue;
775                }
776
777                source_vec.set_scratch_vector_size(m.len());
778                target_vec.set_scratch_vector_size(m[0].len());
779
780                for (i, row) in m.into_iter().enumerate() {
781                    source_vec.set_to_zero();
782                    source_vec.set_entry(i, 1);
783                    target_vec.copy_from_slice(&row);
784
785                    let source = BidegreeElement::new(b, source_vec);
786                    sseq.add_differential(2, &source, target_vec.as_slice());
787
788                    source_vec = source.into_vec();
789                }
790            }
791        }
792
793        for b in self.underlying.iter_stem() {
794            if sseq.invalid(b) {
795                sseq.update_degree(b);
796            }
797        }
798        sseq
799    }
800}
801
802// Rustdoc ICE's when trying to document this struct. See
803// https://github.com/rust-lang/rust/issues/91380
804#[doc(hidden)]
805pub struct SecondaryResolutionHomomorphism<
806    CC1: FreeChainComplex,
807    CC2: FreeChainComplex<Algebra = CC1::Algebra>,
808> where
809    CC1::Algebra: PairAlgebra,
810{
811    source: Arc<SecondaryResolution<CC1>>,
812    target: Arc<SecondaryResolution<CC2>>,
813    underlying: Arc<ResolutionHomomorphism<CC1, CC2>>,
814    /// input s -> homotopy
815    homotopies: OnceBiVec<SecondaryHomotopy<CC1::Algebra>>,
816    intermediates: DashMap<BidegreeGenerator, FpVector>,
817}
818
819impl<CC1: FreeChainComplex, CC2: FreeChainComplex<Algebra = CC1::Algebra>> SecondaryLift
820    for SecondaryResolutionHomomorphism<CC1, CC2>
821where
822    CC1::Algebra: PairAlgebra,
823{
824    type Algebra = CC1::Algebra;
825    type Source = CC1;
826    type Target = CC2;
827    type Underlying = ResolutionHomomorphism<CC1, CC2>;
828
829    fn underlying(&self) -> Arc<Self::Underlying> {
830        Arc::clone(&self.underlying)
831    }
832
833    fn algebra(&self) -> Arc<Self::Algebra> {
834        self.source.algebra()
835    }
836
837    fn source(&self) -> Arc<Self::Source> {
838        Arc::clone(&self.source.underlying)
839    }
840
841    fn target(&self) -> Arc<Self::Target> {
842        Arc::clone(&self.target.underlying)
843    }
844
845    fn shift(&self) -> Bidegree {
846        self.underlying.shift + Bidegree::s_t(1, 0)
847    }
848
849    fn max(&self) -> BidegreeRange<'_, Self> {
850        BidegreeRange::new(
851            self,
852            self.underlying.next_homological_degree(),
853            &|selff, s| {
854                std::cmp::min(
855                    selff.underlying.get_map(s).next_degree(),
856                    std::cmp::min(
857                        selff.source.homotopies[s].homotopies.next_degree(),
858                        if s == selff.shift().s() {
859                            i32::MAX
860                        } else {
861                            selff.target.homotopies[s + 1 - selff.shift().s()]
862                                .composites
863                                .max_degree()
864                                + selff.shift().t()
865                                + 1
866                        },
867                    ),
868                )
869            },
870        )
871    }
872
873    fn homotopies(&self) -> &OnceBiVec<SecondaryHomotopy<Self::Algebra>> {
874        &self.homotopies
875    }
876
877    fn intermediates(&self) -> &DashMap<BidegreeGenerator, FpVector> {
878        &self.intermediates
879    }
880
881    fn save_dir(&self) -> &SaveDirectory {
882        self.underlying.save_dir()
883    }
884
885    fn composite(&self, s: i32) -> CompositeData<Self::Algebra> {
886        let p = self.prime();
887        // This is -1 mod p^2
888        let neg_1 = p * p - 1;
889
890        let d_source = self.source.underlying.differential(s);
891        let d_target = self
892            .target
893            .underlying
894            .differential(s + 1 - self.shift().s());
895
896        let c1 = self.underlying.get_map(s);
897        let c0 = self.underlying.get_map(s - 1);
898
899        vec![(neg_1, d_source, c0), (1, c1, d_target)]
900    }
901
902    fn compute_intermediate(&self, g: BidegreeGenerator) -> FpVector {
903        let p = self.prime();
904        let neg_1 = p - 1;
905        let shifted_b = g.degree() - self.shift();
906        let target = self.target().module(shifted_b.s() - 1);
907
908        let mut result = FpVector::new(p, target.dimension(shifted_b.t() - 1));
909        let d = self.source().differential(g.s());
910
911        self.homotopies[g.s() - 1].act(
912            result.as_slice_mut(),
913            neg_1,
914            g.t(),
915            d.output(g.t(), g.idx()).as_slice(),
916            false,
917        );
918        self.target.homotopy(shifted_b.s() + 1).act(
919            result.as_slice_mut(),
920            neg_1,
921            shifted_b.t(),
922            self.underlying
923                .get_map(g.s())
924                .output(g.t(), g.idx())
925                .as_slice(),
926            true,
927        );
928        self.underlying.get_map(g.s() - 2).apply(
929            result.as_slice_mut(),
930            1,
931            g.t() - 1,
932            self.source
933                .homotopy(g.s())
934                .homotopies
935                .output(g.t(), g.idx())
936                .as_slice(),
937        );
938
939        result
940    }
941}
942
943impl<CC1: FreeChainComplex, CC2: FreeChainComplex<Algebra = CC1::Algebra>>
944    SecondaryResolutionHomomorphism<CC1, CC2>
945where
946    CC1::Algebra: PairAlgebra,
947{
948    pub fn new(
949        source: Arc<SecondaryResolution<CC1>>,
950        target: Arc<SecondaryResolution<CC2>>,
951        underlying: Arc<ResolutionHomomorphism<CC1, CC2>>,
952    ) -> Self {
953        assert!(Arc::ptr_eq(&underlying.source, &source.underlying));
954        assert!(Arc::ptr_eq(&underlying.target, &target.underlying));
955
956        if let Some(p) = underlying.save_dir().write() {
957            for subdir in SaveKind::secondary_data() {
958                subdir.create_dir(p).unwrap();
959            }
960        }
961
962        Self {
963            source,
964            target,
965            homotopies: OnceBiVec::new(underlying.shift.s() + 1),
966            underlying,
967            intermediates: DashMap::new(),
968        }
969    }
970
971    pub fn name(&self) -> String {
972        let name = self.underlying.name();
973        if name.starts_with('[') || name.starts_with('λ') {
974            name.to_owned()
975        } else {
976            format!("[{name}]")
977        }
978    }
979
980    pub fn homotopy(&self, s: i32) -> &SecondaryHomotopy<CC1::Algebra> {
981        &self.homotopies[s]
982    }
983
984    /// A version of [`hom_k`] but with a non-trivial λ part.
985    pub fn hom_k_with<'a>(
986        &self,
987        lambda_part: Option<&ResolutionHomomorphism<CC1, CC2>>,
988        sseq: Option<&sseq::Sseq<2, sseq::Adams>>,
989        b: Bidegree,
990        inputs: impl Iterator<Item = FpSlice<'a>>,
991        outputs: impl Iterator<Item = FpSliceMut<'a>>,
992    ) {
993        let source = b + self.shift() - Bidegree::s_t(1, 0);
994        let lambda_source = source + LAMBDA_BIDEGREE;
995
996        let p = self.prime();
997        let h_0 = self.algebra().p_tilde();
998
999        let source_num_gens = self.source().number_of_gens_in_bidegree(source);
1000        let lambda_num_gens = self.source().number_of_gens_in_bidegree(lambda_source);
1001
1002        let m0 = self.underlying.get_map(source.s()).hom_k(b.t());
1003        let mut m1 = Matrix::from_vec(p, &self.homotopy(lambda_source.s()).homotopies.hom_k(b.t()));
1004        if let Some(lambda_part) = lambda_part {
1005            m1 += &Matrix::from_vec(p, &lambda_part.get_map(lambda_source.s()).hom_k(b.t()));
1006        }
1007
1008        // The multiplication by p map
1009        let mp = Matrix::from_vec(
1010            p,
1011            &self
1012                .source()
1013                .filtration_one_product(1, h_0, source)
1014                .unwrap(),
1015        );
1016
1017        let sign = if (self.underlying.shift.s() * b.t()) % 2 == 1 {
1018            p * p - 1
1019        } else {
1020            1
1021        };
1022        let filtration_one_sign = if (b.t() % 2) == 1 { p - 1 } else { 1 };
1023
1024        let page_data = sseq.map(|sseq| {
1025            let d = sseq.page_data(lambda_source);
1026            &d[std::cmp::min(3, d.len() - 1)]
1027        });
1028
1029        let mut scratch0: Vec<u32> = Vec::new();
1030        for (input, mut out) in inputs.zip_eq(outputs) {
1031            scratch0.clear();
1032            scratch0.resize(source_num_gens, 0);
1033            for (i, v) in input.iter_nonzero() {
1034                scratch0
1035                    .iter_mut()
1036                    .zip_eq(&m0[i])
1037                    .for_each(|(a, b)| *a += v * b * sign);
1038                out.slice_mut(source_num_gens, source_num_gens + lambda_num_gens)
1039                    .add(m1.row(i), (v * sign) % p);
1040            }
1041            for (i, v) in scratch0.iter().enumerate() {
1042                out.add_basis_element(i, *v % p);
1043
1044                let extra = *v / p;
1045                out.slice_mut(source_num_gens, source_num_gens + lambda_num_gens)
1046                    .add(mp.row(i), (extra * filtration_one_sign) % p);
1047            }
1048            if let Some(page_data) = page_data {
1049                page_data.reduce_by_quotient(
1050                    out.slice_mut(source_num_gens, source_num_gens + lambda_num_gens),
1051                );
1052            }
1053        }
1054    }
1055
1056    /// Compute the induced map on Mod_{C\lambda^2} homotopy groups. This only computes it on
1057    /// standard lifts on elements in Ext. `outputs` is an iterator of `FpSliceMut`s whose lengths
1058    /// are equal to the total dimension of `(s + shift_s, t + shift_t)` and `(s + shift_s + 1, t +
1059    /// shift_t + 1)`. The first chunk records the Ext part of the result, and the second chunk
1060    /// records the λ part of the result.
1061    ///
1062    /// This reduces the λ part of the result by the image of d₂.
1063    ///
1064    /// # Arguments
1065    /// - `sseq`: A sseq object that records the $d_2$ differentials. If present, reduce the value
1066    ///   of the map by the image of $d_2$.
1067    pub fn hom_k<'a>(
1068        &self,
1069        sseq: Option<&sseq::Sseq<2, sseq::Adams>>,
1070        b: Bidegree,
1071        inputs: impl Iterator<Item = FpSlice<'a>>,
1072        outputs: impl Iterator<Item = FpSliceMut<'a>>,
1073    ) {
1074        self.hom_k_with(None, sseq, b, inputs, outputs);
1075    }
1076
1077    /// Given an element b whose product with this is null, find the element whose $d_2$ hits the
1078    /// λ part of the composition.
1079    ///
1080    /// # Arguments:
1081    /// - `sseq`: spectral sequence object of the source
1082    pub fn product_nullhomotopy(
1083        &self,
1084        lambda_part: Option<&ResolutionHomomorphism<CC1, CC2>>,
1085        sseq: &sseq::Sseq<2, sseq::Adams>,
1086        b: Bidegree,
1087        class: FpSlice,
1088    ) -> FpVector {
1089        let p = self.prime();
1090        let shift = self.underlying.shift;
1091
1092        let result_num_gens = self
1093            .source()
1094            .number_of_gens_in_bidegree(shift + b - Bidegree::s_t(1, 0));
1095
1096        let lambda_num_gens = self
1097            .source()
1098            .number_of_gens_in_bidegree(b + shift + LAMBDA_BIDEGREE);
1099
1100        let lower_num_gens = self.source().number_of_gens_in_bidegree(b + shift);
1101
1102        let target_num_gens = self.target().number_of_gens_in_bidegree(b);
1103        let target_lambda_num_gens = self
1104            .target()
1105            .number_of_gens_in_bidegree(b + LAMBDA_BIDEGREE);
1106
1107        let mut output_class = FpVector::new(p, result_num_gens);
1108        if result_num_gens == 0 || lambda_num_gens == 0 {
1109            return output_class;
1110        }
1111
1112        let mut prod_value = FpVector::new(p, lower_num_gens + lambda_num_gens);
1113        self.hom_k_with(
1114            lambda_part,
1115            None,
1116            b,
1117            [class.restrict(0, target_num_gens)].into_iter(),
1118            [prod_value.as_slice_mut()].into_iter(),
1119        );
1120        assert!(prod_value.slice(0, lower_num_gens).is_zero());
1121
1122        let matrix = Matrix::from_vec(
1123            p,
1124            &self
1125                .underlying
1126                .get_map((b + shift + LAMBDA_BIDEGREE).s())
1127                .hom_k((b + LAMBDA_BIDEGREE).t()),
1128        );
1129        matrix.apply(
1130            prod_value.slice_mut(lower_num_gens, lower_num_gens + lambda_num_gens),
1131            1,
1132            class.restrict(target_num_gens, target_num_gens + target_lambda_num_gens),
1133        );
1134
1135        let diff_source = b + shift - Bidegree::n_s(-1, 1);
1136        sseq.differentials(diff_source)[2].quasi_inverse(
1137            output_class.as_slice_mut(),
1138            prod_value.slice(lower_num_gens, lower_num_gens + lambda_num_gens),
1139        );
1140
1141        output_class
1142    }
1143}
1144
1145#[doc(hidden)]
1146pub struct SecondaryChainHomotopy<
1147    S: FreeChainComplex,
1148    T: FreeChainComplex<Algebra = S::Algebra> + Sync,
1149    U: FreeChainComplex<Algebra = S::Algebra> + Sync,
1150> where
1151    S::Algebra: PairAlgebra,
1152{
1153    underlying: Arc<ChainHomotopy<S, T, U>>,
1154    left: Arc<SecondaryResolutionHomomorphism<S, T>>,
1155    right: Arc<SecondaryResolutionHomomorphism<T, U>>,
1156    left_lambda: Option<Arc<ResolutionHomomorphism<S, T>>>,
1157    right_lambda: Option<Arc<ResolutionHomomorphism<T, U>>>,
1158    homotopies: OnceBiVec<SecondaryHomotopy<S::Algebra>>,
1159    intermediates: DashMap<BidegreeGenerator, FpVector>,
1160}
1161
1162impl<
1163    S: FreeChainComplex,
1164    T: FreeChainComplex<Algebra = S::Algebra> + Sync,
1165    U: FreeChainComplex<Algebra = S::Algebra> + Sync,
1166> SecondaryLift for SecondaryChainHomotopy<S, T, U>
1167where
1168    S::Algebra: PairAlgebra,
1169{
1170    type Algebra = S::Algebra;
1171    type Source = S;
1172    type Target = U;
1173    type Underlying = ChainHomotopy<S, T, U>;
1174
1175    const HIT_GENERATOR: bool = true;
1176
1177    fn underlying(&self) -> Arc<Self::Underlying> {
1178        Arc::clone(&self.underlying)
1179    }
1180
1181    fn algebra(&self) -> Arc<Self::Algebra> {
1182        self.left.algebra()
1183    }
1184
1185    fn source(&self) -> Arc<Self::Source> {
1186        self.left.source()
1187    }
1188
1189    fn target(&self) -> Arc<Self::Target> {
1190        self.right.target()
1191    }
1192
1193    fn shift(&self) -> Bidegree {
1194        Bidegree::s_t(
1195            self.underlying.shift().s(),
1196            self.left.shift().t() + self.right.shift().t(),
1197        )
1198    }
1199
1200    fn max(&self) -> BidegreeRange<'_, Self> {
1201        BidegreeRange::new(
1202            self,
1203            std::cmp::min(
1204                self.right.target.max().s() + self.shift().s() - 1,
1205                self.left.source.max().s(),
1206            ),
1207            &|selff, s| {
1208                std::cmp::min(
1209                    selff.left.source.max().t(s),
1210                    if s == selff.shift().s() {
1211                        i32::MAX
1212                    } else {
1213                        selff.right.target.max().t(s - selff.shift().s() + 1) + selff.shift().t()
1214                    },
1215                )
1216            },
1217        )
1218    }
1219
1220    fn homotopies(&self) -> &OnceBiVec<SecondaryHomotopy<S::Algebra>> {
1221        &self.homotopies
1222    }
1223
1224    fn intermediates(&self) -> &DashMap<BidegreeGenerator, FpVector> {
1225        &self.intermediates
1226    }
1227
1228    fn save_dir(&self) -> &SaveDirectory {
1229        self.underlying.save_dir()
1230    }
1231
1232    fn compute_intermediate(&self, g: BidegreeGenerator) -> FpVector {
1233        let p = self.prime();
1234        let neg_1 = p - 1;
1235        let shifted_b = g.degree() - self.shift();
1236
1237        let target = self.target().module(shifted_b.s() - 1);
1238
1239        let mut result = FpVector::new(p, target.dimension(shifted_b.t() - 1));
1240
1241        self.homotopies[g.s() - 1].act(
1242            result.as_slice_mut(),
1243            1,
1244            g.t(),
1245            self.source()
1246                .differential(g.s())
1247                .output(g.t(), g.idx())
1248                .as_slice(),
1249            false,
1250        );
1251
1252        self.right.target.homotopies()[shifted_b.s() + 1].act(
1253            result.as_slice_mut(),
1254            1,
1255            shifted_b.t(),
1256            self.underlying
1257                .homotopy(g.s())
1258                .output(g.t(), g.idx())
1259                .as_slice(),
1260            true,
1261        );
1262
1263        self.underlying.homotopy(g.s() - 2).apply(
1264            result.as_slice_mut(),
1265            neg_1,
1266            g.t() - 1,
1267            self.left.source.homotopies()[g.s()]
1268                .homotopies
1269                .output(g.t(), g.idx())
1270                .as_slice(),
1271        );
1272
1273        let left_shifted_b = g.degree() - self.left.underlying.shift;
1274        self.right.homotopies()[left_shifted_b.s()].act(
1275            result.as_slice_mut(),
1276            neg_1,
1277            left_shifted_b.t(),
1278            self.left
1279                .underlying
1280                .get_map(g.s())
1281                .output(g.t(), g.idx())
1282                .as_slice(),
1283            true,
1284        );
1285
1286        // This is inefficient if both right_lambda and right are non-zero, but this is not needed atm
1287        // and the change would not be user-facing.
1288        if let Some(right_lambda) = &self.right_lambda {
1289            right_lambda.get_map(left_shifted_b.s()).apply(
1290                result.as_slice_mut(),
1291                neg_1,
1292                left_shifted_b.t(),
1293                self.left
1294                    .underlying
1295                    .get_map(g.s())
1296                    .output(g.t(), g.idx())
1297                    .as_slice(),
1298            );
1299        }
1300
1301        self.right.underlying.get_map(left_shifted_b.s() - 1).apply(
1302            result.as_slice_mut(),
1303            neg_1,
1304            left_shifted_b.t() - 1,
1305            self.left.homotopies()[g.s()]
1306                .homotopies
1307                .output(g.t(), g.idx())
1308                .as_slice(),
1309        );
1310
1311        if let Some(left_lambda) = &self.left_lambda {
1312            self.right.underlying.get_map(left_shifted_b.s() - 1).apply(
1313                result.as_slice_mut(),
1314                neg_1,
1315                left_shifted_b.t() - 1,
1316                left_lambda.get_map(g.s()).output(g.t(), g.idx()).as_slice(),
1317            );
1318        }
1319        result
1320    }
1321
1322    fn composite(&self, s: i32) -> CompositeData<S::Algebra> {
1323        let p = self.prime();
1324        // This is -1 mod p^2
1325        let neg_1 = p * p - 1;
1326
1327        vec![
1328            (
1329                neg_1,
1330                self.underlying.left().get_map(s),
1331                self.underlying
1332                    .right()
1333                    .get_map(s - self.left.underlying.shift.s()),
1334            ),
1335            (
1336                1,
1337                self.underlying.homotopy(s),
1338                self.target().differential(s - self.shift().s() + 1),
1339            ),
1340            (
1341                1,
1342                self.source().differential(s),
1343                self.underlying.homotopy(s - 1),
1344            ),
1345        ]
1346    }
1347}
1348
1349impl<
1350    S: FreeChainComplex,
1351    T: FreeChainComplex<Algebra = S::Algebra> + Sync,
1352    U: FreeChainComplex<Algebra = S::Algebra> + Sync,
1353> SecondaryChainHomotopy<S, T, U>
1354where
1355    S::Algebra: PairAlgebra,
1356{
1357    pub fn new(
1358        left: Arc<SecondaryResolutionHomomorphism<S, T>>,
1359        right: Arc<SecondaryResolutionHomomorphism<T, U>>,
1360        left_lambda: Option<Arc<ResolutionHomomorphism<S, T>>>,
1361        right_lambda: Option<Arc<ResolutionHomomorphism<T, U>>>,
1362        underlying: Arc<ChainHomotopy<S, T, U>>,
1363    ) -> Self {
1364        assert!(Arc::ptr_eq(&underlying.left(), &left.underlying));
1365        assert!(Arc::ptr_eq(&underlying.right(), &right.underlying));
1366
1367        if let Some(left_lambda) = &left_lambda {
1368            assert!(Arc::ptr_eq(&left_lambda.source, &underlying.left().source));
1369            assert!(Arc::ptr_eq(&left_lambda.target, &underlying.left().target));
1370
1371            assert_eq!(left_lambda.shift, underlying.left().shift + LAMBDA_BIDEGREE);
1372        }
1373
1374        if let Some(right_lambda) = &right_lambda {
1375            assert!(Arc::ptr_eq(
1376                &right_lambda.source,
1377                &underlying.right().source
1378            ));
1379            assert!(Arc::ptr_eq(
1380                &right_lambda.target,
1381                &underlying.right().target
1382            ));
1383
1384            assert_eq!(
1385                right_lambda.shift,
1386                underlying.right().shift + LAMBDA_BIDEGREE
1387            );
1388        }
1389
1390        if let Some(p) = underlying.save_dir().write() {
1391            for subdir in SaveKind::secondary_data() {
1392                subdir.create_dir(p).unwrap();
1393            }
1394        }
1395
1396        Self {
1397            left,
1398            right,
1399            left_lambda,
1400            right_lambda,
1401            homotopies: OnceBiVec::new(underlying.shift().s()),
1402            underlying,
1403            intermediates: DashMap::new(),
1404        }
1405    }
1406}
1407
1408#[cfg(test)]
1409mod tests {
1410    use serde_json::json;
1411
1412    use super::*;
1413    use crate::utils;
1414
1415    #[test]
1416    #[should_panic(
1417        expected = "secondary: Failed to lift at (14, 3). This likely indicates an invalid input."
1418    )]
1419    fn cofib_h4() {
1420        let module = json!({
1421            "type": "finite dimensional module",
1422            "p": 2,
1423            "gens": {
1424                "x0": 0,
1425                "x16": 16,
1426            },
1427            "actions": ["Sq16 x0 = x16"]
1428        });
1429        let resolution = utils::construct((module, algebra::AlgebraType::Milnor), None).unwrap();
1430        resolution.compute_through_stem(Bidegree::n_s(20, 5));
1431        let lift = SecondaryResolution::new(Arc::new(resolution));
1432        lift.extend_all();
1433    }
1434}