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