1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

use crate::markoff::*;
use crate::numbers::*;
use crate::streams::*;

use rayon::iter::plumbing::*;
use rayon::iter::*;

/// A stream which can be run either in sequence or parallel, yielding Markoff numbers modulo `P`.
#[derive(Clone)]
pub struct CoordStream<'a, S, const L_HYPER: usize, const L_ELLIP: usize, const P: u128>
where
    FpNum<P>: SylowDecomposable<S>,
    QuadNum<P>: SylowDecomposable<S>,
{
    hyper_stream: Option<SylowStream<S, L_HYPER, FpNum<P>, ()>>,
    ellip_stream: Option<SylowStream<S, L_ELLIP, QuadNum<P>, ()>>,
    hyper_decomp: &'a SylowDecomp<S, L_HYPER, FpNum<P>>,
    ellip_decomp: &'a SylowDecomp<S, L_ELLIP, QuadNum<P>>,
}

impl<'a, S, const L_HYPER: usize, const L_ELLIP: usize, const P: u128>
    CoordStream<'a, S, L_HYPER, L_ELLIP, P>
where
    FpNum<P>: SylowDecomposable<S>,
    QuadNum<P>: SylowDecomposable<S>,
{
    /// Creates a new `CoordStream` with orders up to `limit`.
    pub fn new(
        hyper_decomp: &'a SylowDecomp<S, L_HYPER, FpNum<P>>,
        ellip_decomp: &'a SylowDecomp<S, L_ELLIP, QuadNum<P>>,
        hyper_lim: u128,
        ellip_lim: u128,
    ) -> CoordStream<'a, S, L_HYPER, L_ELLIP, P> {
        let hyper_stream = DivisorStream::new(FpNum::FACTORS.factors(), hyper_lim, true)
            .map(|v| v.try_into().unwrap())
            .fold(
                SylowStreamBuilder::<S, L_HYPER, FpNum<P>, ()>::new()
                    .add_flag(flags::NO_PARABOLIC)
                    .add_flag(flags::NO_UPPER_HALF)
                    .add_flag(flags::LEQ),
                |b, x| b.add_target(&x),
            )
            .into_iter();
        let ellip_stream = DivisorStream::new(QuadNum::FACTORS.factors(), ellip_lim, true)
            .map(|v| v.try_into().unwrap())
            .fold(
                SylowStreamBuilder::<S, L_ELLIP, QuadNum<P>, ()>::new()
                    .add_flag(flags::NO_PARABOLIC)
                    .add_flag(flags::NO_UPPER_HALF)
                    .add_flag(flags::LEQ),
                |b, x| b.add_target(&x),
            )
            .into_iter();
        CoordStream {
            hyper_stream: Some(hyper_stream),
            ellip_stream: Some(ellip_stream),
            hyper_decomp,
            ellip_decomp,
        }
    }

    /// Returns an iterator yielding pairs of coordinates without repeats up to permutation.
    pub fn upper_triangle(self) -> impl ParallelIterator<Item = (Coord<P>, Coord<P>)> + 'a
    where
        S: Clone + Send + Sync,
    {
        #[derive(Clone)]
        struct Helper<I>(I);
        impl<I: Iterator + Clone> Iterator for Helper<I> {
            type Item = (I, I::Item);
            fn next(&mut self) -> Option<Self::Item> {
                let iter = self.0.clone();
                let Some(a) = self.0.next() else {
                    return None;
                };
                Some((iter, a))
            }
        }
        Helper(self)
            .par_bridge()
            .flat_map(|(iter, a)| ParallelIterator::map(iter, move |b| (a, b)))
    }
}

impl<'a, S, const L_HYPER: usize, const L_ELLIP: usize, const P: u128> Iterator
    for CoordStream<'a, S, L_HYPER, L_ELLIP, P>
where
    FpNum<P>: SylowDecomposable<S>,
    QuadNum<P>: SylowDecomposable<S>,
{
    type Item = Coord<P>;

    fn next(&mut self) -> Option<Self::Item> {
        if let Some(stream) = self.hyper_stream.as_mut() {
            if let Some((a, _)) = stream.next() {
                return Some(Coord(FpNum::from_chi(&a, self.hyper_decomp)));
            }
            self.hyper_stream = None;
        }
        if let Some(stream) = self.ellip_stream.as_mut() {
            if let Some((a, _)) = stream.next() {
                return Some(Coord(QuadNum::from_chi(&a, self.ellip_decomp)));
            }
            self.ellip_stream = None;
        }
        None
    }
}

impl<'a, S, const L_HYPER: usize, const L_ELLIP: usize, const P: u128> ParallelIterator
    for CoordStream<'a, S, L_HYPER, L_ELLIP, P>
where
    S: Send + Sync,
    FpNum<P>: SylowDecomposable<S>,
    QuadNum<P>: SylowDecomposable<S>,
{
    type Item = Coord<P>;

    fn drive_unindexed<C>(self, consumer: C) -> C::Result
    where
        C: UnindexedConsumer<Self::Item>,
    {
        let left = self.hyper_stream.map(|stream| {
            stream
                .parallelize()
                .map(|(x, _)| Coord(FpNum::from_chi(&x, self.hyper_decomp)))
                .drive_unindexed(consumer.split_off_left())
        });
        let right = self.ellip_stream.map(|stream| {
            stream
                .parallelize()
                .map(|(x, _)| Coord(QuadNum::from_chi(&x, self.ellip_decomp)))
                .drive_unindexed(consumer.split_off_left())
        });
        match (left, right) {
            (Some(l), Some(r)) => consumer.to_reducer().reduce(l, r),
            (None, Some(r)) => r,
            (None, None) => consumer.into_folder().complete(),
            (Some(l), None) => l,
        }
    }
}