devela/num/float/wrapper/shared.rs
1// devela::num::float::wrapper::shared
2//
3//! Defines all the shared, cross-platform public methods for `Float`.
4//
5
6use crate::{Float, FloatCategory, Sign, cfor, concat as cc, is, stringify as sfy};
7
8/// Implements methods independently of any features
9///
10/// $f: the floating-point type.
11/// $uf: unsigned int type with the same bit-size.
12/// $ie: signed int type used for integer exponentiation.
13/// $ue: unsigned int type used for integer exponentiation and number of terms (u32).
14macro_rules! impl_float_shared {
15 () => {
16 impl_float_shared![ (f32:u32, i32, u32), (f64:u64, i32, u32)];
17 // #[cfg(feature = "nightly_float")] // TODO
18 // impl_float_shared![ (f16:u16, u32), (f128:u128, u32)];
19 };
20
21 ($( ($f:ty:$uf:ty, $ie:ty, $ue:ty)),+) => {
22 $( impl_float_shared![@$f:$uf, $ie, $ue]; )+
23 };
24 (@$f:ty:$uf:ty, $ie:ty, $ue:ty) => {
25 /// # *Common implementations with or without `std` or `libm`*.
26 impl Float<$f> {
27 /// The largest integer less than or equal to itself.
28 /// # Formulation
29 #[doc = crate::_FLOAT_FORMULA_FLOOR!()]
30 pub const fn const_floor(self) -> Float<$f> {
31 let mut result = self.const_trunc().0;
32 if self.0.is_sign_negative() && Float(self.0 - result).abs().0 > <$f>::EPSILON {
33 result -= 1.0;
34 }
35 Float(result)
36 }
37
38 /// The smallest integer greater than or equal to itself.
39 /// # Formulation
40 #[doc = crate::_FLOAT_FORMULA_CEIL!()]
41 pub const fn const_ceil(self) -> Float<$f> {
42 let mut result = self.const_trunc().0;
43 if self.0.is_sign_positive() && Float(self.0 - result).abs().0 > <$f>::EPSILON {
44 result += 1.0;
45 }
46 Float(result)
47 }
48
49 /// The nearest integer to itself, default rounding
50 ///
51 /// This is the default [`round_ties_away`][Self::round_ties_away] implementation.
52 pub const fn const_round(self) -> Float<$f> { self.const_round_ties_away() }
53
54 /// The nearest integer to itself, rounding ties away from `0.0`.
55 ///
56 /// This is the default `round` method implementation,
57 /// without `std` nor `dep_libm`.
58 ///
59 /// # Formulation
60 #[doc = crate::_FLOAT_FORMULA_ROUND_TIES_AWAY!()]
61 pub const fn const_round_ties_away(self) -> Float<$f> {
62 Float(self.0 +
63 Float(0.5 - 0.25 * <$f>::EPSILON).copysign(self.0).0)
64 .const_trunc()
65 }
66
67 /// Returns the nearest integer to `x`, rounding ties to the nearest even integer.
68 /// # Formulation
69 #[doc = crate::_FLOAT_FORMULA_ROUND_TIES_EVEN!()]
70 pub const fn const_round_ties_even(self) -> Float<$f> {
71 let r = self.const_round_ties_away();
72 if r.0 % 2.0 == 0.0 {
73 r
74 } else {
75 #[allow(clippy::float_cmp, reason = "IMPROVE")]
76 if Float(self.0 - r.0).abs().0 == 0.5 { // -0.5 < error_margin
77 Float(r.0 - self.signum().0)
78 } else {
79 r
80 }
81 }
82 }
83
84 /// The integral part.
85 /// This means that non-integer numbers are always truncated towards zero.
86 ///
87 /// # Formulation
88 #[doc = crate::_FLOAT_FORMULA_TRUNC!()]
89 ///
90 /// This implementation uses bitwise manipulation to remove the fractional part
91 /// of the floating-point number. The exponent is extracted, and a mask is
92 /// created to remove the fractional part. The new bits are then used to create
93 /// the truncated floating-point number.
94 pub const fn const_trunc(self) -> Float<$f> {
95 let bits = self.0.to_bits();
96 const BIAS: $ie = Float::<$f>::EXPONENT_BIAS as $ie;
97 const SIG_BITS: $ie = Float::<$f>::SIGNIFICAND_BITS as $ie;
98 const EXP_MASK: $uf = (1 << Float::<$f>::EXPONENT_BITS) - 1;
99
100 #[allow(clippy::cast_possible_wrap)]
101 let exponent = (((bits >> SIG_BITS) & EXP_MASK) as $ie) - BIAS;
102 if exponent < 0 {
103 is![self.0.is_sign_positive(); Float(0.0); Float(-0.0)]
104 } else if exponent < SIG_BITS {
105 let mask = !(((1 as $uf) << (SIG_BITS - exponent)) - 1);
106 let new_bits = bits & mask;
107 Float(<$f>::from_bits(new_bits))
108 } else {
109 self
110 }
111 }
112
113 /// Returns the nearest integer, rounding ties to the nearest odd integer.
114 /// # Formulation
115 #[doc = crate::_FLOAT_FORMULA_ROUND_TIES_ODD!()]
116 #[allow(clippy::float_cmp, reason = "RETHINK (self.0 - r.0).abs() == 0.5")]
117 pub const fn round_ties_odd(self) -> Float<$f> {
118 let r = self.const_round_ties_away();
119 is![r.0 % 2.0 != 0.0; r ;
120 is![(self.0 - r.0).abs() == 0.5; Float(r.0 + self.0.signum()); r]]
121 }
122
123 /// Returns the floating point category of the number.
124 pub const fn classify(self) -> FloatCategory { self.0.classify() }
125
126 /// The fractional part.
127 /// # Formulation
128 #[doc = crate::_FLOAT_FORMULA_FRACT!()]
129 pub const fn const_fract(self) -> Float<$f> { Float(self.0 - self.const_trunc().0) }
130
131 /// The integral and fractional parts.
132 /// # Formulation
133 #[doc = crate::_FLOAT_FORMULA_SPLIT!()]
134 pub const fn const_split(self) -> (Float<$f>, Float<$f>) {
135 let trunc = self.const_trunc();
136 (trunc, Float(self.0 - trunc.0))
137 }
138
139 /// A number that represents its sign, propagating `NaN`.
140 pub const fn signum(self) -> Float<$f> { Float(self.0.signum()) }
141 #[allow(missing_docs)]
142 #[deprecated(since = "0.23.0", note = "use `signum()` instead")]
143 pub const fn const_signum(self) -> Float<$f> { self.signum() }
144 // if self.0.is_nan() { Float(<$f>::NAN) } else { Self::ONE.copysign(self.0) }
145
146 /// A number composed of the magnitude of itself and the `sign` of other.
147 pub const fn copysign(self, sign: $f) -> Float<$f> { Float(self.0.copysign(sign)) }
148 #[allow(missing_docs)]
149 #[deprecated(since = "0.23.0", note = "use `copysign()` instead")]
150 pub const fn const_copysign(self, sign: $f) -> Float<$f> { self.copysign(sign) }
151 // const SIGN_MASK: $uf = <$uf>::MAX / 2 + 1;
152 // const VALUE_MASK: $uf = <$uf>::MAX / 2;
153 // let sign_bit = sign.to_bits() & SIGN_MASK;
154 // let value_bits = self.0.to_bits() & VALUE_MASK;
155 // Float(<$f>::from_bits(value_bits | sign_bit))
156
157 /// Returns the [`Sign`].
158 pub const fn sign(self) -> Sign {
159 if self.is_sign_positive() { Sign::Positive } else { Sign::Negative }
160 }
161
162 /// Returns the [`Sign`], returning [`None`][Sign::None] for zero
163 pub const fn sign_nonzero(self) -> Sign {
164 if self.is_zero() {
165 Sign::None
166 } else if self.is_sign_positive() {
167 Sign::Positive
168 } else {
169 Sign::Negative
170 }
171 }
172
173 /// Returns `true` if `self` has a positive sign.
174 #[must_use]
175 pub const fn is_sign_positive(self) -> bool { self.0.is_sign_positive() }
176
177 /// Returns `true` if `self` has a negative sign.
178 #[must_use]
179 pub const fn is_sign_negative(self) -> bool { self.0.is_sign_negative() }
180
181 /// Returns `true` if `self` is 0.0 or -0.0.
182 #[must_use]
183 pub const fn is_zero(self) -> bool {
184 let non_sign_bits_mask = !(<$uf>::MAX / 2 + 1);
185 (self.0.to_bits() & non_sign_bits_mask) == 0
186 }
187
188 /// Returns `true` if `self` has a positive sign and is not zero.
189 #[must_use]
190 pub const fn is_sign_positive_nonzero(self) -> bool {
191 !self.is_zero() && self.is_sign_positive()
192 }
193
194 /// Returns `true` if `self` has a negative sign and is not zero.
195 #[must_use]
196 pub const fn is_sign_negative_nonzero(self) -> bool {
197 !self.is_zero() && self.is_sign_negative()
198 }
199
200 /// Computes `(x * mul + add)` normally.
201 pub const fn mul_add_fallback(self, mul: $f, add: $f) -> Float<$f> {
202 Float(self.0 * mul + add)
203 }
204
205 /// The euclidean division.
206 // NOTE: [incorrect computations](https://github.com/rust-lang/rust/issues/107904)
207 pub const fn div_euclid(self, other: $f) -> Float<$f> {
208 let q = Float(self.0 / other).const_trunc().0;
209 if self.0 % other < 0.0 {
210 is![other > 0.0; Float(q - 1.0); Float(q + 1.0)]
211 } else {
212 Float(q)
213 }
214 }
215
216 /// The least non-negative remainder of `self` % `other`.
217 // NOTE: [yield inconsistent results](https://github.com/rust-lang/rust/issues/111405)
218 pub const fn rem_euclid(self, other: $f) -> Float<$f> {
219 let r = self.0 % other;
220 is![r < 0.0; Float(r + Float(other).abs().0); Float(r)]
221 }
222
223 /// Returns `self` between `[min..=max]` scaled to a new range `[u..=v]`.
224 ///
225 /// Also known as *remap*.
226 ///
227 /// Values of `self` outside of `[min..=max]` are not clamped
228 /// and will result in extrapolation.
229 ///
230 /// # Formulation
231 #[doc = crate::_FLOAT_FORMULA_SCALE!()]
232 /// # Examples
233 /// ```
234 /// # use devela::Float;
235 #[doc = cc!["assert_eq![Float(45_", sfy![$f], ").scale(0., 360., 0., 1.), 0.125];"]]
236 #[doc = cc!["assert_eq![Float(45_", sfy![$f], ").scale(0., 360., -1., 1.), -0.75];"]]
237 #[doc = cc!["assert_eq![Float(0.125_", sfy![$f], ").scale(0., 1., 0., 360.), 45.];"]]
238 #[doc = cc!["assert_eq![Float(-0.75_", sfy![$f], ").scale(-1., 1., 0., 360.), 45.];"]]
239 /// ```
240 pub const fn scale(self, min: $f, max: $f, u: $f, v: $f) -> Float<$f> {
241 Float((v - u) * (self.0 - min) / (max - min) + u)
242 }
243
244 /// Calculates a linearly interpolated value between `u..=v`
245 /// based on `self` as a percentage between `[0..=1]`.
246 ///
247 /// Values of `self` outside `[0..=1]` are not clamped
248 /// and will result in extrapolation.
249 ///
250 /// # Formulation
251 #[doc = crate::_FLOAT_FORMULA_LERP!()]
252 /// # Example
253 /// ```
254 /// # use devela::Float;
255 #[doc = cc!["assert_eq![Float(0.5_", sfy![$f], ").lerp(40., 80.), 60.];"]]
256 // TODO more examples extrapolated
257 /// ```
258 pub const fn lerp(self, u: $f, v: $f) -> Float<$f> {
259 Float((1.0 - self.0) * u + self.0 * v)
260 }
261
262 /// $ 1 / \sqrt{x} $ the
263 /// [fast inverse square root algorithm](https://en.wikipedia.org/wiki/Fast_inverse_square_root).
264 pub const fn fisr(self) -> Float<$f> {
265 let (mut i, three_halfs, x2) = (self.0.to_bits(), 1.5, self.0 * 0.5);
266 i = Self::FISR_MAGIC - (i >> 1);
267 let y = <$f>::from_bits(i);
268 Float(y * (three_halfs - (x2 * y * y)))
269 }
270
271 /// $ \sqrt{x} $ The square root calculated calling both [`sqrt_fisr`] and [`sqrt_nr`].
272 ///
273 /// [`sqrt_fisr`]: Self::sqrt_fisr
274 /// [`sqrt_nr`]: Self::sqrt_nr
275 pub const fn sqrt_hybrid(self) -> Float<$f> {
276 is![self.0 < 0.0; return Self::NAN; is![self.0 == 0.0; return Self::ZERO]];
277 let y = self.fisr().0; // fast path using fisr + newton refinement
278 let mut x = self.0 * y; // initial estimate: x ~= sqrt(n)
279 // 1 newton iteration some times can be enough:
280 // Float(0.5 * (x + self.0 / x))
281 // But 2 iterations gives much better precision:
282 cfor! { _ in 0..2 => { x = 0.5 * (x + self.0 / x); }}
283 Float(x)
284 }
285
286 /// $ \sqrt{x} $ The square root calculated using the
287 /// [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method).
288 pub const fn sqrt_nr(self) -> Float<$f> {
289 if self.0 < 0.0 {
290 Self::NAN
291 } else if self.0 == 0.0 {
292 Self::ZERO
293 } else {
294 let mut guess = self.0;
295 let mut guess_next = 0.5 * (guess + self.0 / guess);
296 while Self(guess - guess_next).abs().0 > Self::NR_TOLERANCE {
297 guess = guess_next;
298 guess_next = 0.5 * (guess + self.0 / guess);
299 }
300 Float(guess_next)
301 }
302 }
303
304 /// $ \sqrt{x} $ the square root calculated using the
305 /// [fast inverse square root algorithm](https://en.wikipedia.org/wiki/Fast_inverse_square_root).
306 pub const fn sqrt_fisr(self) -> Float<$f> { Float(1.0 / self.fisr().0) }
307
308 /// The hypothenuse (the euclidean distance) using the
309 /// [fast inverse square root algorithm](https://en.wikipedia.org/wiki/Fast_inverse_square_root).
310 ///
311 /// # Formulation
312 #[doc = crate::_FLOAT_FORMULA_HYPOT_FISR!()]
313 pub const fn hypot_fisr(self, y: $f) -> Float<$f> {
314 Float(self.0 * self.0 + y * y).sqrt_fisr()
315 }
316
317 /// The hypothenuse (the euclidean distance) using the
318 /// [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method).
319 ///
320 /// # Formulation
321 #[doc = crate::_FLOAT_FORMULA_HYPOT_NR!()]
322 pub const fn hypot_nr(self, y: $f) -> Float<$f> {
323 Float(self.0 * self.0 + y * y).sqrt_nr()
324 }
325
326 /// $ \sqrt\[3\]{x} $ The cubic root calculated using the
327 /// [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method).
328 pub const fn cbrt_nr(self) -> Float<$f> {
329 is![self.0 == 0.0; return self];
330 let mut guess = self.0;
331 loop {
332 let next_guess = (2.0 * guess + self.0 / (guess * guess)) / 3.0;
333 if Float(next_guess - guess).abs().0 < Self::NR_TOLERANCE {
334 break Float(next_guess);
335 }
336 guess = next_guess;
337 }
338 }
339
340 /// The factorial of the integer value `x`.
341 ///
342 /// The maximum values with a representable result are:
343 /// 34 for `f32` and 170 for `f64`.
344 ///
345 /// Note that precision is poor for large values.
346 pub const fn factorial(x: $ue) -> Float<$f> {
347 let mut result = Self::ONE.0;
348 cfor! { i in 1..{x+1} => { result *= i as $f; }}
349 Float(result)
350 }
351
352 /// The absolute value of `self`.
353 pub const fn abs(self) -> Float<$f> {
354 // let mask = <$uf>::MAX / 2;
355 // Float(<$f>::from_bits(self.0.to_bits() & mask))
356 Self(self.0.abs())
357 }
358
359 /// The negative absolute value of `self` (sets its sign to be negative).
360 pub const fn neg_abs(self) -> Float<$f> {
361 if self.is_sign_negative() { self } else { self.flip_sign() }
362 }
363
364 /// Flips its sign.
365 pub const fn flip_sign(self) -> Float<$f> {
366 let sign_bit_mask = <$uf>::MAX / 2 + 1;
367 Float(<$f>::from_bits(self.0.to_bits() ^ sign_bit_mask))
368 }
369
370 /// Returns the least number smaller than self.
371 pub const fn next_down(self) -> Float<$f> { Self(self.0.next_down()) }
372
373 /// Returns the least number greater than self.
374 pub const fn next_up(self) -> Float<$f> { Self(self.0.next_up()) }
375
376 /// Takes the reciprocal (inverse), $1/x$.
377 pub const fn recip(self) -> Float<$f> { Self(self.0.recip()) }
378
379 /// Converts radians to degrees.
380 pub const fn to_degrees(self) -> Float<$f> { Self(self.0.to_degrees()) }
381
382 /// Converts degrees to radians.
383 pub const fn to_radians(self) -> Float<$f> { Self(self.0.to_radians()) }
384
385 /// Returns itself clamped between `min` and `max`, ignoring `NaN`.
386 ///
387 /// # Example
388 /// ```
389 /// # use devela::Float;
390 #[doc = cc!["assert_eq![Float(50.0_", sfy![$f], ").clamp(40., 80.), 50.];"]]
391 #[doc = cc!["assert_eq![Float(100.0_", sfy![$f], ").clamp(40., 80.), 80.];"]]
392 #[doc = cc!["assert_eq![Float(10.0_", sfy![$f], ").clamp(40., 80.), 40.];"]]
393 /// ```
394 /// See also: [`clamp_nan`][Self::clamp_nan], [`clamp_total`][Self::clamp_total].
395 pub const fn clamp(self, min: $f, max: $f) -> Float<$f> {
396 // self.max(min).min(max)
397 Self(self.0.clamp(min, max))
398 }
399
400 /// The maximum between itself and `other`, ignoring `NaN`.
401 pub const fn max(self, other: $f) -> Float<$f> {
402 // if self.0.is_nan() || self.0 < other { Float(other) } else { self }
403 Self(self.0.max(other))
404 }
405
406 /// The minimum between itself and other, ignoring `NaN`.
407 pub const fn min(self, other: $f) -> Float<$f> {
408 // if other.is_nan() || self.0 < other { self } else { Float(other) }
409 Self(self.0.min(other))
410 }
411 ///
412 #[deprecated(since = "0.23.0", note = "use `clamp()` instead")]
413 pub const fn const_clamp(self, min: $f, max: $f) -> Float<$f> { self.clamp(min, max) }
414 ///
415 #[deprecated(since = "0.23.0", note = "use `max()` instead")]
416 pub const fn const_max(self, other: $f) -> Float<$f> { self.max(other) }
417 ///
418 #[deprecated(since = "0.23.0", note = "use `min()` instead")]
419 pub const fn const_min(self, other: $f) -> Float<$f> { self.min(other) }
420
421 /// Returns itself clamped between `min` and `max`, using total order.
422 ///
423 /// # Example
424 /// ```
425 /// # use devela::Float;
426 #[doc = cc!["assert_eq![Float(50.0_", sfy![$f], ").clamp_total(40., 80.), 50.];"]]
427 #[doc = cc!["assert_eq![Float(100.0_", sfy![$f], ").clamp_total(40., 80.), 80.];"]]
428 #[doc = cc!["assert_eq![Float(10.0_", sfy![$f], ").clamp_total(40., 80.), 40.];"]]
429 /// ```
430 /// See also: [`clamp`][Self::clamp], [`clamp_nan`][Self::clamp_nan].
431 pub const fn clamp_total(self, min: $f, max: $f) -> Float<$f> {
432 Float(crate::Compare(self.0).clamp(min, max))
433 }
434
435 /// Returns the maximum between itself and `other`, using total order.
436 ///
437 /// See also: [`max_nan`][Self::max_nan].
438 pub const fn max_total(self, other: $f) -> Float<$f> {
439 Float(crate::Compare(self.0).max(other))
440 }
441
442 /// Returns the minimum between itself and `other`, using total order.
443 ///
444 /// See also: [`min_nan`][Self::min_nan].
445 pub const fn min_total(self, other: $f) -> Float<$f> {
446 Float(crate::Compare(self.0).min(other))
447 }
448
449 /// Returns itself clamped between `min` and `max`, propagating `NaN`.
450 ///
451 /// # Example
452 /// ```
453 /// # use devela::Float;
454 #[doc = cc!["assert_eq![Float(50.0_", sfy![$f], ").clamp_nan(40., 80.), 50.];"]]
455 #[doc = cc!["assert_eq![Float(100.0_", sfy![$f], ").clamp_nan(40., 80.), 80.];"]]
456 #[doc = cc!["assert_eq![Float(10.0_", sfy![$f], ").clamp_nan(40., 80.), 40.];"]]
457 /// ```
458 /// See also: [`clamp`][Self::clamp], [`clamp_total`][Self::clamp_total].
459 pub const fn clamp_nan(self, min: $f, max: $f) -> Float<$f> {
460 self.max_nan(min).min_nan(max)
461 }
462
463 /// Returns the maximum between itself and `other`, propagating `Nan`.
464 ///
465 /// # Example
466 /// ```
467 /// # use devela::Float;
468 #[doc = cc!["assert_eq![Float(50.0_", sfy![$f], ").max_nan(80.), 80.];"]]
469 #[doc = cc!["assert_eq![Float(100.0_", sfy![$f], ").max_nan(80.), 100.];"]]
470 /// ```
471 /// See also: [`max_total`][Self::max_total].
472 // WAIT: [float_minimum_maximum](https://github.com/rust-lang/rust/issues/91079)
473 #[expect(clippy::float_cmp, reason = "TODO:CHECK:IMPROVE?")]
474 pub const fn max_nan(self, other: $f) -> Float<$f> {
475 if self.0 > other {
476 self
477 } else if self.0 < other {
478 Float(other)
479 } else if self.0 == other {
480 is![self.is_sign_positive() && other.is_sign_negative(); self; Float(other)]
481 } else {
482 // At least one input is NaN. Use `+` to perform NaN propagation and quieting.
483 Float(self.0 + other)
484 }
485 }
486
487 /// Returns the minimum between itself and `other`, propagating `Nan`.
488 ///
489 /// # Example
490 /// ```
491 /// # use devela::Float;
492 #[doc = cc!["assert_eq![Float(50.0_", sfy![$f], ").min_nan(80.), 50.];"]]
493 #[doc = cc!["assert_eq![Float(100.0_", sfy![$f], ").min_nan(80.), 80.];"]]
494 /// ```
495 /// See also: [`min_total`][Self::min_total].
496 // WAIT: [float_minimum_maximum](https://github.com/rust-lang/rust/issues/91079)
497 #[expect(clippy::float_cmp, reason = "TODO:CHECK:IMPROVE?")]
498 pub const fn min_nan(self, other: $f) -> Float<$f> {
499 if self.0 < other {
500 self
501 } else if self.0 > other {
502 Float(other)
503 } else if self.0 == other {
504 is![self.is_sign_negative() && other.is_sign_positive(); self; Float(other)]
505 } else {
506 // At least one input is NaN. Use `+` to perform NaN propagation and quieting.
507 Float(self.0 + other)
508 }
509 }
510
511 /// Calculates the middle point of `self` and `other`.
512 ///
513 /// This returns `NaN` when either argument is `NaN`,
514 /// or if a combination of `+inf` and `-inf` is provided as arguments.
515 pub const fn midpoint(self, other: $f) -> Float<$f> {
516 Self(self.0.midpoint(other))
517 }
518
519 /// Raises itself to the `p` integer power.
520 pub const fn const_powi(self, p: $ie) -> Float<$f> {
521 match p {
522 0 => Self::ONE,
523 1.. => {
524 let mut result = self.0;
525 cfor! { _ in 1..p => { result *= self.0; }}
526 Float(result)
527 }
528 _ => {
529 let mut result = self.0;
530 let abs_p = p.abs();
531 cfor! { _ in 1..abs_p => { result /= self.0; }}
532 Float(result)
533 }
534 }
535 }
536
537 /// Evaluates a polynomial at the `self` point value, using [Horner's method].
538 ///
539 /// Expects a slice of `coefficients` $[a_n, a_{n-1}, ..., a_1, a_0]$
540 /// representing the polynomial $ a_n * x^n + a_{n-1} * x^{(n-1)} + ... + a_1 * x + a_0 $.
541 ///
542 /// # Examples
543 /// ```
544 /// # use devela::Float;
545 /// let coefficients = [2.0, -6.0, 2.0, -1.0];
546 #[doc = cc!["assert_eq![Float(3.0_", sfy![$f], ").eval_poly(&coefficients), 5.0];"]]
547 #[doc = cc!["assert_eq![Float(3.0_", sfy![$f], ").eval_poly(&[]), 0.0];"]]
548 /// ```
549 ///
550 /// [Horner's method]: https://en.wikipedia.org/wiki/Horner%27s_method#Polynomial_evaluation_and_long_division
551 // WAIT: [for-loops in constants](https://github.com/rust-lang/rust/issues/87575)
552 pub const fn eval_poly(self, coefficients: &[$f]) -> Float<$f> {
553 let coef = coefficients;
554 match coef.len() {
555 0 => Float(0.0),
556 1 => Float(coef[0]),
557 _ => {
558 let mut result = coef[0];
559 // for &c in &coef[1..] { result = result * self.0 + c; }
560 cfor! { i in 1..coef.len() => { result = result * self.0 + coef[i]; }}
561 Float(result)
562 }
563 }
564 }
565
566 /// Approximates the derivative of the 1D function `f` at `x` point using step size `h`.
567 ///
568 /// Uses the [finite difference method].
569 ///
570 /// # Formulation
571 #[doc = crate::_FLOAT_FORMULA_DERIVATIVE!()]
572 ///
573 // FLAG_DISABLED:nightly_autodiff
574 // See also the [`autodiff`] attr macro, enabled with the `nightly_autodiff` cfg flag.
575 //
576 /// [finite difference method]: https://en.wikipedia.org/wiki/Finite_difference_method
577 // [`autodiff`]: crate::autodiff
578 pub fn derivative<F>(f: F, x: $f, h: $f) -> Float<$f>
579 where
580 F: Fn($f) -> $f,
581 {
582 Float((f(x + h) - f(x)) / h)
583 }
584
585 /// Approximates the integral of the 1D function `f` over the range `[x, y]`
586 /// using `n` subdivisions.
587 ///
588 /// Uses the [Riemann Sum](https://en.wikipedia.org/wiki/Riemann_sum).
589 ///
590 /// # Formulation
591 #[doc = crate::_FLOAT_FORMULA_INTEGRATE!()]
592 pub fn integrate<F>(f: F, x: $f, y: $f, n: usize) -> Float<$f>
593 where
594 F: Fn($f) -> $f,
595 {
596 let dx = (y - x) / n as $f;
597 Float(
598 (0..n)
599 .map(|i| { let x = x + i as $f * dx; f(x) * dx })
600 .sum()
601 )
602 }
603
604 /// Approximates the partial derivative of the 2D function `f` at point (`x`, `y`)
605 /// with step size `h`, differentiating over `x`.
606 ///
607 /// # Formulation
608 #[doc = crate::_FLOAT_FORMULA_PARTIAL_DERIVATIVE_X!()]
609 pub fn partial_derivative_x<F>(f: F, x: $f, y: $f, h: $f) -> Float<$f>
610 where
611 F: Fn($f, $f) -> $f,
612 {
613 Float((f(x + h, y) - f(x, y)) / h)
614 }
615
616 /// Approximates the partial derivative of the 2D function `f` at point (`x`, `y`)
617 /// with step size `h`, differentiating over `x`.
618 ///
619 /// # Formulation
620 #[doc = crate::_FLOAT_FORMULA_PARTIAL_DERIVATIVE_Y!()]
621 pub fn partial_derivative_y<F>(f: F, x: $f, y: $f, h: $f) -> Float<$f>
622 where
623 F: Fn($f, $f) -> $f,
624 {
625 Float((f(x, y + h) - f(x, y)) / h)
626 }
627 }
628 };
629}
630impl_float_shared!();