Quaternion.php 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337
  1. <?php
  2. namespace MathPHP\Number;
  3. use MathPHP\Exception;
  4. /**
  5. * Quaternionic Numbers
  6. *
  7. * A quaternion is a number that can be expressed in the form a + bi + cj + dk,
  8. * where a, b, c, andd are real numbers and i, j, and k are the basic quaternions, satisfying the
  9. * equation i² = j² = k² = ijk = −1.
  10. * https://en.wikipedia.org/wiki/Quaternion
  11. */
  12. class Quaternion implements ObjectArithmetic
  13. {
  14. /** @var int|float Real part of the quaternionic number */
  15. protected $r;
  16. /** @var int|float First Imaginary part of the quaternionic number */
  17. protected $i;
  18. /** @var int|float Second Imaginary part of the quaternionic number */
  19. protected $j;
  20. /** @var int|float Third Imaginary part of the quaternionic number */
  21. protected $k;
  22. /** Floating-point range near zero to consider insignificant */
  23. const EPSILON = 1e-6;
  24. /**
  25. * @param int|float $r Real part
  26. * @param int|float $i Imaginary part
  27. * @param int|float $j Imaginary part
  28. * @param int|float $k Imaginary part
  29. */
  30. public function __construct($r, $i, $j, $k)
  31. {
  32. if (!\is_numeric($r) || !\is_numeric($i) || !\is_numeric($j) || !\is_numeric($k)) {
  33. throw new Exception\BadDataException('Values must be real numbers.');
  34. }
  35. $this->r = $r;
  36. $this->i = $i;
  37. $this->j = $j;
  38. $this->k = $k;
  39. }
  40. /**
  41. * Creates 0 + 0i
  42. *
  43. * @return Quaternion
  44. */
  45. public static function createZeroValue(): ObjectArithmetic
  46. {
  47. return new Quaternion(0, 0, 0, 0);
  48. }
  49. /**
  50. * String representation of a complex number
  51. * a + bi + cj + dk, a - bi - cj - dk, etc.
  52. *
  53. * @return string
  54. */
  55. public function __toString(): string
  56. {
  57. if ($this->r == 0 & $this->i == 0 & $this->j == 0 & $this->k == 0) {
  58. return '0';
  59. }
  60. $string = self::stringifyNumberPart($this->r);
  61. $string = self::stringifyNumberPart($this->i, 'i', $string);
  62. $string = self::stringifyNumberPart($this->j, 'j', $string);
  63. return self::stringifyNumberPart($this->k, 'k', $string);
  64. }
  65. /**
  66. * Get r or i or j or k
  67. *
  68. * @param string $part
  69. *
  70. * @return int|float
  71. *
  72. * @throws Exception\BadParameterException if something other than r or i is attempted
  73. */
  74. public function __get(string $part)
  75. {
  76. switch ($part) {
  77. case 'r':
  78. case 'i':
  79. case 'j':
  80. case 'k':
  81. return $this->$part;
  82. default:
  83. throw new Exception\BadParameterException("The $part property does not exist in Quaternion");
  84. }
  85. }
  86. /**************************************************************************
  87. * UNARY FUNCTIONS
  88. **************************************************************************/
  89. /**
  90. * The conjugate of a quaternion
  91. * https://en.wikipedia.org/wiki/Quaternion#Conjugation.2C_the_norm.2C_and_reciprocal
  92. *
  93. * q* = a - bi - cj -dk
  94. *
  95. * @return Quaternion
  96. */
  97. public function complexConjugate(): Quaternion
  98. {
  99. return new Quaternion($this->r, -$this->i, -$this->j, -$this->k);
  100. }
  101. /**
  102. * The absolute value (magnitude) of a quaternion or norm
  103. * https://en.wikipedia.org/wiki/Quaternion#Conjugation.2C_the_norm.2C_and_reciprocal
  104. *
  105. * If z = a + bi + cj + dk
  106. * _________________
  107. * |z| = √a² + b² + c² + d²
  108. *
  109. * @return int|float
  110. */
  111. public function abs()
  112. {
  113. return sqrt($this->r**2 + $this->i**2 + $this->j**2 + $this->k**2);
  114. }
  115. /**
  116. * The inverse of a quaternion (reciprocal)
  117. *
  118. * https://en.wikipedia.org/wiki/Quaternion#Conjugation.2C_the_norm.2C_and_reciprocal
  119. *
  120. * 1
  121. * (a + bi + cj + dk)⁻¹ = ----------------- (a - bi - cj -dk)
  122. * a² + b² + c² + d²
  123. *
  124. * @return Quaternion
  125. *
  126. * @throws Exception\BadDataException if = to 0 + 0i
  127. */
  128. public function inverse(): Quaternion
  129. {
  130. if ($this->r == 0 && $this->i == 0 && $this->j == 0 && $this->k == 0) {
  131. throw new Exception\BadDataException('Cannot take inverse of 0 + 0i');
  132. }
  133. return $this->complexConjugate()->divide($this->abs() ** 2);
  134. }
  135. /**
  136. * Negate the quaternion
  137. * Switches the signs of both the real and imaginary parts.
  138. *
  139. * @return Quaternion
  140. */
  141. public function negate(): Quaternion
  142. {
  143. return new Quaternion(-$this->r, -$this->i, -$this->j, -$this->k);
  144. }
  145. /**************************************************************************
  146. * BINARY FUNCTIONS
  147. **************************************************************************/
  148. /**
  149. * Quaternion addition
  150. *
  151. *
  152. * (a + bi + cj + dk) - (e + fi + gj + hk) = (a + e) + (b + f)i + (c + g)j + (d + h)k
  153. *
  154. * @param int|float|Quaternion $q
  155. *
  156. * @return Quaternion
  157. *
  158. * @throws Exception\IncorrectTypeException if the argument is not numeric or Complex.
  159. */
  160. public function add($q): Quaternion
  161. {
  162. if (!\is_numeric($q) && ! $q instanceof Quaternion) {
  163. throw new Exception\IncorrectTypeException('Argument must be real or quaternion' . \print_r($q, true));
  164. }
  165. if (\is_numeric($q)) {
  166. $r = $this->r + $q;
  167. return new Quaternion($r, $this->i, $this->j, $this->k);
  168. }
  169. $r = $this->r + $q->r;
  170. $i = $this->i + $q->i;
  171. $j = $this->j + $q->j;
  172. $k = $this->k + $q->k;
  173. return new Quaternion($r, $i, $j, $k);
  174. }
  175. /**
  176. * Quaternion subtraction
  177. *
  178. *
  179. * (a + bi + cj + dk) - (e + fi + gj + hk) = (a - e) + (b - f)i + (c - g)j + (d - h)k
  180. *
  181. * @param int|float|Quaternion $q
  182. *
  183. * @return Quaternion
  184. *
  185. * @throws Exception\IncorrectTypeException if the argument is not numeric or Complex.
  186. */
  187. public function subtract($q): Quaternion
  188. {
  189. if (!is_numeric($q) && ! $q instanceof Quaternion) {
  190. throw new Exception\IncorrectTypeException('Argument must be real or quaternion' . print_r($q, true));
  191. }
  192. if (\is_numeric($q)) {
  193. $r = $this->r - $q;
  194. return new Quaternion($r, $this->i, $this->j, $this->k);
  195. }
  196. $r = $this->r - $q->r;
  197. $i = $this->i - $q->i;
  198. $j = $this->j - $q->j;
  199. $k = $this->k - $q->k;
  200. return new Quaternion($r, $i, $j, $k);
  201. }
  202. /**
  203. * Quaternion multiplication (Hamilton product)
  204. *
  205. * (a₁ + b₁i - c₁j - d₁k)(a₂ + b₂i + c₂j + d₂k)
  206. *
  207. * a₁a₂ - b₁b₂ - c₁c₂ - d₁d₂
  208. * + (b₁a₂ + a₁b₂ + c₁d₂ - d₁c₂)i
  209. * + (a₁c₂ - b₁d₂ + c₁a₂ + d₁b₂)k
  210. * + (a₁d₂ + b₁c₂ - c₁b₂ + d₁a₂)k
  211. *
  212. * Note: Quaternion multiplication is not commutative.
  213. *
  214. * @param int|float|Quaternion $q
  215. *
  216. * @return Quaternion
  217. *
  218. * @throws Exception\IncorrectTypeException if the argument is not numeric or Complex.
  219. */
  220. public function multiply($q): Quaternion
  221. {
  222. if (!is_numeric($q) && ! $q instanceof Quaternion) {
  223. throw new Exception\IncorrectTypeException('Argument must be real or quaternion' . print_r($q, true));
  224. }
  225. if (\is_numeric($q)) {
  226. return new Quaternion($this->r * $q, $this->i * $q, $this->j * $q, $this->k * $q);
  227. }
  228. [$a₁, $b₁, $c₁, $d₁] = [$this->r, $this->i, $this->j, $this->k];
  229. [$a₂, $b₂, $c₂, $d₂] = [$q->r, $q->i, $q->j, $q->k];
  230. return new Quaternion(
  231. $a₁*$a₂ - $b₁*$b₂ - $c₁*$c₂ - $d₁*$d₂,
  232. $b₁*$a₂ + $a₁*$b₂ + $c₁*$d₂ - $d₁*$c₂,
  233. $a₁*$c₂ - $b₁*$d₂ + $c₁*$a₂ + $d₁*$b₂,
  234. $a₁*$d₂ + $b₁*$c₂ - $c₁*$b₂ + $d₁*$a₂
  235. );
  236. }
  237. /**
  238. * Quaternion division
  239. * Dividing two quaternions is accomplished by multiplying the first by the inverse of the second
  240. * This is not commutative!
  241. *
  242. * @param int|float|Quaternion $q
  243. *
  244. * @return Quaternion
  245. *
  246. * @throws Exception\IncorrectTypeException if the argument is not numeric or Complex.
  247. */
  248. public function divide($q): Quaternion
  249. {
  250. if (!is_numeric($q) && ! $q instanceof Quaternion) {
  251. throw new Exception\IncorrectTypeException('Argument must be real or quaternion' . print_r($q, true));
  252. }
  253. if (\is_numeric($q)) {
  254. $r = $this->r / $q;
  255. $i = $this->i / $q;
  256. $j = $this->j / $q;
  257. $k = $this->k / $q;
  258. return new Quaternion($r, $i, $j, $k);
  259. }
  260. return $this->multiply($q->inverse());
  261. }
  262. /**************************************************************************
  263. * COMPARISON FUNCTIONS
  264. **************************************************************************/
  265. /**
  266. * Test for equality
  267. * Two quaternions are equal if and only if both their real and imaginary parts are equal.
  268. *
  269. *
  270. * @param Quaternion $q
  271. *
  272. * @return bool
  273. */
  274. public function equals(Quaternion $q): bool
  275. {
  276. return \abs($this->r - $q->r) < self::EPSILON && \abs($this->i - $q->i) < self::EPSILON
  277. && \abs($this->j - $q->j) < self::EPSILON && \abs($this->k - $q->k) < self::EPSILON;
  278. }
  279. /**************************************************************************
  280. * PRIVATE FUNCTIONS
  281. **************************************************************************/
  282. /**
  283. * Stringify an additional part of the quaternion
  284. *
  285. * @param int|float $q
  286. * @param string $unit
  287. * @param string $string
  288. * @return string
  289. */
  290. private static function stringifyNumberPart($q, string $unit = '', string $string = ''): string
  291. {
  292. if ($q == 0) {
  293. return $string;
  294. }
  295. if ($q > 0) {
  296. $plus = $string == '' ? '' : ' + ';
  297. return $string . $plus . "$q" . $unit;
  298. }
  299. $minus = $string == '' ? '-' : ' - ';
  300. return $string . $minus . (string) \abs($q) . $unit;
  301. }
  302. }