1use crate::astronomy::Location;
11use crate::rata_die::{Moment, RataDie};
12#[allow(unused_imports)]
13use core_maths::*;
14
15pub(crate) trait IntegerRoundings {
16 fn div_ceil(self, rhs: Self) -> Self;
17}
18
19impl IntegerRoundings for i64 {
20 fn div_ceil(self, rhs: Self) -> Self {
22 let d = self / rhs;
23 let r = self % rhs;
24 if (r > 0 && rhs > 0) || (r < 0 && rhs < 0) {
25 d + 1
26 } else {
27 d
28 }
29 }
30}
31
32#[test]
33fn test_div_ceil() {
34 assert_eq!(IntegerRoundings::div_ceil(i64::MIN, 1), i64::MIN);
35 assert_eq!(
36 IntegerRoundings::div_ceil(i64::MIN, 2),
37 -4611686018427387904
38 );
39 assert_eq!(
40 IntegerRoundings::div_ceil(i64::MIN, 3),
41 -3074457345618258602
42 );
43
44 assert_eq!(IntegerRoundings::div_ceil(-10, 1), -10);
45 assert_eq!(IntegerRoundings::div_ceil(-10, 2), -5);
46 assert_eq!(IntegerRoundings::div_ceil(-10, 3), -3);
47
48 assert_eq!(IntegerRoundings::div_ceil(-9, 1), -9);
49 assert_eq!(IntegerRoundings::div_ceil(-9, 2), -4);
50 assert_eq!(IntegerRoundings::div_ceil(-9, 3), -3);
51
52 assert_eq!(IntegerRoundings::div_ceil(-8, 1), -8);
53 assert_eq!(IntegerRoundings::div_ceil(-8, 2), -4);
54 assert_eq!(IntegerRoundings::div_ceil(-8, 3), -2);
55
56 assert_eq!(IntegerRoundings::div_ceil(-2, 1), -2);
57 assert_eq!(IntegerRoundings::div_ceil(-2, 2), -1);
58 assert_eq!(IntegerRoundings::div_ceil(-2, 3), 0);
59
60 assert_eq!(IntegerRoundings::div_ceil(-1, 1), -1);
61 assert_eq!(IntegerRoundings::div_ceil(-1, 2), 0);
62 assert_eq!(IntegerRoundings::div_ceil(-1, 3), 0);
63
64 assert_eq!(IntegerRoundings::div_ceil(0, 1), 0);
65 assert_eq!(IntegerRoundings::div_ceil(0, 2), 0);
66 assert_eq!(IntegerRoundings::div_ceil(0, 3), 0);
67
68 assert_eq!(IntegerRoundings::div_ceil(1, 1), 1);
69 assert_eq!(IntegerRoundings::div_ceil(1, 2), 1);
70 assert_eq!(IntegerRoundings::div_ceil(1, 3), 1);
71
72 assert_eq!(IntegerRoundings::div_ceil(2, 1), 2);
73 assert_eq!(IntegerRoundings::div_ceil(2, 2), 1);
74 assert_eq!(IntegerRoundings::div_ceil(2, 3), 1);
75
76 assert_eq!(IntegerRoundings::div_ceil(8, 1), 8);
77 assert_eq!(IntegerRoundings::div_ceil(8, 2), 4);
78 assert_eq!(IntegerRoundings::div_ceil(8, 3), 3);
79
80 assert_eq!(IntegerRoundings::div_ceil(9, 1), 9);
81 assert_eq!(IntegerRoundings::div_ceil(9, 2), 5);
82 assert_eq!(IntegerRoundings::div_ceil(9, 3), 3);
83
84 assert_eq!(IntegerRoundings::div_ceil(10, 1), 10);
85 assert_eq!(IntegerRoundings::div_ceil(10, 2), 5);
86 assert_eq!(IntegerRoundings::div_ceil(10, 3), 4);
87
88 assert_eq!(IntegerRoundings::div_ceil(i64::MAX, 1), 9223372036854775807);
89 assert_eq!(IntegerRoundings::div_ceil(i64::MAX, 2), 4611686018427387904);
90 assert_eq!(IntegerRoundings::div_ceil(i64::MAX, 3), 3074457345618258603);
91
92 for n in -100..100 {
93 for d in 1..5 {
94 let x1 = IntegerRoundings::div_ceil(n, d);
95 let x2 = (n as f64 / d as f64).ceil();
96 assert_eq!(x1, x2 as i64);
97 }
98 }
99}
100
101pub(crate) fn poly(x: f64, coeffs: &[f64]) -> f64 {
103 coeffs.iter().rev().fold(0.0, |a, c| a * x + c)
104}
105
106pub(crate) fn binary_search(
109 mut l: f64,
110 mut h: f64,
111 test: impl Fn(f64) -> bool,
112 epsilon: f64,
113) -> f64 {
114 debug_assert!(l < h);
115
116 loop {
117 let mid = l + (h - l) / 2.0;
118
119 (l, h) = if test(mid) { (l, mid) } else { (mid, h) };
121
122 if (h - l) < epsilon {
124 return mid;
125 }
126 }
127}
128
129pub(crate) fn next_moment<F>(mut index: Moment, location: Location, condition: F) -> RataDie
130where
131 F: Fn(Moment, Location) -> bool,
132{
133 loop {
134 if condition(index, location) {
135 return index.as_rata_die();
136 }
137 index += 1.0;
138 }
139}
140
141pub(crate) fn next<F>(mut index: RataDie, condition: F) -> RataDie
142where
143 F: Fn(RataDie) -> bool,
144{
145 loop {
146 if condition(index) {
147 return index;
148 }
149 index += 1;
150 }
151}
152
153pub(crate) fn next_u8<F>(mut index: u8, condition: F) -> u8
154where
155 F: Fn(u8) -> bool,
156{
157 loop {
158 if condition(index) {
159 return index;
160 }
161 index += 1;
162 }
163}
164
165pub(crate) fn final_func<F>(mut index: i32, condition: F) -> i32
167where
168 F: Fn(i32) -> bool,
169{
170 while condition(index) {
171 index += 1;
172 }
173 index - 1
174}
175
176#[test]
177fn test_binary_search() {
178 struct TestCase {
179 test_fn: fn(f64) -> bool,
180 range: (f64, f64),
181 expected: f64,
182 }
183
184 let test_cases = [
185 TestCase {
186 test_fn: |x: f64| x >= 4.0,
187 range: (0.0, 10.0),
188 expected: 4.0,
189 },
190 TestCase {
191 test_fn: |x: f64| x * x >= 2.0,
192 range: (0.0, 2.0),
193 expected: 2.0f64.sqrt(),
194 },
195 TestCase {
196 test_fn: |x: f64| x >= -4.0,
197 range: (-10.0, 0.0),
198 expected: -4.0,
199 },
200 TestCase {
201 test_fn: |x: f64| x >= 0.0,
202 range: (0.0, 10.0),
203 expected: 0.0,
204 },
205 TestCase {
206 test_fn: |x: f64| x > 10.0,
207 range: (0.0, 10.0),
208 expected: 10.0,
209 },
210 ];
211
212 for case in test_cases {
213 let result = binary_search(case.range.0, case.range.1, case.test_fn, 1e-4);
214 assert!((result - case.expected).abs() < 0.0001);
215 }
216}
217
218pub(crate) fn invert_angular<F: Fn(f64) -> f64>(f: F, y: f64, r: (f64, f64)) -> f64 {
219 binary_search(r.0, r.1, |x| (f(x) - y).rem_euclid(360.0) < 180.0, 1e-5)
220}
221
222#[test]
223fn test_invert_angular() {
224 struct TestCase {
225 f: fn(f64) -> f64,
226 y: f64,
227 r: (f64, f64),
228 expected: f64,
229 }
230
231 fn f1(x: f64) -> f64 {
232 (2.0 * x).rem_euclid(360.0)
233 }
234
235 fn f2(x: f64) -> f64 {
236 (3.0 * x).rem_euclid(360.0)
237 }
238
239 fn f3(x: f64) -> f64 {
240 (x).rem_euclid(360.0)
241 }
242 let tolerance = 1e-5;
244
245 let test_cases = [
246 TestCase {
247 f: f1,
248 y: 4.0,
249 r: (0.0, 10.0),
250 expected: 4.0,
251 },
252 TestCase {
253 f: f2,
254 y: 6.0,
255 r: (0.0, 20.0),
256 expected: 6.0,
257 },
258 TestCase {
259 f: f3,
260 y: 400.0,
261 r: (0.0, 10.0),
262 expected: 10.0,
263 },
264 TestCase {
265 f: f3,
266 y: 0.0,
267 r: (0.0, 10.0),
268 expected: 0.0,
269 },
270 TestCase {
271 f: f3,
272 y: 10.0,
273 r: (0.0, 10.0),
274 expected: 10.0,
275 },
276 ];
277
278 for case in test_cases {
279 let x = invert_angular(case.f, case.y, case.r);
280 assert!((((case.f)(x)).rem_euclid(360.0) - case.expected).abs() < tolerance);
281 }
282}
283
284#[derive(Copy, Clone, Debug, displaydoc::Display)]
286#[allow(clippy::exhaustive_enums)] pub enum I32CastError {
288 BelowMin,
290 AboveMax,
292}
293
294impl core::error::Error for I32CastError {}
295
296impl I32CastError {
297 pub(crate) const fn saturate(self) -> i32 {
298 match self {
299 I32CastError::BelowMin => i32::MIN,
300 I32CastError::AboveMax => i32::MAX,
301 }
302 }
303}
304
305#[inline]
307pub const fn i64_to_i32(input: i64) -> Result<i32, I32CastError> {
308 if input < i32::MIN as i64 {
309 Err(I32CastError::BelowMin)
310 } else if input > i32::MAX as i64 {
311 Err(I32CastError::AboveMax)
312 } else {
313 Ok(input as i32)
314 }
315}
316
317#[inline]
319pub(crate) fn i64_to_saturated_i32(input: i64) -> i32 {
320 i64_to_i32(input).unwrap_or_else(|i| i.saturate())
321}
322
323#[test]
324fn test_i64_to_saturated_i32() {
325 assert_eq!(i64_to_saturated_i32(i64::MIN), i32::MIN);
326 assert_eq!(i64_to_saturated_i32(-2147483649), -2147483648);
327 assert_eq!(i64_to_saturated_i32(-2147483648), -2147483648);
328 assert_eq!(i64_to_saturated_i32(-2147483647), -2147483647);
329 assert_eq!(i64_to_saturated_i32(-2147483646), -2147483646);
330 assert_eq!(i64_to_saturated_i32(-100), -100);
331 assert_eq!(i64_to_saturated_i32(0), 0);
332 assert_eq!(i64_to_saturated_i32(100), 100);
333 assert_eq!(i64_to_saturated_i32(2147483646), 2147483646);
334 assert_eq!(i64_to_saturated_i32(2147483647), 2147483647);
335 assert_eq!(i64_to_saturated_i32(2147483648), 2147483647);
336 assert_eq!(i64_to_saturated_i32(2147483649), 2147483647);
337 assert_eq!(i64_to_saturated_i32(i64::MAX), i32::MAX);
338}