lib.rs 45 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401
  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. #![doc(html_logo_url = "https://rust-num.github.io/num/rust-logo-128x128-blk-v2.png",
  12. html_favicon_url = "https://rust-num.github.io/num/favicon.ico",
  13. html_root_url = "https://rust-num.github.io/num/",
  14. html_playground_url = "http://play.integer32.com/")]
  15. extern crate num_traits as traits;
  16. #[cfg(feature = "rustc-serialize")]
  17. extern crate rustc_serialize;
  18. #[cfg(feature = "serde")]
  19. extern crate serde;
  20. use std::fmt;
  21. #[cfg(test)]
  22. use std::hash;
  23. use std::ops::{Add, Div, Mul, Neg, Sub};
  24. use traits::{Zero, One, Num, Float};
  25. // FIXME #1284: handle complex NaN & infinity etc. This
  26. // probably doesn't map to C's _Complex correctly.
  27. /// A complex number in Cartesian form.
  28. ///
  29. /// ## Representation and Foreign Function Interface Compatibility
  30. ///
  31. /// `Complex<T>` is memory layout compatible with an array `[T; 2]`.
  32. ///
  33. /// Note that `Complex<F>` where F is a floating point type is **only** memory
  34. /// layout compatible with C's complex types, **not** necessarily calling
  35. /// convention compatible. This means that for FFI you can only pass
  36. /// `Complex<F>` behind a pointer, not as a value.
  37. ///
  38. /// ## Examples
  39. ///
  40. /// Example of extern function declaration.
  41. ///
  42. /// ```
  43. /// use num_complex::Complex;
  44. /// use std::os::raw::c_int;
  45. ///
  46. /// extern "C" {
  47. /// fn zaxpy_(n: *const c_int, alpha: *const Complex<f64>,
  48. /// x: *const Complex<f64>, incx: *const c_int,
  49. /// y: *mut Complex<f64>, incy: *const c_int);
  50. /// }
  51. /// ```
  52. #[derive(PartialEq, Eq, Copy, Clone, Hash, Debug, Default)]
  53. #[cfg_attr(feature = "rustc-serialize", derive(RustcEncodable, RustcDecodable))]
  54. #[repr(C)]
  55. pub struct Complex<T> {
  56. /// Real portion of the complex number
  57. pub re: T,
  58. /// Imaginary portion of the complex number
  59. pub im: T
  60. }
  61. pub type Complex32 = Complex<f32>;
  62. pub type Complex64 = Complex<f64>;
  63. impl<T: Clone + Num> Complex<T> {
  64. /// Create a new Complex
  65. #[inline]
  66. pub fn new(re: T, im: T) -> Complex<T> {
  67. Complex { re: re, im: im }
  68. }
  69. /// Returns imaginary unit
  70. #[inline]
  71. pub fn i() -> Complex<T> {
  72. Self::new(T::zero(), T::one())
  73. }
  74. /// Returns the square of the norm (since `T` doesn't necessarily
  75. /// have a sqrt function), i.e. `re^2 + im^2`.
  76. #[inline]
  77. pub fn norm_sqr(&self) -> T {
  78. self.re.clone() * self.re.clone() + self.im.clone() * self.im.clone()
  79. }
  80. /// Multiplies `self` by the scalar `t`.
  81. #[inline]
  82. pub fn scale(&self, t: T) -> Complex<T> {
  83. Complex::new(self.re.clone() * t.clone(), self.im.clone() * t)
  84. }
  85. /// Divides `self` by the scalar `t`.
  86. #[inline]
  87. pub fn unscale(&self, t: T) -> Complex<T> {
  88. Complex::new(self.re.clone() / t.clone(), self.im.clone() / t)
  89. }
  90. }
  91. impl<T: Clone + Num + Neg<Output = T>> Complex<T> {
  92. /// Returns the complex conjugate. i.e. `re - i im`
  93. #[inline]
  94. pub fn conj(&self) -> Complex<T> {
  95. Complex::new(self.re.clone(), -self.im.clone())
  96. }
  97. /// Returns `1/self`
  98. #[inline]
  99. pub fn inv(&self) -> Complex<T> {
  100. let norm_sqr = self.norm_sqr();
  101. Complex::new(self.re.clone() / norm_sqr.clone(),
  102. -self.im.clone() / norm_sqr)
  103. }
  104. }
  105. impl<T: Clone + Float> Complex<T> {
  106. /// Calculate |self|
  107. #[inline]
  108. pub fn norm(&self) -> T {
  109. self.re.hypot(self.im)
  110. }
  111. /// Calculate the principal Arg of self.
  112. #[inline]
  113. pub fn arg(&self) -> T {
  114. self.im.atan2(self.re)
  115. }
  116. /// Convert to polar form (r, theta), such that `self = r * exp(i
  117. /// * theta)`
  118. #[inline]
  119. pub fn to_polar(&self) -> (T, T) {
  120. (self.norm(), self.arg())
  121. }
  122. /// Convert a polar representation into a complex number.
  123. #[inline]
  124. pub fn from_polar(r: &T, theta: &T) -> Complex<T> {
  125. Complex::new(*r * theta.cos(), *r * theta.sin())
  126. }
  127. /// Computes `e^(self)`, where `e` is the base of the natural logarithm.
  128. #[inline]
  129. pub fn exp(&self) -> Complex<T> {
  130. // formula: e^(a + bi) = e^a (cos(b) + i*sin(b))
  131. // = from_polar(e^a, b)
  132. Complex::from_polar(&self.re.exp(), &self.im)
  133. }
  134. /// Computes the principal value of natural logarithm of `self`.
  135. ///
  136. /// This function has one branch cut:
  137. ///
  138. /// * `(-∞, 0]`, continuous from above.
  139. ///
  140. /// The branch satisfies `-π ≤ arg(ln(z)) ≤ π`.
  141. #[inline]
  142. pub fn ln(&self) -> Complex<T> {
  143. // formula: ln(z) = ln|z| + i*arg(z)
  144. let (r, theta) = self.to_polar();
  145. Complex::new(r.ln(), theta)
  146. }
  147. /// Computes the principal value of the square root of `self`.
  148. ///
  149. /// This function has one branch cut:
  150. ///
  151. /// * `(-∞, 0)`, continuous from above.
  152. ///
  153. /// The branch satisfies `-π/2 ≤ arg(sqrt(z)) ≤ π/2`.
  154. #[inline]
  155. pub fn sqrt(&self) -> Complex<T> {
  156. // formula: sqrt(r e^(it)) = sqrt(r) e^(it/2)
  157. let two = T::one() + T::one();
  158. let (r, theta) = self.to_polar();
  159. Complex::from_polar(&(r.sqrt()), &(theta/two))
  160. }
  161. /// Raises `self` to a floating point power.
  162. #[inline]
  163. pub fn powf(&self, exp: T) -> Complex<T> {
  164. // formula: x^y = (ρ e^(i θ))^y = ρ^y e^(i θ y)
  165. // = from_polar(ρ^y, θ y)
  166. let (r, theta) = self.to_polar();
  167. Complex::from_polar(&r.powf(exp), &(theta*exp))
  168. }
  169. /// Returns the logarithm of `self` with respect to an arbitrary base.
  170. #[inline]
  171. pub fn log(&self, base: T) -> Complex<T> {
  172. // formula: log_y(x) = log_y(ρ e^(i θ))
  173. // = log_y(ρ) + log_y(e^(i θ)) = log_y(ρ) + ln(e^(i θ)) / ln(y)
  174. // = log_y(ρ) + i θ / ln(y)
  175. let (r, theta) = self.to_polar();
  176. Complex::new(r.log(base), theta / base.ln())
  177. }
  178. /// Raises `self` to a complex power.
  179. #[inline]
  180. pub fn powc(&self, exp: Complex<T>) -> Complex<T> {
  181. // formula: x^y = (a + i b)^(c + i d)
  182. // = (ρ e^(i θ))^c (ρ e^(i θ))^(i d)
  183. // where ρ=|x| and θ=arg(x)
  184. // = ρ^c e^(−d θ) e^(i c θ) ρ^(i d)
  185. // = p^c e^(−d θ) (cos(c θ)
  186. // + i sin(c θ)) (cos(d ln(ρ)) + i sin(d ln(ρ)))
  187. // = p^c e^(−d θ) (
  188. // cos(c θ) cos(d ln(ρ)) − sin(c θ) sin(d ln(ρ))
  189. // + i(cos(c θ) sin(d ln(ρ)) + sin(c θ) cos(d ln(ρ))))
  190. // = p^c e^(−d θ) (cos(c θ + d ln(ρ)) + i sin(c θ + d ln(ρ)))
  191. // = from_polar(p^c e^(−d θ), c θ + d ln(ρ))
  192. let (r, theta) = self.to_polar();
  193. Complex::from_polar(
  194. &(r.powf(exp.re) * (-exp.im * theta).exp()),
  195. &(exp.re * theta + exp.im * r.ln()))
  196. }
  197. /// Raises a floating point number to the complex power `self`.
  198. #[inline]
  199. pub fn expf(&self, base: T) -> Complex<T> {
  200. // formula: x^(a+bi) = x^a x^bi = x^a e^(b ln(x) i)
  201. // = from_polar(x^a, b ln(x))
  202. Complex::from_polar(&base.powf(self.re), &(self.im * base.ln()))
  203. }
  204. /// Computes the sine of `self`.
  205. #[inline]
  206. pub fn sin(&self) -> Complex<T> {
  207. // formula: sin(a + bi) = sin(a)cosh(b) + i*cos(a)sinh(b)
  208. Complex::new(self.re.sin() * self.im.cosh(), self.re.cos() * self.im.sinh())
  209. }
  210. /// Computes the cosine of `self`.
  211. #[inline]
  212. pub fn cos(&self) -> Complex<T> {
  213. // formula: cos(a + bi) = cos(a)cosh(b) - i*sin(a)sinh(b)
  214. Complex::new(self.re.cos() * self.im.cosh(), -self.re.sin() * self.im.sinh())
  215. }
  216. /// Computes the tangent of `self`.
  217. #[inline]
  218. pub fn tan(&self) -> Complex<T> {
  219. // formula: tan(a + bi) = (sin(2a) + i*sinh(2b))/(cos(2a) + cosh(2b))
  220. let (two_re, two_im) = (self.re + self.re, self.im + self.im);
  221. Complex::new(two_re.sin(), two_im.sinh()).unscale(two_re.cos() + two_im.cosh())
  222. }
  223. /// Computes the principal value of the inverse sine of `self`.
  224. ///
  225. /// This function has two branch cuts:
  226. ///
  227. /// * `(-∞, -1)`, continuous from above.
  228. /// * `(1, ∞)`, continuous from below.
  229. ///
  230. /// The branch satisfies `-π/2 ≤ Re(asin(z)) ≤ π/2`.
  231. #[inline]
  232. pub fn asin(&self) -> Complex<T> {
  233. // formula: arcsin(z) = -i ln(sqrt(1-z^2) + iz)
  234. let i = Complex::i();
  235. -i*((Complex::one() - self*self).sqrt() + i*self).ln()
  236. }
  237. /// Computes the principal value of the inverse cosine of `self`.
  238. ///
  239. /// This function has two branch cuts:
  240. ///
  241. /// * `(-∞, -1)`, continuous from above.
  242. /// * `(1, ∞)`, continuous from below.
  243. ///
  244. /// The branch satisfies `0 ≤ Re(acos(z)) ≤ π`.
  245. #[inline]
  246. pub fn acos(&self) -> Complex<T> {
  247. // formula: arccos(z) = -i ln(i sqrt(1-z^2) + z)
  248. let i = Complex::i();
  249. -i*(i*(Complex::one() - self*self).sqrt() + self).ln()
  250. }
  251. /// Computes the principal value of the inverse tangent of `self`.
  252. ///
  253. /// This function has two branch cuts:
  254. ///
  255. /// * `(-∞i, -i]`, continuous from the left.
  256. /// * `[i, ∞i)`, continuous from the right.
  257. ///
  258. /// The branch satisfies `-π/2 ≤ Re(atan(z)) ≤ π/2`.
  259. #[inline]
  260. pub fn atan(&self) -> Complex<T> {
  261. // formula: arctan(z) = (ln(1+iz) - ln(1-iz))/(2i)
  262. let i = Complex::i();
  263. let one = Complex::one();
  264. let two = one + one;
  265. if *self == i {
  266. return Complex::new(T::zero(), T::infinity());
  267. }
  268. else if *self == -i {
  269. return Complex::new(T::zero(), -T::infinity());
  270. }
  271. ((one + i * self).ln() - (one - i * self).ln()) / (two * i)
  272. }
  273. /// Computes the hyperbolic sine of `self`.
  274. #[inline]
  275. pub fn sinh(&self) -> Complex<T> {
  276. // formula: sinh(a + bi) = sinh(a)cos(b) + i*cosh(a)sin(b)
  277. Complex::new(self.re.sinh() * self.im.cos(), self.re.cosh() * self.im.sin())
  278. }
  279. /// Computes the hyperbolic cosine of `self`.
  280. #[inline]
  281. pub fn cosh(&self) -> Complex<T> {
  282. // formula: cosh(a + bi) = cosh(a)cos(b) + i*sinh(a)sin(b)
  283. Complex::new(self.re.cosh() * self.im.cos(), self.re.sinh() * self.im.sin())
  284. }
  285. /// Computes the hyperbolic tangent of `self`.
  286. #[inline]
  287. pub fn tanh(&self) -> Complex<T> {
  288. // formula: tanh(a + bi) = (sinh(2a) + i*sin(2b))/(cosh(2a) + cos(2b))
  289. let (two_re, two_im) = (self.re + self.re, self.im + self.im);
  290. Complex::new(two_re.sinh(), two_im.sin()).unscale(two_re.cosh() + two_im.cos())
  291. }
  292. /// Computes the principal value of inverse hyperbolic sine of `self`.
  293. ///
  294. /// This function has two branch cuts:
  295. ///
  296. /// * `(-∞i, -i)`, continuous from the left.
  297. /// * `(i, ∞i)`, continuous from the right.
  298. ///
  299. /// The branch satisfies `-π/2 ≤ Im(asinh(z)) ≤ π/2`.
  300. #[inline]
  301. pub fn asinh(&self) -> Complex<T> {
  302. // formula: arcsinh(z) = ln(z + sqrt(1+z^2))
  303. let one = Complex::one();
  304. (self + (one + self * self).sqrt()).ln()
  305. }
  306. /// Computes the principal value of inverse hyperbolic cosine of `self`.
  307. ///
  308. /// This function has one branch cut:
  309. ///
  310. /// * `(-∞, 1)`, continuous from above.
  311. ///
  312. /// The branch satisfies `-π ≤ Im(acosh(z)) ≤ π` and `0 ≤ Re(acosh(z)) < ∞`.
  313. #[inline]
  314. pub fn acosh(&self) -> Complex<T> {
  315. // formula: arccosh(z) = 2 ln(sqrt((z+1)/2) + sqrt((z-1)/2))
  316. let one = Complex::one();
  317. let two = one + one;
  318. two * (((self + one)/two).sqrt() + ((self - one)/two).sqrt()).ln()
  319. }
  320. /// Computes the principal value of inverse hyperbolic tangent of `self`.
  321. ///
  322. /// This function has two branch cuts:
  323. ///
  324. /// * `(-∞, -1]`, continuous from above.
  325. /// * `[1, ∞)`, continuous from below.
  326. ///
  327. /// The branch satisfies `-π/2 ≤ Im(atanh(z)) ≤ π/2`.
  328. #[inline]
  329. pub fn atanh(&self) -> Complex<T> {
  330. // formula: arctanh(z) = (ln(1+z) - ln(1-z))/2
  331. let one = Complex::one();
  332. let two = one + one;
  333. if *self == one {
  334. return Complex::new(T::infinity(), T::zero());
  335. }
  336. else if *self == -one {
  337. return Complex::new(-T::infinity(), T::zero());
  338. }
  339. ((one + self).ln() - (one - self).ln()) / two
  340. }
  341. /// Checks if the given complex number is NaN
  342. #[inline]
  343. pub fn is_nan(self) -> bool {
  344. self.re.is_nan() || self.im.is_nan()
  345. }
  346. /// Checks if the given complex number is infinite
  347. #[inline]
  348. pub fn is_infinite(self) -> bool {
  349. !self.is_nan() && (self.re.is_infinite() || self.im.is_infinite())
  350. }
  351. /// Checks if the given complex number is finite
  352. #[inline]
  353. pub fn is_finite(self) -> bool {
  354. self.re.is_finite() && self.im.is_finite()
  355. }
  356. /// Checks if the given complex number is normal
  357. #[inline]
  358. pub fn is_normal(self) -> bool {
  359. self.re.is_normal() && self.im.is_normal()
  360. }
  361. }
  362. impl<T: Clone + Num> From<T> for Complex<T> {
  363. #[inline]
  364. fn from(re: T) -> Complex<T> {
  365. Complex { re: re, im: T::zero() }
  366. }
  367. }
  368. impl<'a, T: Clone + Num> From<&'a T> for Complex<T> {
  369. #[inline]
  370. fn from(re: &T) -> Complex<T> {
  371. From::from(re.clone())
  372. }
  373. }
  374. macro_rules! forward_ref_ref_binop {
  375. (impl $imp:ident, $method:ident) => {
  376. impl<'a, 'b, T: Clone + Num> $imp<&'b Complex<T>> for &'a Complex<T> {
  377. type Output = Complex<T>;
  378. #[inline]
  379. fn $method(self, other: &Complex<T>) -> Complex<T> {
  380. self.clone().$method(other.clone())
  381. }
  382. }
  383. }
  384. }
  385. macro_rules! forward_ref_val_binop {
  386. (impl $imp:ident, $method:ident) => {
  387. impl<'a, T: Clone + Num> $imp<Complex<T>> for &'a Complex<T> {
  388. type Output = Complex<T>;
  389. #[inline]
  390. fn $method(self, other: Complex<T>) -> Complex<T> {
  391. self.clone().$method(other)
  392. }
  393. }
  394. }
  395. }
  396. macro_rules! forward_val_ref_binop {
  397. (impl $imp:ident, $method:ident) => {
  398. impl<'a, T: Clone + Num> $imp<&'a Complex<T>> for Complex<T> {
  399. type Output = Complex<T>;
  400. #[inline]
  401. fn $method(self, other: &Complex<T>) -> Complex<T> {
  402. self.$method(other.clone())
  403. }
  404. }
  405. }
  406. }
  407. macro_rules! forward_all_binop {
  408. (impl $imp:ident, $method:ident) => {
  409. forward_ref_ref_binop!(impl $imp, $method);
  410. forward_ref_val_binop!(impl $imp, $method);
  411. forward_val_ref_binop!(impl $imp, $method);
  412. };
  413. }
  414. /* arithmetic */
  415. forward_all_binop!(impl Add, add);
  416. // (a + i b) + (c + i d) == (a + c) + i (b + d)
  417. impl<T: Clone + Num> Add<Complex<T>> for Complex<T> {
  418. type Output = Complex<T>;
  419. #[inline]
  420. fn add(self, other: Complex<T>) -> Complex<T> {
  421. Complex::new(self.re + other.re, self.im + other.im)
  422. }
  423. }
  424. forward_all_binop!(impl Sub, sub);
  425. // (a + i b) - (c + i d) == (a - c) + i (b - d)
  426. impl<T: Clone + Num> Sub<Complex<T>> for Complex<T> {
  427. type Output = Complex<T>;
  428. #[inline]
  429. fn sub(self, other: Complex<T>) -> Complex<T> {
  430. Complex::new(self.re - other.re, self.im - other.im)
  431. }
  432. }
  433. forward_all_binop!(impl Mul, mul);
  434. // (a + i b) * (c + i d) == (a*c - b*d) + i (a*d + b*c)
  435. impl<T: Clone + Num> Mul<Complex<T>> for Complex<T> {
  436. type Output = Complex<T>;
  437. #[inline]
  438. fn mul(self, other: Complex<T>) -> Complex<T> {
  439. let re = self.re.clone() * other.re.clone() - self.im.clone() * other.im.clone();
  440. let im = self.re * other.im + self.im * other.re;
  441. Complex::new(re, im)
  442. }
  443. }
  444. forward_all_binop!(impl Div, div);
  445. // (a + i b) / (c + i d) == [(a + i b) * (c - i d)] / (c*c + d*d)
  446. // == [(a*c + b*d) / (c*c + d*d)] + i [(b*c - a*d) / (c*c + d*d)]
  447. impl<T: Clone + Num> Div<Complex<T>> for Complex<T> {
  448. type Output = Complex<T>;
  449. #[inline]
  450. fn div(self, other: Complex<T>) -> Complex<T> {
  451. let norm_sqr = other.norm_sqr();
  452. let re = self.re.clone() * other.re.clone() + self.im.clone() * other.im.clone();
  453. let im = self.im * other.re - self.re * other.im;
  454. Complex::new(re / norm_sqr.clone(), im / norm_sqr)
  455. }
  456. }
  457. impl<T: Clone + Num + Neg<Output = T>> Neg for Complex<T> {
  458. type Output = Complex<T>;
  459. #[inline]
  460. fn neg(self) -> Complex<T> {
  461. Complex::new(-self.re, -self.im)
  462. }
  463. }
  464. impl<'a, T: Clone + Num + Neg<Output = T>> Neg for &'a Complex<T> {
  465. type Output = Complex<T>;
  466. #[inline]
  467. fn neg(self) -> Complex<T> {
  468. -self.clone()
  469. }
  470. }
  471. macro_rules! real_arithmetic {
  472. (@forward $imp:ident::$method:ident for $($real:ident),*) => (
  473. impl<'a, T: Clone + Num> $imp<&'a T> for Complex<T> {
  474. type Output = Complex<T>;
  475. #[inline]
  476. fn $method(self, other: &T) -> Complex<T> {
  477. self.$method(other.clone())
  478. }
  479. }
  480. impl<'a, T: Clone + Num> $imp<T> for &'a Complex<T> {
  481. type Output = Complex<T>;
  482. #[inline]
  483. fn $method(self, other: T) -> Complex<T> {
  484. self.clone().$method(other)
  485. }
  486. }
  487. impl<'a, 'b, T: Clone + Num> $imp<&'a T> for &'b Complex<T> {
  488. type Output = Complex<T>;
  489. #[inline]
  490. fn $method(self, other: &T) -> Complex<T> {
  491. self.clone().$method(other.clone())
  492. }
  493. }
  494. $(
  495. impl<'a> $imp<&'a Complex<$real>> for $real {
  496. type Output = Complex<$real>;
  497. #[inline]
  498. fn $method(self, other: &Complex<$real>) -> Complex<$real> {
  499. self.$method(other.clone())
  500. }
  501. }
  502. impl<'a> $imp<Complex<$real>> for &'a $real {
  503. type Output = Complex<$real>;
  504. #[inline]
  505. fn $method(self, other: Complex<$real>) -> Complex<$real> {
  506. self.clone().$method(other)
  507. }
  508. }
  509. impl<'a, 'b> $imp<&'a Complex<$real>> for &'b $real {
  510. type Output = Complex<$real>;
  511. #[inline]
  512. fn $method(self, other: &Complex<$real>) -> Complex<$real> {
  513. self.clone().$method(other.clone())
  514. }
  515. }
  516. )*
  517. );
  518. (@implement $imp:ident::$method:ident for $($real:ident),*) => (
  519. impl<T: Clone + Num> $imp<T> for Complex<T> {
  520. type Output = Complex<T>;
  521. #[inline]
  522. fn $method(self, other: T) -> Complex<T> {
  523. self.$method(Complex::from(other))
  524. }
  525. }
  526. $(
  527. impl $imp<Complex<$real>> for $real {
  528. type Output = Complex<$real>;
  529. #[inline]
  530. fn $method(self, other: Complex<$real>) -> Complex<$real> {
  531. Complex::from(self).$method(other)
  532. }
  533. }
  534. )*
  535. );
  536. ($($real:ident),*) => (
  537. real_arithmetic!(@forward Add::add for $($real),*);
  538. real_arithmetic!(@forward Sub::sub for $($real),*);
  539. real_arithmetic!(@forward Mul::mul for $($real),*);
  540. real_arithmetic!(@forward Div::div for $($real),*);
  541. real_arithmetic!(@implement Add::add for $($real),*);
  542. real_arithmetic!(@implement Sub::sub for $($real),*);
  543. real_arithmetic!(@implement Mul::mul for $($real),*);
  544. real_arithmetic!(@implement Div::div for $($real),*);
  545. );
  546. }
  547. real_arithmetic!(usize, u8, u16, u32, u64, isize, i8, i16, i32, i64, f32, f64);
  548. /* constants */
  549. impl<T: Clone + Num> Zero for Complex<T> {
  550. #[inline]
  551. fn zero() -> Complex<T> {
  552. Complex::new(Zero::zero(), Zero::zero())
  553. }
  554. #[inline]
  555. fn is_zero(&self) -> bool {
  556. self.re.is_zero() && self.im.is_zero()
  557. }
  558. }
  559. impl<T: Clone + Num> One for Complex<T> {
  560. #[inline]
  561. fn one() -> Complex<T> {
  562. Complex::new(One::one(), Zero::zero())
  563. }
  564. }
  565. macro_rules! write_complex {
  566. ($f:ident, $t:expr, $prefix:expr, $re:expr, $im:expr, $T:ident) => {{
  567. let abs_re = if $re < Zero::zero() { $T::zero() - $re.clone() } else { $re.clone() };
  568. let abs_im = if $im < Zero::zero() { $T::zero() - $im.clone() } else { $im.clone() };
  569. let real: String;
  570. let imag: String;
  571. if let Some(prec) = $f.precision() {
  572. real = format!(concat!("{:.1$", $t, "}"), abs_re, prec);
  573. imag = format!(concat!("{:.1$", $t, "}"), abs_im, prec);
  574. }
  575. else {
  576. real = format!(concat!("{:", $t, "}"), abs_re);
  577. imag = format!(concat!("{:", $t, "}"), abs_im);
  578. }
  579. let prefix = if $f.alternate() { $prefix } else { "" };
  580. let sign = if $re < Zero::zero() {
  581. "-"
  582. } else if $f.sign_plus() {
  583. "+"
  584. } else {
  585. ""
  586. };
  587. let complex = if $im < Zero::zero() {
  588. format!("{}{pre}{re}-{pre}{im}i", sign, re=real, im=imag, pre=prefix)
  589. }
  590. else {
  591. format!("{}{pre}{re}+{pre}{im}i", sign, re=real, im=imag, pre=prefix)
  592. };
  593. if let Some(width) = $f.width() {
  594. write!($f, "{0: >1$}", complex, width)
  595. }
  596. else {
  597. write!($f, "{}", complex)
  598. }
  599. }}
  600. }
  601. /* string conversions */
  602. impl<T> fmt::Display for Complex<T> where
  603. T: fmt::Display + Num + PartialOrd + Clone
  604. {
  605. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  606. write_complex!(f, "", "", self.re, self.im, T)
  607. }
  608. }
  609. impl<T> fmt::LowerExp for Complex<T> where
  610. T: fmt::LowerExp + Num + PartialOrd + Clone
  611. {
  612. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  613. write_complex!(f, "e", "", self.re, self.im, T)
  614. }
  615. }
  616. impl<T> fmt::UpperExp for Complex<T> where
  617. T: fmt::UpperExp + Num + PartialOrd + Clone
  618. {
  619. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  620. write_complex!(f, "E", "", self.re, self.im, T)
  621. }
  622. }
  623. impl<T> fmt::LowerHex for Complex<T> where
  624. T: fmt::LowerHex + Num + PartialOrd + Clone
  625. {
  626. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  627. write_complex!(f, "x", "0x", self.re, self.im, T)
  628. }
  629. }
  630. impl<T> fmt::UpperHex for Complex<T> where
  631. T: fmt::UpperHex + Num + PartialOrd + Clone
  632. {
  633. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  634. write_complex!(f, "X", "0x", self.re, self.im, T)
  635. }
  636. }
  637. impl<T> fmt::Octal for Complex<T> where
  638. T: fmt::Octal + Num + PartialOrd + Clone
  639. {
  640. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  641. write_complex!(f, "o", "0o", self.re, self.im, T)
  642. }
  643. }
  644. impl<T> fmt::Binary for Complex<T> where
  645. T: fmt::Binary + Num + PartialOrd + Clone
  646. {
  647. fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
  648. write_complex!(f, "b", "0b", self.re, self.im, T)
  649. }
  650. }
  651. #[cfg(feature = "serde")]
  652. impl<T> serde::Serialize for Complex<T>
  653. where T: serde::Serialize
  654. {
  655. fn serialize<S>(&self, serializer: &mut S) -> Result<(), S::Error> where
  656. S: serde::Serializer
  657. {
  658. (&self.re, &self.im).serialize(serializer)
  659. }
  660. }
  661. #[cfg(feature = "serde")]
  662. impl<T> serde::Deserialize for Complex<T> where
  663. T: serde::Deserialize + Num + Clone
  664. {
  665. fn deserialize<D>(deserializer: &mut D) -> Result<Self, D::Error> where
  666. D: serde::Deserializer,
  667. {
  668. let (re, im) = try!(serde::Deserialize::deserialize(deserializer));
  669. Ok(Complex::new(re, im))
  670. }
  671. }
  672. #[cfg(test)]
  673. fn hash<T: hash::Hash>(x: &T) -> u64 {
  674. use std::hash::Hasher;
  675. let mut hasher = hash::SipHasher::new();
  676. x.hash(&mut hasher);
  677. hasher.finish()
  678. }
  679. #[cfg(test)]
  680. mod test {
  681. #![allow(non_upper_case_globals)]
  682. use super::{Complex64, Complex};
  683. use std::f64;
  684. use traits::{Zero, One, Float};
  685. pub const _0_0i : Complex64 = Complex { re: 0.0, im: 0.0 };
  686. pub const _1_0i : Complex64 = Complex { re: 1.0, im: 0.0 };
  687. pub const _1_1i : Complex64 = Complex { re: 1.0, im: 1.0 };
  688. pub const _0_1i : Complex64 = Complex { re: 0.0, im: 1.0 };
  689. pub const _neg1_1i : Complex64 = Complex { re: -1.0, im: 1.0 };
  690. pub const _05_05i : Complex64 = Complex { re: 0.5, im: 0.5 };
  691. pub const all_consts : [Complex64; 5] = [_0_0i, _1_0i, _1_1i, _neg1_1i, _05_05i];
  692. #[test]
  693. fn test_consts() {
  694. // check our constants are what Complex::new creates
  695. fn test(c : Complex64, r : f64, i: f64) {
  696. assert_eq!(c, Complex::new(r,i));
  697. }
  698. test(_0_0i, 0.0, 0.0);
  699. test(_1_0i, 1.0, 0.0);
  700. test(_1_1i, 1.0, 1.0);
  701. test(_neg1_1i, -1.0, 1.0);
  702. test(_05_05i, 0.5, 0.5);
  703. assert_eq!(_0_0i, Zero::zero());
  704. assert_eq!(_1_0i, One::one());
  705. }
  706. #[test]
  707. #[cfg_attr(target_arch = "x86", ignore)]
  708. // FIXME #7158: (maybe?) currently failing on x86.
  709. fn test_norm() {
  710. fn test(c: Complex64, ns: f64) {
  711. assert_eq!(c.norm_sqr(), ns);
  712. assert_eq!(c.norm(), ns.sqrt())
  713. }
  714. test(_0_0i, 0.0);
  715. test(_1_0i, 1.0);
  716. test(_1_1i, 2.0);
  717. test(_neg1_1i, 2.0);
  718. test(_05_05i, 0.5);
  719. }
  720. #[test]
  721. fn test_scale_unscale() {
  722. assert_eq!(_05_05i.scale(2.0), _1_1i);
  723. assert_eq!(_1_1i.unscale(2.0), _05_05i);
  724. for &c in all_consts.iter() {
  725. assert_eq!(c.scale(2.0).unscale(2.0), c);
  726. }
  727. }
  728. #[test]
  729. fn test_conj() {
  730. for &c in all_consts.iter() {
  731. assert_eq!(c.conj(), Complex::new(c.re, -c.im));
  732. assert_eq!(c.conj().conj(), c);
  733. }
  734. }
  735. #[test]
  736. fn test_inv() {
  737. assert_eq!(_1_1i.inv(), _05_05i.conj());
  738. assert_eq!(_1_0i.inv(), _1_0i.inv());
  739. }
  740. #[test]
  741. #[should_panic]
  742. fn test_divide_by_zero_natural() {
  743. let n = Complex::new(2, 3);
  744. let d = Complex::new(0, 0);
  745. let _x = n / d;
  746. }
  747. #[test]
  748. fn test_inv_zero() {
  749. // FIXME #20: should this really fail, or just NaN?
  750. assert!(_0_0i.inv().is_nan());
  751. }
  752. #[test]
  753. fn test_arg() {
  754. fn test(c: Complex64, arg: f64) {
  755. assert!((c.arg() - arg).abs() < 1.0e-6)
  756. }
  757. test(_1_0i, 0.0);
  758. test(_1_1i, 0.25 * f64::consts::PI);
  759. test(_neg1_1i, 0.75 * f64::consts::PI);
  760. test(_05_05i, 0.25 * f64::consts::PI);
  761. }
  762. #[test]
  763. fn test_polar_conv() {
  764. fn test(c: Complex64) {
  765. let (r, theta) = c.to_polar();
  766. assert!((c - Complex::from_polar(&r, &theta)).norm() < 1e-6);
  767. }
  768. for &c in all_consts.iter() { test(c); }
  769. }
  770. fn close(a: Complex64, b: Complex64) -> bool {
  771. close_to_tol(a, b, 1e-10)
  772. }
  773. fn close_to_tol(a: Complex64, b: Complex64, tol: f64) -> bool {
  774. // returns true if a and b are reasonably close
  775. (a == b) || (a-b).norm() < tol
  776. }
  777. #[test]
  778. fn test_exp() {
  779. assert!(close(_1_0i.exp(), _1_0i.scale(f64::consts::E)));
  780. assert!(close(_0_0i.exp(), _1_0i));
  781. assert!(close(_0_1i.exp(), Complex::new(1.0.cos(), 1.0.sin())));
  782. assert!(close(_05_05i.exp()*_05_05i.exp(), _1_1i.exp()));
  783. assert!(close(_0_1i.scale(-f64::consts::PI).exp(), _1_0i.scale(-1.0)));
  784. for &c in all_consts.iter() {
  785. // e^conj(z) = conj(e^z)
  786. assert!(close(c.conj().exp(), c.exp().conj()));
  787. // e^(z + 2 pi i) = e^z
  788. assert!(close(c.exp(), (c + _0_1i.scale(f64::consts::PI*2.0)).exp()));
  789. }
  790. }
  791. #[test]
  792. fn test_ln() {
  793. assert!(close(_1_0i.ln(), _0_0i));
  794. assert!(close(_0_1i.ln(), _0_1i.scale(f64::consts::PI/2.0)));
  795. assert!(close(_0_0i.ln(), Complex::new(f64::neg_infinity(), 0.0)));
  796. assert!(close((_neg1_1i * _05_05i).ln(), _neg1_1i.ln() + _05_05i.ln()));
  797. for &c in all_consts.iter() {
  798. // ln(conj(z() = conj(ln(z))
  799. assert!(close(c.conj().ln(), c.ln().conj()));
  800. // for this branch, -pi <= arg(ln(z)) <= pi
  801. assert!(-f64::consts::PI <= c.ln().arg() && c.ln().arg() <= f64::consts::PI);
  802. }
  803. }
  804. #[test]
  805. fn test_powc()
  806. {
  807. let a = Complex::new(2.0, -3.0);
  808. let b = Complex::new(3.0, 0.0);
  809. assert!(close(a.powc(b), a.powf(b.re)));
  810. assert!(close(b.powc(a), a.expf(b.re)));
  811. let c = Complex::new(1.0 / 3.0, 0.1);
  812. assert!(close_to_tol(a.powc(c), Complex::new(1.65826, -0.33502), 1e-5));
  813. }
  814. #[test]
  815. fn test_powf()
  816. {
  817. let c = Complex::new(2.0, -1.0);
  818. let r = c.powf(3.5);
  819. assert!(close_to_tol(r, Complex::new(-0.8684746, -16.695934), 1e-5));
  820. }
  821. #[test]
  822. fn test_log()
  823. {
  824. let c = Complex::new(2.0, -1.0);
  825. let r = c.log(10.0);
  826. assert!(close_to_tol(r, Complex::new(0.349485, -0.20135958), 1e-5));
  827. }
  828. #[test]
  829. fn test_some_expf_cases()
  830. {
  831. let c = Complex::new(2.0, -1.0);
  832. let r = c.expf(10.0);
  833. assert!(close_to_tol(r, Complex::new(-66.82015, -74.39803), 1e-5));
  834. let c = Complex::new(5.0, -2.0);
  835. let r = c.expf(3.4);
  836. assert!(close_to_tol(r, Complex::new(-349.25, -290.63), 1e-2));
  837. let c = Complex::new(-1.5, 2.0 / 3.0);
  838. let r = c.expf(1.0 / 3.0);
  839. assert!(close_to_tol(r, Complex::new(3.8637, -3.4745), 1e-2));
  840. }
  841. #[test]
  842. fn test_sqrt() {
  843. assert!(close(_0_0i.sqrt(), _0_0i));
  844. assert!(close(_1_0i.sqrt(), _1_0i));
  845. assert!(close(Complex::new(-1.0, 0.0).sqrt(), _0_1i));
  846. assert!(close(Complex::new(-1.0, -0.0).sqrt(), _0_1i.scale(-1.0)));
  847. assert!(close(_0_1i.sqrt(), _05_05i.scale(2.0.sqrt())));
  848. for &c in all_consts.iter() {
  849. // sqrt(conj(z() = conj(sqrt(z))
  850. assert!(close(c.conj().sqrt(), c.sqrt().conj()));
  851. // for this branch, -pi/2 <= arg(sqrt(z)) <= pi/2
  852. assert!(-f64::consts::PI/2.0 <= c.sqrt().arg() && c.sqrt().arg() <= f64::consts::PI/2.0);
  853. // sqrt(z) * sqrt(z) = z
  854. assert!(close(c.sqrt()*c.sqrt(), c));
  855. }
  856. }
  857. #[test]
  858. fn test_sin() {
  859. assert!(close(_0_0i.sin(), _0_0i));
  860. assert!(close(_1_0i.scale(f64::consts::PI*2.0).sin(), _0_0i));
  861. assert!(close(_0_1i.sin(), _0_1i.scale(1.0.sinh())));
  862. for &c in all_consts.iter() {
  863. // sin(conj(z)) = conj(sin(z))
  864. assert!(close(c.conj().sin(), c.sin().conj()));
  865. // sin(-z) = -sin(z)
  866. assert!(close(c.scale(-1.0).sin(), c.sin().scale(-1.0)));
  867. }
  868. }
  869. #[test]
  870. fn test_cos() {
  871. assert!(close(_0_0i.cos(), _1_0i));
  872. assert!(close(_1_0i.scale(f64::consts::PI*2.0).cos(), _1_0i));
  873. assert!(close(_0_1i.cos(), _1_0i.scale(1.0.cosh())));
  874. for &c in all_consts.iter() {
  875. // cos(conj(z)) = conj(cos(z))
  876. assert!(close(c.conj().cos(), c.cos().conj()));
  877. // cos(-z) = cos(z)
  878. assert!(close(c.scale(-1.0).cos(), c.cos()));
  879. }
  880. }
  881. #[test]
  882. fn test_tan() {
  883. assert!(close(_0_0i.tan(), _0_0i));
  884. assert!(close(_1_0i.scale(f64::consts::PI/4.0).tan(), _1_0i));
  885. assert!(close(_1_0i.scale(f64::consts::PI).tan(), _0_0i));
  886. for &c in all_consts.iter() {
  887. // tan(conj(z)) = conj(tan(z))
  888. assert!(close(c.conj().tan(), c.tan().conj()));
  889. // tan(-z) = -tan(z)
  890. assert!(close(c.scale(-1.0).tan(), c.tan().scale(-1.0)));
  891. }
  892. }
  893. #[test]
  894. fn test_asin() {
  895. assert!(close(_0_0i.asin(), _0_0i));
  896. assert!(close(_1_0i.asin(), _1_0i.scale(f64::consts::PI/2.0)));
  897. assert!(close(_1_0i.scale(-1.0).asin(), _1_0i.scale(-f64::consts::PI/2.0)));
  898. assert!(close(_0_1i.asin(), _0_1i.scale((1.0 + 2.0.sqrt()).ln())));
  899. for &c in all_consts.iter() {
  900. // asin(conj(z)) = conj(asin(z))
  901. assert!(close(c.conj().asin(), c.asin().conj()));
  902. // asin(-z) = -asin(z)
  903. assert!(close(c.scale(-1.0).asin(), c.asin().scale(-1.0)));
  904. // for this branch, -pi/2 <= asin(z).re <= pi/2
  905. assert!(-f64::consts::PI/2.0 <= c.asin().re && c.asin().re <= f64::consts::PI/2.0);
  906. }
  907. }
  908. #[test]
  909. fn test_acos() {
  910. assert!(close(_0_0i.acos(), _1_0i.scale(f64::consts::PI/2.0)));
  911. assert!(close(_1_0i.acos(), _0_0i));
  912. assert!(close(_1_0i.scale(-1.0).acos(), _1_0i.scale(f64::consts::PI)));
  913. assert!(close(_0_1i.acos(), Complex::new(f64::consts::PI/2.0, (2.0.sqrt() - 1.0).ln())));
  914. for &c in all_consts.iter() {
  915. // acos(conj(z)) = conj(acos(z))
  916. assert!(close(c.conj().acos(), c.acos().conj()));
  917. // for this branch, 0 <= acos(z).re <= pi
  918. assert!(0.0 <= c.acos().re && c.acos().re <= f64::consts::PI);
  919. }
  920. }
  921. #[test]
  922. fn test_atan() {
  923. assert!(close(_0_0i.atan(), _0_0i));
  924. assert!(close(_1_0i.atan(), _1_0i.scale(f64::consts::PI/4.0)));
  925. assert!(close(_1_0i.scale(-1.0).atan(), _1_0i.scale(-f64::consts::PI/4.0)));
  926. assert!(close(_0_1i.atan(), Complex::new(0.0, f64::infinity())));
  927. for &c in all_consts.iter() {
  928. // atan(conj(z)) = conj(atan(z))
  929. assert!(close(c.conj().atan(), c.atan().conj()));
  930. // atan(-z) = -atan(z)
  931. assert!(close(c.scale(-1.0).atan(), c.atan().scale(-1.0)));
  932. // for this branch, -pi/2 <= atan(z).re <= pi/2
  933. assert!(-f64::consts::PI/2.0 <= c.atan().re && c.atan().re <= f64::consts::PI/2.0);
  934. }
  935. }
  936. #[test]
  937. fn test_sinh() {
  938. assert!(close(_0_0i.sinh(), _0_0i));
  939. assert!(close(_1_0i.sinh(), _1_0i.scale((f64::consts::E - 1.0/f64::consts::E)/2.0)));
  940. assert!(close(_0_1i.sinh(), _0_1i.scale(1.0.sin())));
  941. for &c in all_consts.iter() {
  942. // sinh(conj(z)) = conj(sinh(z))
  943. assert!(close(c.conj().sinh(), c.sinh().conj()));
  944. // sinh(-z) = -sinh(z)
  945. assert!(close(c.scale(-1.0).sinh(), c.sinh().scale(-1.0)));
  946. }
  947. }
  948. #[test]
  949. fn test_cosh() {
  950. assert!(close(_0_0i.cosh(), _1_0i));
  951. assert!(close(_1_0i.cosh(), _1_0i.scale((f64::consts::E + 1.0/f64::consts::E)/2.0)));
  952. assert!(close(_0_1i.cosh(), _1_0i.scale(1.0.cos())));
  953. for &c in all_consts.iter() {
  954. // cosh(conj(z)) = conj(cosh(z))
  955. assert!(close(c.conj().cosh(), c.cosh().conj()));
  956. // cosh(-z) = cosh(z)
  957. assert!(close(c.scale(-1.0).cosh(), c.cosh()));
  958. }
  959. }
  960. #[test]
  961. fn test_tanh() {
  962. assert!(close(_0_0i.tanh(), _0_0i));
  963. assert!(close(_1_0i.tanh(), _1_0i.scale((f64::consts::E.powi(2) - 1.0)/(f64::consts::E.powi(2) + 1.0))));
  964. assert!(close(_0_1i.tanh(), _0_1i.scale(1.0.tan())));
  965. for &c in all_consts.iter() {
  966. // tanh(conj(z)) = conj(tanh(z))
  967. assert!(close(c.conj().tanh(), c.conj().tanh()));
  968. // tanh(-z) = -tanh(z)
  969. assert!(close(c.scale(-1.0).tanh(), c.tanh().scale(-1.0)));
  970. }
  971. }
  972. #[test]
  973. fn test_asinh() {
  974. assert!(close(_0_0i.asinh(), _0_0i));
  975. assert!(close(_1_0i.asinh(), _1_0i.scale(1.0 + 2.0.sqrt()).ln()));
  976. assert!(close(_0_1i.asinh(), _0_1i.scale(f64::consts::PI/2.0)));
  977. assert!(close(_0_1i.asinh().scale(-1.0), _0_1i.scale(-f64::consts::PI/2.0)));
  978. for &c in all_consts.iter() {
  979. // asinh(conj(z)) = conj(asinh(z))
  980. assert!(close(c.conj().asinh(), c.conj().asinh()));
  981. // asinh(-z) = -asinh(z)
  982. assert!(close(c.scale(-1.0).asinh(), c.asinh().scale(-1.0)));
  983. // for this branch, -pi/2 <= asinh(z).im <= pi/2
  984. assert!(-f64::consts::PI/2.0 <= c.asinh().im && c.asinh().im <= f64::consts::PI/2.0);
  985. }
  986. }
  987. #[test]
  988. fn test_acosh() {
  989. assert!(close(_0_0i.acosh(), _0_1i.scale(f64::consts::PI/2.0)));
  990. assert!(close(_1_0i.acosh(), _0_0i));
  991. assert!(close(_1_0i.scale(-1.0).acosh(), _0_1i.scale(f64::consts::PI)));
  992. for &c in all_consts.iter() {
  993. // acosh(conj(z)) = conj(acosh(z))
  994. assert!(close(c.conj().acosh(), c.conj().acosh()));
  995. // for this branch, -pi <= acosh(z).im <= pi and 0 <= acosh(z).re
  996. assert!(-f64::consts::PI <= c.acosh().im && c.acosh().im <= f64::consts::PI && 0.0 <= c.cosh().re);
  997. }
  998. }
  999. #[test]
  1000. fn test_atanh() {
  1001. assert!(close(_0_0i.atanh(), _0_0i));
  1002. assert!(close(_0_1i.atanh(), _0_1i.scale(f64::consts::PI/4.0)));
  1003. assert!(close(_1_0i.atanh(), Complex::new(f64::infinity(), 0.0)));
  1004. for &c in all_consts.iter() {
  1005. // atanh(conj(z)) = conj(atanh(z))
  1006. assert!(close(c.conj().atanh(), c.conj().atanh()));
  1007. // atanh(-z) = -atanh(z)
  1008. assert!(close(c.scale(-1.0).atanh(), c.atanh().scale(-1.0)));
  1009. // for this branch, -pi/2 <= atanh(z).im <= pi/2
  1010. assert!(-f64::consts::PI/2.0 <= c.atanh().im && c.atanh().im <= f64::consts::PI/2.0);
  1011. }
  1012. }
  1013. #[test]
  1014. fn test_exp_ln() {
  1015. for &c in all_consts.iter() {
  1016. // e^ln(z) = z
  1017. assert!(close(c.ln().exp(), c));
  1018. }
  1019. }
  1020. #[test]
  1021. fn test_trig_to_hyperbolic() {
  1022. for &c in all_consts.iter() {
  1023. // sin(iz) = i sinh(z)
  1024. assert!(close((_0_1i * c).sin(), _0_1i * c.sinh()));
  1025. // cos(iz) = cosh(z)
  1026. assert!(close((_0_1i * c).cos(), c.cosh()));
  1027. // tan(iz) = i tanh(z)
  1028. assert!(close((_0_1i * c).tan(), _0_1i * c.tanh()));
  1029. }
  1030. }
  1031. #[test]
  1032. fn test_trig_identities() {
  1033. for &c in all_consts.iter() {
  1034. // tan(z) = sin(z)/cos(z)
  1035. assert!(close(c.tan(), c.sin()/c.cos()));
  1036. // sin(z)^2 + cos(z)^2 = 1
  1037. assert!(close(c.sin()*c.sin() + c.cos()*c.cos(), _1_0i));
  1038. // sin(asin(z)) = z
  1039. assert!(close(c.asin().sin(), c));
  1040. // cos(acos(z)) = z
  1041. assert!(close(c.acos().cos(), c));
  1042. // tan(atan(z)) = z
  1043. // i and -i are branch points
  1044. if c != _0_1i && c != _0_1i.scale(-1.0) {
  1045. assert!(close(c.atan().tan(), c));
  1046. }
  1047. // sin(z) = (e^(iz) - e^(-iz))/(2i)
  1048. assert!(close(((_0_1i*c).exp() - (_0_1i*c).exp().inv())/_0_1i.scale(2.0), c.sin()));
  1049. // cos(z) = (e^(iz) + e^(-iz))/2
  1050. assert!(close(((_0_1i*c).exp() + (_0_1i*c).exp().inv()).unscale(2.0), c.cos()));
  1051. // tan(z) = i (1 - e^(2iz))/(1 + e^(2iz))
  1052. assert!(close(_0_1i * (_1_0i - (_0_1i*c).scale(2.0).exp())/(_1_0i + (_0_1i*c).scale(2.0).exp()), c.tan()));
  1053. }
  1054. }
  1055. #[test]
  1056. fn test_hyperbolic_identites() {
  1057. for &c in all_consts.iter() {
  1058. // tanh(z) = sinh(z)/cosh(z)
  1059. assert!(close(c.tanh(), c.sinh()/c.cosh()));
  1060. // cosh(z)^2 - sinh(z)^2 = 1
  1061. assert!(close(c.cosh()*c.cosh() - c.sinh()*c.sinh(), _1_0i));
  1062. // sinh(asinh(z)) = z
  1063. assert!(close(c.asinh().sinh(), c));
  1064. // cosh(acosh(z)) = z
  1065. assert!(close(c.acosh().cosh(), c));
  1066. // tanh(atanh(z)) = z
  1067. // 1 and -1 are branch points
  1068. if c != _1_0i && c != _1_0i.scale(-1.0) {
  1069. assert!(close(c.atanh().tanh(), c));
  1070. }
  1071. // sinh(z) = (e^z - e^(-z))/2
  1072. assert!(close((c.exp() - c.exp().inv()).unscale(2.0), c.sinh()));
  1073. // cosh(z) = (e^z + e^(-z))/2
  1074. assert!(close((c.exp() + c.exp().inv()).unscale(2.0), c.cosh()));
  1075. // tanh(z) = ( e^(2z) - 1)/(e^(2z) + 1)
  1076. assert!(close((c.scale(2.0).exp() - _1_0i)/(c.scale(2.0).exp() + _1_0i), c.tanh()));
  1077. }
  1078. }
  1079. mod complex_arithmetic {
  1080. use super::{_0_0i, _1_0i, _1_1i, _0_1i, _neg1_1i, _05_05i, all_consts};
  1081. use traits::Zero;
  1082. #[test]
  1083. fn test_add() {
  1084. assert_eq!(_05_05i + _05_05i, _1_1i);
  1085. assert_eq!(_0_1i + _1_0i, _1_1i);
  1086. assert_eq!(_1_0i + _neg1_1i, _0_1i);
  1087. for &c in all_consts.iter() {
  1088. assert_eq!(_0_0i + c, c);
  1089. assert_eq!(c + _0_0i, c);
  1090. }
  1091. }
  1092. #[test]
  1093. fn test_sub() {
  1094. assert_eq!(_05_05i - _05_05i, _0_0i);
  1095. assert_eq!(_0_1i - _1_0i, _neg1_1i);
  1096. assert_eq!(_0_1i - _neg1_1i, _1_0i);
  1097. for &c in all_consts.iter() {
  1098. assert_eq!(c - _0_0i, c);
  1099. assert_eq!(c - c, _0_0i);
  1100. }
  1101. }
  1102. #[test]
  1103. fn test_mul() {
  1104. assert_eq!(_05_05i * _05_05i, _0_1i.unscale(2.0));
  1105. assert_eq!(_1_1i * _0_1i, _neg1_1i);
  1106. // i^2 & i^4
  1107. assert_eq!(_0_1i * _0_1i, -_1_0i);
  1108. assert_eq!(_0_1i * _0_1i * _0_1i * _0_1i, _1_0i);
  1109. for &c in all_consts.iter() {
  1110. assert_eq!(c * _1_0i, c);
  1111. assert_eq!(_1_0i * c, c);
  1112. }
  1113. }
  1114. #[test]
  1115. fn test_div() {
  1116. assert_eq!(_neg1_1i / _0_1i, _1_1i);
  1117. for &c in all_consts.iter() {
  1118. if c != Zero::zero() {
  1119. assert_eq!(c / c, _1_0i);
  1120. }
  1121. }
  1122. }
  1123. #[test]
  1124. fn test_neg() {
  1125. assert_eq!(-_1_0i + _0_1i, _neg1_1i);
  1126. assert_eq!((-_0_1i) * _0_1i, _1_0i);
  1127. for &c in all_consts.iter() {
  1128. assert_eq!(-(-c), c);
  1129. }
  1130. }
  1131. }
  1132. mod real_arithmetic {
  1133. use super::super::Complex;
  1134. #[test]
  1135. fn test_add() {
  1136. assert_eq!(Complex::new(4.0, 2.0) + 0.5, Complex::new(4.5, 2.0));
  1137. assert_eq!(0.5 + Complex::new(4.0, 2.0), Complex::new(4.5, 2.0));
  1138. }
  1139. #[test]
  1140. fn test_sub() {
  1141. assert_eq!(Complex::new(4.0, 2.0) - 0.5, Complex::new(3.5, 2.0));
  1142. assert_eq!(0.5 - Complex::new(4.0, 2.0), Complex::new(-3.5, -2.0));
  1143. }
  1144. #[test]
  1145. fn test_mul() {
  1146. assert_eq!(Complex::new(4.0, 2.0) * 0.5, Complex::new(2.0, 1.0));
  1147. assert_eq!(0.5 * Complex::new(4.0, 2.0), Complex::new(2.0, 1.0));
  1148. }
  1149. #[test]
  1150. fn test_div() {
  1151. assert_eq!(Complex::new(4.0, 2.0) / 0.5, Complex::new(8.0, 4.0));
  1152. assert_eq!(0.5 / Complex::new(4.0, 2.0), Complex::new(0.1, -0.05));
  1153. }
  1154. }
  1155. #[test]
  1156. fn test_to_string() {
  1157. fn test(c : Complex64, s: String) {
  1158. assert_eq!(c.to_string(), s);
  1159. }
  1160. test(_0_0i, "0+0i".to_string());
  1161. test(_1_0i, "1+0i".to_string());
  1162. test(_0_1i, "0+1i".to_string());
  1163. test(_1_1i, "1+1i".to_string());
  1164. test(_neg1_1i, "-1+1i".to_string());
  1165. test(-_neg1_1i, "1-1i".to_string());
  1166. test(_05_05i, "0.5+0.5i".to_string());
  1167. }
  1168. #[test]
  1169. fn test_string_formatting() {
  1170. let a = Complex::new(1.23456, 123.456);
  1171. assert_eq!(format!("{}", a), "1.23456+123.456i");
  1172. assert_eq!(format!("{:.2}", a), "1.23+123.46i");
  1173. assert_eq!(format!("{:.2e}", a), "1.23e0+1.23e2i");
  1174. assert_eq!(format!("{:+20.2E}", a), " +1.23E0+1.23E2i");
  1175. let b = Complex::new(0x80, 0xff);
  1176. assert_eq!(format!("{:X}", b), "80+FFi");
  1177. assert_eq!(format!("{:#x}", b), "0x80+0xffi");
  1178. assert_eq!(format!("{:+#b}", b), "+0b10000000+0b11111111i");
  1179. assert_eq!(format!("{:+#16o}", b), " +0o200+0o377i");
  1180. let c = Complex::new(-10, -10000);
  1181. assert_eq!(format!("{}", c), "-10-10000i");
  1182. assert_eq!(format!("{:16}", c), " -10-10000i");
  1183. }
  1184. #[test]
  1185. fn test_hash() {
  1186. let a = Complex::new(0i32, 0i32);
  1187. let b = Complex::new(1i32, 0i32);
  1188. let c = Complex::new(0i32, 1i32);
  1189. assert!(::hash(&a) != ::hash(&b));
  1190. assert!(::hash(&b) != ::hash(&c));
  1191. assert!(::hash(&c) != ::hash(&a));
  1192. }
  1193. #[test]
  1194. fn test_hashset() {
  1195. use std::collections::HashSet;
  1196. let a = Complex::new(0i32, 0i32);
  1197. let b = Complex::new(1i32, 0i32);
  1198. let c = Complex::new(0i32, 1i32);
  1199. let set: HashSet<_> = [a, b, c].iter().cloned().collect();
  1200. assert!(set.contains(&a));
  1201. assert!(set.contains(&b));
  1202. assert!(set.contains(&c));
  1203. assert!(!set.contains(&(a + b + c)));
  1204. }
  1205. #[test]
  1206. fn test_is_nan() {
  1207. assert!(!_1_1i.is_nan());
  1208. let a = Complex::new(f64::NAN, f64::NAN);
  1209. assert!(a.is_nan());
  1210. }
  1211. #[test]
  1212. fn test_is_nan_special_cases() {
  1213. let a = Complex::new(0f64, f64::NAN);
  1214. let b = Complex::new(f64::NAN, 0f64);
  1215. assert!(a.is_nan());
  1216. assert!(b.is_nan());
  1217. }
  1218. #[test]
  1219. fn test_is_infinite() {
  1220. let a = Complex::new(2f64, f64::INFINITY);
  1221. assert!(a.is_infinite());
  1222. }
  1223. #[test]
  1224. fn test_is_finite() {
  1225. assert!(_1_1i.is_finite())
  1226. }
  1227. #[test]
  1228. fn test_is_normal() {
  1229. let a = Complex::new(0f64, f64::NAN);
  1230. let b = Complex::new(2f64, f64::INFINITY);
  1231. assert!(!a.is_normal());
  1232. assert!(!b.is_normal());
  1233. assert!(_1_1i.is_normal());
  1234. }
  1235. }