1use super::value::Ival;
2use crate::mpfr::{
3 mpfr_add, mpfr_div, mpfr_fmod, mpfr_max, mpfr_min, mpfr_mul, mpfr_neg, mpfr_remainder,
4 mpfr_round, mpfr_sign, mpfr_sub, mpfr_trunc, zero,
5};
6use rug::{Assign, Float, float::Round};
7
8impl Ival {
9 pub fn fmod_assign(&mut self, x: &Ival, y: &Ival) {
11 let y_lo = y.lo.as_float();
12 let y_hi = y.hi.as_float();
13 let y_straddles_zero = (y_lo.is_sign_negative() || y_lo.is_zero())
14 && (y_hi.is_sign_positive() || y_hi.is_zero());
15
16 self.err = x.err.union(&y.err);
17 if y_straddles_zero {
18 self.err.partial = true;
19 }
20 if y_lo.is_zero() && y_hi.is_zero() {
21 self.err.total = true;
22 }
23
24 let mut y_abs = Ival::zero(y.max_prec());
25 y_abs.exact_fabs_assign(y);
26
27 let x_hi = x.hi.as_float();
28 let x_lo = x.lo.as_float();
29 let hi_is_neg = mpfr_sign(x_hi) == -1 && !x_hi.is_zero();
30 let lo_is_pos = mpfr_sign(x_lo) == 1 || x_lo.is_zero();
31
32 if hi_is_neg {
33 let mut neg_x = Ival::zero(x.max_prec());
34 neg_x.exact_neg_assign(x);
35 self.fmod_pos_assign(&neg_x, &y_abs);
36 self.neg_inplace();
37 } else if lo_is_pos {
38 self.fmod_pos_assign(x, &y_abs);
39 } else {
40 let zero_val = zero(x.prec());
41 let (neg, pos) = x.split_at(&zero_val);
42
43 let mut neg_x = Ival::zero(neg.max_prec());
44 neg_x.exact_neg_assign(&neg);
45
46 let mut neg_result = Ival::zero(self.prec());
47 let mut pos_result = Ival::zero(self.prec());
48
49 neg_result.fmod_pos_assign(&neg_x, &y_abs);
50 pos_result.fmod_pos_assign(&pos, &y_abs);
51
52 neg_result.neg_inplace();
53
54 self.assign_from(&pos_result);
55 self.union_assign(neg_result);
56 }
57 }
58
59 pub fn remainder_assign(&mut self, x: &Ival, y: &Ival) {
61 let y_lo = y.lo.as_float();
62 let y_hi = y.hi.as_float();
63 let y_straddles_zero = (y_lo.is_sign_negative() || y_lo.is_zero())
64 && (y_hi.is_sign_positive() || y_hi.is_zero());
65
66 self.err = x.err.union(&y.err);
67 if y_straddles_zero {
68 self.err.partial = true;
69 }
70 if y_lo.is_zero() && y_hi.is_zero() {
71 self.err.total = true;
72 }
73
74 let mut y_abs = Ival::zero(y.max_prec());
75 y_abs.exact_fabs_assign(y);
76
77 let x_hi = x.hi.as_float();
78 let x_lo = x.lo.as_float();
79 let hi_is_neg = mpfr_sign(x_hi) == -1 && !x_hi.is_zero();
80 let lo_is_pos = mpfr_sign(x_lo) == 1 || x_lo.is_zero();
81
82 if hi_is_neg {
83 let mut neg_x = Ival::zero(x.max_prec());
84 neg_x.exact_neg_assign(x);
85 self.remainder_pos_assign(&neg_x, &y_abs);
86 self.neg_inplace();
87 } else if lo_is_pos {
88 self.remainder_pos_assign(x, &y_abs);
89 } else {
90 let zero_val = zero(x.prec());
91 let (neg, pos) = x.split_at(&zero_val);
92
93 let mut neg_x = Ival::zero(neg.max_prec());
94 neg_x.exact_neg_assign(&neg);
95
96 let mut neg_result = Ival::zero(self.prec());
97 let mut pos_result = Ival::zero(self.prec());
98
99 neg_result.remainder_pos_assign(&neg_x, &y_abs);
100 pos_result.remainder_pos_assign(&pos, &y_abs);
101
102 neg_result.neg_inplace();
103
104 self.assign_from(&pos_result);
105 self.union_assign(neg_result);
106 }
107 }
108
109 fn fmod_pos_assign(&mut self, x: &Ival, y: &Ival) {
110 let prec = self.prec();
111 let x_lo = x.lo.as_float();
112 let x_hi = x.hi.as_float();
113 let y_lo = y.lo.as_float();
114 let y_hi = y.hi.as_float();
115
116 let mut tmp_div = zero(prec);
117 mpfr_div(x_lo, y_hi, &mut tmp_div, Round::Down);
118 let mut a = zero(prec);
119 mpfr_trunc(&tmp_div, &mut a, Round::Down);
120
121 mpfr_div(x_hi, y_hi, &mut tmp_div, Round::Up);
122 let mut b = zero(prec);
123 mpfr_trunc(&tmp_div, &mut b, Round::Up);
124
125 if a == b {
126 mpfr_div(x_hi, y_hi, &mut tmp_div, Round::Down);
127 let mut c = zero(prec);
128 mpfr_trunc(&tmp_div, &mut c, Round::Down);
129
130 mpfr_div(x_hi, y_lo, &mut tmp_div, Round::Up);
131 let mut d = zero(prec);
132 mpfr_trunc(&tmp_div, &mut d, Round::Up);
133
134 if c == d {
135 let mut lo = zero(prec);
136 let mut hi = zero(prec);
137
138 mpfr_fmod(x_lo, y_hi, &mut lo, Round::Down);
139 mpfr_fmod(x_hi, y_lo, &mut hi, Round::Up);
140
141 self.lo.as_float_mut().assign(&lo);
142 self.lo.immovable = false;
143 self.hi.as_float_mut().assign(&hi);
144 self.hi.immovable = false;
145 } else {
146 let mut c_plus_1 = zero(prec);
147 let one = Float::with_val(prec, 1);
148 mpfr_add(&c, &one, &mut c_plus_1, Round::Down);
149
150 let mut hi = zero(prec);
151 mpfr_div(x_hi, &c_plus_1, &mut hi, Round::Up);
152
153 self.lo.as_float_mut().assign(0);
154 self.lo.immovable = false;
155 self.hi.as_float_mut().assign(&hi);
156 self.hi.immovable = false;
157 }
158 } else {
159 self.lo.as_float_mut().assign(0);
160 self.lo.immovable = false;
161 self.hi.as_float_mut().assign(y_hi);
162 self.hi.immovable = false;
163 }
164 }
165
166 fn remainder_pos_assign(&mut self, x: &Ival, y: &Ival) {
167 let prec = self.prec();
168 let x_lo = x.lo.as_float();
169 let x_hi = x.hi.as_float();
170 let y_lo = y.lo.as_float();
171 let y_hi = y.hi.as_float();
172
173 let mut tmp_div = zero(prec);
174 mpfr_div(x_lo, y_hi, &mut tmp_div, Round::Down);
175 let mut a = zero(prec);
176 mpfr_round(&tmp_div, &mut a, Round::Down);
177
178 mpfr_div(x_hi, y_hi, &mut tmp_div, Round::Up);
179 let mut b = zero(prec);
180 mpfr_round(&tmp_div, &mut b, Round::Up);
181
182 if a == b {
183 mpfr_div(x_hi, y_hi, &mut tmp_div, Round::Down);
184 let mut c = zero(prec);
185 mpfr_round(&tmp_div, &mut c, Round::Down);
186
187 mpfr_div(x_hi, y_lo, &mut tmp_div, Round::Up);
188 let mut d = zero(prec);
189 mpfr_round(&tmp_div, &mut d, Round::Up);
190
191 if c == d {
192 let mut lo = zero(prec);
193 let mut hi = zero(prec);
194
195 mpfr_remainder(x_lo, y_hi, &mut lo, Round::Down);
196 mpfr_remainder(x_hi, y_lo, &mut hi, Round::Up);
197
198 self.lo.as_float_mut().assign(&lo);
199 self.lo.immovable = false;
200 self.hi.as_float_mut().assign(&hi);
201 self.hi.immovable = false;
202 } else {
203 let half = Float::with_val(prec, 0.5);
204 let mut c_plus_half = zero(prec);
205 mpfr_add(&c, &half, &mut c_plus_half, Round::Down);
206
207 let mut tmp1 = zero(prec);
208 mpfr_div(x_hi, &c_plus_half, &mut tmp1, Round::Up);
209
210 let two = Float::with_val(prec, 2);
211 let mut y_star_hi = zero(prec);
212 mpfr_div(&tmp1, &two, &mut y_star_hi, Round::Up);
213
214 let mut c_times_yhi = zero(prec);
215 mpfr_mul(&c, y_hi, &mut c_times_yhi, Round::Up);
216
217 let mut x_lo_minus_c_yhi = zero(prec);
218 mpfr_sub(x_lo, &c_times_yhi, &mut x_lo_minus_c_yhi, Round::Down);
219
220 let mut yhi_half = zero(prec);
221 mpfr_div(y_hi, &two, &mut yhi_half, Round::Up);
222 let mut neg_yhi_half = zero(prec);
223 mpfr_neg(&yhi_half, &mut neg_yhi_half, Round::Down);
224
225 let mut y_star_lo = zero(prec);
226 mpfr_max(
227 &x_lo_minus_c_yhi,
228 &neg_yhi_half,
229 &mut y_star_lo,
230 Round::Down,
231 );
232
233 let mut neg_y_star_hi = zero(prec);
234 mpfr_neg(&y_star_hi, &mut neg_y_star_hi, Round::Down);
235
236 let mut final_lo = zero(prec);
237 mpfr_min(&y_star_lo, &neg_y_star_hi, &mut final_lo, Round::Down);
238
239 self.lo.as_float_mut().assign(&final_lo);
240 self.lo.immovable = false;
241 self.hi.as_float_mut().assign(&y_star_hi);
242 self.hi.immovable = false;
243 }
244 } else {
245 let two = Float::with_val(prec, 2);
246 let mut y_half = zero(prec);
247 mpfr_div(y_hi, &two, &mut y_half, Round::Up);
248
249 let mut neg_y_half = zero(prec);
250 mpfr_neg(&y_half, &mut neg_y_half, Round::Down);
251
252 self.lo.as_float_mut().assign(&neg_y_half);
253 self.lo.immovable = false;
254 self.hi.as_float_mut().assign(&y_half);
255 self.hi.immovable = false;
256 }
257 }
258}