hoomd_manifold/biquaternion.rs
1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Implement [`Biquaternion`] and a four-dimensional matrix representation
5//! of SO(3,1).
6
7use num::complex::Complex;
8use rand::{
9 Rng,
10 distr::{Distribution, StandardUniform, Uniform},
11};
12use serde::{Deserialize, Serialize};
13use std::{
14 array, fmt,
15 iter::zip,
16 ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Sub, SubAssign},
17};
18
19use crate::{Error, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski};
20#[expect(unused_imports, reason = "Needed for doc link")]
21use hoomd_vector::Quaternion;
22
23/// A quaternion with complex coefficients.
24///
25/// Biquaternions are the set of numbers $`a + b\mathbf{i} + c\mathbf{j} + d\mathbf{k}`$
26/// where $`a,b,c,d`$ are complex numbers and $`\{1,\mathbf{i},\mathbf{j},\mathbf{k}\}`$
27/// are the quaternion algebra. Biquaternions can be thought of as a
28/// generalization of quaternions which allow for complex coefficients.
29/// Analogous to quaternions and SO(3), biquaternions furnish a representation
30/// of SO(3,1)
31///
32/// ## Construction of Biquaternions
33///
34/// Create a biquaternion from an array of four complex numbers. Note that
35/// components are in the order $`[\mathbf{i},\mathbf{j},\mathbf{k},1]`$
36/// (i.e., the scalar component is at the end)
37/// ```
38/// use hoomd_manifold::Biquaternion;
39/// use num::complex::Complex;
40///
41/// let q = Biquaternion::from([
42/// Complex::new(1.0, 4.0),
43/// Complex::new(2.0, 3.0),
44/// Complex::new(3.0, 2.0),
45/// Complex::new(4.0, 1.0),
46/// ]);
47/// assert_eq!(4.0, q.components[0].im);
48/// ```
49///
50/// ## Operations with Biquaternions.
51///
52/// Similar to [`Quaternion`], biquaternions support vector operations
53/// (addition, multiplication by a scalar, etc.):
54/// ```
55/// use hoomd_manifold::Biquaternion;
56/// use num::complex::Complex;
57///
58/// let mut a = Biquaternion::from([
59/// Complex::new(1.0, 0.0),
60/// Complex::new(2.0, 0.0),
61/// Complex::new(3.0, 0.0),
62/// Complex::new(0.0, 1.0),
63/// ]);
64/// let mut b = Biquaternion::from([
65/// Complex::new(0.0, 4.0),
66/// Complex::new(0.0, 3.0),
67/// Complex::new(0.0, 2.0),
68/// Complex::new(1.0, 0.0),
69/// ]);
70/// b /= 2.0;
71/// let mut c = a + b;
72/// assert_eq!(Complex::new(1.0, 2.0), c.components[0]);
73/// ```
74///
75/// Biquaternions also support the following operations:
76///
77/// Hamiltonian conjugate/ biconjugate:
78/// Denoted by the method "bar", the Hamiltonian conjugate multiplies the
79/// vector part of the biquaternion by -1.0.
80/// ```
81/// use hoomd_manifold::Biquaternion;
82/// use num::complex::Complex;
83///
84/// let q = Biquaternion::from([
85/// Complex::new(-1.0, 0.0),
86/// Complex::new(-1.0, 2.0),
87/// Complex::new(1.0, 0.0),
88/// Complex::new(1.0, 0.0),
89/// ]);
90/// let p = Biquaternion::from([
91/// Complex::new(1.0, 0.0),
92/// Complex::new(1.0, -2.0),
93/// Complex::new(-1.0, 0.0),
94/// Complex::new(1.0, 0.0),
95/// ]);
96///
97/// assert_eq!(p, q.bar());
98/// ```
99///
100/// Complex conjugation:
101/// Denoted by method "conj", takes the complex conjugate of all components of
102/// the biquaternion
103/// ```
104/// use hoomd_manifold::Biquaternion;
105/// use num::complex::Complex;
106///
107/// let q = Biquaternion::from([
108/// Complex::new(1.0, 8.0),
109/// Complex::new(2.0, 7.0),
110/// Complex::new(3.0, 6.0),
111/// Complex::new(4.0, 5.0),
112/// ]);
113/// let p = Biquaternion::from([
114/// Complex::new(1.0, -8.0),
115/// Complex::new(2.0, -7.0),
116/// Complex::new(3.0, -6.0),
117/// Complex::new(4.0, -5.0),
118/// ]);
119///
120/// assert_eq!(p, q.conj());
121/// ```
122///
123/// Biquaternion Product:
124/// The biquaternion product takes two biquaternions and outputs another
125/// biquaternion.
126/// ```
127/// use hoomd_manifold::Biquaternion;
128/// use num::complex::Complex;
129///
130/// let q = Biquaternion::from([
131/// Complex::new(2.0, 0.0),
132/// Complex::new(0.0, 1.0),
133/// Complex::new(1.0, 0.0),
134/// Complex::new(1.0, 0.0),
135/// ]);
136/// let p = Biquaternion::from([
137/// Complex::new(3.0, 0.0),
138/// Complex::new(2.0, 0.0),
139/// Complex::new(1.0, 0.0),
140/// Complex::new(0.0, 1.0),
141/// ]);
142/// let c = Biquaternion::from([
143/// Complex::new(1.0, 3.0),
144/// Complex::new(2.0, 0.0),
145/// Complex::new(5.0, -2.0),
146/// Complex::new(-7.0, -1.0),
147/// ]);
148/// assert_eq!(c, q.dot(&p));
149/// ```
150///
151/// Scalar Product:
152/// The scalar product takes two biquaternions and outputs a complex number
153/// according to
154/// ```math
155/// \frac{1}{2}(a\overline{b} + b\overline{a})
156/// ```
157/// ```
158/// use hoomd_manifold::Biquaternion;
159/// use num::complex::Complex;
160///
161/// let q = Biquaternion::from([
162/// Complex::new(2.0, 0.0),
163/// Complex::new(0.0, 1.0),
164/// Complex::new(1.0, 0.0),
165/// Complex::new(1.0, 0.0),
166/// ]);
167/// let p = Biquaternion::from([
168/// Complex::new(3.0, 0.0),
169/// Complex::new(2.0, 0.0),
170/// Complex::new(1.0, 0.0),
171/// Complex::new(0.0, 1.0),
172/// ]);
173/// assert_eq!(Complex::new(7.0, 3.0), q.scalar_product(&p));
174/// ```
175///
176/// Biquaternion Norm:
177/// The scalar product furnishes a "norm" for the biquaternion.
178/// ```
179/// use hoomd_manifold::Biquaternion;
180/// use num::complex::Complex;
181///
182/// let q = Biquaternion::from([
183/// Complex::new(3.0, 0.0),
184/// Complex::new(0.0, 1.0),
185/// Complex::new(4.0, 0.0),
186/// Complex::new(0.0, 2.0),
187/// ]);
188/// assert_eq!(Complex::new(20.0_f64, 0.0).sqrt(), q.norm());
189/// ```
190#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
191pub struct Biquaternion {
192 /// Components of the biquaternion, in the order `[i,j,k,1]`.
193 pub components: [Complex<f64>; 4],
194}
195
196impl Biquaternion {
197 /// Compute the Hamiltonian conjugate or biconjugate of a biquaternion.
198 ///
199 /// # Example
200 /// ```
201 /// use hoomd_manifold::Biquaternion;
202 /// use num::complex::Complex;
203 ///
204 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
205 /// let q = Biquaternion::from([
206 /// Complex::new(-1.0, 0.0),
207 /// Complex::new(0.0, 1.0),
208 /// Complex::new(1.0, 0.0),
209 /// Complex::new(1.0, 0.0),
210 /// ]);
211 /// let p = Biquaternion::from([
212 /// Complex::new(1.0, 0.0),
213 /// Complex::new(0.0, -1.0),
214 /// Complex::new(-1.0, 0.0),
215 /// Complex::new(1.0, 0.0),
216 /// ]);
217 /// assert_eq!(p, q.bar());
218 /// # Ok(())
219 /// # }
220 /// ```
221 #[inline]
222 #[must_use]
223 pub fn bar(&self) -> Self {
224 Biquaternion::from([
225 (self.components[0]).scale(-1.0),
226 (self.components[1]).scale(-1.0),
227 (self.components[2]).scale(-1.0),
228 (self.components[3]),
229 ])
230 }
231 /// Compute the complex conjugate of a biquaternion.
232 ///
233 /// # Example
234 /// ```
235 /// use hoomd_manifold::Biquaternion;
236 /// use num::complex::Complex;
237 ///
238 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
239 /// let q = Biquaternion::from([
240 /// Complex::new(1.0, 0.0),
241 /// Complex::new(0.0, 1.0),
242 /// Complex::new(1.0, 0.0),
243 /// Complex::new(1.0, 2.0),
244 /// ]);
245 /// let p = Biquaternion::from([
246 /// Complex::new(1.0, 0.0),
247 /// Complex::new(0.0, -1.0),
248 /// Complex::new(1.0, 0.0),
249 /// Complex::new(1.0, -2.0),
250 /// ]);
251 /// assert_eq!(p, q.conj());
252 /// # Ok(())
253 /// # }
254 /// ```
255 #[inline]
256 #[must_use]
257 pub fn conj(&self) -> Self {
258 Biquaternion::from([
259 (self.components[0]).conj(),
260 (self.components[1]).conj(),
261 (self.components[2]).conj(),
262 (self.components[3]).conj(),
263 ])
264 }
265 /// Compute the squared norm of a biquaternion.
266 ///
267 /// # Example
268 /// ```
269 /// use hoomd_manifold::Biquaternion;
270 /// use num::complex::Complex;
271 ///
272 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
273 /// let q = Biquaternion::from([
274 /// Complex::new(1.0, 0.0),
275 /// Complex::new(0.0, 1.0),
276 /// Complex::new(1.0, 0.0),
277 /// Complex::new(1.0, 0.0),
278 /// ]);
279 /// assert_eq!(2.0, q.norm_squared().re);
280 /// # Ok(())
281 /// # }
282 /// ```
283 #[inline]
284 #[must_use]
285 pub fn norm_squared(&self) -> Complex<f64> {
286 self.scalar_product(self)
287 }
288 /// the norm of a biquaternion
289 ///
290 /// # Example
291 /// ```
292 /// use hoomd_manifold::Biquaternion;
293 /// use num::complex::Complex;
294 ///
295 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
296 /// let q = Biquaternion::from([
297 /// Complex::new(3.0, 0.0),
298 /// Complex::new(0.0, 1.0),
299 /// Complex::new(4.0, 0.0),
300 /// Complex::new(1.0, 0.0),
301 /// ]);
302 /// assert_eq!(Complex::new(5.0, 0.0), q.norm());
303 /// # Ok(())
304 /// # }
305 /// ```
306 #[inline]
307 #[must_use]
308 pub fn norm(&self) -> Complex<f64> {
309 self.norm_squared().sqrt()
310 }
311 /// Compute the quaternion product of two biquaternions.
312 ///
313 /// # Example
314 /// ```
315 /// use hoomd_manifold::Biquaternion;
316 /// use num::complex::Complex;
317 ///
318 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
319 /// let q = Biquaternion::from([
320 /// Complex::new(2.0, 0.0),
321 /// Complex::new(0.0, 1.0),
322 /// Complex::new(1.0, 0.0),
323 /// Complex::new(1.0, 0.0),
324 /// ]);
325 /// let p = Biquaternion::from([
326 /// Complex::new(3.0, 0.0),
327 /// Complex::new(2.0, 0.0),
328 /// Complex::new(1.0, 0.0),
329 /// Complex::new(0.0, 1.0),
330 /// ]);
331 /// let c = Biquaternion::from([
332 /// Complex::new(1.0, 3.0),
333 /// Complex::new(2.0, 0.0),
334 /// Complex::new(5.0, -2.0),
335 /// Complex::new(-7.0, -1.0),
336 /// ]);
337 /// assert_eq!(c, q.dot(&p));
338 /// # Ok(())
339 /// # }
340 /// ```
341 #[inline]
342 #[must_use]
343 pub fn dot(&self, other: &Self) -> Self {
344 Biquaternion::from([
345 self.components[3] * other.components[0]
346 + other.components[3] * self.components[0]
347 + self.components[1] * other.components[2]
348 - other.components[1] * self.components[2],
349 self.components[3] * other.components[1]
350 + other.components[3] * self.components[1]
351 + self.components[2] * other.components[0]
352 - other.components[2] * self.components[0],
353 self.components[3] * other.components[2]
354 + other.components[3] * self.components[2]
355 + self.components[0] * other.components[1]
356 - other.components[0] * self.components[1],
357 self.components[3] * other.components[3]
358 - self.components[0] * other.components[0]
359 - self.components[1] * other.components[1]
360 - self.components[2] * other.components[2],
361 ])
362 }
363 /// Compute the scalar product of two biquaternions.
364 ///
365 /// # Example
366 /// ```
367 /// use hoomd_manifold::Biquaternion;
368 /// use num::complex::Complex;
369 ///
370 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
371 /// let q = Biquaternion::from([
372 /// Complex::new(2.0, 0.0),
373 /// Complex::new(0.0, 1.0),
374 /// Complex::new(1.0, 0.0),
375 /// Complex::new(1.0, 0.0),
376 /// ]);
377 /// let p = Biquaternion::from([
378 /// Complex::new(3.0, 0.0),
379 /// Complex::new(2.0, 0.0),
380 /// Complex::new(1.0, 0.0),
381 /// Complex::new(0.0, 1.0),
382 /// ]);
383 /// assert_eq!(Complex::new(7.0, 3.0), q.scalar_product(&p));
384 /// # Ok(())
385 /// # }
386 /// ```
387 #[inline]
388 #[must_use]
389 pub fn scalar_product(&self, other: &Self) -> Complex<f64> {
390 zip(self.components.iter(), other.components.iter())
391 .fold(Complex::new(0.0, 0.0), |product, x| product + x.0 * x.1)
392 }
393
394 /// Create a [`UnitBiquaternion`] by normalizing the given biquaternion.
395 ///
396 /// # Errors
397 ///
398 /// [`Error::InvalidBiquaternionMagnitude`] when `self` is the 0 quaternion.
399 #[inline]
400 pub fn to_unit(self) -> Result<UnitBiquaternion, Error> {
401 let mag = self.norm();
402 if mag == Complex::new(0.0, 0.0) {
403 return Err(Error::InvalidBiquaternionMagnitude);
404 }
405 Ok(UnitBiquaternion(self / mag))
406 }
407
408 /// Create a [`UnitBiquaternion`] by normalizing the given biquaternion.
409 #[inline]
410 #[must_use]
411 pub fn to_unit_unchecked(self) -> UnitBiquaternion {
412 let mag = self.norm();
413 UnitBiquaternion(self / mag)
414 }
415}
416
417impl Default for Biquaternion {
418 /// Create a biquaternion with all zeros.
419 #[inline]
420 fn default() -> Self {
421 Self {
422 components: array::from_fn(|_| Complex::new(0.0, 0.0)),
423 }
424 }
425}
426
427impl From<[Complex<f64>; 4]> for Biquaternion {
428 /// Construct a [`Biquaternion`] from 4 complex values.
429 ///
430 /// # Example
431 /// ```
432 /// use hoomd_manifold::Biquaternion;
433 /// use num::complex::Complex;
434 ///
435 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
436 /// let q = Biquaternion::from([
437 /// Complex::new(1.0, 0.0),
438 /// Complex::new(0.0, 0.1),
439 /// Complex::new(1.0, 0.0),
440 /// Complex::new(1.0, 1.0),
441 /// ]);
442 /// assert_eq!(
443 /// q.components,
444 /// [
445 /// Complex::new(1.0, 0.0),
446 /// Complex::new(0.0, 0.1),
447 /// Complex::new(1.0, 0.0),
448 /// Complex::new(1.0, 1.0)
449 /// ]
450 /// );
451 /// # Ok(())
452 /// # }
453 /// ```
454 #[inline]
455 fn from(value: [Complex<f64>; 4]) -> Self {
456 Self { components: value }
457 }
458}
459
460impl fmt::Display for Biquaternion {
461 #[inline]
462 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
463 write!(
464 f,
465 "[{}, {}, {}, {}]",
466 self.components[0], self.components[1], self.components[2], self.components[3]
467 )
468 }
469}
470
471impl Add for Biquaternion {
472 type Output = Self;
473
474 #[inline]
475 fn add(self, rhs: Self) -> Self {
476 Self {
477 components: array::from_fn(|i| self.components[i] + rhs.components[i]),
478 }
479 }
480}
481
482impl AddAssign for Biquaternion {
483 #[inline]
484 fn add_assign(&mut self, rhs: Self) {
485 for n in 0..4 {
486 self.components[n] += rhs.components[n];
487 }
488 }
489}
490
491impl Sub for Biquaternion {
492 type Output = Self;
493
494 #[inline]
495 fn sub(self, rhs: Self) -> Self {
496 Self {
497 components: array::from_fn(|i| self.components[i] - rhs.components[i]),
498 }
499 }
500}
501
502impl SubAssign for Biquaternion {
503 #[inline]
504 fn sub_assign(&mut self, rhs: Self) {
505 for n in 0..4 {
506 self.components[n] -= rhs.components[n];
507 }
508 }
509}
510
511impl Mul<f64> for Biquaternion {
512 type Output = Self;
513
514 #[inline]
515 fn mul(self, rhs: f64) -> Self {
516 Self {
517 components: array::from_fn(|i| (self.components[i]).scale(rhs)),
518 }
519 }
520}
521
522impl Mul<Complex<f64>> for Biquaternion {
523 type Output = Self;
524
525 #[inline]
526 fn mul(self, rhs: Complex<f64>) -> Self {
527 Self {
528 components: array::from_fn(|i| self.components[i] * rhs),
529 }
530 }
531}
532
533impl MulAssign<f64> for Biquaternion {
534 #[inline]
535 fn mul_assign(&mut self, rhs: f64) {
536 for n in 0..4 {
537 self.components[n] *= Complex::new(rhs, 0.0);
538 }
539 }
540}
541impl Div<f64> for Biquaternion {
542 type Output = Self;
543
544 #[inline]
545 fn div(self, rhs: f64) -> Self {
546 Self {
547 components: array::from_fn(|i| (self.components[i]).scale(1.0 / rhs)),
548 }
549 }
550}
551
552impl Div<Complex<f64>> for Biquaternion {
553 type Output = Self;
554
555 #[inline]
556 fn div(self, rhs: Complex<f64>) -> Self {
557 Self {
558 components: array::from_fn(|i| self.components[i] / rhs),
559 }
560 }
561}
562
563impl DivAssign<f64> for Biquaternion {
564 #[inline]
565 fn div_assign(&mut self, rhs: f64) {
566 for n in 0..4 {
567 self.components[n] /= Complex::new(rhs, 0.0);
568 }
569 }
570}
571
572#[derive(Clone, Copy, Debug, PartialEq)]
573/// Represent SO(3,1) with a normalized biquaternion.
574///
575/// Unit-norm Biquaternions furnish a representation of SO(3,1), analogous to
576/// quaternions and SO(3). If $`\vec{x} = (x_1, x_2, x_3, x_4)`$ is a vector
577/// in Minkowski space, then $`\vec{x}`$ can be mapped to a biquaternion
578/// ```math
579/// \vec{x} \mapsto X = [x_1, x_2, x_3,h x_4]
580/// ```
581/// (where h is the imaginary number) whose squared norm is
582/// ```math
583/// |X|^2 = x_1^2 + x_2^2 + x_3^2 - x_4^2
584/// ```
585/// It can be shown that, for
586/// a unit biquaternion $`q`$, the transformation
587/// ```math
588/// q^* X \overline{q} = X'
589/// ```
590/// preserves the norm, i.e.,
591/// ```math
592/// |X|^2 = |X'|^2
593/// ```
594/// We therefore have that this action by unit biquaternions
595/// produces a representation of SO(3,1). The biquaternion algebra can be used
596/// directly to transform Minkowski 4-vectors, or unit biquaternions can be
597/// represented as matrices using [`HyperbolicRotationMatrix<4>`].
598///
599/// Like quaternions, the unit biquaternion
600/// ```math
601/// q = \cos(\theta/2) + \bf{i}\sin(\theta/2)
602/// ```
603/// generates a rotation about the $` \mathbf{i} `$ axis by angle $`\theta`$:
604/// ```
605/// use approxim::assert_relative_eq;
606/// use hoomd_manifold::{
607/// Biquaternion, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
608/// UnitBiquaternion,
609/// };
610/// use num::complex::Complex;
611/// use std::f64::consts::PI;
612///
613/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
614/// let q = Biquaternion::from([
615/// Complex::new((PI / 4.0).sin(), 0.0),
616/// Complex::new(0.0, 0.0),
617/// Complex::new(0.0, 0.0),
618/// Complex::new((PI / 4.0).cos(), 0.0),
619/// ]);
620/// let v = q.to_unit()?;
621/// let x = Minkowski::from([1.0, 1.0, 1.0, 1.0]);
622/// let rotation = HyperbolicRotationMatrix::from(v);
623/// let rotated = rotation.hyperbolic_rotate(&x);
624/// assert_relative_eq!(rotated, [1.0, -1.0, 1.0, 1.0].into(), epsilon = 1e-12);
625/// # Ok(())
626/// # }
627/// ```
628///
629/// However, biquaternions also generate boosts via
630/// ```math
631/// q = \cosh(v) + \mathbf{i}h\sinh(v)
632/// ```
633/// which represents a boost of rapidity $`v`$ in the $`\mathbf{i}`$ direction:
634/// ```
635/// use approxim::assert_relative_eq;
636/// use hoomd_manifold::{
637/// Biquaternion, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
638/// UnitBiquaternion,
639/// };
640/// use num::complex::Complex;
641/// use std::f64::consts::PI;
642///
643/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
644/// let q = Biquaternion::from([
645/// Complex::new(0.0, (0.2_f64).sinh()),
646/// Complex::new(0.0, 0.0),
647/// Complex::new(0.0, 0.0),
648/// Complex::new((0.2_f64).cosh(), 0.0),
649/// ]);
650/// let v = q.to_unit()?;
651/// let x = Minkowski::from([0.0, 0.0, 0.0, 1.0]);
652/// let boost = HyperbolicRotationMatrix::from(v);
653/// let boosted = boost.hyperbolic_rotate(&x);
654/// assert_relative_eq!(
655/// boosted,
656/// [(0.4_f64).sinh(), 0.0, 0.0, (0.4_f64).cosh()].into(),
657/// epsilon = 1e-12
658/// );
659/// # Ok(())
660/// # }
661/// ```
662pub struct UnitBiquaternion(Biquaternion);
663
664impl UnitBiquaternion {
665 /// Normalize a biquaternion.
666 #[inline]
667 #[must_use]
668 pub fn normalized(self) -> Self {
669 let UnitBiquaternion(q) = self;
670 let f = 1.0 / q.norm();
671 Self(q * f)
672 }
673 /// Compute the square of the norm of a biquaternion.
674 #[inline]
675 #[must_use]
676 pub fn norm_squared(self) -> Complex<f64> {
677 let UnitBiquaternion(q) = self;
678 q.norm_squared()
679 }
680}
681
682impl Distribution<UnitBiquaternion> for StandardUniform {
683 /// Sample a random [`UnitBiquaternion`]
684 ///
685 /// # Example
686 ///
687 /// ```
688 /// use approxim::assert_relative_eq;
689 /// use hoomd_manifold::{Biquaternion, UnitBiquaternion};
690 /// use num::complex::Complex;
691 /// use rand::{RngExt, SeedableRng, rngs::StdRng};
692 ///
693 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
694 /// let mut rng = StdRng::seed_from_u64(1);
695 /// let v: UnitBiquaternion = rng.random();
696 /// assert_relative_eq!(v.norm_squared().re, 1.0, epsilon = 1e-12);
697 /// # Ok(())
698 /// # }
699 /// ```
700 #[inline]
701 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> UnitBiquaternion {
702 let uniform = Uniform::new(-1.0, 1.0).expect("hard-coded distribution should be valid");
703
704 let array_re: [f64; 4] = array::from_fn(|_| uniform.sample(rng));
705 let array_im: [f64; 4] = array::from_fn(|_| uniform.sample(rng));
706 let mut scale =
707 zip(array_re.iter(), array_im.iter()).fold(Complex::new(0.0, 0.0), |product, x| {
708 product + Complex::new((x.0).powi(2) - (x.1).powi(2), 2.0_f64 * (x.1) * (x.0))
709 });
710 scale = scale.sqrt();
711 UnitBiquaternion(Biquaternion {
712 components: array::from_fn(|i| Complex::new(array_re[i], array_im[i]) / scale),
713 })
714 }
715}
716
717impl From<UnitBiquaternion> for HyperbolicRotationMatrix<4> {
718 #[inline]
719 #[expect(clippy::many_single_char_names, reason = "dummy variables")]
720 fn from(q: UnitBiquaternion) -> HyperbolicRotationMatrix<4> {
721 let UnitBiquaternion(biquaternion) = q;
722 let [a, b, c, d]: [Complex<f64>; 4] = array::from_fn(|i| biquaternion.components[i]);
723
724 HyperbolicRotationMatrix {
725 rows: [
726 [
727 (d * d.conj() + a * a.conj() - b * b.conj() - c * c.conj()).re,
728 (a * b.conj() + b * a.conj() - c * d.conj() - d * c.conj()).re,
729 (a * c.conj() + c * a.conj() + b * d.conj() + d * b.conj()).re,
730 -(d * a.conj() - a * d.conj() + b * c.conj() - c * b.conj()).im,
731 ]
732 .into(),
733 [
734 (b * a.conj() + a * b.conj() + c * d.conj() + d * c.conj()).re,
735 (d * d.conj() - a * a.conj() + b * b.conj() - c * c.conj()).re,
736 (b * c.conj() + c * b.conj() - a * d.conj() - d * a.conj()).re,
737 -(d * b.conj() - b * d.conj() + c * a.conj() - a * c.conj()).im,
738 ]
739 .into(),
740 [
741 (c * a.conj() + a * c.conj() - b * d.conj() - d * b.conj()).re,
742 (c * b.conj() + b * c.conj() + a * d.conj() + d * a.conj()).re,
743 (d * d.conj() - a * a.conj() - b * b.conj() + c * c.conj()).re,
744 -(d * c.conj() - c * d.conj() + a * b.conj() - b * a.conj()).im,
745 ]
746 .into(),
747 [
748 (a * d.conj() - d * a.conj() + b * c.conj() - c * b.conj()).im,
749 (b * d.conj() - d * b.conj() + c * a.conj() - a * c.conj()).im,
750 (c * d.conj() - d * c.conj() + a * b.conj() - b * a.conj()).im,
751 (a * a.conj() + b * b.conj() + c * c.conj() + d * d.conj()).re,
752 ]
753 .into(),
754 ],
755 }
756 }
757}
758
759impl HyperbolicRotate<Minkowski<4>> for UnitBiquaternion {
760 type Matrix = HyperbolicRotationMatrix<4>;
761
762 /// Transform a [`Minkowski<4>`] by a [`UnitBiquaternion`].
763 ///
764 /// ```math
765 /// \overline{\mathbf{q}} \vec{a} \mathbf{q}^*
766 /// ```
767 ///
768 /// # Examples
769 /// Rotation about z axis:
770 /// ```
771 /// use approxim::assert_relative_eq;
772 /// use hoomd_manifold::{
773 /// Biquaternion, HyperbolicRotate, Minkowski, UnitBiquaternion,
774 /// };
775 /// use num::complex::Complex;
776 /// use std::f64::consts::PI;
777 ///
778 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
779 /// let x = Minkowski::from([1.0, 0.0, 0.0, 1.0]);
780 /// let q = Biquaternion::from([
781 /// Complex::new(0.0, 0.0),
782 /// Complex::new(0.0, 0.0),
783 /// Complex::new((PI / 4.0).sin(), 0.0),
784 /// Complex::new((PI / 4.0).cos(), 0.0),
785 /// ]);
786 /// let v = q.to_unit_unchecked();
787 /// let rotated = v.hyperbolic_rotate(&x);
788 /// assert_relative_eq!(rotated, [0.0, 1.0, 0.0, 1.0].into(), epsilon = 1e-12);
789 /// # Ok(())
790 /// # }
791 /// ```
792 ///
793 /// Boost in x direction:
794 /// ```
795 /// use approxim::assert_relative_eq;
796 /// use hoomd_manifold::{
797 /// Biquaternion, HyperbolicRotate, Minkowski, UnitBiquaternion,
798 /// };
799 /// use num::complex::Complex;
800 /// use std::f64::consts::PI;
801 ///
802 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
803 /// let x = Minkowski::from([0.0, 0.0, 0.0, 1.0]);
804 /// let q = Biquaternion::from([
805 /// Complex::new(0.0, PI / 4.0).sin(),
806 /// Complex::new(0.0, 0.0),
807 /// Complex::new(0.0, 0.0),
808 /// Complex::new(0.0, PI / 4.0).cos(),
809 /// ]);
810 /// let v = q.to_unit_unchecked();
811 /// let boosted = v.hyperbolic_rotate(&x);
812 /// assert_relative_eq!(
813 /// boosted,
814 /// [(PI / 2.0).sinh(), 0.0, 0.0, (PI / 2.0).cosh()].into(),
815 /// epsilon = 1e-12
816 /// );
817 /// # Ok(())
818 /// # }
819 /// ```
820 #[inline]
821 fn hyperbolic_rotate(&self, vector: &Minkowski<4>) -> Minkowski<4> {
822 let UnitBiquaternion(biquaternion) = self;
823 let x = Biquaternion::from([
824 Complex::new(vector[0], 0.0),
825 Complex::new(vector[1], 0.0),
826 Complex::new(vector[2], 0.0),
827 Complex::new(0.0, vector[3]),
828 ]);
829 let x_transformed = (biquaternion.conj()).dot(&x.dot(&(biquaternion.bar())));
830 Minkowski::from([
831 x_transformed.components[0].re,
832 x_transformed.components[1].re,
833 x_transformed.components[2].re,
834 x_transformed.components[3].im,
835 ])
836 }
837}
838
839#[cfg(test)]
840mod tests {
841 use super::*;
842 use approxim::assert_relative_eq;
843 use num::complex::Complex;
844 use rstest::*;
845 use std::f64::consts::PI;
846
847 #[test]
848 fn from_array() {
849 let q = Biquaternion::from([
850 Complex::new(-1.0, 0.0),
851 Complex::new(0.0, 1.0),
852 Complex::new(1.0, 0.0),
853 Complex::new(1.0, 0.0),
854 ]);
855 assert!(q.components[0] == Complex::new(-1.0, 0.0));
856 assert!(q.components[1] == Complex::new(0.0, 1.0));
857 assert!(q.components[2] == Complex::new(1.0, 0.0));
858 assert!(q.components[3] == Complex::new(1.0, 0.0));
859 }
860
861 #[test]
862 fn bar() {
863 let q = Biquaternion::from([
864 Complex::new(-2.0, 0.0),
865 Complex::new(-1.0, 1.0),
866 Complex::new(1.0, 0.0),
867 Complex::new(1.0, 0.0),
868 ]);
869 let p = Biquaternion::from([
870 Complex::new(2.0, 0.0),
871 Complex::new(1.0, -1.0),
872 Complex::new(-1.0, 0.0),
873 Complex::new(1.0, 0.0),
874 ]);
875 assert_eq!(p, q.bar());
876 }
877
878 #[test]
879 fn conjugate() {
880 let q = Biquaternion::from([
881 Complex::new(1.0, 8.0),
882 Complex::new(2.0, 7.0),
883 Complex::new(3.0, 6.0),
884 Complex::new(4.0, 5.0),
885 ]);
886 let p = Biquaternion::from([
887 Complex::new(1.0, -8.0),
888 Complex::new(2.0, -7.0),
889 Complex::new(3.0, -6.0),
890 Complex::new(4.0, -5.0),
891 ]);
892 assert_eq!(p, q.conj());
893 }
894
895 #[test]
896 fn biquat_product() {
897 let q = Biquaternion::from([
898 Complex::new(2.0, 0.0),
899 Complex::new(0.0, 1.0),
900 Complex::new(1.0, 0.0),
901 Complex::new(1.0, 0.0),
902 ]);
903 let p = Biquaternion::from([
904 Complex::new(3.0, 0.0),
905 Complex::new(2.0, 0.0),
906 Complex::new(1.0, 0.0),
907 Complex::new(0.0, 1.0),
908 ]);
909 let c = Biquaternion::from([
910 Complex::new(1.0, 3.0),
911 Complex::new(2.0, 0.0),
912 Complex::new(5.0, -2.0),
913 Complex::new(-7.0, -1.0),
914 ]);
915 assert_eq!(c, q.dot(&p));
916 }
917
918 #[test]
919 fn scalar_product() {
920 let q = Biquaternion::from([
921 Complex::new(2.0, 0.0),
922 Complex::new(0.0, 1.0),
923 Complex::new(1.0, 0.0),
924 Complex::new(1.0, 0.0),
925 ]);
926 let p = Biquaternion::from([
927 Complex::new(3.0, 0.0),
928 Complex::new(2.0, 0.0),
929 Complex::new(1.0, 0.0),
930 Complex::new(0.0, 1.0),
931 ]);
932 assert_eq!(Complex::new(7.0, 3.0), q.scalar_product(&p));
933 }
934
935 #[test]
936 fn norm() {
937 let q = Biquaternion::from([
938 Complex::new(3.0, 0.0),
939 Complex::new(0.0, 1.0),
940 Complex::new(4.0, 0.0),
941 Complex::new(0.0, 2.0),
942 ]);
943 assert_eq!(Complex::new(20.0_f64, 0.0).sqrt(), q.norm());
944 }
945
946 #[test]
947 fn ops() {
948 let a = Biquaternion::from([
949 Complex::new(1.0, 0.0),
950 Complex::new(2.0, 1.0),
951 Complex::new(-3.0, 0.0),
952 Complex::new(4.0, 2.0),
953 ]);
954 let b = Biquaternion::from([
955 Complex::new(4.0, 0.0),
956 Complex::new(0.0, 3.0),
957 Complex::new(-2.0, 0.0),
958 Complex::new(1.0, 0.0),
959 ]);
960
961 // +, +=
962 assert_eq!(
963 a + b,
964 Biquaternion::from([
965 Complex::new(5.0, 0.0),
966 Complex::new(2.0, 4.0),
967 Complex::new(-5.0, 0.0),
968 Complex::new(5.0, 2.0)
969 ])
970 );
971 let mut c = a;
972 c += b;
973 assert_eq!(
974 c,
975 Biquaternion::from([
976 Complex::new(5.0, 0.0),
977 Complex::new(2.0, 4.0),
978 Complex::new(-5.0, 0.0),
979 Complex::new(5.0, 2.0)
980 ])
981 );
982
983 // -, -=
984 assert_eq!(
985 a - b,
986 Biquaternion::from([
987 Complex::new(-3.0, 0.0),
988 Complex::new(2.0, -2.0),
989 Complex::new(-1.0, 0.0),
990 Complex::new(3.0, 2.0)
991 ])
992 );
993 let mut c = a;
994 c -= b;
995 assert_eq!(
996 c,
997 Biquaternion::from([
998 Complex::new(-3.0, 0.0),
999 Complex::new(2.0, -2.0),
1000 Complex::new(-1.0, 0.0),
1001 Complex::new(3.0, 2.0)
1002 ])
1003 );
1004
1005 // Scalar * and /
1006 assert_eq!(
1007 a * 2.0,
1008 Biquaternion::from([
1009 Complex::new(2.0, 0.0),
1010 Complex::new(4.0, 2.0),
1011 Complex::new(-6.0, 0.0),
1012 Complex::new(8.0, 4.0)
1013 ])
1014 );
1015 let mut c = a;
1016 c *= 2.0;
1017 assert_eq!(
1018 c,
1019 Biquaternion::from([
1020 Complex::new(2.0, 0.0),
1021 Complex::new(4.0, 2.0),
1022 Complex::new(-6.0, 0.0),
1023 Complex::new(8.0, 4.0)
1024 ])
1025 );
1026
1027 assert_eq!(
1028 a / 2.0,
1029 Biquaternion::from([
1030 Complex::new(0.5, 0.0),
1031 Complex::new(1.0, 0.5),
1032 Complex::new(-1.5, 0.0),
1033 Complex::new(2.0, 1.0)
1034 ])
1035 );
1036 let mut c = a;
1037 c /= 2.0;
1038 assert_eq!(
1039 c,
1040 Biquaternion::from([
1041 Complex::new(0.5, 0.0),
1042 Complex::new(1.0, 0.5),
1043 Complex::new(-1.5, 0.0),
1044 Complex::new(2.0, 1.0)
1045 ])
1046 );
1047 }
1048
1049 #[test]
1050 fn display() {
1051 let a = Biquaternion::from([
1052 Complex::new(0.5, -3.1),
1053 Complex::new(1.2, 5.2),
1054 Complex::new(-1.5, 1.5),
1055 Complex::new(2.1, 1.5),
1056 ]);
1057 let s = format!("{a}");
1058 assert_eq!(s, "[0.5-3.1i, 1.2+5.2i, -1.5+1.5i, 2.1+1.5i]");
1059 }
1060
1061 #[test]
1062 fn to_unit() {
1063 let a = Biquaternion::from([
1064 Complex::new(0.5, -3.0),
1065 Complex::new(1.0, 5.0),
1066 Complex::new(-1.5, 1.5),
1067 Complex::new(2.0, 1.0),
1068 ]);
1069 let a_unit = a.to_unit().expect("hard-coded to be nonzero");
1070 assert_relative_eq!(a_unit.norm_squared().re, 1.0, epsilon = 1e-12);
1071 }
1072
1073 #[test]
1074 fn to_unit_unchecked() {
1075 let a = Biquaternion::from([
1076 Complex::new(2.8, 0.5),
1077 Complex::new(-1.0, 2.0),
1078 Complex::new(0.5, 1.3),
1079 Complex::new(2.3, 3.4),
1080 ]);
1081 let a_unit = a.to_unit_unchecked();
1082 assert_relative_eq!(a_unit.norm_squared().re, 1.0, epsilon = 1e-12);
1083 }
1084
1085 // Test named cases of three input values (rotation biquaternion, Minkowski input, and answer)
1086 #[rstest]
1087 #[case::y_rotate_pi_halves([Complex::new(0.0,0.0),
1088 Complex::new((PI/4.0).sin(),0.0),
1089 Complex::new(0.0, 0.0),
1090 Complex::new((PI/4.0).cos(), 0.0)],
1091 [1.0,0.0,0.0,1.0],
1092 [0.0,0.0,-1.0,1.0])]
1093 #[case::x_boost_half([Complex::new(0.0, 0.25).sin(),
1094 Complex::new(0.0,0.0),
1095 Complex::new(0.0, 0.0),
1096 Complex::new(0.0, 0.25).cos()],
1097 [0.0,0.0,0.0,1.0],
1098 [(0.5_f64).sinh(),0.0,0.0,(0.5_f64).cosh()])]
1099 #[case::z_boost_one([Complex::new(0.0, 0.0),
1100 Complex::new(0.0,0.0),
1101 Complex::new(0.0, 0.5).sin(),
1102 Complex::new(0.0, 0.5).cos()],
1103 [0.0,0.0,0.0,1.0],
1104 [0.0,0.0,(1.0_f64).sinh(),(1.0_f64).cosh()])]
1105 #[case::z_rotate_pi([Complex::new(0.0, 0.0),
1106 Complex::new(0.0,0.0),
1107 Complex::new((PI/2.0).sin(),0.0),
1108 Complex::new((PI/2.0).cos(),0.0)],
1109 [1.0,0.0,0.0,1.0],
1110 [-1.0,0.0,0.0,1.0])]
1111 fn rotate(#[case] biquat: [Complex<f64>; 4], #[case] vec: [f64; 4], #[case] ans: [f64; 4]) {
1112 let q = Biquaternion::from(biquat);
1113 let v = q.to_unit().expect("hard-coded to be nonzero");
1114 let x = Minkowski::from(vec);
1115
1116 // using matrix representation
1117 let matrix = HyperbolicRotationMatrix::from(v);
1118 let from_matrix = matrix.hyperbolic_rotate(&x);
1119 assert_relative_eq!(from_matrix, ans.into(), epsilon = 1e-12);
1120
1121 // using biquaternion algebra
1122 let from_algebra = v.hyperbolic_rotate(&x);
1123 assert_relative_eq!(from_algebra, ans.into(), epsilon = 1e-12);
1124 }
1125}