unitforge/
impl_macros.rs

1#[macro_export]
2macro_rules! impl_quantity {
3    ($type:ident, $unit:ty, $display_unit:expr) => {
4        #[cfg(feature = "pyo3")]
5        use pyo3::{
6            basic::CompareOp, exceptions::PyValueError, pymethods, types::PyType, Bound, IntoPy,
7            Py, PyAny, PyObject, PyRef, PyResult,
8        };
9        use $crate::small_linalg::{Matrix3, Vector3};
10        #[cfg(feature = "pyo3")]
11        use $crate::Quantity;
12
13        #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
14        #[cfg(feature = "pyo3")]
15        #[pyclass]
16        #[derive(Copy, Clone)]
17        pub struct $type {
18            pub(crate) multiplier: f64,
19            pub(crate) power: i32,
20        }
21        #[cfg(not(feature = "pyo3"))]
22        #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
23        #[derive(Copy, Clone)]
24        pub struct $type {
25            pub(crate) multiplier: f64,
26            pub(crate) power: i32,
27        }
28
29        impl PhysicsQuantity for $type {
30            fn as_f64(&self) -> f64 {
31                self.multiplier * 10_f64.powi(self.power)
32            }
33
34            type Unit = $unit;
35
36            fn new(value: f64, unit: Self::Unit) -> $type {
37                if value.is_infinite() {
38                    if value.is_sign_negative() {
39                        return Self::NEG_INFINITY;
40                    } else {
41                        return Self::INFINITY;
42                    }
43                }
44                if value.is_zero() {
45                    return $type {
46                        multiplier: 0.0,
47                        power: 0,
48                    };
49                }
50                let (unit_multiplier, unit_power) = unit.base_per_x();
51                let (multiplier, power) = Self::split_value(value);
52                let r = $type {
53                    multiplier: multiplier * unit_multiplier,
54                    power: power + unit_power,
55                };
56                r
57            }
58
59            fn split_value(v: f64) -> (f64, i32) {
60                if v.is_zero() {
61                    (0.0, 0)
62                } else {
63                    let power = v.abs().log10().floor() as i32;
64                    let multiplier = v / 10f64.powi(power);
65                    (multiplier, power)
66                }
67            }
68
69            fn get_value(&self) -> f64 {
70                self.multiplier * 10_f64.powi(self.power)
71            }
72
73            fn get_power(&self) -> i32 {
74                self.power
75            }
76
77            fn get_multiplier(&self) -> f64 {
78                self.multiplier
79            }
80
81            fn get_tuple(&self) -> (f64, i32) {
82                (self.multiplier, self.power)
83            }
84
85            fn to(&self, unit: Self::Unit) -> f64 {
86                let (unit_multiplier, unit_power) = unit.base_per_x();
87                self.multiplier / unit_multiplier * 10_f64.powi(self.power - unit_power)
88            }
89
90            fn abs(self) -> Self {
91                Self {
92                    multiplier: self.multiplier.abs(),
93                    power: self.power,
94                }
95            }
96
97            fn is_nan(&self) -> bool {
98                self.multiplier.is_nan()
99            }
100
101            fn from_raw(value: f64) -> Self {
102                if value.is_infinite() {
103                    if value.is_sign_negative() {
104                        return Self::NEG_INFINITY;
105                    } else {
106                        return Self::INFINITY;
107                    }
108                }
109                let (multiplier, power) = Self::split_value(value);
110                Self { multiplier, power }
111            }
112
113            fn from_exponential(multiplier: f64, power: i32) -> Self {
114                Self { multiplier, power }
115            }
116
117            fn min(self, other: Self) -> Self {
118                if self < other {
119                    self
120                } else {
121                    other
122                }
123            }
124
125            fn max(self, other: Self) -> Self {
126                if self > other {
127                    self
128                } else {
129                    other
130                }
131            }
132
133            fn is_close(&self, other: &Self, tolerance: &Self) -> bool {
134                (self.as_f64() - other.as_f64()).abs() <= tolerance.as_f64().abs()
135            }
136
137            fn optimize(&mut self) {
138                if self.multiplier.abs() < f64::EPSILON {
139                    return;
140                }
141                let power_on_multiplier = self.multiplier.abs().log10().round() as i32;
142                self.multiplier /= 10_f64.powi(power_on_multiplier);
143                self.power += power_on_multiplier;
144            }
145
146            const INFINITY: Self = Self {
147                multiplier: f64::INFINITY,
148                power: 0,
149            };
150
151            const NEG_INFINITY: Self = Self {
152                multiplier: f64::NEG_INFINITY,
153                power: 0,
154            };
155        }
156
157        impl From<f64> for $type {
158            fn from(value: f64) -> Self {
159                if value.is_infinite() {
160                    if value.is_sign_negative() {
161                        return Self::NEG_INFINITY;
162                    } else {
163                        return Self::INFINITY;
164                    }
165                }
166                let (multiplier, power) = Self::split_value(value);
167                Self { multiplier, power }
168            }
169        }
170
171        #[cfg(feature = "pyo3")]
172        #[pymethods]
173        impl $type {
174            fn __mul__(lhs: PyRef<Self>, rhs: Py<PyAny>) -> PyResult<PyObject> {
175                let py = lhs.py();
176                let rhs_ref = rhs.bind(py);
177                let lhs_py = lhs.into_py(py);
178                let lhs_ref = lhs_py.bind(py);
179                let rhs_quantity = Quantity::from_py_any(rhs_ref)
180                    .map_err(|e| PyValueError::new_err(e.to_string()))?;
181                let lhs_quantity = Quantity::from_py_any(lhs_ref)
182                    .map_err(|e| PyValueError::new_err(e.to_string()))?;
183                match lhs_quantity * rhs_quantity {
184                    Ok(value) => Ok(value.to_pyobject(py)),
185                    Err(_) => Err(PyValueError::new_err(
186                        "Multiplication of given objects is not possible.",
187                    )),
188                }
189            }
190
191            fn __truediv__(lhs: PyRef<Self>, rhs: Py<PyAny>) -> PyResult<PyObject> {
192                let py = lhs.py();
193                let rhs_ref = rhs.bind(py);
194                let lhs_py = lhs.into_py(py);
195                let lhs_ref = lhs_py.bind(py);
196                let rhs_quantity = Quantity::from_py_any(rhs_ref)
197                    .map_err(|e| PyValueError::new_err(e.to_string()))?;
198                let lhs_quantity = Quantity::from_py_any(lhs_ref)
199                    .map_err(|e| PyValueError::new_err(e.to_string()))?;
200                match lhs_quantity / rhs_quantity {
201                    Ok(value) => Ok(value.to_pyobject(py)),
202                    Err(_) => Err(PyValueError::new_err(
203                        "Division of given objects is not possible.",
204                    )),
205                }
206            }
207
208            fn __rmul__(&self, rhs: f64) -> PyResult<Self> {
209                Ok(*self * rhs)
210            }
211        }
212
213        impl fmt::Display for $type {
214            fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
215                let base_unit = $display_unit;
216                let unit_name = if base_unit.name() == "°" {
217                    "°".to_string()
218                } else {
219                    format!(" {}", base_unit.name())
220                };
221                let value = self.to(base_unit);
222                let rounded = if value.is_zero() {
223                    value
224                } else {
225                    let didgits = value.log10().ceil() as i32;
226                    let rounding_multiplier = 10_f64.powi(3 - didgits);
227                    (value * rounding_multiplier).round() as f64 / rounding_multiplier as f64
228                };
229
230                write!(f, "{}{}", rounded, unit_name)
231            }
232        }
233
234        impl Neg for $type {
235            type Output = Self;
236
237            fn neg(self) -> Self::Output {
238                Self {
239                    multiplier: -self.multiplier,
240                    power: self.power,
241                }
242            }
243        }
244
245        impl PartialEq<Self> for $type
246        where
247            $type: PhysicsQuantity,
248        {
249            fn eq(&self, other: &Self) -> bool {
250                self.is_close(
251                    other,
252                    &Self {
253                        multiplier: self.multiplier,
254                        power: self.power.clone() - 9,
255                    },
256                )
257            }
258        }
259
260        impl PartialOrd for $type
261        where
262            $type: PhysicsQuantity,
263        {
264            fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
265                self.as_f64().partial_cmp(&other.as_f64())
266            }
267        }
268
269        impl FromPrimitive for $type {
270            fn from_i64(n: i64) -> Option<Self> {
271                Some(Self::from_raw(n as f64))
272            }
273
274            fn from_u64(n: u64) -> Option<Self> {
275                Some(Self::from_raw(n as f64))
276            }
277
278            fn from_f64(n: f64) -> Option<Self> {
279                Some(Self::from_raw(n))
280            }
281        }
282
283        impl Add for $type {
284            type Output = Self;
285
286            fn add(self, other: Self) -> Self {
287                let common_power = self.power.max(other.power);
288                let multiplier = self.multiplier * 10_f64.powi(self.power - common_power)
289                    + other.multiplier * 10_f64.powi(other.power - common_power);
290                let mut res = Self {
291                    multiplier,
292                    power: common_power,
293                };
294                res.optimize();
295                res
296            }
297        }
298
299        impl Sub for $type {
300            type Output = Self;
301
302            fn sub(self, other: Self) -> Self {
303                let common_power = (self.power + other.power) / 2;
304                let multiplier = self.multiplier * 10_f64.powi(self.power - common_power)
305                    - other.multiplier * 10_f64.powi(other.power - common_power);
306
307                let mut res = Self {
308                    multiplier,
309                    power: common_power,
310                };
311                res.optimize();
312                res
313            }
314        }
315
316        impl Div<f64> for $type {
317            type Output = Self;
318
319            fn div(self, rhs: f64) -> Self::Output {
320                let (rhs_multiplier, rhs_power) = Self::split_value(rhs);
321                Self {
322                    multiplier: self.multiplier / rhs_multiplier,
323                    power: self.power - rhs_power,
324                }
325            }
326        }
327
328        impl Mul<f64> for $type {
329            type Output = Self;
330
331            fn mul(self, rhs: f64) -> Self::Output {
332                let (rhs_multiplier, rhs_power) = Self::split_value(rhs);
333                Self {
334                    multiplier: self.multiplier * rhs_multiplier,
335                    power: self.power + rhs_power,
336                }
337            }
338        }
339
340        impl Mul<$type> for f64 {
341            type Output = $type;
342
343            fn mul(self, rhs: $type) -> Self::Output {
344                rhs * self
345            }
346        }
347
348        impl AddAssign for $type {
349            fn add_assign(&mut self, other: Self) {
350                let common_power = (self.power + other.power) / 2;
351                self.multiplier = self.multiplier * 10_f64.powi(self.power - common_power)
352                    + other.multiplier * 10_f64.powi(other.power - common_power);
353                self.power = common_power;
354                self.optimize();
355            }
356        }
357
358        impl SubAssign for $type {
359            fn sub_assign(&mut self, other: Self) {
360                let common_power = (self.power + other.power) / 2;
361                self.multiplier = self.multiplier * 10_f64.powi(self.power - common_power)
362                    - other.multiplier * 10_f64.powi(other.power - common_power);
363                self.power = common_power;
364                self.optimize();
365            }
366        }
367
368        impl MulAssign<f64> for $type {
369            fn mul_assign(&mut self, rhs: f64) {
370                let (rhs_multiplier, rhs_power) = Self::split_value(rhs);
371                self.multiplier *= rhs_multiplier;
372                self.power += rhs_power;
373            }
374        }
375
376        impl DivAssign<f64> for $type {
377            fn div_assign(&mut self, rhs: f64) {
378                let (rhs_multiplier, rhs_power) = Self::split_value(rhs);
379                self.multiplier /= rhs_multiplier;
380                self.power -= rhs_power;
381            }
382        }
383
384        impl Mul<Vector3<f64>> for $type {
385            type Output = Vector3<$type>;
386
387            fn mul(self, rhs: Vector3<f64>) -> Self::Output {
388                Vector3::new([self * rhs[0], self * rhs[1], self * rhs[2]])
389            }
390        }
391
392        impl Mul<Matrix3<f64>> for $type {
393            type Output = Matrix3<$type>;
394
395            fn mul(self, rhs: Matrix3<f64>) -> Self::Output {
396                Matrix3::new([
397                    [self * rhs[(0, 0)], self * rhs[(0, 1)], self * rhs[(0, 2)]],
398                    [self * rhs[(1, 0)], self * rhs[(1, 1)], self * rhs[(1, 2)]],
399                    [self * rhs[(2, 0)], self * rhs[(2, 1)], self * rhs[(2, 2)]],
400                ])
401            }
402        }
403
404        impl Zero for $type {
405            fn zero() -> Self {
406                Self {
407                    multiplier: 0.0,
408                    power: 0,
409                }
410            }
411
412            fn is_zero(&self) -> bool {
413                self.multiplier == 0.0
414            }
415        }
416
417        impl fmt::Debug for $type {
418            fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
419                fmt::Display::fmt(self, f)
420            }
421        }
422
423        impl QuantityArray2<$unit> for Array2<$type> {
424            fn from_raw(raw: ArrayView2<f64>, unit: $unit) -> Self {
425                let mut res = Array2::zeros(raw.dim());
426                for i in 0..raw.dim().0 {
427                    for j in 0..raw.dim().1 {
428                        res[[i, j]] = <$type>::new(raw[[i, j]], unit);
429                    }
430                }
431                res
432            }
433
434            fn to_raw(&self) -> Array2<f64> {
435                let mut res = Array2::zeros(self.dim());
436                for i in 0..self.dim().0 {
437                    for j in 0..self.dim().1 {
438                        res[[i, j]] = self[[i, j]].as_f64();
439                    }
440                }
441                res
442            }
443
444            fn to(&self, unit: $unit) -> Array2<f64> {
445                let mut res = Array2::zeros(self.dim());
446                for i in 0..self.dim().0 {
447                    for j in 0..self.dim().1 {
448                        res[[i, j]] = self[[i, j]].to(unit);
449                    }
450                }
451                res
452            }
453        }
454
455        impl QuantityArray1<$unit> for Array1<$type> {
456            fn from_raw(raw: ArrayView1<f64>, unit: $unit) -> Self {
457                let mut res = Array1::zeros(raw.dim());
458                for i in 0..raw.dim() {
459                    res[i] = <$type>::new(raw[i], unit);
460                }
461                res
462            }
463
464            fn to_raw(&self) -> Array1<f64> {
465                let mut res = Array1::zeros(self.dim());
466                for i in 0..self.dim() {
467                    res[i] = self[i].as_f64();
468                }
469                res
470            }
471
472            fn to(&self, unit: $unit) -> Array1<f64> {
473                let mut res = Array1::zeros(self.dim());
474                for i in 0..self.dim() {
475                    res[i] = self[i].to(unit);
476                }
477                res
478            }
479        }
480
481        #[cfg(feature = "pyo3")]
482        #[pymethods]
483        impl $type {
484            #[classmethod]
485            #[pyo3(name = "zero")]
486            fn zero_py(_cls: &Bound<'_, PyType>) -> Self {
487                Self::zero()
488            }
489
490            #[new]
491            fn new_py(value: f64, unit: $unit) -> PyResult<Self> {
492                Ok(Self::new(value, unit))
493            }
494
495            fn close_abs(&self, other: PyRef<Self>, threshold: Self) -> bool {
496                (self.clone() - other.clone()).abs() <= threshold
497            }
498
499            fn close_rel(&self, other: PyRef<Self>, threshold: f64) -> bool {
500                let mean = (self.clone() + other.clone()) / 2.;
501                (self.clone() - other.clone()).abs() <= mean * threshold
502            }
503
504            fn __richcmp__(&self, other: PyRef<Self>, op: CompareOp) -> Py<PyAny> {
505                let py = other.py();
506                match op {
507                    CompareOp::Lt => (self.clone() < other.clone()).into_py(py),
508                    CompareOp::Le => (self.clone() <= other.clone()).into_py(py),
509                    CompareOp::Eq => (self.clone() == other.clone()).into_py(py),
510                    CompareOp::Ne => (self.clone() != other.clone()).into_py(py),
511                    CompareOp::Gt => (self.clone() > other.clone()).into_py(py),
512                    CompareOp::Ge => (self.clone() >= other.clone()).into_py(py),
513                }
514            }
515
516            fn __repr__(&self) -> PyResult<String> {
517                Ok(format!("{:?}", self))
518            }
519
520            fn __str__(&self) -> PyResult<String> {
521                Ok(format!("{:?}", self))
522            }
523
524            fn __add__(lhs: PyRef<Self>, rhs: PyRef<Self>) -> PyResult<Self> {
525                Ok(lhs.clone() + rhs.clone())
526            }
527
528            fn __sub__(lhs: PyRef<Self>, rhs: PyRef<Self>) -> PyResult<Self> {
529                Ok(lhs.clone() - rhs.clone())
530            }
531
532            fn __abs__(&self) -> Self {
533                self.abs()
534            }
535
536            #[allow(clippy::wrong_self_convention)]
537            #[pyo3(name = "to")]
538            fn to_py(&self, unit: $unit) -> f64 {
539                self.to(unit)
540            }
541        }
542    };
543}
544
545#[macro_export]
546macro_rules! impl_const {
547    ($type:ident, $name:ident, $multiplier:expr, $power:expr) => {
548        #[cfg(feature = "pyo3")]
549        #[pymethods]
550        impl $type {
551            #[staticmethod]
552            pub fn $name() -> Self {
553                Self {
554                    multiplier: $multiplier,
555                    power: $power,
556                }
557            }
558        }
559        #[cfg(not(feature = "pyo3"))]
560        impl $type {
561            pub fn $name() -> Self {
562                Self {
563                    multiplier: $multiplier,
564                    power: $power,
565                }
566            }
567        }
568    };
569}
570
571#[macro_export]
572macro_rules! impl_div_with_self_to_f64 {
573    ($lhs:ty) => {
574        impl Div<$lhs> for $lhs {
575            type Output = f64;
576
577            fn div(self, rhs: Self) -> Self::Output {
578                (self.multiplier / rhs.multiplier) * 10_f64.powi(self.power - rhs.power)
579            }
580        }
581    };
582}
583
584#[macro_export]
585macro_rules! impl_mul {
586    ($lhs:ty, $rhs:ty, $result:ty) => {
587        impl std::ops::Mul<$rhs> for $lhs {
588            type Output = $result;
589
590            fn mul(self, rhs: $rhs) -> Self::Output {
591                <$result>::from_exponential(
592                    self.multiplier * rhs.multiplier,
593                    self.power + rhs.power,
594                )
595            }
596        }
597
598        impl std::ops::Mul<$lhs> for $rhs {
599            type Output = $result;
600
601            fn mul(self, rhs: $lhs) -> Self::Output {
602                <$result>::from_exponential(
603                    self.multiplier * rhs.multiplier,
604                    self.power + rhs.power,
605                )
606            }
607        }
608
609        impl MulArray1<$rhs> for Array1<$lhs> {
610            type Output = Array1<$result>;
611
612            fn mul_array1(self, rhs: Array1<$rhs>) -> Array1<$result> {
613                self.into_iter()
614                    .zip(rhs.into_iter())
615                    .map(|(force, distance)| force * distance)
616                    .collect()
617            }
618        }
619
620        impl MulArray2<$rhs> for Array2<$lhs> {
621            type Output = Array2<$result>;
622
623            fn mul_array2(self, rhs: Array2<$rhs>) -> Result<Array2<$result>, String> {
624                let mut results = Vec::new();
625
626                for (lhs_row, rhs_row) in self.outer_iter().zip(rhs.outer_iter()) {
627                    let result_row: Array1<$result> = lhs_row
628                        .iter()
629                        .zip(rhs_row.iter())
630                        .map(|(&lhs, &rhs)| lhs * rhs)
631                        .collect();
632                    results.push(result_row);
633                }
634
635                let nrows = results.len();
636                let ncols = if nrows > 0 { results[0].len() } else { 0 };
637                let data: Vec<$result> = results
638                    .into_iter()
639                    .flat_map(|r| {
640                        let (raw_vec, _) = r.into_raw_vec_and_offset();
641                        raw_vec
642                    })
643                    .collect();
644
645                Array2::from_shape_vec((nrows, ncols), data)
646                    .map_err(|_| "Shape mismatch".to_string())
647            }
648        }
649
650        impl MulArray1<$lhs> for Array1<$rhs> {
651            type Output = Array1<$result>;
652
653            fn mul_array1(self, rhs: Array1<$lhs>) -> Array1<$result> {
654                self.into_iter()
655                    .zip(rhs.into_iter())
656                    .map(|(force, distance)| force * distance)
657                    .collect()
658            }
659        }
660
661        impl MulArray2<$lhs> for Array2<$rhs> {
662            type Output = Array2<$result>;
663
664            fn mul_array2(self, rhs: Array2<$lhs>) -> Result<Array2<$result>, String> {
665                let mut results = Vec::new();
666
667                for (force_row, distance_row) in self.outer_iter().zip(rhs.outer_iter()) {
668                    let result_row: Array1<$result> = force_row
669                        .iter()
670                        .zip(distance_row.iter())
671                        .map(|(&force, &distance)| force * distance)
672                        .collect();
673                    results.push(result_row);
674                }
675
676                let nrows = results.len();
677                let ncols = if nrows > 0 { results[0].len() } else { 0 };
678                let data: Vec<$result> = results
679                    .into_iter()
680                    .flat_map(|r| {
681                        let (raw_vec, _) = r.into_raw_vec_and_offset();
682                        raw_vec
683                    })
684                    .collect();
685
686                Array2::from_shape_vec((nrows, ncols), data)
687                    .map_err(|_| "Shape mismatch".to_string())
688            }
689        }
690    };
691}
692
693#[macro_export]
694macro_rules! impl_mul_with_self {
695    ($lhs:ty,$result:ty) => {
696        impl std::ops::Mul<$lhs> for $lhs {
697            type Output = $result;
698
699            fn mul(self, rhs: $lhs) -> Self::Output {
700                <$result>::from_exponential(
701                    self.multiplier * rhs.multiplier,
702                    self.power + rhs.power,
703                )
704            }
705        }
706
707        impl MulArray1<$lhs> for Array1<$lhs> {
708            type Output = Array1<$result>;
709
710            fn mul_array1(self, rhs: Array1<$lhs>) -> Array1<$result> {
711                self.into_iter()
712                    .zip(rhs.into_iter())
713                    .map(|(force, distance)| force * distance)
714                    .collect()
715            }
716        }
717
718        impl MulArray2<$lhs> for Array2<$lhs> {
719            type Output = Array2<$result>;
720
721            fn mul_array2(self, rhs: Array2<$lhs>) -> Result<Array2<$result>, String> {
722                let mut results = Vec::new();
723
724                for (force_row, distance_row) in self.outer_iter().zip(rhs.outer_iter()) {
725                    let result_row: Array1<$result> = force_row
726                        .iter()
727                        .zip(distance_row.iter())
728                        .map(|(&force, &distance)| force * distance)
729                        .collect();
730                    results.push(result_row);
731                }
732
733                let nrows = results.len();
734                let ncols = if nrows > 0 { results[0].len() } else { 0 };
735                let data: Vec<$result> = results
736                    .into_iter()
737                    .flat_map(|r| {
738                        let (raw_vec, _) = r.into_raw_vec_and_offset();
739                        raw_vec
740                    })
741                    .collect();
742
743                Array2::from_shape_vec((nrows, ncols), data)
744                    .map_err(|_| "Shape mismatch".to_string())
745            }
746        }
747    };
748}
749
750#[macro_export]
751macro_rules! impl_div {
752    ($lhs:ty, $rhs:ty, $result:ty) => {
753        impl std::ops::Div<$rhs> for $lhs {
754            type Output = $result;
755
756            fn div(self, rhs: $rhs) -> Self::Output {
757                <$result>::from_exponential(
758                    self.multiplier / rhs.multiplier,
759                    self.power - rhs.power,
760                )
761            }
762        }
763
764        impl DivArray1<$rhs> for Array1<$lhs> {
765            type Output = Array1<$result>;
766
767            fn div_array1(self, rhs: Array1<$rhs>) -> Array1<$result> {
768                self.into_iter()
769                    .zip(rhs.into_iter())
770                    .map(|(force, distance)| force / distance)
771                    .collect()
772            }
773        }
774
775        impl DivArray2<$rhs> for Array2<$lhs> {
776            type Output = Array2<$result>;
777
778            fn div_array2(self, rhs: Array2<$rhs>) -> Result<Array2<$result>, String> {
779                let mut results = Vec::new();
780
781                for (force_row, distance_row) in self.outer_iter().zip(rhs.outer_iter()) {
782                    let result_row: Array1<$result> = force_row
783                        .iter()
784                        .zip(distance_row.iter())
785                        .map(|(&force, &distance)| force / distance)
786                        .collect();
787                    results.push(result_row);
788                }
789
790                let nrows = results.len();
791                let ncols = if nrows > 0 { results[0].len() } else { 0 };
792                let data: Vec<$result> = results
793                    .into_iter()
794                    .flat_map(|r| {
795                        let (raw_vec, _) = r.into_raw_vec_and_offset();
796                        raw_vec
797                    })
798                    .collect();
799
800                Array2::from_shape_vec((nrows, ncols), data)
801                    .map_err(|_| "Shape mismatch".to_string())
802            }
803        }
804    };
805}
806
807#[macro_export]
808macro_rules! impl_sqrt {
809    ($lhs:ty, $res:ty) => {
810        impl Sqrt<$res> for $lhs {
811            fn sqrt(self) -> $res {
812                if self.power % 2 == 0 {
813                    <$res>::from_exponential(self.multiplier.sqrt(), self.power / 2)
814                } else {
815                    <$res>::from_exponential(
816                        self.multiplier.sqrt() * 10_f64.sqrt(),
817                        (self.power - 1) / 2,
818                    )
819                }
820            }
821        }
822        #[cfg(feature = "pyo3")]
823        #[pymethods]
824        impl $lhs {
825            #[pyo3(name = "sqrt")]
826            fn sqrt_py(&self) -> $res {
827                self.sqrt()
828            }
829        }
830    };
831}
832
833#[macro_export]
834macro_rules! impl_mul_relation_with_self {
835    ($lhs:ty, $res:ty) => {
836        impl_mul_with_self!($lhs, $res);
837        impl_sqrt!($res, $lhs);
838        impl_div!($res, $lhs, $lhs);
839    };
840}
841
842#[macro_export]
843macro_rules! impl_mul_relation_with_other {
844    ($lhs:ty, $rhs:ty, $res:ty) => {
845        impl_mul!($lhs, $rhs, $res);
846        impl_div!($res, $lhs, $rhs);
847        impl_div!($res, $rhs, $lhs);
848    };
849}
850
851pub mod macros {
852    pub use crate::impl_const;
853    pub use crate::impl_div;
854    pub use crate::impl_div_with_self_to_f64;
855    pub use crate::impl_mul;
856    pub use crate::impl_mul_relation_with_other;
857    pub use crate::impl_mul_relation_with_self;
858    pub use crate::impl_mul_with_self;
859    pub use crate::impl_quantity;
860    pub use crate::impl_sqrt;
861}