//! Pollard rho: Teske r-adding / mixed walks, Floyd / Brent single-walk cycle detection,
//! van Oorschot–Wiener distinguished-point parallel collision search, negation map and
//! order-3 automorphism classes with fruitless-cycle handling.
//!
//! Iteration function (see README.md for the full specification and references):
//!   j = h(x(P)) ∈ [0, r),  f(P) = ‖P + M_j‖  (adding partition)  or  ‖2P‖ (doubling partition, --mixed),
//!   where M_j = u_j G + v_j Q and ‖·‖ is the canonical representative of the equivalence class
//!   (identity, {P, −P} with --neg, {P, ζP, ζ²P} with --aut3, the order-6 orbit with both).
//! With --neg the look-ahead rule of Bos–Kleinjung–Lenstra avoids fruitless 2-cycles, and every
//! walk remembers one point every `cycle_check` steps to detect longer fruitless cycles, escaped by
//! doubling the minimum point of the cycle (Bernstein–Lange–Schwabe 2011).

use crate::curve::{CombTable, Counters, Pt};
use crate::instance::{is_distinguished, partition, Instance};
use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha20Rng;
use serde::{Deserialize, Serialize};
use std::sync::atomic::{AtomicBool, Ordering};
use std::sync::Mutex;
use std::time::Instant;

#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]
#[serde(rename_all = "lowercase")]
pub enum RhoMode {
    Floyd,
    Brent,
    Dp,
}

#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct RhoConfig {
    pub r: usize,
    /// 0 = pure r-adding walk; k > 0 = partitions j with j % k == k−1 are doubling partitions
    pub mixed: usize,
    pub neg: bool,
    pub aut3: bool,
    /// remember one point every `cycle_check` steps (detects fruitless cycles of length ≤ cycle_check)
    pub cycle_check: u64,
    pub dp_bits: u32,
    /// walks per thread (batch size of the simultaneous inversion)
    pub walks: usize,
    pub threads: usize,
    pub mode: RhoMode,
    pub seed: u64,
    /// walks longer than this (steps since last distinguished point) are abandoned and restarted
    pub max_walk_len: u64,
}

impl RhoConfig {
    pub fn default_dp_bits(bits: u32) -> u32 {
        (bits / 4).saturating_sub(2).max(1)
    }
}

#[derive(Clone, Copy, Debug)]
pub struct Mult {
    pub pt: Pt,
    pub u: u64,
    pub v: u64,
}

#[derive(Clone, Debug, Default, Serialize, Deserialize)]
pub struct WalkStats {
    pub ops: Counters,
    /// steps of the iteration function actually taken (advances), all walks
    pub iterations: u64,
    pub fruitless_cycles: u64,
    /// fruitless cycles by length (index = length, up to 64)
    pub fruitless_by_length: Vec<u64>,
    /// look-ahead rejections (2-cycle avoidance), each costs one wasted addition
    pub lookahead_rejections: u64,
    pub doubling_fallbacks: u64,
    pub abandoned_walks: u64,
    pub restarts: u64,
    pub dp_found: u64,
    pub useless_collisions: u64,
    pub cpu_s: f64,
}

impl WalkStats {
    pub fn merge(&mut self, o: &WalkStats) {
        self.ops.merge(&o.ops);
        self.iterations += o.iterations;
        self.fruitless_cycles += o.fruitless_cycles;
        if self.fruitless_by_length.len() < o.fruitless_by_length.len() {
            self.fruitless_by_length.resize(o.fruitless_by_length.len(), 0);
        }
        for (i, v) in o.fruitless_by_length.iter().enumerate() {
            self.fruitless_by_length[i] += v;
        }
        self.lookahead_rejections += o.lookahead_rejections;
        self.doubling_fallbacks += o.doubling_fallbacks;
        self.abandoned_walks += o.abandoned_walks;
        self.restarts += o.restarts;
        self.dp_found += o.dp_found;
        self.useless_collisions += o.useless_collisions;
        self.cpu_s += o.cpu_s;
    }
}

#[derive(Clone, Copy, Debug)]
pub struct Walk {
    pub p: Pt,
    pub u: u64,
    pub v: u64,
    /// partition override for the look-ahead retry (u32::MAX = none, FORCE_DBL = doubling fallback)
    ovr: u32,
    tries: u32,
    /// steps since the last (re)start
    pub steps: u64,
    saved: Pt,
    saved_step: u64,
    last_escape_x: u64,
}

const NO_OVR: u32 = u32::MAX;
const FORCE_DBL: u32 = u32::MAX - 1;

/// Data shared by all threads of one rho run.
pub struct Shared<'a> {
    pub inst: &'a Instance,
    pub cfg: &'a RhoConfig,
    pub mults: Vec<Mult>,
    pub comb_g: CombTable,
    pub comb_q: CombTable,
    pub beta2: u64,
    pub lam2: u64,
    pub dp_mask: u64,
    pub setup_ops: Counters,
}

impl<'a> Shared<'a> {
    pub fn new(inst: &'a Instance, cfg: &'a RhoConfig) -> Shared<'a> {
        let mut c = Counters::default();
        let mut rng = ChaCha20Rng::seed_from_u64(cfg.seed ^ 0x5EED_0000_0000_0001);
        let comb_g = CombTable::new(&inst.curve, &inst.g, &mut c);
        let comb_q = CombTable::new(&inst.curve, &inst.q, &mut c);
        let n = inst.n.n;
        let mut mults = Vec::with_capacity(cfg.r);
        while mults.len() < cfg.r {
            let u = rng.random_range(1..n);
            let v = rng.random_range(1..n);
            let pt = inst.curve.add(&comb_g.mul(&inst.curve, u, &mut c), &comb_q.mul(&inst.curve, v, &mut c), &mut c);
            if pt.inf {
                continue; // u + v·d ≡ 0 — astronomically unlikely; resample
            }
            mults.push(Mult { pt, u, v });
        }
        let (beta2, lam2) = match (inst.beta, inst.lambda) {
            (Some(b), Some(l)) => (inst.curve.f.sqr(b), inst.n.mul(l, l)),
            _ => (0, 0),
        };
        let dp_mask = if cfg.dp_bits >= 64 { u64::MAX } else { (1u64 << cfg.dp_bits) - 1 };
        Shared { inst, cfg, mults, comb_g, comb_q, beta2, lam2, dp_mask, setup_ops: c }
    }
}

pub enum Outcome {
    Advanced,
    Rejected,
    Solved(u64),
}

/// Per-thread walk context.
pub struct Ctx<'a> {
    pub sh: &'a Shared<'a>,
    pub rng: ChaCha20Rng,
    pub st: WalkStats,
    dens: Vec<u64>,
    scratch: Vec<u64>,
    js: Vec<u32>,
    dbls: Vec<bool>,
    special: Vec<bool>,
}

#[derive(Clone, Copy, Debug)]
pub struct DpRecord {
    pub x: u64,
    pub y: u64,
    pub u: u64,
    pub v: u64,
}

impl<'a> Ctx<'a> {
    pub fn new(sh: &'a Shared<'a>, thread_id: u64) -> Ctx<'a> {
        let mut seed = [0u8; 32];
        seed[..8].copy_from_slice(&sh.cfg.seed.to_le_bytes());
        seed[8..16].copy_from_slice(&thread_id.to_le_bytes());
        seed[16..24].copy_from_slice(b"rho_walk");
        let w = sh.cfg.walks.max(1);
        Ctx {
            sh,
            rng: ChaCha20Rng::from_seed(seed),
            st: WalkStats { fruitless_by_length: vec![0; 65], ..Default::default() },
            dens: vec![0; w],
            scratch: vec![0; w],
            js: vec![0; w],
            dbls: vec![false; w],
            special: vec![false; w],
        }
    }

    #[inline(always)]
    fn partition_of(&self, x: u64) -> u32 {
        partition(x, self.sh.cfg.r as u64)
    }

    #[inline(always)]
    fn is_dbl_partition(&self, j: u32) -> bool {
        let k = self.sh.cfg.mixed;
        k > 0 && (j as usize + 1) % k == 0
    }

    /// Canonical representative of the class of (P; u, v) and the adjusted scalars.
    #[inline(always)]
    pub fn canon(&mut self, mut p: Pt, mut u: u64, mut v: u64) -> (Pt, u64, u64) {
        let inst = self.sh.inst;
        let f = &inst.curve.f;
        let cfg = self.sh.cfg;
        if cfg.aut3 {
            // candidates x, βx, β²x  ↔  P, λP, λ²P ; pick the minimum standard-form x
            let beta = inst.beta.unwrap();
            let lam = inst.lambda.unwrap();
            let x1 = f.mul(beta, p.x);
            let x2 = f.mul(self.sh.beta2, p.x);
            let s0 = f.from_mont(p.x);
            let s1 = f.from_mont(x1);
            let s2 = f.from_mont(x2);
            self.st.ops.mul += 5;
            if s1 < s0 && s1 <= s2 {
                p.x = x1;
                u = inst.n.mul(u, lam);
                v = inst.n.mul(v, lam);
            } else if s2 < s0 && s2 < s1 {
                p.x = x2;
                u = inst.n.mul(u, self.sh.lam2);
                v = inst.n.mul(v, self.sh.lam2);
            }
        }
        if cfg.neg {
            // keep the representative with y_std < p − y_std, i.e. 2·y_std < p
            let ys = f.from_mont(p.y);
            self.st.ops.mul += 1;
            if 2 * ys > f.p {
                p.y = f.neg(p.y);
                u = inst.n.neg(u);
                v = inst.n.neg(v);
            }
        }
        (p, u, v)
    }

    /// d = num / den mod n, verified against Q. None if den ≡ 0 or verification fails.
    fn try_solve(&mut self, num: u64, den: u64) -> Option<u64> {
        let inst = self.sh.inst;
        let d = inst.n.solve(den, num)?;
        if inst.verify(d, &mut self.st.ops) {
            Some(d)
        } else {
            None
        }
    }

    /// Fresh random walk start ‖u G + v Q‖. Err(d) if the start happens to be O (then d = −u/v).
    pub fn random_walk(&mut self) -> Result<Walk, u64> {
        let inst = self.sh.inst;
        let n = inst.n.n;
        loop {
            let u = self.rng.random_range(0..n);
            let v = self.rng.random_range(1..n);
            let p = inst.curve.add(
                &self.sh.comb_g.mul(&inst.curve, u, &mut self.st.ops),
                &self.sh.comb_q.mul(&inst.curve, v, &mut self.st.ops),
                &mut self.st.ops,
            );
            self.st.restarts += 1;
            if p.inf {
                if let Some(d) = self.try_solve(inst.n.neg(u), v) {
                    return Err(d);
                }
                continue;
            }
            let (p, u, v) = self.canon(p, u, v);
            return Ok(Walk {
                p,
                u,
                v,
                ovr: NO_OVR,
                tries: 0,
                steps: 0,
                saved: Pt::INF,
                saved_step: 0,
                last_escape_x: u64::MAX,
            });
        }
    }

    #[inline(always)]
    fn choose(&self, w: &Walk) -> (u32, bool) {
        if w.ovr == NO_OVR {
            let j = self.partition_of(w.p.x);
            (j, self.is_dbl_partition(j))
        } else if w.ovr == FORCE_DBL {
            (0, true)
        } else {
            (w.ovr, self.is_dbl_partition(w.ovr))
        }
    }

    /// Denominator of the slope for the chosen step; 0 marks a special case.
    #[inline(always)]
    fn denominator(&self, w: &Walk, j: u32, is_dbl: bool) -> u64 {
        let f = &self.sh.inst.curve.f;
        if is_dbl {
            f.dbl(w.p.y)
        } else {
            f.sub(self.sh.mults[j as usize].pt.x, w.p.x)
        }
    }

    /// Apply one step given the inverse of the denominator. Does not touch cycle bookkeeping.
    #[inline(always)]
    fn apply(&mut self, w: &mut Walk, j: u32, is_dbl: bool, inv: u64) -> Outcome {
        let inst = self.sh.inst;
        let f = &inst.curve.f;
        let (r, u, v) = if is_dbl {
            let x2 = f.sqr(w.p.x);
            let num = f.add(f.add(f.add(x2, x2), x2), inst.curve.a);
            let lam = f.mul(num, inv);
            let r = inst.curve.finish_dbl(&w.p, lam);
            self.st.ops.dbl += 1;
            self.st.ops.mul += 4;
            (r, inst.n.dbl(w.u), inst.n.dbl(w.v))
        } else {
            let m = &self.sh.mults[j as usize];
            let lam = f.mul(f.sub(m.pt.y, w.p.y), inv);
            let r = inst.curve.finish_add(&w.p, m.pt.x, lam);
            self.st.ops.add += 1;
            self.st.ops.mul += 3;
            (r, inst.n.add(w.u, m.u), inst.n.add(w.v, m.v))
        };
        let (r, u, v) = self.canon(r, u, v);
        if self.sh.cfg.neg && !is_dbl && self.partition_of(r.x) == j {
            // look-ahead: ‖R + M_j‖ could be ‖P‖ (fruitless 2-cycle) → use the next partition
            self.st.lookahead_rejections += 1;
            w.tries += 1;
            if w.tries as usize >= self.sh.cfg.r {
                w.ovr = FORCE_DBL;
                self.st.doubling_fallbacks += 1;
            } else {
                w.ovr = (j + 1) % self.sh.cfg.r as u32;
            }
            return Outcome::Rejected;
        }
        w.p = r;
        w.u = u;
        w.v = v;
        w.ovr = NO_OVR;
        w.tries = 0;
        w.steps += 1;
        self.st.iterations += 1;
        Outcome::Advanced
    }

    /// Slow path for a zero denominator (P = ±M_j, or y = 0 for a doubling).
    fn special(&mut self, w: &mut Walk, j: u32, is_dbl: bool) -> Option<u64> {
        let inst = self.sh.inst;
        if is_dbl {
            // y = 0 ⇒ 2P = O: impossible in a prime-order group of odd order; restart defensively.
            return self.restart(w);
        }
        let m = self.sh.mults[j as usize];
        if w.p.y == m.pt.y {
            // P = M_j : u + v d = u_j + v_j d
            if let Some(d) = self.try_solve(inst.n.sub(m.u, w.u), inst.n.sub(w.v, m.v)) {
                return Some(d);
            }
            // same scalars: genuine doubling
            let r = inst.curve.dbl(&w.p, &mut self.st.ops);
            if r.inf {
                return self.restart(w);
            }
            let (r, u, v) = self.canon(r, inst.n.dbl(w.u), inst.n.dbl(w.v));
            w.p = r;
            w.u = u;
            w.v = v;
            w.ovr = NO_OVR;
            w.tries = 0;
            w.steps += 1;
            self.st.iterations += 1;
            None
        } else {
            // P = −M_j : u + v d = −u_j − v_j d  ⇒ (v + v_j) d = −(u + u_j)
            if let Some(d) = self.try_solve(inst.n.neg(inst.n.add(w.u, m.u)), inst.n.add(w.v, m.v)) {
                return Some(d);
            }
            self.st.useless_collisions += 1;
            self.restart(w)
        }
    }

    fn restart(&mut self, w: &mut Walk) -> Option<u64> {
        match self.random_walk() {
            Ok(nw) => {
                *w = nw;
                None
            }
            Err(d) => Some(d),
        }
    }

    /// One full step of the iteration function with a single (non-batched) inversion.
    pub fn step_single(&mut self, w: &mut Walk) -> Option<u64> {
        loop {
            let (j, is_dbl) = self.choose(w);
            let den = self.denominator(w, j, is_dbl);
            if den == 0 {
                return self.special(w, j, is_dbl);
            }
            let inv = self.sh.inst.curve.f.inv(den);
            self.st.ops.inv += 1;
            self.st.ops.mul += 1;
            match self.apply(w, j, is_dbl, inv) {
                Outcome::Advanced => return None,
                Outcome::Rejected => continue,
                Outcome::Solved(d) => return Some(d),
            }
        }
    }

    /// Cycle bookkeeping after an advance: detect a fruitless cycle (point equal to the remembered
    /// point) and escape it; otherwise remember the current point every `cycle_check` steps.
    #[inline(always)]
    fn after_advance(&mut self, w: &mut Walk) -> Option<u64> {
        if !w.saved.inf && w.p.x == w.saved.x && w.p.y == w.saved.y {
            let len = w.steps - w.saved_step;
            return self.escape_cycle(w, len);
        }
        if w.steps % self.sh.cfg.cycle_check == 0 {
            w.saved = w.p;
            w.saved_step = w.steps;
        }
        None
    }

    /// Escape a fruitless cycle of length `len` containing the current point: walk the cycle once
    /// to find its minimum point (smallest Montgomery x), then continue from ‖2·P_min‖.
    pub fn escape_cycle(&mut self, w: &mut Walk, len: u64) -> Option<u64> {
        let inst = self.sh.inst;
        self.st.fruitless_cycles += 1;
        let idx = (len as usize).min(64);
        self.st.fruitless_by_length[idx] += 1;
        let start = w.p;
        let mut best = *w;
        let mut cur = *w;
        for _ in 0..len.max(1) {
            if let Some(d) = self.step_single(&mut cur) {
                return Some(d);
            }
            if cur.p.x < best.p.x || (cur.p.x == best.p.x && cur.p.y < best.p.y) {
                best = cur;
            }
            if cur.p == start {
                break;
            }
        }
        if best.p.x == w.last_escape_x {
            // the escape led back into the same cycle: give up on this trajectory
            return self.restart(w);
        }
        let r = inst.curve.dbl(&best.p, &mut self.st.ops);
        if r.inf {
            return self.restart(w);
        }
        let (r, u, v) = self.canon(r, inst.n.dbl(best.u), inst.n.dbl(best.v));
        w.p = r;
        w.u = u;
        w.v = v;
        w.ovr = NO_OVR;
        w.tries = 0;
        w.steps += 1;
        w.saved = Pt::INF;
        w.last_escape_x = best.p.x;
        self.st.iterations += 1;
        None
    }

    /// Full single step including cycle bookkeeping (used by Floyd / Brent).
    pub fn step(&mut self, w: &mut Walk) -> Option<u64> {
        if let Some(d) = self.step_single(w) {
            return Some(d);
        }
        self.after_advance(w)
    }

    /// Batched step of all walks (one simultaneous inversion). Distinguished points are appended to
    /// `dps` and the corresponding walks restarted. Returns Some(d) if the DLP was solved on the way.
    pub fn batch_step(&mut self, walks: &mut [Walk], dps: &mut Vec<DpRecord>) -> Option<u64> {
        let k = walks.len();
        for i in 0..k {
            let (j, is_dbl) = self.choose(&walks[i]);
            let den = self.denominator(&walks[i], j, is_dbl);
            self.js[i] = j;
            self.dbls[i] = is_dbl;
            self.special[i] = den == 0;
            self.dens[i] = if den == 0 { self.sh.inst.curve.f.one } else { den };
        }
        let (m, inv) = self.sh.inst.curve.f.batch_inv(&mut self.dens[..k], &mut self.scratch[..k]);
        self.st.ops.mul += m;
        self.st.ops.inv += inv;
        for i in 0..k {
            let w = &mut walks[i];
            if self.special[i] {
                let (j, is_dbl) = (self.js[i], self.dbls[i]);
                if let Some(d) = self.special(w, j, is_dbl) {
                    return Some(d);
                }
                continue;
            }
            let inv = self.dens[i];
            match self.apply(w, self.js[i], self.dbls[i], inv) {
                Outcome::Rejected => continue,
                Outcome::Solved(d) => return Some(d),
                Outcome::Advanced => {}
            }
            if let Some(d) = self.after_advance(w) {
                return Some(d);
            }
            if is_distinguished(w.p.x, self.sh.dp_mask) {
                dps.push(DpRecord { x: w.p.x, y: w.p.y, u: w.u, v: w.v });
                self.st.dp_found += 1;
                if let Some(d) = self.restart(w) {
                    return Some(d);
                }
            } else if w.steps > self.sh.cfg.max_walk_len {
                self.st.abandoned_walks += 1;
                if let Some(d) = self.restart(w) {
                    return Some(d);
                }
            }
        }
        None
    }

    /// Resolve a collision of two walks at the same point (or at P and −P).
    fn collide(&mut self, a: &Walk, b: &Walk) -> Option<u64> {
        let inst = self.sh.inst;
        if a.p.y == b.p.y {
            if a.u == b.u && a.v == b.v {
                return None;
            }
            self.try_solve(inst.n.sub(a.u, b.u), inst.n.sub(b.v, a.v))
        } else {
            self.try_solve(inst.n.neg(inst.n.add(a.u, b.u)), inst.n.add(a.v, b.v))
        }
    }

    /// Floyd's cycle finding on a single walk: x_i vs x_{2i}.
    pub fn floyd(&mut self) -> u64 {
        let mut t = match self.random_walk() {
            Ok(w) => w,
            Err(d) => return d,
        };
        let mut h = t;
        loop {
            if let Some(d) = self.step(&mut t) {
                return d;
            }
            if let Some(d) = self.step(&mut h) {
                return d;
            }
            if let Some(d) = self.step(&mut h) {
                return d;
            }
            if t.p.x == h.p.x {
                if let Some(d) = self.collide(&t, &h) {
                    return d;
                }
                // same scalars (fruitless cycle) or degenerate: escape and continue
                let len = if t.p == h.p && t.u == h.u && t.v == h.v { self.cycle_length(&h) } else { None };
                if let Some(len) = len {
                    if let Some(d) = self.escape_cycle(&mut h, len) {
                        return d;
                    }
                    t = h;
                } else {
                    self.st.useless_collisions += 1;
                    match self.random_walk() {
                        Ok(w) => {
                            t = w;
                            h = t;
                        }
                        Err(d) => return d,
                    }
                }
            }
        }
    }

    /// Brent's cycle finding on a single walk (power-of-two teleporting tortoise).
    pub fn brent(&mut self) -> u64 {
        let mut t = match self.random_walk() {
            Ok(w) => w,
            Err(d) => return d,
        };
        let mut h = t;
        let mut power: u64 = 1;
        let mut lam: u64 = 0;
        loop {
            if let Some(d) = self.step(&mut h) {
                return d;
            }
            lam += 1;
            if t.p.x == h.p.x {
                if let Some(d) = self.collide(&t, &h) {
                    return d;
                }
                let len = if t.p == h.p && t.u == h.u && t.v == h.v { self.cycle_length(&h) } else { None };
                if let Some(len) = len {
                    if let Some(d) = self.escape_cycle(&mut h, len) {
                        return d;
                    }
                    t = h;
                } else {
                    self.st.useless_collisions += 1;
                    match self.random_walk() {
                        Ok(w) => {
                            t = w;
                            h = t;
                        }
                        Err(d) => return d,
                    }
                }
                power = 1;
                lam = 0;
                continue;
            }
            if power == lam {
                t = h;
                power *= 2;
                lam = 0;
            }
        }
    }

    /// Length of the pure-f cycle through w.p, if w.p lies on one of length ≤ 4·cycle_check (+1024).
    /// Floyd/Brent can also meet an earlier point of the *augmented* trajectory (after an escape),
    /// which is not an f-cycle: then None is returned and the caller restarts the walk.
    fn cycle_length(&mut self, w: &Walk) -> Option<u64> {
        let bound = 4 * self.sh.cfg.cycle_check + 1024;
        let mut cur = *w;
        let mut len = 0u64;
        loop {
            self.step_single(&mut cur)?;
            len += 1;
            if cur.p == w.p {
                return Some(len);
            }
            if len > bound {
                return None;
            }
        }
    }
}

#[derive(Clone, Copy, Debug)]
pub struct DpEntry {
    pub y: u64,
    pub u: u64,
    pub v: u64,
}

#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct RhoResult {
    pub d: Option<u64>,
    pub stats: WalkStats,
    pub dp_count: u64,
    pub dp_bytes: u64,
    pub table_entries: u64,
    pub threads: usize,
    pub wall_s: f64,
}

/// Expected number of iterations of a random mapping on a set of size n: √(πn/2).
pub fn expected_iterations(n: u64, class_size: u64) -> f64 {
    (std::f64::consts::PI * (n as f64) / (2.0 * class_size as f64)).sqrt()
}

pub fn solve_rho(inst: &Instance, cfg: &RhoConfig) -> RhoResult {
    let t0 = Instant::now();
    let sh = Shared::new(inst, cfg);
    let mut total = WalkStats { fruitless_by_length: vec![0; 65], ..Default::default() };
    total.ops.merge(&sh.setup_ops);

    // trivial cases
    if inst.q.inf {
        return RhoResult { d: Some(0), stats: total, dp_count: 0, dp_bytes: 0, table_entries: 0, threads: 1, wall_s: t0.elapsed().as_secs_f64() };
    }

    match cfg.mode {
        RhoMode::Floyd | RhoMode::Brent => {
            let cpu0 = crate::stats::thread_cpu_s();
            let mut ctx = Ctx::new(&sh, 0);
            let d = if cfg.mode == RhoMode::Floyd { ctx.floyd() } else { ctx.brent() };
            ctx.st.cpu_s = crate::stats::thread_cpu_s() - cpu0;
            total.merge(&ctx.st);
            RhoResult { d: Some(d), stats: total, dp_count: 0, dp_bytes: 0, table_entries: 0, threads: 1, wall_s: t0.elapsed().as_secs_f64() }
        }
        RhoMode::Dp => {
            let threads = cfg.threads.max(1);
            let class = (if cfg.neg { 2 } else { 1 }) * (if cfg.aut3 { 3 } else { 1 });
            let expected_dps = expected_iterations(inst.n.n, class) / (1u64 << cfg.dp_bits.min(63)) as f64;
            let table: Mutex<hashbrown::HashMap<u64, DpEntry, crate::stats::FxBuild>> =
                Mutex::new(hashbrown::HashMap::with_capacity_and_hasher((expected_dps * 2.0) as usize + 1024, Default::default()));
            let found: Mutex<Option<u64>> = Mutex::new(None);
            let stop = AtomicBool::new(false);
            let stats: Mutex<Vec<WalkStats>> = Mutex::new(Vec::new());

            std::thread::scope(|s| {
                for t in 0..threads {
                    let sh = &sh;
                    let table = &table;
                    let found = &found;
                    let stop = &stop;
                    let stats = &stats;
                    s.spawn(move || {
                        let cpu0 = crate::stats::thread_cpu_s();
                        let mut ctx = Ctx::new(sh, t as u64);
                        let mut walks: Vec<Walk> = Vec::with_capacity(sh.cfg.walks);
                        let mut solved: Option<u64> = None;
                        for _ in 0..sh.cfg.walks.max(1) {
                            match ctx.random_walk() {
                                Ok(w) => walks.push(w),
                                Err(d) => {
                                    solved = Some(d);
                                    break;
                                }
                            }
                        }
                        let mut dps: Vec<DpRecord> = Vec::new();
                        while solved.is_none() && !stop.load(Ordering::Relaxed) {
                            if let Some(d) = ctx.batch_step(&mut walks, &mut dps) {
                                solved = Some(d);
                                break;
                            }
                            if !dps.is_empty() {
                                let mut tab = table.lock().unwrap();
                                for rec in dps.drain(..) {
                                    match tab.get(&rec.x) {
                                        None => {
                                            tab.insert(rec.x, DpEntry { y: rec.y, u: rec.u, v: rec.v });
                                        }
                                        Some(e) => {
                                            let e = *e;
                                            let a = Walk { p: Pt::new(rec.x, e.y), u: e.u, v: e.v, ovr: NO_OVR, tries: 0, steps: 0, saved: Pt::INF, saved_step: 0, last_escape_x: 0 };
                                            let b = Walk { p: Pt::new(rec.x, rec.y), u: rec.u, v: rec.v, ovr: NO_OVR, tries: 0, steps: 0, saved: Pt::INF, saved_step: 0, last_escape_x: 0 };
                                            if a.u == b.u && a.v == b.v && a.p.y == b.p.y {
                                                continue; // identical trajectory data, nothing to learn
                                            }
                                            if let Some(d) = ctx.collide(&a, &b) {
                                                solved = Some(d);
                                                break;
                                            }
                                            ctx.st.useless_collisions += 1;
                                        }
                                    }
                                }
                            }
                        }
                        if let Some(d) = solved {
                            let mut f = found.lock().unwrap();
                            if f.is_none() {
                                *f = Some(d);
                            }
                            stop.store(true, Ordering::Relaxed);
                        }
                        ctx.st.cpu_s = crate::stats::thread_cpu_s() - cpu0;
                        stats.lock().unwrap().push(ctx.st);
                    });
                }
            });
            for st in stats.lock().unwrap().iter() {
                total.merge(st);
            }
            let tab = table.lock().unwrap();
            let entries = tab.len() as u64;
            let bytes = (tab.capacity() as u64) * (std::mem::size_of::<(u64, DpEntry)>() as u64 + 1);
            let d = *found.lock().unwrap();
            RhoResult { d, stats: total, dp_count: entries, dp_bytes: bytes, table_entries: entries, threads, wall_s: t0.elapsed().as_secs_f64() }
        }
    }
}