Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 16 additions & 17 deletions src/odgp.rs
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ fn solve_odgp_internal(i: &Interval, j: &Interval) -> Vec<ZRootTwo> {
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()
}

Expand All @@ -133,31 +133,30 @@ 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<DRootTwo> {
pub fn solve_scaled_odgp(i: &Interval, j: &Interval, k: i64) -> impl Iterator<Item = DRootTwo> {
let scale = pow_sqrt2(k);
let neg_scale = -scale.clone();
let scaled_j = if k & 1 == 0 {
j.scale(&scale)
} 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<DRootTwo> {
if k == 0 {
let base = beta.renew_denomexp(0);
return solve_odgp_with_parity(i, j, &base)
.map(DRootTwo::from_zroottwo)
.collect();
}
) -> impl Iterator<Item = DRootTwo> {
// 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 {
Expand All @@ -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)]
Expand All @@ -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());
}
}
75 changes: 35 additions & 40 deletions src/tdgp.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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};

Expand All @@ -28,84 +28,79 @@ pub fn solve_tdgp(
) -> Vec<DOmega> {
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() {
continue;
}
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));
}
Expand Down