1//! Pollard rho: Teske r-adding / mixed walks, Floyd / Brent single-walk cycle detection,2//! van Oorschot–Wiener distinguished-point parallel collision search, negation map and3//! order-3 automorphism classes with fruitless-cycle handling.4//!5//! Iteration function (see README.md for the full specification and references):6//! j = h(x(P)) ∈ [0, r), f(P) = ‖P + M_j‖ (adding partition) or ‖2P‖ (doubling partition, --mixed),7//! where M_j = u_j G + v_j Q and ‖·‖ is the canonical representative of the equivalence class8//! (identity, {P, −P} with --neg, {P, ζP, ζ²P} with --aut3, the order-6 orbit with both).9//! With --neg the look-ahead rule of Bos–Kleinjung–Lenstra avoids fruitless 2-cycles, and every10//! walk remembers one point every `cycle_check` steps to detect longer fruitless cycles, escaped by11//! doubling the minimum point of the cycle (Bernstein–Lange–Schwabe 2011).1213use crate::curve::{CombTable, Counters, Pt};14use crate::instance::{is_distinguished, partition, Instance};15use rand::{Rng, SeedableRng};16use rand_chacha::ChaCha20Rng;17use serde::{Deserialize, Serialize};18use std::sync::atomic::{AtomicBool, Ordering};19use std::sync::Mutex;20use std::time::Instant;2122#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]23#[serde(rename_all = "lowercase")]24pub enum RhoMode {25 Floyd,26 Brent,27 Dp,28}2930#[derive(Clone, Debug, Serialize, Deserialize)]31pub struct RhoConfig {32 pub r: usize,33 /// 0 = pure r-adding walk; k > 0 = partitions j with j % k == k−1 are doubling partitions34 pub mixed: usize,35 pub neg: bool,36 pub aut3: bool,37 /// remember one point every `cycle_check` steps (detects fruitless cycles of length ≤ cycle_check)38 pub cycle_check: u64,39 pub dp_bits: u32,40 /// walks per thread (batch size of the simultaneous inversion)41 pub walks: usize,42 pub threads: usize,43 pub mode: RhoMode,44 pub seed: u64,45 /// walks longer than this (steps since last distinguished point) are abandoned and restarted46 pub max_walk_len: u64,47}4849impl RhoConfig {50 pub fn default_dp_bits(bits: u32) -> u32 {51 (bits / 4).saturating_sub(2).max(1)52 }53}5455#[derive(Clone, Copy, Debug)]56pub struct Mult {57 pub pt: Pt,58 pub u: u64,59 pub v: u64,60}6162#[derive(Clone, Debug, Default, Serialize, Deserialize)]63pub struct WalkStats {64 pub ops: Counters,65 /// steps of the iteration function actually taken (advances), all walks66 pub iterations: u64,67 pub fruitless_cycles: u64,68 /// fruitless cycles by length (index = length, up to 64)69 pub fruitless_by_length: Vec<u64>,70 /// look-ahead rejections (2-cycle avoidance), each costs one wasted addition71 pub lookahead_rejections: u64,72 pub doubling_fallbacks: u64,73 pub abandoned_walks: u64,74 pub restarts: u64,75 pub dp_found: u64,76 pub useless_collisions: u64,77 pub cpu_s: f64,78}7980impl WalkStats {81 pub fn merge(&mut self, o: &WalkStats) {82 self.ops.merge(&o.ops);83 self.iterations += o.iterations;84 self.fruitless_cycles += o.fruitless_cycles;85 if self.fruitless_by_length.len() < o.fruitless_by_length.len() {86 self.fruitless_by_length.resize(o.fruitless_by_length.len(), 0);87 }88 for (i, v) in o.fruitless_by_length.iter().enumerate() {89 self.fruitless_by_length[i] += v;90 }91 self.lookahead_rejections += o.lookahead_rejections;92 self.doubling_fallbacks += o.doubling_fallbacks;93 self.abandoned_walks += o.abandoned_walks;94 self.restarts += o.restarts;95 self.dp_found += o.dp_found;96 self.useless_collisions += o.useless_collisions;97 self.cpu_s += o.cpu_s;98 }99}100101#[derive(Clone, Copy, Debug)]102pub struct Walk {103 pub p: Pt,104 pub u: u64,105 pub v: u64,106 /// partition override for the look-ahead retry (u32::MAX = none, FORCE_DBL = doubling fallback)107 ovr: u32,108 tries: u32,109 /// steps since the last (re)start110 pub steps: u64,111 saved: Pt,112 saved_step: u64,113 last_escape_x: u64,114}115116const NO_OVR: u32 = u32::MAX;117const FORCE_DBL: u32 = u32::MAX - 1;118119/// Data shared by all threads of one rho run.120pub struct Shared<'a> {121 pub inst: &'a Instance,122 pub cfg: &'a RhoConfig,123 pub mults: Vec<Mult>,124 pub comb_g: CombTable,125 pub comb_q: CombTable,126 pub beta2: u64,127 pub lam2: u64,128 pub dp_mask: u64,129 pub setup_ops: Counters,130}131132impl<'a> Shared<'a> {133 pub fn new(inst: &'a Instance, cfg: &'a RhoConfig) -> Shared<'a> {134 let mut c = Counters::default();135 let mut rng = ChaCha20Rng::seed_from_u64(cfg.seed ^ 0x5EED_0000_0000_0001);136 let comb_g = CombTable::new(&inst.curve, &inst.g, &mut c);137 let comb_q = CombTable::new(&inst.curve, &inst.q, &mut c);138 let n = inst.n.n;139 let mut mults = Vec::with_capacity(cfg.r);140 while mults.len() < cfg.r {141 let u = rng.random_range(1..n);142 let v = rng.random_range(1..n);143 let pt = inst.curve.add(&comb_g.mul(&inst.curve, u, &mut c), &comb_q.mul(&inst.curve, v, &mut c), &mut c);144 if pt.inf {145 continue; // u + v·d ≡ 0 — astronomically unlikely; resample146 }147 mults.push(Mult { pt, u, v });148 }149 let (beta2, lam2) = match (inst.beta, inst.lambda) {150 (Some(b), Some(l)) => (inst.curve.f.sqr(b), inst.n.mul(l, l)),151 _ => (0, 0),152 };153 let dp_mask = if cfg.dp_bits >= 64 { u64::MAX } else { (1u64 << cfg.dp_bits) - 1 };154 Shared { inst, cfg, mults, comb_g, comb_q, beta2, lam2, dp_mask, setup_ops: c }155 }156}157158pub enum Outcome {159 Advanced,160 Rejected,161 Solved(u64),162}163164/// Per-thread walk context.165pub struct Ctx<'a> {166 pub sh: &'a Shared<'a>,167 pub rng: ChaCha20Rng,168 pub st: WalkStats,169 dens: Vec<u64>,170 scratch: Vec<u64>,171 js: Vec<u32>,172 dbls: Vec<bool>,173 special: Vec<bool>,174}175176#[derive(Clone, Copy, Debug)]177pub struct DpRecord {178 pub x: u64,179 pub y: u64,180 pub u: u64,181 pub v: u64,182}183184impl<'a> Ctx<'a> {185 pub fn new(sh: &'a Shared<'a>, thread_id: u64) -> Ctx<'a> {186 let mut seed = [0u8; 32];187 seed[..8].copy_from_slice(&sh.cfg.seed.to_le_bytes());188 seed[8..16].copy_from_slice(&thread_id.to_le_bytes());189 seed[16..24].copy_from_slice(b"rho_walk");190 let w = sh.cfg.walks.max(1);191 Ctx {192 sh,193 rng: ChaCha20Rng::from_seed(seed),194 st: WalkStats { fruitless_by_length: vec![0; 65], ..Default::default() },195 dens: vec![0; w],196 scratch: vec![0; w],197 js: vec![0; w],198 dbls: vec![false; w],199 special: vec![false; w],200 }201 }202203 #[inline(always)]204 fn partition_of(&self, x: u64) -> u32 {205 partition(x, self.sh.cfg.r as u64)206 }207208 #[inline(always)]209 fn is_dbl_partition(&self, j: u32) -> bool {210 let k = self.sh.cfg.mixed;211 k > 0 && (j as usize + 1) % k == 0212 }213214 /// Canonical representative of the class of (P; u, v) and the adjusted scalars.215 #[inline(always)]216 pub fn canon(&mut self, mut p: Pt, mut u: u64, mut v: u64) -> (Pt, u64, u64) {217 let inst = self.sh.inst;218 let f = &inst.curve.f;219 let cfg = self.sh.cfg;220 if cfg.aut3 {221 // candidates x, βx, β²x ↔ P, λP, λ²P ; pick the minimum standard-form x222 let beta = inst.beta.unwrap();223 let lam = inst.lambda.unwrap();224 let x1 = f.mul(beta, p.x);225 let x2 = f.mul(self.sh.beta2, p.x);226 let s0 = f.from_mont(p.x);227 let s1 = f.from_mont(x1);228 let s2 = f.from_mont(x2);229 self.st.ops.mul += 5;230 if s1 < s0 && s1 <= s2 {231 p.x = x1;232 u = inst.n.mul(u, lam);233 v = inst.n.mul(v, lam);234 } else if s2 < s0 && s2 < s1 {235 p.x = x2;236 u = inst.n.mul(u, self.sh.lam2);237 v = inst.n.mul(v, self.sh.lam2);238 }239 }240 if cfg.neg {241 // keep the representative with y_std < p − y_std, i.e. 2·y_std < p242 let ys = f.from_mont(p.y);243 self.st.ops.mul += 1;244 if 2 * ys > f.p {245 p.y = f.neg(p.y);246 u = inst.n.neg(u);247 v = inst.n.neg(v);248 }249 }250 (p, u, v)251 }252253 /// d = num / den mod n, verified against Q. None if den ≡ 0 or verification fails.254 fn try_solve(&mut self, num: u64, den: u64) -> Option<u64> {255 let inst = self.sh.inst;256 let d = inst.n.solve(den, num)?;257 if inst.verify(d, &mut self.st.ops) {258 Some(d)259 } else {260 None261 }262 }263264 /// Fresh random walk start ‖u G + v Q‖. Err(d) if the start happens to be O (then d = −u/v).265 pub fn random_walk(&mut self) -> Result<Walk, u64> {266 let inst = self.sh.inst;267 let n = inst.n.n;268 loop {269 let u = self.rng.random_range(0..n);270 let v = self.rng.random_range(1..n);271 let p = inst.curve.add(272 &self.sh.comb_g.mul(&inst.curve, u, &mut self.st.ops),273 &self.sh.comb_q.mul(&inst.curve, v, &mut self.st.ops),274 &mut self.st.ops,275 );276 self.st.restarts += 1;277 if p.inf {278 if let Some(d) = self.try_solve(inst.n.neg(u), v) {279 return Err(d);280 }281 continue;282 }283 let (p, u, v) = self.canon(p, u, v);284 return Ok(Walk {285 p,286 u,287 v,288 ovr: NO_OVR,289 tries: 0,290 steps: 0,291 saved: Pt::INF,292 saved_step: 0,293 last_escape_x: u64::MAX,294 });295 }296 }297298 #[inline(always)]299 fn choose(&self, w: &Walk) -> (u32, bool) {300 if w.ovr == NO_OVR {301 let j = self.partition_of(w.p.x);302 (j, self.is_dbl_partition(j))303 } else if w.ovr == FORCE_DBL {304 (0, true)305 } else {306 (w.ovr, self.is_dbl_partition(w.ovr))307 }308 }309310 /// Denominator of the slope for the chosen step; 0 marks a special case.311 #[inline(always)]312 fn denominator(&self, w: &Walk, j: u32, is_dbl: bool) -> u64 {313 let f = &self.sh.inst.curve.f;314 if is_dbl {315 f.dbl(w.p.y)316 } else {317 f.sub(self.sh.mults[j as usize].pt.x, w.p.x)318 }319 }320321 /// Apply one step given the inverse of the denominator. Does not touch cycle bookkeeping.322 #[inline(always)]323 fn apply(&mut self, w: &mut Walk, j: u32, is_dbl: bool, inv: u64) -> Outcome {324 let inst = self.sh.inst;325 let f = &inst.curve.f;326 let (r, u, v) = if is_dbl {327 let x2 = f.sqr(w.p.x);328 let num = f.add(f.add(f.add(x2, x2), x2), inst.curve.a);329 let lam = f.mul(num, inv);330 let r = inst.curve.finish_dbl(&w.p, lam);331 self.st.ops.dbl += 1;332 self.st.ops.mul += 4;333 (r, inst.n.dbl(w.u), inst.n.dbl(w.v))334 } else {335 let m = &self.sh.mults[j as usize];336 let lam = f.mul(f.sub(m.pt.y, w.p.y), inv);337 let r = inst.curve.finish_add(&w.p, m.pt.x, lam);338 self.st.ops.add += 1;339 self.st.ops.mul += 3;340 (r, inst.n.add(w.u, m.u), inst.n.add(w.v, m.v))341 };342 let (r, u, v) = self.canon(r, u, v);343 if self.sh.cfg.neg && !is_dbl && self.partition_of(r.x) == j {344 // look-ahead: ‖R + M_j‖ could be ‖P‖ (fruitless 2-cycle) → use the next partition345 self.st.lookahead_rejections += 1;346 w.tries += 1;347 if w.tries as usize >= self.sh.cfg.r {348 w.ovr = FORCE_DBL;349 self.st.doubling_fallbacks += 1;350 } else {351 w.ovr = (j + 1) % self.sh.cfg.r as u32;352 }353 return Outcome::Rejected;354 }355 w.p = r;356 w.u = u;357 w.v = v;358 w.ovr = NO_OVR;359 w.tries = 0;360 w.steps += 1;361 self.st.iterations += 1;362 Outcome::Advanced363 }364365 /// Slow path for a zero denominator (P = ±M_j, or y = 0 for a doubling).366 fn special(&mut self, w: &mut Walk, j: u32, is_dbl: bool) -> Option<u64> {367 let inst = self.sh.inst;368 if is_dbl {369 // y = 0 ⇒ 2P = O: impossible in a prime-order group of odd order; restart defensively.370 return self.restart(w);371 }372 let m = self.sh.mults[j as usize];373 if w.p.y == m.pt.y {374 // P = M_j : u + v d = u_j + v_j d375 if let Some(d) = self.try_solve(inst.n.sub(m.u, w.u), inst.n.sub(w.v, m.v)) {376 return Some(d);377 }378 // same scalars: genuine doubling379 let r = inst.curve.dbl(&w.p, &mut self.st.ops);380 if r.inf {381 return self.restart(w);382 }383 let (r, u, v) = self.canon(r, inst.n.dbl(w.u), inst.n.dbl(w.v));384 w.p = r;385 w.u = u;386 w.v = v;387 w.ovr = NO_OVR;388 w.tries = 0;389 w.steps += 1;390 self.st.iterations += 1;391 None392 } else {393 // P = −M_j : u + v d = −u_j − v_j d ⇒ (v + v_j) d = −(u + u_j)394 if let Some(d) = self.try_solve(inst.n.neg(inst.n.add(w.u, m.u)), inst.n.add(w.v, m.v)) {395 return Some(d);396 }397 self.st.useless_collisions += 1;398 self.restart(w)399 }400 }401402 fn restart(&mut self, w: &mut Walk) -> Option<u64> {403 match self.random_walk() {404 Ok(nw) => {405 *w = nw;406 None407 }408 Err(d) => Some(d),409 }410 }411412 /// One full step of the iteration function with a single (non-batched) inversion.413 pub fn step_single(&mut self, w: &mut Walk) -> Option<u64> {414 loop {415 let (j, is_dbl) = self.choose(w);416 let den = self.denominator(w, j, is_dbl);417 if den == 0 {418 return self.special(w, j, is_dbl);419 }420 let inv = self.sh.inst.curve.f.inv(den);421 self.st.ops.inv += 1;422 self.st.ops.mul += 1;423 match self.apply(w, j, is_dbl, inv) {424 Outcome::Advanced => return None,425 Outcome::Rejected => continue,426 Outcome::Solved(d) => return Some(d),427 }428 }429 }430431 /// Cycle bookkeeping after an advance: detect a fruitless cycle (point equal to the remembered432 /// point) and escape it; otherwise remember the current point every `cycle_check` steps.433 #[inline(always)]434 fn after_advance(&mut self, w: &mut Walk) -> Option<u64> {435 if !w.saved.inf && w.p.x == w.saved.x && w.p.y == w.saved.y {436 let len = w.steps - w.saved_step;437 return self.escape_cycle(w, len);438 }439 if w.steps % self.sh.cfg.cycle_check == 0 {440 w.saved = w.p;441 w.saved_step = w.steps;442 }443 None444 }445446 /// Escape a fruitless cycle of length `len` containing the current point: walk the cycle once447 /// to find its minimum point (smallest Montgomery x), then continue from ‖2·P_min‖.448 pub fn escape_cycle(&mut self, w: &mut Walk, len: u64) -> Option<u64> {449 let inst = self.sh.inst;450 self.st.fruitless_cycles += 1;451 let idx = (len as usize).min(64);452 self.st.fruitless_by_length[idx] += 1;453 let start = w.p;454 let mut best = *w;455 let mut cur = *w;456 for _ in 0..len.max(1) {457 if let Some(d) = self.step_single(&mut cur) {458 return Some(d);459 }460 if cur.p.x < best.p.x || (cur.p.x == best.p.x && cur.p.y < best.p.y) {461 best = cur;462 }463 if cur.p == start {464 break;465 }466 }467 if best.p.x == w.last_escape_x {468 // the escape led back into the same cycle: give up on this trajectory469 return self.restart(w);470 }471 let r = inst.curve.dbl(&best.p, &mut self.st.ops);472 if r.inf {473 return self.restart(w);474 }475 let (r, u, v) = self.canon(r, inst.n.dbl(best.u), inst.n.dbl(best.v));476 w.p = r;477 w.u = u;478 w.v = v;479 w.ovr = NO_OVR;480 w.tries = 0;481 w.steps += 1;482 w.saved = Pt::INF;483 w.last_escape_x = best.p.x;484 self.st.iterations += 1;485 None486 }487488 /// Full single step including cycle bookkeeping (used by Floyd / Brent).489 pub fn step(&mut self, w: &mut Walk) -> Option<u64> {490 if let Some(d) = self.step_single(w) {491 return Some(d);492 }493 self.after_advance(w)494 }495496 /// Batched step of all walks (one simultaneous inversion). Distinguished points are appended to497 /// `dps` and the corresponding walks restarted. Returns Some(d) if the DLP was solved on the way.498 pub fn batch_step(&mut self, walks: &mut [Walk], dps: &mut Vec<DpRecord>) -> Option<u64> {499 let k = walks.len();500 for i in 0..k {501 let (j, is_dbl) = self.choose(&walks[i]);502 let den = self.denominator(&walks[i], j, is_dbl);503 self.js[i] = j;504 self.dbls[i] = is_dbl;505 self.special[i] = den == 0;506 self.dens[i] = if den == 0 { self.sh.inst.curve.f.one } else { den };507 }508 let (m, inv) = self.sh.inst.curve.f.batch_inv(&mut self.dens[..k], &mut self.scratch[..k]);509 self.st.ops.mul += m;510 self.st.ops.inv += inv;511 for i in 0..k {512 let w = &mut walks[i];513 if self.special[i] {514 let (j, is_dbl) = (self.js[i], self.dbls[i]);515 if let Some(d) = self.special(w, j, is_dbl) {516 return Some(d);517 }518 continue;519 }520 let inv = self.dens[i];521 match self.apply(w, self.js[i], self.dbls[i], inv) {522 Outcome::Rejected => continue,523 Outcome::Solved(d) => return Some(d),524 Outcome::Advanced => {}525 }526 if let Some(d) = self.after_advance(w) {527 return Some(d);528 }529 if is_distinguished(w.p.x, self.sh.dp_mask) {530 dps.push(DpRecord { x: w.p.x, y: w.p.y, u: w.u, v: w.v });531 self.st.dp_found += 1;532 if let Some(d) = self.restart(w) {533 return Some(d);534 }535 } else if w.steps > self.sh.cfg.max_walk_len {536 self.st.abandoned_walks += 1;537 if let Some(d) = self.restart(w) {538 return Some(d);539 }540 }541 }542 None543 }544545 /// Resolve a collision of two walks at the same point (or at P and −P).546 fn collide(&mut self, a: &Walk, b: &Walk) -> Option<u64> {547 let inst = self.sh.inst;548 if a.p.y == b.p.y {549 if a.u == b.u && a.v == b.v {550 return None;551 }552 self.try_solve(inst.n.sub(a.u, b.u), inst.n.sub(b.v, a.v))553 } else {554 self.try_solve(inst.n.neg(inst.n.add(a.u, b.u)), inst.n.add(a.v, b.v))555 }556 }557558 /// Floyd's cycle finding on a single walk: x_i vs x_{2i}.559 pub fn floyd(&mut self) -> u64 {560 let mut t = match self.random_walk() {561 Ok(w) => w,562 Err(d) => return d,563 };564 let mut h = t;565 loop {566 if let Some(d) = self.step(&mut t) {567 return d;568 }569 if let Some(d) = self.step(&mut h) {570 return d;571 }572 if let Some(d) = self.step(&mut h) {573 return d;574 }575 if t.p.x == h.p.x {576 if let Some(d) = self.collide(&t, &h) {577 return d;578 }579 // same scalars (fruitless cycle) or degenerate: escape and continue580 let len = if t.p == h.p && t.u == h.u && t.v == h.v { self.cycle_length(&h) } else { None };581 if let Some(len) = len {582 if let Some(d) = self.escape_cycle(&mut h, len) {583 return d;584 }585 t = h;586 } else {587 self.st.useless_collisions += 1;588 match self.random_walk() {589 Ok(w) => {590 t = w;591 h = t;592 }593 Err(d) => return d,594 }595 }596 }597 }598 }599600 /// Brent's cycle finding on a single walk (power-of-two teleporting tortoise).601 pub fn brent(&mut self) -> u64 {602 let mut t = match self.random_walk() {603 Ok(w) => w,604 Err(d) => return d,605 };606 let mut h = t;607 let mut power: u64 = 1;608 let mut lam: u64 = 0;609 loop {610 if let Some(d) = self.step(&mut h) {611 return d;612 }613 lam += 1;614 if t.p.x == h.p.x {615 if let Some(d) = self.collide(&t, &h) {616 return d;617 }618 let len = if t.p == h.p && t.u == h.u && t.v == h.v { self.cycle_length(&h) } else { None };619 if let Some(len) = len {620 if let Some(d) = self.escape_cycle(&mut h, len) {621 return d;622 }623 t = h;624 } else {625 self.st.useless_collisions += 1;626 match self.random_walk() {627 Ok(w) => {628 t = w;629 h = t;630 }631 Err(d) => return d,632 }633 }634 power = 1;635 lam = 0;636 continue;637 }638 if power == lam {639 t = h;640 power *= 2;641 lam = 0;642 }643 }644 }645646 /// Length of the pure-f cycle through w.p, if w.p lies on one of length ≤ 4·cycle_check (+1024).647 /// Floyd/Brent can also meet an earlier point of the *augmented* trajectory (after an escape),648 /// which is not an f-cycle: then None is returned and the caller restarts the walk.649 fn cycle_length(&mut self, w: &Walk) -> Option<u64> {650 let bound = 4 * self.sh.cfg.cycle_check + 1024;651 let mut cur = *w;652 let mut len = 0u64;653 loop {654 self.step_single(&mut cur)?;655 len += 1;656 if cur.p == w.p {657 return Some(len);658 }659 if len > bound {660 return None;661 }662 }663 }664}665666#[derive(Clone, Copy, Debug)]667pub struct DpEntry {668 pub y: u64,669 pub u: u64,670 pub v: u64,671}672673#[derive(Clone, Debug, Serialize, Deserialize)]674pub struct RhoResult {675 pub d: Option<u64>,676 pub stats: WalkStats,677 pub dp_count: u64,678 pub dp_bytes: u64,679 pub table_entries: u64,680 pub threads: usize,681 pub wall_s: f64,682}683684/// Expected number of iterations of a random mapping on a set of size n: √(πn/2).685pub fn expected_iterations(n: u64, class_size: u64) -> f64 {686 (std::f64::consts::PI * (n as f64) / (2.0 * class_size as f64)).sqrt()687}688689pub fn solve_rho(inst: &Instance, cfg: &RhoConfig) -> RhoResult {690 let t0 = Instant::now();691 let sh = Shared::new(inst, cfg);692 let mut total = WalkStats { fruitless_by_length: vec![0; 65], ..Default::default() };693 total.ops.merge(&sh.setup_ops);694695 // trivial cases696 if inst.q.inf {697 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() };698 }699700 match cfg.mode {701 RhoMode::Floyd | RhoMode::Brent => {702 let cpu0 = crate::stats::thread_cpu_s();703 let mut ctx = Ctx::new(&sh, 0);704 let d = if cfg.mode == RhoMode::Floyd { ctx.floyd() } else { ctx.brent() };705 ctx.st.cpu_s = crate::stats::thread_cpu_s() - cpu0;706 total.merge(&ctx.st);707 RhoResult { d: Some(d), stats: total, dp_count: 0, dp_bytes: 0, table_entries: 0, threads: 1, wall_s: t0.elapsed().as_secs_f64() }708 }709 RhoMode::Dp => {710 let threads = cfg.threads.max(1);711 let class = (if cfg.neg { 2 } else { 1 }) * (if cfg.aut3 { 3 } else { 1 });712 let expected_dps = expected_iterations(inst.n.n, class) / (1u64 << cfg.dp_bits.min(63)) as f64;713 let table: Mutex<hashbrown::HashMap<u64, DpEntry, crate::stats::FxBuild>> =714 Mutex::new(hashbrown::HashMap::with_capacity_and_hasher((expected_dps * 2.0) as usize + 1024, Default::default()));715 let found: Mutex<Option<u64>> = Mutex::new(None);716 let stop = AtomicBool::new(false);717 let stats: Mutex<Vec<WalkStats>> = Mutex::new(Vec::new());718719 std::thread::scope(|s| {720 for t in 0..threads {721 let sh = &sh;722 let table = &table;723 let found = &found;724 let stop = &stop;725 let stats = &stats;726 s.spawn(move || {727 let cpu0 = crate::stats::thread_cpu_s();728 let mut ctx = Ctx::new(sh, t as u64);729 let mut walks: Vec<Walk> = Vec::with_capacity(sh.cfg.walks);730 let mut solved: Option<u64> = None;731 for _ in 0..sh.cfg.walks.max(1) {732 match ctx.random_walk() {733 Ok(w) => walks.push(w),734 Err(d) => {735 solved = Some(d);736 break;737 }738 }739 }740 let mut dps: Vec<DpRecord> = Vec::new();741 while solved.is_none() && !stop.load(Ordering::Relaxed) {742 if let Some(d) = ctx.batch_step(&mut walks, &mut dps) {743 solved = Some(d);744 break;745 }746 if !dps.is_empty() {747 let mut tab = table.lock().unwrap();748 for rec in dps.drain(..) {749 match tab.get(&rec.x) {750 None => {751 tab.insert(rec.x, DpEntry { y: rec.y, u: rec.u, v: rec.v });752 }753 Some(e) => {754 let e = *e;755 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 };756 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 };757 if a.u == b.u && a.v == b.v && a.p.y == b.p.y {758 continue; // identical trajectory data, nothing to learn759 }760 if let Some(d) = ctx.collide(&a, &b) {761 solved = Some(d);762 break;763 }764 ctx.st.useless_collisions += 1;765 }766 }767 }768 }769 }770 if let Some(d) = solved {771 let mut f = found.lock().unwrap();772 if f.is_none() {773 *f = Some(d);774 }775 stop.store(true, Ordering::Relaxed);776 }777 ctx.st.cpu_s = crate::stats::thread_cpu_s() - cpu0;778 stats.lock().unwrap().push(ctx.st);779 });780 }781 });782 for st in stats.lock().unwrap().iter() {783 total.merge(st);784 }785 let tab = table.lock().unwrap();786 let entries = tab.len() as u64;787 let bytes = (tab.capacity() as u64) * (std::mem::size_of::<(u64, DpEntry)>() as u64 + 1);788 let d = *found.lock().unwrap();789 RhoResult { d, stats: total, dp_count: entries, dp_bytes: bytes, table_entries: entries, threads, wall_s: t0.elapsed().as_secs_f64() }790 }791 }792}793