@@ -7,7 +7,7 @@ use dashu_int::IBig;
77
88use crate :: common:: { fb_with_prec, ib_to_bf_prec} ;
99use crate :: grid_op:: GridOp ;
10- use crate :: odgp:: { solve_scaled_odgp, solve_scaled_odgp_with_parity } ;
10+ use crate :: odgp:: { solve_scaled_odgp, solve_scaled_odgp_with_parity_k_ne_0 } ;
1111use crate :: region:: { Ellipse , Interval , Rectangle } ;
1212use crate :: ring:: { DOmega , DRootTwo } ;
1313
@@ -28,84 +28,79 @@ pub fn solve_tdgp(
2828) -> Vec < DOmega > {
2929 let mut sol_sufficient = Vec :: with_capacity ( 100 ) ; // Pre-allocate reasonable capacity
3030
31- let sol_x = solve_scaled_odgp ( bbox_a. x . clone ( ) , bbox_b. x . clone ( ) , k + 1 ) ;
32- if sol_x. is_empty ( ) {
33- return vec ! [ ] ;
34- }
31+ let mut sol_x = solve_scaled_odgp ( & bbox_a. x , & bbox_b. x , k + 1 ) ;
3532
36- let sol_y = solve_scaled_odgp (
37- bbox_a
38- . y
39- . fatten ( & ( bbox_a. y . width ( ) / ib_to_bf_prec ( IBig :: from ( 10000 ) ) ) ) ,
40- bbox_b
41- . y
42- . fatten ( & ( bbox_b. y . width ( ) / ib_to_bf_prec ( IBig :: from ( 10000 ) ) ) ) ,
43- k + 1 ,
44- ) ;
45-
46- if sol_y. is_empty ( ) {
47- return vec ! [ ] ;
48- }
33+ let alpha0 = match sol_x. next ( ) {
34+ Some ( val) => val,
35+ None => return vec ! [ ] ,
36+ } ;
4937
50- let alpha0 = & sol_x[ 0 ] ;
5138 let droot_zero = DRootTwo :: from_int ( IBig :: ZERO ) ;
5239 let _k_ibig = IBig :: from ( k) ;
5340 let dx = DRootTwo :: power_of_inv_sqrt2 ( k) ;
5441 let op_g_inv_result = op_g. inv ( ) ;
55- let op_g_inv = op_g_inv_result. as_ref ( ) . unwrap ( ) ;
42+
43+ let op_g_inv = op_g_inv_result. unwrap ( ) ;
5644 let zero_droottwo = DRootTwo :: from_int ( IBig :: ZERO ) ;
45+ let v = op_g_inv * DOmega :: from_droottwo_vector ( & dx, & zero_droottwo, k) ;
5746
58- let v = op_g_inv. clone ( ) * DOmega :: from_droottwo_vector ( & dx, & zero_droottwo, k) ;
5947 let v_conj_sq2 = v. conj_sq2 ( ) ;
48+
49+ let bbox_a_new = bbox_a
50+ . y
51+ . fatten ( & ( bbox_a. y . width ( ) / ib_to_bf_prec ( IBig :: from ( 10000 ) ) ) ) ;
52+ let bbox_b_new = bbox_b
53+ . y
54+ . fatten ( & ( bbox_b. y . width ( ) / ib_to_bf_prec ( IBig :: from ( 10000 ) ) ) ) ;
55+ let sol_y = solve_scaled_odgp ( & bbox_a_new, & bbox_b_new, k + 1 ) ;
56+
6057 for beta in sol_y {
6158 let dx = DRootTwo :: power_of_inv_sqrt2 ( k) ;
62- let z0 = op_g. inv ( ) . as_ref ( ) . unwrap ( ) . clone ( )
63- * DOmega :: from_droottwo_vector ( alpha0, & beta, k + 1 ) ;
64- let v = op_g. inv ( ) . as_ref ( ) . unwrap ( ) . clone ( )
65- * DOmega :: from_droottwo_vector ( & dx, & droot_zero, k) ;
59+ let z0 = op_g. inv ( ) . unwrap ( ) * DOmega :: from_droottwo_vector ( & alpha0, & beta, k + 1 ) ;
60+ let v = op_g. inv ( ) . unwrap ( ) * DOmega :: from_droottwo_vector ( & dx, & droot_zero, k) ;
61+
6662 let t_a = set_a. intersect ( & z0, & v) ;
6763 let t_b = set_b. intersect ( z0. conj_sq2 ( ) , v_conj_sq2) ;
6864 if t_a. is_none ( ) || t_b. is_none ( ) {
6965 continue ;
7066 }
7167 let ( t_a, t_b) = ( t_a. unwrap ( ) , t_b. unwrap ( ) ) ;
7268
73- let parity = ( & beta - alpha0) . mul_by_sqrt2_power_renewing_denomexp ( k) ;
69+ let parity = ( & beta - & alpha0) . mul_by_sqrt2_power_renewing_denomexp ( k) ;
7470 let ( mut int_a, mut int_b) = ( Interval :: new ( t_a. 0 , t_a. 1 ) , Interval :: new ( t_b. 0 , t_b. 1 ) ) ;
7571 let dt_a = {
7672 let ten = ib_to_bf_prec ( IBig :: from ( 10 ) ) ;
73+
74+ let shift_k = IBig :: ONE << ( k as usize ) ;
75+ let width_product = shift_k * int_b. width ( ) ;
7776 let max_val = {
78- let shift_k = IBig :: ONE << ( k as usize ) ;
79- let width_product = shift_k * int_b. width ( ) ;
8077 if ten > width_product {
81- ten. clone ( )
78+ & ten
8279 } else {
83- width_product
80+ & width_product
8481 }
8582 } ;
86- fb_with_prec ( & ten / & max_val)
83+ fb_with_prec ( & ten / max_val)
8784 } ;
8885 let dt_b = {
8986 let ten = ib_to_bf_prec ( IBig :: from ( 10 ) ) ;
87+ let shift_k = IBig :: from ( 1 ) << ( k as usize ) ;
88+ let width_product = shift_k * int_a. width ( ) ;
9089 let max_val = {
91- let shift_k = IBig :: from ( 1 ) << ( k as usize ) ;
92- let width_product = shift_k * int_a. width ( ) ;
9390 if ten > width_product {
94- ten. clone ( )
91+ & ten
9592 } else {
96- width_product
93+ & width_product
9794 }
9895 } ;
99- fb_with_prec ( & ten / & max_val)
96+ fb_with_prec ( & ten / max_val)
10097 } ;
10198
10299 int_a = int_a. fatten ( & dt_a) ;
103100 int_b = int_b. fatten ( & dt_b) ;
104101
105- let sol_t = solve_scaled_odgp_with_parity ( int_a, int_b, 1 , & parity) ;
106- let sol_x = sol_t
107- . into_iter ( )
108- . map ( |alpha| alpha * dx. clone ( ) + alpha0. clone ( ) ) ;
102+ let sol_t = solve_scaled_odgp_with_parity_k_ne_0 ( int_a, int_b, 1 , & parity) ;
103+ let sol_x = sol_t. map ( |alpha| alpha * dx. clone ( ) + alpha0. clone ( ) ) ;
109104 for alpha in sol_x {
110105 sol_sufficient. push ( DOmega :: from_droottwo_vector ( & alpha, & beta, k) ) ;
111106 }
0 commit comments