complex.rs 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448
  1. // Copyright 2013 The Rust Project Developers. See the COPYRIGHT
  2. // file at the top-level directory of this distribution and at
  3. // http://rust-lang.org/COPYRIGHT.
  4. //
  5. // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
  6. // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
  7. // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
  8. // option. This file may not be copied, modified, or distributed
  9. // except according to those terms.
  10. //! Complex numbers.
  11. use std::fmt;
  12. use std::num::Float;
  13. use std::ops::{Add, Div, Mul, Neg, Sub};
  14. use {Zero, One, Num};
  15. // FIXME #1284: handle complex NaN & infinity etc. This
  16. // probably doesn't map to C's _Complex correctly.
  17. /// A complex number in Cartesian form.
  18. #[derive(PartialEq, Copy, Clone, Hash, RustcEncodable, RustcDecodable, Debug)]
  19. pub struct Complex<T> {
  20. /// Real portion of the complex number
  21. pub re: T,
  22. /// Imaginary portion of the complex number
  23. pub im: T
  24. }
  25. pub type Complex32 = Complex<f32>;
  26. pub type Complex64 = Complex<f64>;
  27. impl<T: Clone + Num> Complex<T> {
  28. /// Create a new Complex
  29. #[inline]
  30. pub fn new(re: T, im: T) -> Complex<T> {
  31. Complex { re: re, im: im }
  32. }
  33. /// Returns the square of the norm (since `T` doesn't necessarily
  34. /// have a sqrt function), i.e. `re^2 + im^2`.
  35. #[inline]
  36. pub fn norm_sqr(&self) -> T {
  37. self.re.clone() * self.re.clone() + self.im.clone() * self.im.clone()
  38. }
  39. /// Returns the complex conjugate. i.e. `re - i im`
  40. #[inline]
  41. pub fn conj(&self) -> Complex<T> {
  42. Complex::new(self.re.clone(), -self.im.clone())
  43. }
  44. /// Multiplies `self` by the scalar `t`.
  45. #[inline]
  46. pub fn scale(&self, t: T) -> Complex<T> {
  47. Complex::new(self.re.clone() * t.clone(), self.im.clone() * t)
  48. }
  49. /// Divides `self` by the scalar `t`.
  50. #[inline]
  51. pub fn unscale(&self, t: T) -> Complex<T> {
  52. Complex::new(self.re.clone() / t.clone(), self.im.clone() / t)
  53. }
  54. /// Returns `1/self`
  55. #[inline]
  56. pub fn inv(&self) -> Complex<T> {
  57. let norm_sqr = self.norm_sqr();
  58. Complex::new(self.re.clone() / norm_sqr.clone(),
  59. -self.im.clone() / norm_sqr)
  60. }
  61. }
  62. impl<T: Clone + Float> Complex<T> {
  63. /// Calculate |self|
  64. #[inline]
  65. pub fn norm(&self) -> T {
  66. self.re.clone().hypot(self.im.clone())
  67. }
  68. }
  69. impl<T: Clone + Float + Num> Complex<T> {
  70. /// Calculate the principal Arg of self.
  71. #[inline]
  72. pub fn arg(&self) -> T {
  73. self.im.clone().atan2(self.re.clone())
  74. }
  75. /// Convert to polar form (r, theta), such that `self = r * exp(i
  76. /// * theta)`
  77. #[inline]
  78. pub fn to_polar(&self) -> (T, T) {
  79. (self.norm(), self.arg())
  80. }
  81. /// Convert a polar representation into a complex number.
  82. #[inline]
  83. pub fn from_polar(r: &T, theta: &T) -> Complex<T> {
  84. Complex::new(*r * theta.cos(), *r * theta.sin())
  85. }
  86. }
  87. macro_rules! forward_val_val_binop {
  88. (impl $imp:ident, $method:ident) => {
  89. impl<T: Clone + Num> $imp<Complex<T>> for Complex<T> {
  90. type Output = Complex<T>;
  91. #[inline]
  92. fn $method(self, other: Complex<T>) -> Complex<T> {
  93. (&self).$method(&other)
  94. }
  95. }
  96. }
  97. }
  98. macro_rules! forward_ref_val_binop {
  99. (impl $imp:ident, $method:ident) => {
  100. impl<'a, T: Clone + Num> $imp<Complex<T>> for &'a Complex<T> {
  101. type Output = Complex<T>;
  102. #[inline]
  103. fn $method(self, other: Complex<T>) -> Complex<T> {
  104. self.$method(&other)
  105. }
  106. }
  107. }
  108. }
  109. macro_rules! forward_val_ref_binop {
  110. (impl $imp:ident, $method:ident) => {
  111. impl<'a, T: Clone + Num> $imp<&'a Complex<T>> for Complex<T> {
  112. type Output = Complex<T>;
  113. #[inline]
  114. fn $method(self, other: &Complex<T>) -> Complex<T> {
  115. (&self).$method(other)
  116. }
  117. }
  118. }
  119. }
  120. macro_rules! forward_all_binop {
  121. (impl $imp:ident, $method:ident) => {
  122. forward_val_val_binop!(impl $imp, $method);
  123. forward_ref_val_binop!(impl $imp, $method);
  124. forward_val_ref_binop!(impl $imp, $method);
  125. };
  126. }
  127. /* arithmetic */
  128. forward_all_binop!(impl Add, add);
  129. // (a + i b) + (c + i d) == (a + c) + i (b + d)
  130. impl<'a, 'b, T: Clone + Num> Add<&'b Complex<T>> for &'a Complex<T> {
  131. type Output = Complex<T>;
  132. #[inline]
  133. fn add(self, other: &Complex<T>) -> Complex<T> {
  134. Complex::new(self.re.clone() + other.re.clone(),
  135. self.im.clone() + other.im.clone())
  136. }
  137. }
  138. forward_all_binop!(impl Sub, sub);
  139. // (a + i b) - (c + i d) == (a - c) + i (b - d)
  140. impl<'a, 'b, T: Clone + Num> Sub<&'b Complex<T>> for &'a Complex<T> {
  141. type Output = Complex<T>;
  142. #[inline]
  143. fn sub(self, other: &Complex<T>) -> Complex<T> {
  144. Complex::new(self.re.clone() - other.re.clone(),
  145. self.im.clone() - other.im.clone())
  146. }
  147. }
  148. forward_all_binop!(impl Mul, mul);
  149. // (a + i b) * (c + i d) == (a*c - b*d) + i (a*d + b*c)
  150. impl<'a, 'b, T: Clone + Num> Mul<&'b Complex<T>> for &'a Complex<T> {
  151. type Output = Complex<T>;
  152. #[inline]
  153. fn mul(self, other: &Complex<T>) -> Complex<T> {
  154. Complex::new(self.re.clone() * other.re.clone() - self.im.clone() * other.im.clone(),
  155. self.re.clone() * other.im.clone() + self.im.clone() * other.re.clone())
  156. }
  157. }
  158. forward_all_binop!(impl Div, div);
  159. // (a + i b) / (c + i d) == [(a + i b) * (c - i d)] / (c*c + d*d)
  160. // == [(a*c + b*d) / (c*c + d*d)] + i [(b*c - a*d) / (c*c + d*d)]
  161. impl<'a, 'b, T: Clone + Num> Div<&'b Complex<T>> for &'a Complex<T> {
  162. type Output = Complex<T>;
  163. #[inline]
  164. fn div(self, other: &Complex<T>) -> Complex<T> {
  165. let norm_sqr = other.norm_sqr();
  166. Complex::new((self.re.clone() * other.re.clone() + self.im.clone() * other.im.clone()) / norm_sqr.clone(),
  167. (self.im.clone() * other.re.clone() - self.re.clone() * other.im.clone()) / norm_sqr)
  168. }
  169. }
  170. impl<T: Clone + Num> Neg for Complex<T> {
  171. type Output = Complex<T>;
  172. #[inline]
  173. fn neg(self) -> Complex<T> { -&self }
  174. }
  175. impl<'a, T: Clone + Num> Neg for &'a Complex<T> {
  176. type Output = Complex<T>;
  177. #[inline]
  178. fn neg(self) -> Complex<T> {
  179. Complex::new(-self.re.clone(), -self.im.clone())
  180. }
  181. }
  182. /* constants */
  183. impl<T: Clone + Num> Zero for Complex<T> {
  184. #[inline]
  185. fn zero() -> Complex<T> {
  186. Complex::new(Zero::zero(), Zero::zero())
  187. }
  188. #[inline]
  189. fn is_zero(&self) -> bool {
  190. self.re.is_zero() && self.im.is_zero()
  191. }
  192. }
  193. impl<T: Clone + Num> One for Complex<T> {
  194. #[inline]
  195. fn one() -> Complex<T> {
  196. Complex::new(One::one(), Zero::zero())
  197. }
  198. }
  199. /* string conversions */
  200. impl<T: fmt::Display + Num + PartialOrd + Clone> fmt::Display for Complex<T> {
  201. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  202. if self.im < Zero::zero() {
  203. write!(f, "{}-{}i", self.re, -self.im.clone())
  204. } else {
  205. write!(f, "{}+{}i", self.re, self.im)
  206. }
  207. }
  208. }
  209. #[cfg(test)]
  210. mod test {
  211. #![allow(non_upper_case_globals)]
  212. use super::{Complex64, Complex};
  213. use std::f64;
  214. use std::num::Float;
  215. use {Zero, One};
  216. pub const _0_0i : Complex64 = Complex { re: 0.0, im: 0.0 };
  217. pub const _1_0i : Complex64 = Complex { re: 1.0, im: 0.0 };
  218. pub const _1_1i : Complex64 = Complex { re: 1.0, im: 1.0 };
  219. pub const _0_1i : Complex64 = Complex { re: 0.0, im: 1.0 };
  220. pub const _neg1_1i : Complex64 = Complex { re: -1.0, im: 1.0 };
  221. pub const _05_05i : Complex64 = Complex { re: 0.5, im: 0.5 };
  222. pub const all_consts : [Complex64; 5] = [_0_0i, _1_0i, _1_1i, _neg1_1i, _05_05i];
  223. #[test]
  224. fn test_consts() {
  225. // check our constants are what Complex::new creates
  226. fn test(c : Complex64, r : f64, i: f64) {
  227. assert_eq!(c, Complex::new(r,i));
  228. }
  229. test(_0_0i, 0.0, 0.0);
  230. test(_1_0i, 1.0, 0.0);
  231. test(_1_1i, 1.0, 1.0);
  232. test(_neg1_1i, -1.0, 1.0);
  233. test(_05_05i, 0.5, 0.5);
  234. assert_eq!(_0_0i, Zero::zero());
  235. assert_eq!(_1_0i, One::one());
  236. }
  237. #[test]
  238. #[cfg_attr(target_arch = "x86", ignore)]
  239. // FIXME #7158: (maybe?) currently failing on x86.
  240. fn test_norm() {
  241. fn test(c: Complex64, ns: f64) {
  242. assert_eq!(c.norm_sqr(), ns);
  243. assert_eq!(c.norm(), ns.sqrt())
  244. }
  245. test(_0_0i, 0.0);
  246. test(_1_0i, 1.0);
  247. test(_1_1i, 2.0);
  248. test(_neg1_1i, 2.0);
  249. test(_05_05i, 0.5);
  250. }
  251. #[test]
  252. fn test_scale_unscale() {
  253. assert_eq!(_05_05i.scale(2.0), _1_1i);
  254. assert_eq!(_1_1i.unscale(2.0), _05_05i);
  255. for &c in all_consts.iter() {
  256. assert_eq!(c.scale(2.0).unscale(2.0), c);
  257. }
  258. }
  259. #[test]
  260. fn test_conj() {
  261. for &c in all_consts.iter() {
  262. assert_eq!(c.conj(), Complex::new(c.re, -c.im));
  263. assert_eq!(c.conj().conj(), c);
  264. }
  265. }
  266. #[test]
  267. fn test_inv() {
  268. assert_eq!(_1_1i.inv(), _05_05i.conj());
  269. assert_eq!(_1_0i.inv(), _1_0i.inv());
  270. }
  271. #[test]
  272. #[should_fail]
  273. fn test_divide_by_zero_natural() {
  274. let n = Complex::new(2, 3);
  275. let d = Complex::new(0, 0);
  276. let _x = n / d;
  277. }
  278. #[test]
  279. #[should_fail]
  280. #[ignore]
  281. fn test_inv_zero() {
  282. // FIXME #5736: should this really fail, or just NaN?
  283. _0_0i.inv();
  284. }
  285. #[test]
  286. fn test_arg() {
  287. fn test(c: Complex64, arg: f64) {
  288. assert!((c.arg() - arg).abs() < 1.0e-6)
  289. }
  290. test(_1_0i, 0.0);
  291. test(_1_1i, 0.25 * f64::consts::PI);
  292. test(_neg1_1i, 0.75 * f64::consts::PI);
  293. test(_05_05i, 0.25 * f64::consts::PI);
  294. }
  295. #[test]
  296. fn test_polar_conv() {
  297. fn test(c: Complex64) {
  298. let (r, theta) = c.to_polar();
  299. assert!((c - Complex::from_polar(&r, &theta)).norm() < 1e-6);
  300. }
  301. for &c in all_consts.iter() { test(c); }
  302. }
  303. mod arith {
  304. use super::{_0_0i, _1_0i, _1_1i, _0_1i, _neg1_1i, _05_05i, all_consts};
  305. use Zero;
  306. #[test]
  307. fn test_add() {
  308. assert_eq!(_05_05i + _05_05i, _1_1i);
  309. assert_eq!(_0_1i + _1_0i, _1_1i);
  310. assert_eq!(_1_0i + _neg1_1i, _0_1i);
  311. for &c in all_consts.iter() {
  312. assert_eq!(_0_0i + c, c);
  313. assert_eq!(c + _0_0i, c);
  314. }
  315. }
  316. #[test]
  317. fn test_sub() {
  318. assert_eq!(_05_05i - _05_05i, _0_0i);
  319. assert_eq!(_0_1i - _1_0i, _neg1_1i);
  320. assert_eq!(_0_1i - _neg1_1i, _1_0i);
  321. for &c in all_consts.iter() {
  322. assert_eq!(c - _0_0i, c);
  323. assert_eq!(c - c, _0_0i);
  324. }
  325. }
  326. #[test]
  327. fn test_mul() {
  328. assert_eq!(_05_05i * _05_05i, _0_1i.unscale(2.0));
  329. assert_eq!(_1_1i * _0_1i, _neg1_1i);
  330. // i^2 & i^4
  331. assert_eq!(_0_1i * _0_1i, -_1_0i);
  332. assert_eq!(_0_1i * _0_1i * _0_1i * _0_1i, _1_0i);
  333. for &c in all_consts.iter() {
  334. assert_eq!(c * _1_0i, c);
  335. assert_eq!(_1_0i * c, c);
  336. }
  337. }
  338. #[test]
  339. fn test_div() {
  340. assert_eq!(_neg1_1i / _0_1i, _1_1i);
  341. for &c in all_consts.iter() {
  342. if c != Zero::zero() {
  343. assert_eq!(c / c, _1_0i);
  344. }
  345. }
  346. }
  347. #[test]
  348. fn test_neg() {
  349. assert_eq!(-_1_0i + _0_1i, _neg1_1i);
  350. assert_eq!((-_0_1i) * _0_1i, _1_0i);
  351. for &c in all_consts.iter() {
  352. assert_eq!(-(-c), c);
  353. }
  354. }
  355. }
  356. #[test]
  357. fn test_to_string() {
  358. fn test(c : Complex64, s: String) {
  359. assert_eq!(c.to_string(), s);
  360. }
  361. test(_0_0i, "0+0i".to_string());
  362. test(_1_0i, "1+0i".to_string());
  363. test(_0_1i, "0+1i".to_string());
  364. test(_1_1i, "1+1i".to_string());
  365. test(_neg1_1i, "-1+1i".to_string());
  366. test(-_neg1_1i, "1-1i".to_string());
  367. test(_05_05i, "0.5+0.5i".to_string());
  368. }
  369. #[test]
  370. fn test_hash() {
  371. let a = Complex::new(0i32, 0i32);
  372. let b = Complex::new(1i32, 0i32);
  373. let c = Complex::new(0i32, 1i32);
  374. assert!(::hash(&a) != ::hash(&b));
  375. assert!(::hash(&b) != ::hash(&c));
  376. assert!(::hash(&c) != ::hash(&a));
  377. }
  378. }