diff --git a/src/odgp.rs b/src/odgp.rs index 30a8b38..de0ac7f 100644 --- a/src/odgp.rs +++ b/src/odgp.rs @@ -113,7 +113,7 @@ fn solve_odgp_internal(i: &Interval, j: &Interval) -> Vec { let scaled_j = j.scale(&lambda_conj_sq2_n_f); solve_odgp_internal(&scaled_i, &scaled_j) .into_iter() - .map(|beta| beta * lambda_inv_n.clone()) + .map(|beta| &beta * &lambda_inv_n) .collect() } @@ -133,7 +133,7 @@ pub fn solve_odgp_with_parity( .map(move |alpha| (alpha * ZRootTwo::new(IBig::ZERO, IBig::ONE)) + &p) } -pub fn solve_scaled_odgp(i: Interval, j: Interval, k: i64) -> Vec { +pub fn solve_scaled_odgp(i: &Interval, j: &Interval, k: i64) -> impl Iterator { let scale = pow_sqrt2(k); let neg_scale = -scale.clone(); let scaled_j = if k & 1 == 0 { @@ -141,23 +141,22 @@ pub fn solve_scaled_odgp(i: Interval, j: Interval, k: i64) -> Vec { } else { j.scale(&neg_scale) }; - solve_odgp(i.scale(&scale), scaled_j) - .map(|alpha| DRootTwo::new(alpha, k)) - .collect() + solve_odgp(i.scale(&scale), scaled_j).map(move |alpha| DRootTwo::new(alpha, k)) } -pub fn solve_scaled_odgp_with_parity( +pub fn solve_scaled_odgp_with_parity_k_ne_0( i: Interval, j: Interval, k: i64, beta: &DRootTwo, -) -> Vec { - if k == 0 { - let base = beta.renew_denomexp(0); - return solve_odgp_with_parity(i, j, &base) - .map(DRootTwo::from_zroottwo) - .collect(); - } +) -> impl Iterator { + // Can't do this because the iterators are not of the same type. + // But this function is only called with k == 1. So we don't need the k == 0 branch. + // if k == 0 { + // let base = beta.renew_denomexp(0); + // return solve_odgp_with_parity(i, j, &base) + // .map(DRootTwo::from_zroottwo); + // } let p = beta.renew_denomexp(k).parity(); let offset = if p == IBig::ZERO { @@ -168,8 +167,8 @@ pub fn solve_scaled_odgp_with_parity( let sub_i = i - offset.to_real(); let sub_j = j - offset.conj_sq2().to_real(); - let sol = solve_scaled_odgp(sub_i, sub_j, k - 1); - sol.into_iter().map(|a| a + offset.clone()).collect() + let sol = solve_scaled_odgp(&sub_i, &sub_j, k - 1); + sol.map(move |a| a + offset.clone()) } #[cfg(test)] @@ -192,7 +191,7 @@ mod tests { #[test] fn test_use_empty_interval() { let (inti, intj) = create_empty_interval(); - let result = solve_scaled_odgp(inti, intj, 2); - assert!(result.is_empty()); + let mut result = solve_scaled_odgp(&inti, &intj, 2); + assert!(result.next().is_none()); } } diff --git a/src/tdgp.rs b/src/tdgp.rs index 013902f..167cf29 100644 --- a/src/tdgp.rs +++ b/src/tdgp.rs @@ -7,7 +7,7 @@ use dashu_int::IBig; use crate::common::{fb_with_prec, ib_to_bf_prec}; use crate::grid_op::GridOp; -use crate::odgp::{solve_scaled_odgp, solve_scaled_odgp_with_parity}; +use crate::odgp::{solve_scaled_odgp, solve_scaled_odgp_with_parity_k_ne_0}; use crate::region::{Ellipse, Interval, Rectangle}; use crate::ring::{DOmega, DRootTwo}; @@ -28,41 +28,37 @@ pub fn solve_tdgp( ) -> Vec { let mut sol_sufficient = Vec::with_capacity(100); // Pre-allocate reasonable capacity - let sol_x = solve_scaled_odgp(bbox_a.x.clone(), bbox_b.x.clone(), k + 1); - if sol_x.is_empty() { - return vec![]; - } + let mut sol_x = solve_scaled_odgp(&bbox_a.x, &bbox_b.x, k + 1); - let sol_y = solve_scaled_odgp( - bbox_a - .y - .fatten(&(bbox_a.y.width() / ib_to_bf_prec(IBig::from(10000)))), - bbox_b - .y - .fatten(&(bbox_b.y.width() / ib_to_bf_prec(IBig::from(10000)))), - k + 1, - ); - - if sol_y.is_empty() { - return vec![]; - } + let alpha0 = match sol_x.next() { + Some(val) => val, + None => return vec![], + }; - let alpha0 = &sol_x[0]; let droot_zero = DRootTwo::from_int(IBig::ZERO); let _k_ibig = IBig::from(k); let dx = DRootTwo::power_of_inv_sqrt2(k); let op_g_inv_result = op_g.inv(); - let op_g_inv = op_g_inv_result.as_ref().unwrap(); + + let op_g_inv = op_g_inv_result.unwrap(); let zero_droottwo = DRootTwo::from_int(IBig::ZERO); + let v = op_g_inv * DOmega::from_droottwo_vector(&dx, &zero_droottwo, k); - let v = op_g_inv.clone() * DOmega::from_droottwo_vector(&dx, &zero_droottwo, k); let v_conj_sq2 = v.conj_sq2(); + + let bbox_a_new = bbox_a + .y + .fatten(&(bbox_a.y.width() / ib_to_bf_prec(IBig::from(10000)))); + let bbox_b_new = bbox_b + .y + .fatten(&(bbox_b.y.width() / ib_to_bf_prec(IBig::from(10000)))); + let sol_y = solve_scaled_odgp(&bbox_a_new, &bbox_b_new, k + 1); + for beta in sol_y { let dx = DRootTwo::power_of_inv_sqrt2(k); - let z0 = op_g.inv().as_ref().unwrap().clone() - * DOmega::from_droottwo_vector(alpha0, &beta, k + 1); - let v = op_g.inv().as_ref().unwrap().clone() - * DOmega::from_droottwo_vector(&dx, &droot_zero, k); + let z0 = op_g.inv().unwrap() * DOmega::from_droottwo_vector(&alpha0, &beta, k + 1); + let v = op_g.inv().unwrap() * DOmega::from_droottwo_vector(&dx, &droot_zero, k); + let t_a = set_a.intersect(&z0, &v); let t_b = set_b.intersect(z0.conj_sq2(), v_conj_sq2); if t_a.is_none() || t_b.is_none() { @@ -70,42 +66,41 @@ pub fn solve_tdgp( } let (t_a, t_b) = (t_a.unwrap(), t_b.unwrap()); - let parity = (&beta - alpha0).mul_by_sqrt2_power_renewing_denomexp(k); + let parity = (&beta - &alpha0).mul_by_sqrt2_power_renewing_denomexp(k); let (mut int_a, mut int_b) = (Interval::new(t_a.0, t_a.1), Interval::new(t_b.0, t_b.1)); let dt_a = { let ten = ib_to_bf_prec(IBig::from(10)); + + let shift_k = IBig::ONE << (k as usize); + let width_product = shift_k * int_b.width(); let max_val = { - let shift_k = IBig::ONE << (k as usize); - let width_product = shift_k * int_b.width(); if ten > width_product { - ten.clone() + &ten } else { - width_product + &width_product } }; - fb_with_prec(&ten / &max_val) + fb_with_prec(&ten / max_val) }; let dt_b = { let ten = ib_to_bf_prec(IBig::from(10)); + let shift_k = IBig::from(1) << (k as usize); + let width_product = shift_k * int_a.width(); let max_val = { - let shift_k = IBig::from(1) << (k as usize); - let width_product = shift_k * int_a.width(); if ten > width_product { - ten.clone() + &ten } else { - width_product + &width_product } }; - fb_with_prec(&ten / &max_val) + fb_with_prec(&ten / max_val) }; int_a = int_a.fatten(&dt_a); int_b = int_b.fatten(&dt_b); - let sol_t = solve_scaled_odgp_with_parity(int_a, int_b, 1, &parity); - let sol_x = sol_t - .into_iter() - .map(|alpha| alpha * dx.clone() + alpha0.clone()); + let sol_t = solve_scaled_odgp_with_parity_k_ne_0(int_a, int_b, 1, &parity); + let sol_x = sol_t.map(|alpha| alpha * dx.clone() + alpha0.clone()); for alpha in sol_x { sol_sufficient.push(DOmega::from_droottwo_vector(&alpha, &beta, k)); }