use log::debug;
use std::fmt::Debug;
use crate::common::ib_to_bf_prec;
use crate::grid_op::{EllipsePair, GridOp};
use crate::math::{floorsqrt, log};
use crate::region::{Ellipse, Rectangle};
use crate::ring::z_root_two::LAMBDA;
use crate::ring::ZOmega;
use crate::tdgp::Region;
use dashu_int::IBig;
use num_traits::Pow;
fn reduction(
ellipse_pair: EllipsePair,
op_g_l: &GridOp,
op_g_r: &GridOp,
new_op_g: &GridOp,
) -> (EllipsePair, GridOp, GridOp, bool) {
(
new_op_g.clone() * ellipse_pair,
op_g_l.clone(),
new_op_g.clone() * op_g_r.clone(),
false,
)
}
fn shift_ellipse_pair(mut ep: EllipsePair, n: i32) -> EllipsePair {
let lambda_n = LAMBDA.pow(&IBig::from(n));
let lambda_inv_n = LAMBDA.pow(&IBig::from(-n));
let lambda_n_f = lambda_n.to_real();
let lambda_inv_n_f = lambda_inv_n.to_real();
ep.a.d[(0, 0)] = &ep.a.d[(0, 0)] * &lambda_inv_n_f;
ep.a.d[(1, 1)] = &ep.a.d[(1, 1)] * &lambda_n_f;
ep.b.d[(0, 0)] = &ep.b.d[(0, 0)] * &lambda_n_f;
ep.b.d[(1, 1)] = &ep.b.d[(1, 1)] * &lambda_inv_n_f;
if n & 1 == 1 {
let b_val = ep.b.d[(0, 1)].clone();
ep.b.d[(0, 1)] = -b_val;
}
ep
}
pub fn step_lemma(
mut ellipse_pair: EllipsePair,
op_g_l: &GridOp,
op_g_r: &GridOp,
verbose: bool,
) -> (EllipsePair, GridOp, GridOp, bool) {
let a = &ellipse_pair.a;
let b = &ellipse_pair.b;
let check0 = ib_to_bf_prec(IBig::from(33971)) / ib_to_bf_prec(IBig::from(1000)); let check1 = ib_to_bf_prec(IBig::from(29437)) / ib_to_bf_prec(IBig::from(1000000)); let check2 = ib_to_bf_prec(IBig::from(58285)) / ib_to_bf_prec(IBig::from(10000)); let check3 = ib_to_bf_prec(IBig::from(17157)) / ib_to_bf_prec(IBig::from(100000)); let check4 = ib_to_bf_prec(IBig::from(2441)) / ib_to_bf_prec(IBig::from(10000)); let check5 = ib_to_bf_prec(IBig::from(40968)) / ib_to_bf_prec(IBig::from(10000)); let check6 = ib_to_bf_prec(IBig::from(16969)) / ib_to_bf_prec(IBig::from(10000)); if b.b() < &ib_to_bf_prec(IBig::ZERO) {
if verbose {
debug!("Z");
}
let op_z = GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ZERO, IBig::ONE),
ZOmega::new(IBig::ZERO, IBig::NEG_ONE, IBig::ZERO, IBig::ZERO),
);
reduction(ellipse_pair, op_g_l, op_g_r, &op_z)
} else if a.bias() * b.bias() < ib_to_bf_prec(IBig::ONE) {
if verbose {
debug!("X");
}
let op_x = GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ONE, IBig::ZERO, IBig::ZERO),
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ZERO, IBig::ONE),
);
reduction(ellipse_pair, op_g_l, op_g_r, &op_x)
} else if ellipse_pair.bias() > check0 || ellipse_pair.bias() < check1 {
let n: IBig = ((log(ellipse_pair.bias()) / log(LAMBDA.to_real()))
/ ib_to_bf_prec(IBig::from(8)))
.round()
.try_into()
.unwrap();
let op_s = GridOp::new(
ZOmega::new(IBig::NEG_ONE, IBig::ZERO, IBig::ONE, IBig::ONE),
ZOmega::new(IBig::ONE, IBig::NEG_ONE, IBig::ONE, IBig::ZERO),
)
.pow(n.clone());
if verbose {
debug!("S (n={})", n);
}
reduction(ellipse_pair, op_g_l, op_g_r, &op_s)
} else if ellipse_pair.skew() <= ib_to_bf_prec(IBig::from(15)) {
(ellipse_pair, op_g_l.clone(), op_g_r.clone(), true)
} else if ellipse_pair.bias() > check2 || ellipse_pair.bias() < check3 {
let n: i32 = ((log(ellipse_pair.bias()) / log(LAMBDA.to_real()))
/ ib_to_bf_prec(IBig::from(4)))
.round()
.try_into()
.unwrap();
ellipse_pair = shift_ellipse_pair(ellipse_pair, n);
if verbose {
debug!("sigma (n={})", n);
}
let op_sigma_l = if n >= 0 {
GridOp::new(
ZOmega::new(IBig::NEG_ONE, IBig::ZERO, IBig::ONE, IBig::ONE),
ZOmega::new(IBig::ZERO, IBig::ONE, IBig::ZERO, IBig::ZERO),
)
.pow(n)
} else {
GridOp::new(
ZOmega::new(IBig::NEG_ONE, IBig::ZERO, IBig::ONE, IBig::NEG_ONE),
ZOmega::new(IBig::ZERO, IBig::ONE, IBig::ZERO, IBig::ZERO),
)
.pow(-n)
};
let op_sigma_r = if n >= 0 {
GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ZERO, IBig::ONE),
ZOmega::new(IBig::ONE, IBig::NEG_ONE, IBig::ONE, IBig::ZERO),
)
.pow(n)
} else {
GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ZERO, IBig::ONE),
ZOmega::new(IBig::ONE, IBig::ONE, IBig::ONE, IBig::ZERO),
)
.pow(-n)
};
(
ellipse_pair,
op_g_l.clone() * op_sigma_l,
op_sigma_r * op_g_r.clone(),
false,
)
} else if (check4.clone()..=check5.clone()).contains(&a.bias())
&& (check4..=check5).contains(&b.bias())
{
if verbose {
debug!("R");
}
let op_r = GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ONE, IBig::ZERO),
ZOmega::new(IBig::ONE, IBig::ZERO, IBig::ZERO, IBig::ZERO),
);
reduction(ellipse_pair, op_g_l, op_g_r, &op_r)
} else if a.b() >= &ib_to_bf_prec(IBig::ZERO) && a.bias() <= check6.clone() {
if verbose {
debug!("K");
}
let op_k = GridOp::new(
ZOmega::new(IBig::NEG_ONE, IBig::NEG_ONE, IBig::ZERO, IBig::ZERO),
ZOmega::new(IBig::ZERO, IBig::NEG_ONE, IBig::ONE, IBig::ZERO),
);
reduction(ellipse_pair, op_g_l, op_g_r, &op_k)
} else if a.b() >= &ib_to_bf_prec(IBig::ZERO) && b.bias() <= check6 {
if verbose {
debug!("K_conj_sq2");
}
let op_kc = GridOp::new(
ZOmega::new(IBig::ONE, IBig::NEG_ONE, IBig::ZERO, IBig::ZERO),
ZOmega::new(IBig::ZERO, IBig::NEG_ONE, IBig::NEG_ONE, IBig::ZERO),
);
reduction(ellipse_pair, op_g_l, op_g_r, &op_kc)
} else if a.b() >= &ib_to_bf_prec(IBig::ZERO) {
let n: IBig =
floorsqrt((a.bias().min(b.bias())) / ib_to_bf_prec(IBig::from(4))).max(IBig::ONE);
if verbose {
debug!("A (n={})", n);
}
let op_a = GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ZERO, IBig::ONE),
ZOmega::new(IBig::ZERO, IBig::ONE, IBig::ZERO, 2 * n),
);
reduction(ellipse_pair, op_g_l, op_g_r, &op_a)
} else {
let n: IBig =
floorsqrt((a.bias().min(b.bias())) / ib_to_bf_prec(IBig::from(2))).max(IBig::ONE);
if verbose {
debug!("B (n={})", n);
}
let op_b = GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ZERO, IBig::ONE),
ZOmega::new(n.clone(), IBig::ONE, -n, IBig::ZERO),
);
reduction(ellipse_pair, op_g_l, op_g_r, &op_b)
}
}
fn to_upright_ellipse_pair(a: &Ellipse, b: &Ellipse, verbose: bool) -> GridOp {
if verbose {
debug!(
"Converting ellipse pair to upright: A = {:?}, B = {:?}",
a, b
);
}
let mut ep = EllipsePair::new(a.normalize(), b.normalize());
let op_i = GridOp::new(
ZOmega::new(IBig::ZERO, IBig::ZERO, IBig::ZERO, IBig::ONE),
ZOmega::new(IBig::ZERO, IBig::ONE, IBig::ZERO, IBig::ZERO),
);
let mut op_g_l = op_i.clone();
let mut op_g_r = op_i;
loop {
let (next_ep, next_l, next_r, end) = step_lemma(ep, &op_g_l, &op_g_r, verbose);
ep = next_ep;
op_g_l = next_l;
op_g_r = next_r;
if end {
break;
}
}
op_g_l * op_g_r
}
pub fn to_upright_set_pair(
set_a: &(impl Region + Debug),
set_b: &(impl Region + Debug),
verbose: bool,
) -> (GridOp, Ellipse, Ellipse, Rectangle, Rectangle) {
let op_g = to_upright_ellipse_pair(&set_a.ellipse(), &set_b.ellipse(), verbose);
let ellipse_pair = op_g.clone() * EllipsePair::new(set_a.ellipse(), set_b.ellipse());
let ellipse_a_upright = ellipse_pair.a;
let ellipse_b_upright = ellipse_pair.b;
let bbox_a = ellipse_a_upright.bbox();
let bbox_b = ellipse_b_upright.bbox();
(op_g, ellipse_a_upright, ellipse_b_upright, bbox_a, bbox_b)
}