From 9b1489411c1873a6fa2af686311519e125158214 Mon Sep 17 00:00:00 2001 From: DongHun Kwak Date: Tue, 16 May 2023 16:00:01 +0900 Subject: [PATCH 1/1] Import num-complex 0.4.3 --- .cargo_vcs_info.json | 6 + .gitignore | 2 + Cargo.toml | 81 ++ Cargo.toml.orig | 53 + LICENSE-APACHE | 201 ++++ LICENSE-MIT | 25 + README.md | 55 + RELEASES.md | 167 +++ src/cast.rs | 119 +++ src/complex_float.rs | 439 ++++++++ src/crand.rs | 148 +++ src/lib.rs | 2887 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/pow.rs | 186 ++++ 13 files changed, 4369 insertions(+) create mode 100644 .cargo_vcs_info.json create mode 100644 .gitignore create mode 100644 Cargo.toml create mode 100644 Cargo.toml.orig create mode 100644 LICENSE-APACHE create mode 100644 LICENSE-MIT create mode 100644 README.md create mode 100644 RELEASES.md create mode 100644 src/cast.rs create mode 100644 src/complex_float.rs create mode 100644 src/crand.rs create mode 100644 src/lib.rs create mode 100644 src/pow.rs diff --git a/.cargo_vcs_info.json b/.cargo_vcs_info.json new file mode 100644 index 0000000..f8cd311 --- /dev/null +++ b/.cargo_vcs_info.json @@ -0,0 +1,6 @@ +{ + "git": { + "sha1": "23ccbd9362e3b3574b9a8253f67603ae4919f484" + }, + "path_in_vcs": "" +} \ No newline at end of file diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fa8d85a --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +Cargo.lock +target diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..690acba --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,81 @@ +# THIS FILE IS AUTOMATICALLY GENERATED BY CARGO +# +# When uploading crates to the registry Cargo will automatically +# "normalize" Cargo.toml files for maximal compatibility +# with all versions of Cargo and also rewrite `path` dependencies +# to registry (e.g., crates.io) dependencies. +# +# If you are reading this file be aware that the original Cargo.toml +# will likely look very different (and much more reasonable). +# See Cargo.toml.orig for the original contents. + +[package] +edition = "2018" +name = "num-complex" +version = "0.4.3" +authors = ["The Rust Project Developers"] +exclude = [ + "/bors.toml", + "/ci/*", + "/.github/*", +] +description = "Complex numbers implementation for Rust" +homepage = "https://github.com/rust-num/num-complex" +documentation = "https://docs.rs/num-complex" +readme = "README.md" +keywords = [ + "mathematics", + "numerics", +] +categories = [ + "algorithms", + "data-structures", + "science", + "no-std", +] +license = "MIT OR Apache-2.0" +repository = "https://github.com/rust-num/num-complex" + +[package.metadata.docs.rs] +features = [ + "bytemuck", + "std", + "serde", + "rkyv/size_64", + "bytecheck", + "rand", +] + +[dependencies.bytecheck] +version = "0.6" +optional = true +default-features = false + +[dependencies.bytemuck] +version = "1" +optional = true + +[dependencies.num-traits] +version = "0.2.11" +features = ["i128"] +default-features = false + +[dependencies.rand] +version = "0.8" +optional = true +default-features = false + +[dependencies.rkyv] +version = "0.7" +optional = true +default-features = false + +[dependencies.serde] +version = "1.0" +optional = true +default-features = false + +[features] +default = ["std"] +libm = ["num-traits/libm"] +std = ["num-traits/std"] diff --git a/Cargo.toml.orig b/Cargo.toml.orig new file mode 100644 index 0000000..149072c --- /dev/null +++ b/Cargo.toml.orig @@ -0,0 +1,53 @@ +[package] +authors = ["The Rust Project Developers"] +description = "Complex numbers implementation for Rust" +documentation = "https://docs.rs/num-complex" +homepage = "https://github.com/rust-num/num-complex" +keywords = ["mathematics", "numerics"] +categories = ["algorithms", "data-structures", "science", "no-std"] +license = "MIT OR Apache-2.0" +name = "num-complex" +repository = "https://github.com/rust-num/num-complex" +version = "0.4.3" +readme = "README.md" +exclude = ["/bors.toml", "/ci/*", "/.github/*"] +edition = "2018" + +[package.metadata.docs.rs] +features = ["bytemuck", "std", "serde", "rkyv/size_64", "bytecheck", "rand"] + +[dependencies] + +[dependencies.bytemuck] +optional = true +version = "1" + +[dependencies.num-traits] +version = "0.2.11" +default-features = false +features = ["i128"] + +[dependencies.serde] +optional = true +version = "1.0" +default-features = false + +[dependencies.rkyv] +optional = true +version = "0.7" +default-features = false + +[dependencies.bytecheck] +optional = true +version = "0.6" +default-features = false + +[dependencies.rand] +optional = true +version = "0.8" +default-features = false + +[features] +default = ["std"] +std = ["num-traits/std"] +libm = ["num-traits/libm"] diff --git a/LICENSE-APACHE b/LICENSE-APACHE new file mode 100644 index 0000000..16fe87b --- /dev/null +++ b/LICENSE-APACHE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + +TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + +1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + +2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + +3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + +4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + +5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + +6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + +7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + +8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + +9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + +END OF TERMS AND CONDITIONS + +APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + +Copyright [yyyy] [name of copyright owner] + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/LICENSE-MIT b/LICENSE-MIT new file mode 100644 index 0000000..39d4bdb --- /dev/null +++ b/LICENSE-MIT @@ -0,0 +1,25 @@ +Copyright (c) 2014 The Rust Project Developers + +Permission is hereby granted, free of charge, to any +person obtaining a copy of this software and associated +documentation files (the "Software"), to deal in the +Software without restriction, including without +limitation the rights to use, copy, modify, merge, +publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software +is furnished to do so, subject to the following +conditions: + +The above copyright notice and this permission notice +shall be included in all copies or substantial portions +of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF +ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED +TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A +PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT +SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR +IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..fdb1d38 --- /dev/null +++ b/README.md @@ -0,0 +1,55 @@ +# num-complex + +[![crate](https://img.shields.io/crates/v/num-complex.svg)](https://crates.io/crates/num-complex) +[![documentation](https://docs.rs/num-complex/badge.svg)](https://docs.rs/num-complex) +[![minimum rustc 1.31](https://img.shields.io/badge/rustc-1.31+-red.svg)](https://rust-lang.github.io/rfcs/2495-min-rust-version.html) +[![build status](https://github.com/rust-num/num-complex/workflows/master/badge.svg)](https://github.com/rust-num/num-complex/actions) + +`Complex` numbers for Rust. + +## Usage + +Add this to your `Cargo.toml`: + +```toml +[dependencies] +num-complex = "0.4" +``` + +## Features + +This crate can be used without the standard library (`#![no_std]`) by disabling +the default `std` feature. Use this in `Cargo.toml`: + +```toml +[dependencies.num-complex] +version = "0.4" +default-features = false +``` + +Features based on `Float` types are only available when `std` or `libm` is +enabled. Where possible, `FloatCore` is used instead. Formatting complex +numbers only supports format width when `std` is enabled. + +## Releases + +Release notes are available in [RELEASES.md](RELEASES.md). + +## Compatibility + +The `num-complex` crate is tested for rustc 1.31 and greater. + +## License + +Licensed under either of + + * [Apache License, Version 2.0](http://www.apache.org/licenses/LICENSE-2.0) + * [MIT license](http://opensource.org/licenses/MIT) + +at your option. + +### Contribution + +Unless you explicitly state otherwise, any contribution intentionally submitted +for inclusion in the work by you, as defined in the Apache-2.0 license, shall be +dual licensed as above, without any additional terms or conditions. diff --git a/RELEASES.md b/RELEASES.md new file mode 100644 index 0000000..d6244d8 --- /dev/null +++ b/RELEASES.md @@ -0,0 +1,167 @@ +# Release 0.4.3 (2023-01-19) + +- [`Complex` now optionally supports `bytecheck` 0.6 and `rkyv` 0.7][110]. + +**Contributors**: @cuviper, @zyansheep + +# Release 0.4.2 (2022-06-17) + +- [The new `ComplexFloat` trait][95] provides a generic abstraction between + floating-point `T` and `Complex`. +- [`Complex::exp` now handles edge cases with NaN and infinite parts][104]. + +**Contributors**: @cuviper, @JorisDeRidder, @obsgolem, @YakoYakoYokuYoku + +[95]: https://github.com/rust-num/num-complex/pull/95 +[104]: https://github.com/rust-num/num-complex/pull/104 + +# Release 0.4.1 (2022-04-29) + +- [`Complex::from_str_radix` now returns an error for radix > 18][90], because + 'i' and 'j' as digits are ambiguous with _i_ or _j_ imaginary parts. +- [`Complex` now implements `bytemuck` traits when `T` does][100]. +- [`Complex::cis` creates a complex with the given phase][101], _e__i_ θ. + +**Contributors**: @bluss, @bradleyharden, @cuviper, @rayhem + +[90]: https://github.com/rust-num/num-complex/pull/90 +[100]: https://github.com/rust-num/num-complex/pull/100 +[101]: https://github.com/rust-num/num-complex/pull/101 + +# Release 0.4.0 (2021-03-05) + +- `rand` support has been updated to 0.8, requiring Rust 1.36. + +**Contributors**: @cuviper + +# Release 0.3.1 (2020-10-29) + +- Clarify the license specification as "MIT OR Apache-2.0". + +**Contributors**: @cuviper + +# Release 0.3.0 (2020-06-13) + +### Enhancements + +- [The new "libm" feature passes through to `num-traits`][73], enabling `Float` + features on no-`std` builds. + +### Breaking Changes + +- `num-complex` now requires Rust 1.31 or greater. + - The "i128" opt-in feature was removed, now always available. +- [Updated public dependences][65]: + - `rand` support has been updated to 0.7, requiring Rust 1.32. +- [Methods for `T: Float` now take values instead of references][82], most + notably affecting the constructor `from_polar`. + +**Contributors**: @cuviper, @SOF3, @vks + +[65]: https://github.com/rust-num/num-complex/pull/65 +[73]: https://github.com/rust-num/num-complex/pull/73 +[82]: https://github.com/rust-num/num-complex/pull/82 + +# Release 0.2.4 (2020-01-09) + +- [`Complex::new` is now a `const fn` for Rust 1.31 and later][63]. +- [Updated the `autocfg` build dependency to 1.0][68]. + +**Contributors**: @burrbull, @cuviper, @dingelish + +[63]: https://github.com/rust-num/num-complex/pull/63 +[68]: https://github.com/rust-num/num-complex/pull/68 + +# Release 0.2.3 (2019-06-11) + +- [`Complex::sqrt()` is now more accurate for negative reals][60]. +- [`Complex::cbrt()` computes the principal cube root][61]. + +**Contributors**: @cuviper + +[60]: https://github.com/rust-num/num-complex/pull/60 +[61]: https://github.com/rust-num/num-complex/pull/61 + +# Release 0.2.2 (2019-06-10) + +- [`Complex::l1_norm()` computes the Manhattan distance from the origin][43]. +- [`Complex::fdiv()` and `finv()` use floating-point for inversion][41], which + may avoid overflows for some inputs, at the cost of trigonometric rounding. +- [`Complex` now implements `num_traits::MulAdd` and `MulAddAssign`][44]. +- [`Complex` now implements `Zero::set_zero` and `One::set_one`][57]. +- [`Complex` now implements `num_traits::Pow` and adds `powi` and `powu`][56]. + +**Contributors**: @adamnemecek, @cuviper, @ignatenkobrain, @Schultzer + +[41]: https://github.com/rust-num/num-complex/pull/41 +[43]: https://github.com/rust-num/num-complex/pull/43 +[44]: https://github.com/rust-num/num-complex/pull/44 +[56]: https://github.com/rust-num/num-complex/pull/56 +[57]: https://github.com/rust-num/num-complex/pull/57 + +# Release 0.2.1 (2018-10-08) + +- [`Complex` now implements `ToPrimitive`, `FromPrimitive`, `AsPrimitive`, and `NumCast`][33]. + +**Contributors**: @cuviper, @termoshtt + +[33]: https://github.com/rust-num/num-complex/pull/33 + +# Release 0.2.0 (2018-05-24) + +### Enhancements + +- [`Complex` now implements `num_traits::Inv` and `One::is_one`][17]. +- [`Complex` now implements `Sum` and `Product`][11]. +- [`Complex` now supports `i128` and `u128` components][27] with Rust 1.26+. +- [`Complex` now optionally supports `rand` 0.5][28], implementing the + `Standard` distribution and [a generic `ComplexDistribution`][30]. +- [`Rem` with a scalar divisor now avoids `norm_sqr` overflow][25]. + +### Breaking Changes + +- [`num-complex` now requires rustc 1.15 or greater][16]. +- [There is now a `std` feature][22], enabled by default, along with the + implication that building *without* this feature makes this a `#![no_std]` + crate. A few methods now require `FloatCore`, and the remaining methods + based on `Float` are only supported with `std`. +- [The `serde` dependency has been updated to 1.0][7], and `rustc-serialize` + is no longer supported by `num-complex`. + +**Contributors**: @clarcharr, @cuviper, @shingtaklam1324, @termoshtt + +[7]: https://github.com/rust-num/num-complex/pull/7 +[11]: https://github.com/rust-num/num-complex/pull/11 +[16]: https://github.com/rust-num/num-complex/pull/16 +[17]: https://github.com/rust-num/num-complex/pull/17 +[22]: https://github.com/rust-num/num-complex/pull/22 +[25]: https://github.com/rust-num/num-complex/pull/25 +[27]: https://github.com/rust-num/num-complex/pull/27 +[28]: https://github.com/rust-num/num-complex/pull/28 +[30]: https://github.com/rust-num/num-complex/pull/30 + + +# Release 0.1.43 (2018-03-08) + +- [Fix a usage typo in README.md][20]. + +**Contributors**: @shingtaklam1324 + +[20]: https://github.com/rust-num/num-complex/pull/20 + + +# Release 0.1.42 (2018-02-07) + +- [num-complex now has its own source repository][num-356] at [rust-num/num-complex][home]. + +**Contributors**: @cuviper + +[home]: https://github.com/rust-num/num-complex +[num-356]: https://github.com/rust-num/num/pull/356 + + +# Prior releases + +No prior release notes were kept. Thanks all the same to the many +contributors that have made this crate what it is! + diff --git a/src/cast.rs b/src/cast.rs new file mode 100644 index 0000000..e12f5cd --- /dev/null +++ b/src/cast.rs @@ -0,0 +1,119 @@ +use super::Complex; +use num_traits::{AsPrimitive, FromPrimitive, Num, NumCast, ToPrimitive}; + +macro_rules! impl_to_primitive { + ($ty:ty, $to:ident) => { + #[inline] + fn $to(&self) -> Option<$ty> { + if self.im.is_zero() { + self.re.$to() + } else { + None + } + } + }; +} // impl_to_primitive + +// Returns None if Complex part is non-zero +impl ToPrimitive for Complex { + impl_to_primitive!(usize, to_usize); + impl_to_primitive!(isize, to_isize); + impl_to_primitive!(u8, to_u8); + impl_to_primitive!(u16, to_u16); + impl_to_primitive!(u32, to_u32); + impl_to_primitive!(u64, to_u64); + impl_to_primitive!(i8, to_i8); + impl_to_primitive!(i16, to_i16); + impl_to_primitive!(i32, to_i32); + impl_to_primitive!(i64, to_i64); + impl_to_primitive!(u128, to_u128); + impl_to_primitive!(i128, to_i128); + impl_to_primitive!(f32, to_f32); + impl_to_primitive!(f64, to_f64); +} + +macro_rules! impl_from_primitive { + ($ty:ty, $from_xx:ident) => { + #[inline] + fn $from_xx(n: $ty) -> Option { + Some(Complex { + re: T::$from_xx(n)?, + im: T::zero(), + }) + } + }; +} // impl_from_primitive + +impl FromPrimitive for Complex { + impl_from_primitive!(usize, from_usize); + impl_from_primitive!(isize, from_isize); + impl_from_primitive!(u8, from_u8); + impl_from_primitive!(u16, from_u16); + impl_from_primitive!(u32, from_u32); + impl_from_primitive!(u64, from_u64); + impl_from_primitive!(i8, from_i8); + impl_from_primitive!(i16, from_i16); + impl_from_primitive!(i32, from_i32); + impl_from_primitive!(i64, from_i64); + impl_from_primitive!(u128, from_u128); + impl_from_primitive!(i128, from_i128); + impl_from_primitive!(f32, from_f32); + impl_from_primitive!(f64, from_f64); +} + +impl NumCast for Complex { + fn from(n: U) -> Option { + Some(Complex { + re: T::from(n)?, + im: T::zero(), + }) + } +} + +impl AsPrimitive for Complex +where + T: AsPrimitive, + U: 'static + Copy, +{ + fn as_(self) -> U { + self.re.as_() + } +} + +#[cfg(test)] +mod test { + use super::*; + + #[test] + fn test_to_primitive() { + let a: Complex = Complex { re: 3, im: 0 }; + assert_eq!(a.to_i32(), Some(3_i32)); + let b: Complex = Complex { re: 3, im: 1 }; + assert_eq!(b.to_i32(), None); + let x: Complex = Complex { re: 1.0, im: 0.1 }; + assert_eq!(x.to_f32(), None); + let y: Complex = Complex { re: 1.0, im: 0.0 }; + assert_eq!(y.to_f32(), Some(1.0)); + let z: Complex = Complex { re: 1.0, im: 0.0 }; + assert_eq!(z.to_i32(), Some(1)); + } + + #[test] + fn test_from_primitive() { + let a: Complex = FromPrimitive::from_i32(2).unwrap(); + assert_eq!(a, Complex { re: 2.0, im: 0.0 }); + } + + #[test] + fn test_num_cast() { + let a: Complex = NumCast::from(2_i32).unwrap(); + assert_eq!(a, Complex { re: 2.0, im: 0.0 }); + } + + #[test] + fn test_as_primitive() { + let a: Complex = Complex { re: 2.0, im: 0.2 }; + let a_: i32 = a.as_(); + assert_eq!(a_, 2_i32); + } +} diff --git a/src/complex_float.rs b/src/complex_float.rs new file mode 100644 index 0000000..873fe73 --- /dev/null +++ b/src/complex_float.rs @@ -0,0 +1,439 @@ +// Keeps us from accidentally creating a recursive impl rather than a real one. +#![deny(unconditional_recursion)] + +use core::ops::Neg; + +use num_traits::{Float, FloatConst, Num, NumCast}; + +use crate::Complex; + +mod private { + use num_traits::{Float, FloatConst}; + + use crate::Complex; + + pub trait Seal {} + + impl Seal for T where T: Float + FloatConst {} + impl Seal for Complex {} +} + +/// Generic trait for floating point complex numbers. +/// +/// This trait defines methods which are common to complex floating point +/// numbers and regular floating point numbers. +/// +/// This trait is sealed to prevent it from being implemented by anything other +/// than floating point scalars and [Complex] floats. +pub trait ComplexFloat: Num + NumCast + Copy + Neg + private::Seal { + /// The type used to represent the real coefficients of this complex number. + type Real: Float + FloatConst; + + /// Returns `true` if this value is `NaN` and false otherwise. + fn is_nan(self) -> bool; + + /// Returns `true` if this value is positive infinity or negative infinity and + /// false otherwise. + fn is_infinite(self) -> bool; + + /// Returns `true` if this number is neither infinite nor `NaN`. + fn is_finite(self) -> bool; + + /// Returns `true` if the number is neither zero, infinite, + /// [subnormal](http://en.wikipedia.org/wiki/Denormal_number), or `NaN`. + fn is_normal(self) -> bool; + + /// Take the reciprocal (inverse) of a number, `1/x`. See also [Complex::finv]. + fn recip(self) -> Self; + + /// Raises `self` to a signed integer power. + fn powi(self, exp: i32) -> Self; + + /// Raises `self` to a real power. + fn powf(self, exp: Self::Real) -> Self; + + /// Raises `self` to a complex power. + fn powc(self, exp: Complex) -> Complex; + + /// Take the square root of a number. + fn sqrt(self) -> Self; + + /// Returns `e^(self)`, (the exponential function). + fn exp(self) -> Self; + + /// Returns `2^(self)`. + fn exp2(self) -> Self; + + /// Returns `base^(self)`. + fn expf(self, base: Self::Real) -> Self; + + /// Returns the natural logarithm of the number. + fn ln(self) -> Self; + + /// Returns the logarithm of the number with respect to an arbitrary base. + fn log(self, base: Self::Real) -> Self; + + /// Returns the base 2 logarithm of the number. + fn log2(self) -> Self; + + /// Returns the base 10 logarithm of the number. + fn log10(self) -> Self; + + /// Take the cubic root of a number. + fn cbrt(self) -> Self; + + /// Computes the sine of a number (in radians). + fn sin(self) -> Self; + + /// Computes the cosine of a number (in radians). + fn cos(self) -> Self; + + /// Computes the tangent of a number (in radians). + fn tan(self) -> Self; + + /// Computes the arcsine of a number. Return value is in radians in + /// the range [-pi/2, pi/2] or NaN if the number is outside the range + /// [-1, 1]. + fn asin(self) -> Self; + + /// Computes the arccosine of a number. Return value is in radians in + /// the range [0, pi] or NaN if the number is outside the range + /// [-1, 1]. + fn acos(self) -> Self; + + /// Computes the arctangent of a number. Return value is in radians in the + /// range [-pi/2, pi/2]; + fn atan(self) -> Self; + + /// Hyperbolic sine function. + fn sinh(self) -> Self; + + /// Hyperbolic cosine function. + fn cosh(self) -> Self; + + /// Hyperbolic tangent function. + fn tanh(self) -> Self; + + /// Inverse hyperbolic sine function. + fn asinh(self) -> Self; + + /// Inverse hyperbolic cosine function. + fn acosh(self) -> Self; + + /// Inverse hyperbolic tangent function. + fn atanh(self) -> Self; + + /// Returns the real part of the number. + fn re(self) -> Self::Real; + + /// Returns the imaginary part of the number. + fn im(self) -> Self::Real; + + /// Returns the absolute value of the number. See also [Complex::norm] + fn abs(self) -> Self::Real; + + /// Returns the L1 norm `|re| + |im|` -- the [Manhattan distance] from the origin. + /// + /// [Manhattan distance]: https://en.wikipedia.org/wiki/Taxicab_geometry + fn l1_norm(&self) -> Self::Real; + + /// Computes the argument of the number. + fn arg(self) -> Self::Real; + + /// Computes the complex conjugate of the number. + /// + /// Formula: `a+bi -> a-bi` + fn conj(self) -> Self; +} + +macro_rules! forward { + ($( $base:ident :: $method:ident ( self $( , $arg:ident : $ty:ty )* ) -> $ret:ty ; )*) + => {$( + #[inline] + fn $method(self $( , $arg : $ty )* ) -> $ret { + $base::$method(self $( , $arg )* ) + } + )*}; +} + +macro_rules! forward_ref { + ($( Self :: $method:ident ( & self $( , $arg:ident : $ty:ty )* ) -> $ret:ty ; )*) + => {$( + #[inline] + fn $method(self $( , $arg : $ty )* ) -> $ret { + Self::$method(&self $( , $arg )* ) + } + )*}; +} + +impl ComplexFloat for T +where + T: Float + FloatConst, +{ + type Real = T; + + fn re(self) -> Self::Real { + self + } + + fn im(self) -> Self::Real { + T::zero() + } + + fn l1_norm(&self) -> Self::Real { + self.abs() + } + + fn arg(self) -> Self::Real { + if self.is_nan() { + self + } else if self.is_sign_negative() { + T::PI() + } else { + T::zero() + } + } + + fn powc(self, exp: Complex) -> Complex { + Complex::new(self, T::zero()).powc(exp) + } + + fn conj(self) -> Self { + self + } + + fn expf(self, base: Self::Real) -> Self { + base.powf(self) + } + + forward! { + Float::is_normal(self) -> bool; + Float::is_infinite(self) -> bool; + Float::is_finite(self) -> bool; + Float::is_nan(self) -> bool; + Float::recip(self) -> Self; + Float::powi(self, n: i32) -> Self; + Float::powf(self, f: Self) -> Self; + Float::sqrt(self) -> Self; + Float::cbrt(self) -> Self; + Float::exp(self) -> Self; + Float::exp2(self) -> Self; + Float::ln(self) -> Self; + Float::log(self, base: Self) -> Self; + Float::log2(self) -> Self; + Float::log10(self) -> Self; + Float::sin(self) -> Self; + Float::cos(self) -> Self; + Float::tan(self) -> Self; + Float::asin(self) -> Self; + Float::acos(self) -> Self; + Float::atan(self) -> Self; + Float::sinh(self) -> Self; + Float::cosh(self) -> Self; + Float::tanh(self) -> Self; + Float::asinh(self) -> Self; + Float::acosh(self) -> Self; + Float::atanh(self) -> Self; + Float::abs(self) -> Self; + } +} + +impl ComplexFloat for Complex { + type Real = T; + + fn re(self) -> Self::Real { + self.re + } + + fn im(self) -> Self::Real { + self.im + } + + fn abs(self) -> Self::Real { + self.norm() + } + + fn recip(self) -> Self { + self.finv() + } + + // `Complex::l1_norm` uses `Signed::abs` to let it work + // for integers too, but we can just use `Float::abs`. + fn l1_norm(&self) -> Self::Real { + self.re.abs() + self.im.abs() + } + + // `Complex::is_*` methods use `T: FloatCore`, but we + // have `T: Float` that can do them as well. + fn is_nan(self) -> bool { + self.re.is_nan() || self.im.is_nan() + } + + fn is_infinite(self) -> bool { + !self.is_nan() && (self.re.is_infinite() || self.im.is_infinite()) + } + + fn is_finite(self) -> bool { + self.re.is_finite() && self.im.is_finite() + } + + fn is_normal(self) -> bool { + self.re.is_normal() && self.im.is_normal() + } + + forward! { + Complex::arg(self) -> Self::Real; + Complex::powc(self, exp: Complex) -> Complex; + Complex::exp2(self) -> Self; + Complex::log(self, base: Self::Real) -> Self; + Complex::log2(self) -> Self; + Complex::log10(self) -> Self; + Complex::powf(self, f: Self::Real) -> Self; + Complex::sqrt(self) -> Self; + Complex::cbrt(self) -> Self; + Complex::exp(self) -> Self; + Complex::expf(self, base: Self::Real) -> Self; + Complex::ln(self) -> Self; + Complex::sin(self) -> Self; + Complex::cos(self) -> Self; + Complex::tan(self) -> Self; + Complex::asin(self) -> Self; + Complex::acos(self) -> Self; + Complex::atan(self) -> Self; + Complex::sinh(self) -> Self; + Complex::cosh(self) -> Self; + Complex::tanh(self) -> Self; + Complex::asinh(self) -> Self; + Complex::acosh(self) -> Self; + Complex::atanh(self) -> Self; + } + + forward_ref! { + Self::powi(&self, n: i32) -> Self; + Self::conj(&self) -> Self; + } +} + +#[cfg(test)] +mod test { + use crate::{ + complex_float::ComplexFloat, + test::{_0_0i, _0_1i, _1_0i, _1_1i, float::close}, + Complex, + }; + use std::f64; // for constants before Rust 1.43. + + fn closef(a: f64, b: f64) -> bool { + close_to_tolf(a, b, 1e-10) + } + + fn close_to_tolf(a: f64, b: f64, tol: f64) -> bool { + // returns true if a and b are reasonably close + let close = (a == b) || (a - b).abs() < tol; + if !close { + println!("{:?} != {:?}", a, b); + } + close + } + + #[test] + fn test_exp2() { + assert!(close(ComplexFloat::exp2(_0_0i), _1_0i)); + assert!(closef(::exp2(0.), 1.)); + } + + #[test] + fn test_exp() { + assert!(close(ComplexFloat::exp(_0_0i), _1_0i)); + assert!(closef(ComplexFloat::exp(0.), 1.)); + } + + #[test] + fn test_powi() { + assert!(close(ComplexFloat::powi(_0_1i, 4), _1_0i)); + assert!(closef(ComplexFloat::powi(-1., 4), 1.)); + } + + #[test] + fn test_powz() { + assert!(close(ComplexFloat::powc(_1_0i, _0_1i), _1_0i)); + assert!(close(ComplexFloat::powc(1., _0_1i), _1_0i)); + } + + #[test] + fn test_log2() { + assert!(close(ComplexFloat::log2(_1_0i), _0_0i)); + assert!(closef(ComplexFloat::log2(1.), 0.)); + } + + #[test] + fn test_log10() { + assert!(close(ComplexFloat::log10(_1_0i), _0_0i)); + assert!(closef(ComplexFloat::log10(1.), 0.)); + } + + #[test] + fn test_conj() { + assert_eq!(ComplexFloat::conj(_0_1i), Complex::new(0., -1.)); + assert_eq!(ComplexFloat::conj(1.), 1.); + } + + #[test] + fn test_is_nan() { + assert!(!ComplexFloat::is_nan(_1_0i)); + assert!(!ComplexFloat::is_nan(1.)); + + assert!(ComplexFloat::is_nan(Complex::new(f64::NAN, f64::NAN))); + assert!(ComplexFloat::is_nan(f64::NAN)); + } + + #[test] + fn test_is_infinite() { + assert!(!ComplexFloat::is_infinite(_1_0i)); + assert!(!ComplexFloat::is_infinite(1.)); + + assert!(ComplexFloat::is_infinite(Complex::new( + f64::INFINITY, + f64::INFINITY + ))); + assert!(ComplexFloat::is_infinite(f64::INFINITY)); + } + + #[test] + fn test_is_finite() { + assert!(ComplexFloat::is_finite(_1_0i)); + assert!(ComplexFloat::is_finite(1.)); + + assert!(!ComplexFloat::is_finite(Complex::new( + f64::INFINITY, + f64::INFINITY + ))); + assert!(!ComplexFloat::is_finite(f64::INFINITY)); + } + + #[test] + fn test_is_normal() { + assert!(ComplexFloat::is_normal(_1_1i)); + assert!(ComplexFloat::is_normal(1.)); + + assert!(!ComplexFloat::is_normal(Complex::new( + f64::INFINITY, + f64::INFINITY + ))); + assert!(!ComplexFloat::is_normal(f64::INFINITY)); + } + + #[test] + fn test_arg() { + assert!(closef( + ComplexFloat::arg(_0_1i), + core::f64::consts::FRAC_PI_2 + )); + + assert!(closef(ComplexFloat::arg(-1.), core::f64::consts::PI)); + assert!(closef(ComplexFloat::arg(-0.), core::f64::consts::PI)); + assert!(closef(ComplexFloat::arg(0.), 0.)); + assert!(closef(ComplexFloat::arg(1.), 0.)); + assert!(ComplexFloat::arg(f64::NAN).is_nan()); + } +} diff --git a/src/crand.rs b/src/crand.rs new file mode 100644 index 0000000..5edf8a7 --- /dev/null +++ b/src/crand.rs @@ -0,0 +1,148 @@ +//! Rand implementations for complex numbers + +use crate::Complex; +use num_traits::Num; +use rand::distributions::Standard; +use rand::prelude::*; + +impl Distribution> for Standard +where + T: Num + Clone, + Standard: Distribution, +{ + fn sample(&self, rng: &mut R) -> Complex { + Complex::new(self.sample(rng), self.sample(rng)) + } +} + +/// A generic random value distribution for complex numbers. +#[derive(Clone, Copy, Debug)] +pub struct ComplexDistribution { + re: Re, + im: Im, +} + +impl ComplexDistribution { + /// Creates a complex distribution from independent + /// distributions of the real and imaginary parts. + pub fn new(re: Re, im: Im) -> Self { + ComplexDistribution { re, im } + } +} + +impl Distribution> for ComplexDistribution +where + T: Num + Clone, + Re: Distribution, + Im: Distribution, +{ + fn sample(&self, rng: &mut R) -> Complex { + Complex::new(self.re.sample(rng), self.im.sample(rng)) + } +} + +#[cfg(test)] +fn test_rng() -> impl RngCore { + /// Simple `Rng` for testing without additional dependencies + struct XorShiftStar { + a: u64, + } + + impl RngCore for XorShiftStar { + fn next_u32(&mut self) -> u32 { + self.next_u64() as u32 + } + + fn next_u64(&mut self) -> u64 { + // https://en.wikipedia.org/wiki/Xorshift#xorshift* + self.a ^= self.a >> 12; + self.a ^= self.a << 25; + self.a ^= self.a >> 27; + self.a.wrapping_mul(0x2545_F491_4F6C_DD1D) + } + + fn fill_bytes(&mut self, dest: &mut [u8]) { + for chunk in dest.chunks_mut(8) { + let bytes = self.next_u64().to_le_bytes(); + let slice = &bytes[..chunk.len()]; + chunk.copy_from_slice(slice) + } + } + + fn try_fill_bytes(&mut self, dest: &mut [u8]) -> Result<(), rand::Error> { + Ok(self.fill_bytes(dest)) + } + } + + XorShiftStar { + a: 0x0123_4567_89AB_CDEF, + } +} + +#[test] +fn standard_f64() { + let mut rng = test_rng(); + for _ in 0..100 { + let c: Complex = rng.gen(); + assert!(c.re >= 0.0 && c.re < 1.0); + assert!(c.im >= 0.0 && c.im < 1.0); + } +} + +#[test] +fn generic_standard_f64() { + let mut rng = test_rng(); + let dist = ComplexDistribution::new(Standard, Standard); + for _ in 0..100 { + let c: Complex = rng.sample(&dist); + assert!(c.re >= 0.0 && c.re < 1.0); + assert!(c.im >= 0.0 && c.im < 1.0); + } +} + +#[test] +fn generic_uniform_f64() { + use rand::distributions::Uniform; + + let mut rng = test_rng(); + let re = Uniform::new(-100.0, 0.0); + let im = Uniform::new(0.0, 100.0); + let dist = ComplexDistribution::new(re, im); + for _ in 0..100 { + // no type annotation required, since `Uniform` only produces one type. + let c = rng.sample(&dist); + assert!(c.re >= -100.0 && c.re < 0.0); + assert!(c.im >= 0.0 && c.im < 100.0); + } +} + +#[test] +fn generic_mixed_f64() { + use rand::distributions::Uniform; + + let mut rng = test_rng(); + let re = Uniform::new(-100.0, 0.0); + let dist = ComplexDistribution::new(re, Standard); + for _ in 0..100 { + // no type annotation required, since `Uniform` only produces one type. + let c = rng.sample(&dist); + assert!(c.re >= -100.0 && c.re < 0.0); + assert!(c.im >= 0.0 && c.im < 1.0); + } +} + +#[test] +fn generic_uniform_i32() { + use rand::distributions::Uniform; + + let mut rng = test_rng(); + let re = Uniform::new(-100, 0); + let im = Uniform::new(0, 100); + let dist = ComplexDistribution::new(re, im); + for _ in 0..100 { + // no type annotation required, since `Uniform` only produces one type. + let c = rng.sample(&dist); + assert!(c.re >= -100 && c.re < 0); + assert!(c.im >= 0 && c.im < 100); + } +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..1860d28 --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,2887 @@ +// Copyright 2013 The Rust Project Developers. See the COPYRIGHT +// file at the top-level directory of this distribution and at +// http://rust-lang.org/COPYRIGHT. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +//! Complex numbers. +//! +//! ## Compatibility +//! +//! The `num-complex` crate is tested for rustc 1.31 and greater. + +#![doc(html_root_url = "https://docs.rs/num-complex/0.4")] +#![no_std] + +#[cfg(any(test, feature = "std"))] +#[cfg_attr(test, macro_use)] +extern crate std; + +use core::fmt; +#[cfg(test)] +use core::hash; +use core::iter::{Product, Sum}; +use core::ops::{Add, Div, Mul, Neg, Rem, Sub}; +use core::str::FromStr; +#[cfg(feature = "std")] +use std::error::Error; + +use num_traits::{Inv, MulAdd, Num, One, Pow, Signed, Zero}; + +use num_traits::float::FloatCore; +#[cfg(any(feature = "std", feature = "libm"))] +use num_traits::float::{Float, FloatConst}; + +mod cast; +mod pow; + +#[cfg(any(feature = "std", feature = "libm"))] +mod complex_float; +#[cfg(any(feature = "std", feature = "libm"))] +pub use crate::complex_float::ComplexFloat; + +#[cfg(feature = "rand")] +mod crand; +#[cfg(feature = "rand")] +pub use crate::crand::ComplexDistribution; + +// FIXME #1284: handle complex NaN & infinity etc. This +// probably doesn't map to C's _Complex correctly. + +/// A complex number in Cartesian form. +/// +/// ## Representation and Foreign Function Interface Compatibility +/// +/// `Complex` is memory layout compatible with an array `[T; 2]`. +/// +/// Note that `Complex` where F is a floating point type is **only** memory +/// layout compatible with C's complex types, **not** necessarily calling +/// convention compatible. This means that for FFI you can only pass +/// `Complex` behind a pointer, not as a value. +/// +/// ## Examples +/// +/// Example of extern function declaration. +/// +/// ``` +/// use num_complex::Complex; +/// use std::os::raw::c_int; +/// +/// extern "C" { +/// fn zaxpy_(n: *const c_int, alpha: *const Complex, +/// x: *const Complex, incx: *const c_int, +/// y: *mut Complex, incy: *const c_int); +/// } +/// ``` +#[derive(PartialEq, Eq, Copy, Clone, Hash, Debug, Default)] +#[repr(C)] +#[cfg_attr( + feature = "rkyv", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] +#[cfg_attr(feature = "rkyv", archive(as = "Complex"))] +#[cfg_attr(feature = "bytecheck", derive(bytecheck::CheckBytes))] +pub struct Complex { + /// Real portion of the complex number + pub re: T, + /// Imaginary portion of the complex number + pub im: T, +} + +pub type Complex32 = Complex; +pub type Complex64 = Complex; + +impl Complex { + /// Create a new Complex + #[inline] + pub const fn new(re: T, im: T) -> Self { + Complex { re, im } + } +} + +impl Complex { + /// Returns imaginary unit + #[inline] + pub fn i() -> Self { + Self::new(T::zero(), T::one()) + } + + /// Returns the square of the norm (since `T` doesn't necessarily + /// have a sqrt function), i.e. `re^2 + im^2`. + #[inline] + pub fn norm_sqr(&self) -> T { + self.re.clone() * self.re.clone() + self.im.clone() * self.im.clone() + } + + /// Multiplies `self` by the scalar `t`. + #[inline] + pub fn scale(&self, t: T) -> Self { + Self::new(self.re.clone() * t.clone(), self.im.clone() * t) + } + + /// Divides `self` by the scalar `t`. + #[inline] + pub fn unscale(&self, t: T) -> Self { + Self::new(self.re.clone() / t.clone(), self.im.clone() / t) + } + + /// Raises `self` to an unsigned integer power. + #[inline] + pub fn powu(&self, exp: u32) -> Self { + Pow::pow(self, exp) + } +} + +impl> Complex { + /// Returns the complex conjugate. i.e. `re - i im` + #[inline] + pub fn conj(&self) -> Self { + Self::new(self.re.clone(), -self.im.clone()) + } + + /// Returns `1/self` + #[inline] + pub fn inv(&self) -> Self { + let norm_sqr = self.norm_sqr(); + Self::new( + self.re.clone() / norm_sqr.clone(), + -self.im.clone() / norm_sqr, + ) + } + + /// Raises `self` to a signed integer power. + #[inline] + pub fn powi(&self, exp: i32) -> Self { + Pow::pow(self, exp) + } +} + +impl Complex { + /// Returns the L1 norm `|re| + |im|` -- the [Manhattan distance] from the origin. + /// + /// [Manhattan distance]: https://en.wikipedia.org/wiki/Taxicab_geometry + #[inline] + pub fn l1_norm(&self) -> T { + self.re.abs() + self.im.abs() + } +} + +#[cfg(any(feature = "std", feature = "libm"))] +impl Complex { + /// Create a new Complex with a given phase: `exp(i * phase)`. + /// See [cis (mathematics)](https://en.wikipedia.org/wiki/Cis_(mathematics)). + #[inline] + pub fn cis(phase: T) -> Self { + Self::new(phase.cos(), phase.sin()) + } + + /// Calculate |self| + #[inline] + pub fn norm(self) -> T { + self.re.hypot(self.im) + } + /// Calculate the principal Arg of self. + #[inline] + pub fn arg(self) -> T { + self.im.atan2(self.re) + } + /// Convert to polar form (r, theta), such that + /// `self = r * exp(i * theta)` + #[inline] + pub fn to_polar(self) -> (T, T) { + (self.norm(), self.arg()) + } + /// Convert a polar representation into a complex number. + #[inline] + pub fn from_polar(r: T, theta: T) -> Self { + Self::new(r * theta.cos(), r * theta.sin()) + } + + /// Computes `e^(self)`, where `e` is the base of the natural logarithm. + #[inline] + pub fn exp(self) -> Self { + // formula: e^(a + bi) = e^a (cos(b) + i*sin(b)) = from_polar(e^a, b) + + let Complex { re, mut im } = self; + // Treat the corner cases +∞, -∞, and NaN + if re.is_infinite() { + if re < T::zero() { + if !im.is_finite() { + return Self::new(T::zero(), T::zero()); + } + } else { + if im == T::zero() || !im.is_finite() { + if im.is_infinite() { + im = T::nan(); + } + return Self::new(re, im); + } + } + } else if re.is_nan() && im == T::zero() { + return self; + } + + Self::from_polar(re.exp(), im) + } + + /// Computes the principal value of natural logarithm of `self`. + /// + /// This function has one branch cut: + /// + /// * `(-∞, 0]`, continuous from above. + /// + /// The branch satisfies `-π ≤ arg(ln(z)) ≤ π`. + #[inline] + pub fn ln(self) -> Self { + // formula: ln(z) = ln|z| + i*arg(z) + let (r, theta) = self.to_polar(); + Self::new(r.ln(), theta) + } + + /// Computes the principal value of the square root of `self`. + /// + /// This function has one branch cut: + /// + /// * `(-∞, 0)`, continuous from above. + /// + /// The branch satisfies `-π/2 ≤ arg(sqrt(z)) ≤ π/2`. + #[inline] + pub fn sqrt(self) -> Self { + if self.im.is_zero() { + if self.re.is_sign_positive() { + // simple positive real √r, and copy `im` for its sign + Self::new(self.re.sqrt(), self.im) + } else { + // √(r e^(iπ)) = √r e^(iπ/2) = i√r + // √(r e^(-iπ)) = √r e^(-iπ/2) = -i√r + let re = T::zero(); + let im = (-self.re).sqrt(); + if self.im.is_sign_positive() { + Self::new(re, im) + } else { + Self::new(re, -im) + } + } + } else if self.re.is_zero() { + // √(r e^(iπ/2)) = √r e^(iπ/4) = √(r/2) + i√(r/2) + // √(r e^(-iπ/2)) = √r e^(-iπ/4) = √(r/2) - i√(r/2) + let one = T::one(); + let two = one + one; + let x = (self.im.abs() / two).sqrt(); + if self.im.is_sign_positive() { + Self::new(x, x) + } else { + Self::new(x, -x) + } + } else { + // formula: sqrt(r e^(it)) = sqrt(r) e^(it/2) + let one = T::one(); + let two = one + one; + let (r, theta) = self.to_polar(); + Self::from_polar(r.sqrt(), theta / two) + } + } + + /// Computes the principal value of the cube root of `self`. + /// + /// This function has one branch cut: + /// + /// * `(-∞, 0)`, continuous from above. + /// + /// The branch satisfies `-π/3 ≤ arg(cbrt(z)) ≤ π/3`. + /// + /// Note that this does not match the usual result for the cube root of + /// negative real numbers. For example, the real cube root of `-8` is `-2`, + /// but the principal complex cube root of `-8` is `1 + i√3`. + #[inline] + pub fn cbrt(self) -> Self { + if self.im.is_zero() { + if self.re.is_sign_positive() { + // simple positive real ∛r, and copy `im` for its sign + Self::new(self.re.cbrt(), self.im) + } else { + // ∛(r e^(iπ)) = ∛r e^(iπ/3) = ∛r/2 + i∛r√3/2 + // ∛(r e^(-iπ)) = ∛r e^(-iπ/3) = ∛r/2 - i∛r√3/2 + let one = T::one(); + let two = one + one; + let three = two + one; + let re = (-self.re).cbrt() / two; + let im = three.sqrt() * re; + if self.im.is_sign_positive() { + Self::new(re, im) + } else { + Self::new(re, -im) + } + } + } else if self.re.is_zero() { + // ∛(r e^(iπ/2)) = ∛r e^(iπ/6) = ∛r√3/2 + i∛r/2 + // ∛(r e^(-iπ/2)) = ∛r e^(-iπ/6) = ∛r√3/2 - i∛r/2 + let one = T::one(); + let two = one + one; + let three = two + one; + let im = self.im.abs().cbrt() / two; + let re = three.sqrt() * im; + if self.im.is_sign_positive() { + Self::new(re, im) + } else { + Self::new(re, -im) + } + } else { + // formula: cbrt(r e^(it)) = cbrt(r) e^(it/3) + let one = T::one(); + let three = one + one + one; + let (r, theta) = self.to_polar(); + Self::from_polar(r.cbrt(), theta / three) + } + } + + /// Raises `self` to a floating point power. + #[inline] + pub fn powf(self, exp: T) -> Self { + // formula: x^y = (ρ e^(i θ))^y = ρ^y e^(i θ y) + // = from_polar(ρ^y, θ y) + let (r, theta) = self.to_polar(); + Self::from_polar(r.powf(exp), theta * exp) + } + + /// Returns the logarithm of `self` with respect to an arbitrary base. + #[inline] + pub fn log(self, base: T) -> Self { + // formula: log_y(x) = log_y(ρ e^(i θ)) + // = log_y(ρ) + log_y(e^(i θ)) = log_y(ρ) + ln(e^(i θ)) / ln(y) + // = log_y(ρ) + i θ / ln(y) + let (r, theta) = self.to_polar(); + Self::new(r.log(base), theta / base.ln()) + } + + /// Raises `self` to a complex power. + #[inline] + pub fn powc(self, exp: Self) -> Self { + // formula: x^y = (a + i b)^(c + i d) + // = (ρ e^(i θ))^c (ρ e^(i θ))^(i d) + // where ρ=|x| and θ=arg(x) + // = ρ^c e^(−d θ) e^(i c θ) ρ^(i d) + // = p^c e^(−d θ) (cos(c θ) + // + i sin(c θ)) (cos(d ln(ρ)) + i sin(d ln(ρ))) + // = p^c e^(−d θ) ( + // cos(c θ) cos(d ln(ρ)) − sin(c θ) sin(d ln(ρ)) + // + i(cos(c θ) sin(d ln(ρ)) + sin(c θ) cos(d ln(ρ)))) + // = p^c e^(−d θ) (cos(c θ + d ln(ρ)) + i sin(c θ + d ln(ρ))) + // = from_polar(p^c e^(−d θ), c θ + d ln(ρ)) + let (r, theta) = self.to_polar(); + Self::from_polar( + r.powf(exp.re) * (-exp.im * theta).exp(), + exp.re * theta + exp.im * r.ln(), + ) + } + + /// Raises a floating point number to the complex power `self`. + #[inline] + pub fn expf(self, base: T) -> Self { + // formula: x^(a+bi) = x^a x^bi = x^a e^(b ln(x) i) + // = from_polar(x^a, b ln(x)) + Self::from_polar(base.powf(self.re), self.im * base.ln()) + } + + /// Computes the sine of `self`. + #[inline] + pub fn sin(self) -> Self { + // formula: sin(a + bi) = sin(a)cosh(b) + i*cos(a)sinh(b) + Self::new( + self.re.sin() * self.im.cosh(), + self.re.cos() * self.im.sinh(), + ) + } + + /// Computes the cosine of `self`. + #[inline] + pub fn cos(self) -> Self { + // formula: cos(a + bi) = cos(a)cosh(b) - i*sin(a)sinh(b) + Self::new( + self.re.cos() * self.im.cosh(), + -self.re.sin() * self.im.sinh(), + ) + } + + /// Computes the tangent of `self`. + #[inline] + pub fn tan(self) -> Self { + // formula: tan(a + bi) = (sin(2a) + i*sinh(2b))/(cos(2a) + cosh(2b)) + let (two_re, two_im) = (self.re + self.re, self.im + self.im); + Self::new(two_re.sin(), two_im.sinh()).unscale(two_re.cos() + two_im.cosh()) + } + + /// Computes the principal value of the inverse sine of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞, -1)`, continuous from above. + /// * `(1, ∞)`, continuous from below. + /// + /// The branch satisfies `-π/2 ≤ Re(asin(z)) ≤ π/2`. + #[inline] + pub fn asin(self) -> Self { + // formula: arcsin(z) = -i ln(sqrt(1-z^2) + iz) + let i = Self::i(); + -i * ((Self::one() - self * self).sqrt() + i * self).ln() + } + + /// Computes the principal value of the inverse cosine of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞, -1)`, continuous from above. + /// * `(1, ∞)`, continuous from below. + /// + /// The branch satisfies `0 ≤ Re(acos(z)) ≤ π`. + #[inline] + pub fn acos(self) -> Self { + // formula: arccos(z) = -i ln(i sqrt(1-z^2) + z) + let i = Self::i(); + -i * (i * (Self::one() - self * self).sqrt() + self).ln() + } + + /// Computes the principal value of the inverse tangent of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞i, -i]`, continuous from the left. + /// * `[i, ∞i)`, continuous from the right. + /// + /// The branch satisfies `-π/2 ≤ Re(atan(z)) ≤ π/2`. + #[inline] + pub fn atan(self) -> Self { + // formula: arctan(z) = (ln(1+iz) - ln(1-iz))/(2i) + let i = Self::i(); + let one = Self::one(); + let two = one + one; + if self == i { + return Self::new(T::zero(), T::infinity()); + } else if self == -i { + return Self::new(T::zero(), -T::infinity()); + } + ((one + i * self).ln() - (one - i * self).ln()) / (two * i) + } + + /// Computes the hyperbolic sine of `self`. + #[inline] + pub fn sinh(self) -> Self { + // formula: sinh(a + bi) = sinh(a)cos(b) + i*cosh(a)sin(b) + Self::new( + self.re.sinh() * self.im.cos(), + self.re.cosh() * self.im.sin(), + ) + } + + /// Computes the hyperbolic cosine of `self`. + #[inline] + pub fn cosh(self) -> Self { + // formula: cosh(a + bi) = cosh(a)cos(b) + i*sinh(a)sin(b) + Self::new( + self.re.cosh() * self.im.cos(), + self.re.sinh() * self.im.sin(), + ) + } + + /// Computes the hyperbolic tangent of `self`. + #[inline] + pub fn tanh(self) -> Self { + // formula: tanh(a + bi) = (sinh(2a) + i*sin(2b))/(cosh(2a) + cos(2b)) + let (two_re, two_im) = (self.re + self.re, self.im + self.im); + Self::new(two_re.sinh(), two_im.sin()).unscale(two_re.cosh() + two_im.cos()) + } + + /// Computes the principal value of inverse hyperbolic sine of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞i, -i)`, continuous from the left. + /// * `(i, ∞i)`, continuous from the right. + /// + /// The branch satisfies `-π/2 ≤ Im(asinh(z)) ≤ π/2`. + #[inline] + pub fn asinh(self) -> Self { + // formula: arcsinh(z) = ln(z + sqrt(1+z^2)) + let one = Self::one(); + (self + (one + self * self).sqrt()).ln() + } + + /// Computes the principal value of inverse hyperbolic cosine of `self`. + /// + /// This function has one branch cut: + /// + /// * `(-∞, 1)`, continuous from above. + /// + /// The branch satisfies `-π ≤ Im(acosh(z)) ≤ π` and `0 ≤ Re(acosh(z)) < ∞`. + #[inline] + pub fn acosh(self) -> Self { + // formula: arccosh(z) = 2 ln(sqrt((z+1)/2) + sqrt((z-1)/2)) + let one = Self::one(); + let two = one + one; + two * (((self + one) / two).sqrt() + ((self - one) / two).sqrt()).ln() + } + + /// Computes the principal value of inverse hyperbolic tangent of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞, -1]`, continuous from above. + /// * `[1, ∞)`, continuous from below. + /// + /// The branch satisfies `-π/2 ≤ Im(atanh(z)) ≤ π/2`. + #[inline] + pub fn atanh(self) -> Self { + // formula: arctanh(z) = (ln(1+z) - ln(1-z))/2 + let one = Self::one(); + let two = one + one; + if self == one { + return Self::new(T::infinity(), T::zero()); + } else if self == -one { + return Self::new(-T::infinity(), T::zero()); + } + ((one + self).ln() - (one - self).ln()) / two + } + + /// Returns `1/self` using floating-point operations. + /// + /// This may be more accurate than the generic `self.inv()` in cases + /// where `self.norm_sqr()` would overflow to ∞ or underflow to 0. + /// + /// # Examples + /// + /// ``` + /// use num_complex::Complex64; + /// let c = Complex64::new(1e300, 1e300); + /// + /// // The generic `inv()` will overflow. + /// assert!(!c.inv().is_normal()); + /// + /// // But we can do better for `Float` types. + /// let inv = c.finv(); + /// assert!(inv.is_normal()); + /// println!("{:e}", inv); + /// + /// let expected = Complex64::new(5e-301, -5e-301); + /// assert!((inv - expected).norm() < 1e-315); + /// ``` + #[inline] + pub fn finv(self) -> Complex { + let norm = self.norm(); + self.conj() / norm / norm + } + + /// Returns `self/other` using floating-point operations. + /// + /// This may be more accurate than the generic `Div` implementation in cases + /// where `other.norm_sqr()` would overflow to ∞ or underflow to 0. + /// + /// # Examples + /// + /// ``` + /// use num_complex::Complex64; + /// let a = Complex64::new(2.0, 3.0); + /// let b = Complex64::new(1e300, 1e300); + /// + /// // Generic division will overflow. + /// assert!(!(a / b).is_normal()); + /// + /// // But we can do better for `Float` types. + /// let quotient = a.fdiv(b); + /// assert!(quotient.is_normal()); + /// println!("{:e}", quotient); + /// + /// let expected = Complex64::new(2.5e-300, 5e-301); + /// assert!((quotient - expected).norm() < 1e-315); + /// ``` + #[inline] + pub fn fdiv(self, other: Complex) -> Complex { + self * other.finv() + } +} + +#[cfg(any(feature = "std", feature = "libm"))] +impl Complex { + /// Computes `2^(self)`. + #[inline] + pub fn exp2(self) -> Self { + // formula: 2^(a + bi) = 2^a (cos(b*log2) + i*sin(b*log2)) + // = from_polar(2^a, b*log2) + Self::from_polar(self.re.exp2(), self.im * T::LN_2()) + } + + /// Computes the principal value of log base 2 of `self`. + #[inline] + pub fn log2(self) -> Self { + Self::ln(self) / T::LN_2() + } + + /// Computes the principal value of log base 10 of `self`. + #[inline] + pub fn log10(self) -> Self { + Self::ln(self) / T::LN_10() + } +} + +impl Complex { + /// Checks if the given complex number is NaN + #[inline] + pub fn is_nan(self) -> bool { + self.re.is_nan() || self.im.is_nan() + } + + /// Checks if the given complex number is infinite + #[inline] + pub fn is_infinite(self) -> bool { + !self.is_nan() && (self.re.is_infinite() || self.im.is_infinite()) + } + + /// Checks if the given complex number is finite + #[inline] + pub fn is_finite(self) -> bool { + self.re.is_finite() && self.im.is_finite() + } + + /// Checks if the given complex number is normal + #[inline] + pub fn is_normal(self) -> bool { + self.re.is_normal() && self.im.is_normal() + } +} + +// Safety: `Complex` is `repr(C)` and contains only instances of `T`, so we +// can guarantee it contains no *added* padding. Thus, if `T: Zeroable`, +// `Complex` is also `Zeroable` +#[cfg(feature = "bytemuck")] +unsafe impl bytemuck::Zeroable for Complex {} + +// Safety: `Complex` is `repr(C)` and contains only instances of `T`, so we +// can guarantee it contains no *added* padding. Thus, if `T: Pod`, +// `Complex` is also `Pod` +#[cfg(feature = "bytemuck")] +unsafe impl bytemuck::Pod for Complex {} + +impl From for Complex { + #[inline] + fn from(re: T) -> Self { + Self::new(re, T::zero()) + } +} + +impl<'a, T: Clone + Num> From<&'a T> for Complex { + #[inline] + fn from(re: &T) -> Self { + From::from(re.clone()) + } +} + +macro_rules! forward_ref_ref_binop { + (impl $imp:ident, $method:ident) => { + impl<'a, 'b, T: Clone + Num> $imp<&'b Complex> for &'a Complex { + type Output = Complex; + + #[inline] + fn $method(self, other: &Complex) -> Self::Output { + self.clone().$method(other.clone()) + } + } + }; +} + +macro_rules! forward_ref_val_binop { + (impl $imp:ident, $method:ident) => { + impl<'a, T: Clone + Num> $imp> for &'a Complex { + type Output = Complex; + + #[inline] + fn $method(self, other: Complex) -> Self::Output { + self.clone().$method(other) + } + } + }; +} + +macro_rules! forward_val_ref_binop { + (impl $imp:ident, $method:ident) => { + impl<'a, T: Clone + Num> $imp<&'a Complex> for Complex { + type Output = Complex; + + #[inline] + fn $method(self, other: &Complex) -> Self::Output { + self.$method(other.clone()) + } + } + }; +} + +macro_rules! forward_all_binop { + (impl $imp:ident, $method:ident) => { + forward_ref_ref_binop!(impl $imp, $method); + forward_ref_val_binop!(impl $imp, $method); + forward_val_ref_binop!(impl $imp, $method); + }; +} + +// arithmetic +forward_all_binop!(impl Add, add); + +// (a + i b) + (c + i d) == (a + c) + i (b + d) +impl Add> for Complex { + type Output = Self; + + #[inline] + fn add(self, other: Self) -> Self::Output { + Self::Output::new(self.re + other.re, self.im + other.im) + } +} + +forward_all_binop!(impl Sub, sub); + +// (a + i b) - (c + i d) == (a - c) + i (b - d) +impl Sub> for Complex { + type Output = Self; + + #[inline] + fn sub(self, other: Self) -> Self::Output { + Self::Output::new(self.re - other.re, self.im - other.im) + } +} + +forward_all_binop!(impl Mul, mul); + +// (a + i b) * (c + i d) == (a*c - b*d) + i (a*d + b*c) +impl Mul> for Complex { + type Output = Self; + + #[inline] + fn mul(self, other: Self) -> Self::Output { + let re = self.re.clone() * other.re.clone() - self.im.clone() * other.im.clone(); + let im = self.re * other.im + self.im * other.re; + Self::Output::new(re, im) + } +} + +// (a + i b) * (c + i d) + (e + i f) == ((a*c + e) - b*d) + i (a*d + (b*c + f)) +impl> MulAdd> for Complex { + type Output = Complex; + + #[inline] + fn mul_add(self, other: Complex, add: Complex) -> Complex { + let re = self.re.clone().mul_add(other.re.clone(), add.re) + - (self.im.clone() * other.im.clone()); // FIXME: use mulsub when available in rust + let im = self.re.mul_add(other.im, self.im.mul_add(other.re, add.im)); + Complex::new(re, im) + } +} +impl<'a, 'b, T: Clone + Num + MulAdd> MulAdd<&'b Complex> for &'a Complex { + type Output = Complex; + + #[inline] + fn mul_add(self, other: &Complex, add: &Complex) -> Complex { + self.clone().mul_add(other.clone(), add.clone()) + } +} + +forward_all_binop!(impl Div, div); + +// (a + i b) / (c + i d) == [(a + i b) * (c - i d)] / (c*c + d*d) +// == [(a*c + b*d) / (c*c + d*d)] + i [(b*c - a*d) / (c*c + d*d)] +impl Div> for Complex { + type Output = Self; + + #[inline] + fn div(self, other: Self) -> Self::Output { + let norm_sqr = other.norm_sqr(); + let re = self.re.clone() * other.re.clone() + self.im.clone() * other.im.clone(); + let im = self.im * other.re - self.re * other.im; + Self::Output::new(re / norm_sqr.clone(), im / norm_sqr) + } +} + +forward_all_binop!(impl Rem, rem); + +impl Complex { + /// Find the gaussian integer corresponding to the true ratio rounded towards zero. + fn div_trunc(&self, divisor: &Self) -> Self { + let Complex { re, im } = self / divisor; + Complex::new(re.clone() - re % T::one(), im.clone() - im % T::one()) + } +} + +impl Rem> for Complex { + type Output = Self; + + #[inline] + fn rem(self, modulus: Self) -> Self::Output { + let gaussian = self.div_trunc(&modulus); + self - modulus * gaussian + } +} + +// Op Assign + +mod opassign { + use core::ops::{AddAssign, DivAssign, MulAssign, RemAssign, SubAssign}; + + use num_traits::{MulAddAssign, NumAssign}; + + use crate::Complex; + + impl AddAssign for Complex { + fn add_assign(&mut self, other: Self) { + self.re += other.re; + self.im += other.im; + } + } + + impl SubAssign for Complex { + fn sub_assign(&mut self, other: Self) { + self.re -= other.re; + self.im -= other.im; + } + } + + // (a + i b) * (c + i d) == (a*c - b*d) + i (a*d + b*c) + impl MulAssign for Complex { + fn mul_assign(&mut self, other: Self) { + let a = self.re.clone(); + + self.re *= other.re.clone(); + self.re -= self.im.clone() * other.im.clone(); + + self.im *= other.re; + self.im += a * other.im; + } + } + + // (a + i b) * (c + i d) + (e + i f) == ((a*c + e) - b*d) + i (b*c + (a*d + f)) + impl MulAddAssign for Complex { + fn mul_add_assign(&mut self, other: Complex, add: Complex) { + let a = self.re.clone(); + + self.re.mul_add_assign(other.re.clone(), add.re); // (a*c + e) + self.re -= self.im.clone() * other.im.clone(); // ((a*c + e) - b*d) + + let mut adf = a; + adf.mul_add_assign(other.im, add.im); // (a*d + f) + self.im.mul_add_assign(other.re, adf); // (b*c + (a*d + f)) + } + } + + impl<'a, 'b, T: Clone + NumAssign + MulAddAssign> MulAddAssign<&'a Complex, &'b Complex> + for Complex + { + fn mul_add_assign(&mut self, other: &Complex, add: &Complex) { + self.mul_add_assign(other.clone(), add.clone()); + } + } + + // (a + i b) / (c + i d) == [(a + i b) * (c - i d)] / (c*c + d*d) + // == [(a*c + b*d) / (c*c + d*d)] + i [(b*c - a*d) / (c*c + d*d)] + impl DivAssign for Complex { + fn div_assign(&mut self, other: Self) { + let a = self.re.clone(); + let norm_sqr = other.norm_sqr(); + + self.re *= other.re.clone(); + self.re += self.im.clone() * other.im.clone(); + self.re /= norm_sqr.clone(); + + self.im *= other.re; + self.im -= a * other.im; + self.im /= norm_sqr; + } + } + + impl RemAssign for Complex { + fn rem_assign(&mut self, modulus: Self) { + let gaussian = self.div_trunc(&modulus); + *self -= modulus * gaussian; + } + } + + impl AddAssign for Complex { + fn add_assign(&mut self, other: T) { + self.re += other; + } + } + + impl SubAssign for Complex { + fn sub_assign(&mut self, other: T) { + self.re -= other; + } + } + + impl MulAssign for Complex { + fn mul_assign(&mut self, other: T) { + self.re *= other.clone(); + self.im *= other; + } + } + + impl DivAssign for Complex { + fn div_assign(&mut self, other: T) { + self.re /= other.clone(); + self.im /= other; + } + } + + impl RemAssign for Complex { + fn rem_assign(&mut self, other: T) { + self.re %= other.clone(); + self.im %= other; + } + } + + macro_rules! forward_op_assign { + (impl $imp:ident, $method:ident) => { + impl<'a, T: Clone + NumAssign> $imp<&'a Complex> for Complex { + #[inline] + fn $method(&mut self, other: &Self) { + self.$method(other.clone()) + } + } + impl<'a, T: Clone + NumAssign> $imp<&'a T> for Complex { + #[inline] + fn $method(&mut self, other: &T) { + self.$method(other.clone()) + } + } + }; + } + + forward_op_assign!(impl AddAssign, add_assign); + forward_op_assign!(impl SubAssign, sub_assign); + forward_op_assign!(impl MulAssign, mul_assign); + forward_op_assign!(impl DivAssign, div_assign); + forward_op_assign!(impl RemAssign, rem_assign); +} + +impl> Neg for Complex { + type Output = Self; + + #[inline] + fn neg(self) -> Self::Output { + Self::Output::new(-self.re, -self.im) + } +} + +impl<'a, T: Clone + Num + Neg> Neg for &'a Complex { + type Output = Complex; + + #[inline] + fn neg(self) -> Self::Output { + -self.clone() + } +} + +impl> Inv for Complex { + type Output = Self; + + #[inline] + fn inv(self) -> Self::Output { + (&self).inv() + } +} + +impl<'a, T: Clone + Num + Neg> Inv for &'a Complex { + type Output = Complex; + + #[inline] + fn inv(self) -> Self::Output { + self.inv() + } +} + +macro_rules! real_arithmetic { + (@forward $imp:ident::$method:ident for $($real:ident),*) => ( + impl<'a, T: Clone + Num> $imp<&'a T> for Complex { + type Output = Complex; + + #[inline] + fn $method(self, other: &T) -> Self::Output { + self.$method(other.clone()) + } + } + impl<'a, T: Clone + Num> $imp for &'a Complex { + type Output = Complex; + + #[inline] + fn $method(self, other: T) -> Self::Output { + self.clone().$method(other) + } + } + impl<'a, 'b, T: Clone + Num> $imp<&'a T> for &'b Complex { + type Output = Complex; + + #[inline] + fn $method(self, other: &T) -> Self::Output { + self.clone().$method(other.clone()) + } + } + $( + impl<'a> $imp<&'a Complex<$real>> for $real { + type Output = Complex<$real>; + + #[inline] + fn $method(self, other: &Complex<$real>) -> Complex<$real> { + self.$method(other.clone()) + } + } + impl<'a> $imp> for &'a $real { + type Output = Complex<$real>; + + #[inline] + fn $method(self, other: Complex<$real>) -> Complex<$real> { + self.clone().$method(other) + } + } + impl<'a, 'b> $imp<&'a Complex<$real>> for &'b $real { + type Output = Complex<$real>; + + #[inline] + fn $method(self, other: &Complex<$real>) -> Complex<$real> { + self.clone().$method(other.clone()) + } + } + )* + ); + ($($real:ident),*) => ( + real_arithmetic!(@forward Add::add for $($real),*); + real_arithmetic!(@forward Sub::sub for $($real),*); + real_arithmetic!(@forward Mul::mul for $($real),*); + real_arithmetic!(@forward Div::div for $($real),*); + real_arithmetic!(@forward Rem::rem for $($real),*); + + $( + impl Add> for $real { + type Output = Complex<$real>; + + #[inline] + fn add(self, other: Complex<$real>) -> Self::Output { + Self::Output::new(self + other.re, other.im) + } + } + + impl Sub> for $real { + type Output = Complex<$real>; + + #[inline] + fn sub(self, other: Complex<$real>) -> Self::Output { + Self::Output::new(self - other.re, $real::zero() - other.im) + } + } + + impl Mul> for $real { + type Output = Complex<$real>; + + #[inline] + fn mul(self, other: Complex<$real>) -> Self::Output { + Self::Output::new(self * other.re, self * other.im) + } + } + + impl Div> for $real { + type Output = Complex<$real>; + + #[inline] + fn div(self, other: Complex<$real>) -> Self::Output { + // a / (c + i d) == [a * (c - i d)] / (c*c + d*d) + let norm_sqr = other.norm_sqr(); + Self::Output::new(self * other.re / norm_sqr.clone(), + $real::zero() - self * other.im / norm_sqr) + } + } + + impl Rem> for $real { + type Output = Complex<$real>; + + #[inline] + fn rem(self, other: Complex<$real>) -> Self::Output { + Self::Output::new(self, Self::zero()) % other + } + } + )* + ); +} + +impl Add for Complex { + type Output = Complex; + + #[inline] + fn add(self, other: T) -> Self::Output { + Self::Output::new(self.re + other, self.im) + } +} + +impl Sub for Complex { + type Output = Complex; + + #[inline] + fn sub(self, other: T) -> Self::Output { + Self::Output::new(self.re - other, self.im) + } +} + +impl Mul for Complex { + type Output = Complex; + + #[inline] + fn mul(self, other: T) -> Self::Output { + Self::Output::new(self.re * other.clone(), self.im * other) + } +} + +impl Div for Complex { + type Output = Self; + + #[inline] + fn div(self, other: T) -> Self::Output { + Self::Output::new(self.re / other.clone(), self.im / other) + } +} + +impl Rem for Complex { + type Output = Complex; + + #[inline] + fn rem(self, other: T) -> Self::Output { + Self::Output::new(self.re % other.clone(), self.im % other) + } +} + +real_arithmetic!(usize, u8, u16, u32, u64, u128, isize, i8, i16, i32, i64, i128, f32, f64); + +// constants +impl Zero for Complex { + #[inline] + fn zero() -> Self { + Self::new(Zero::zero(), Zero::zero()) + } + + #[inline] + fn is_zero(&self) -> bool { + self.re.is_zero() && self.im.is_zero() + } + + #[inline] + fn set_zero(&mut self) { + self.re.set_zero(); + self.im.set_zero(); + } +} + +impl One for Complex { + #[inline] + fn one() -> Self { + Self::new(One::one(), Zero::zero()) + } + + #[inline] + fn is_one(&self) -> bool { + self.re.is_one() && self.im.is_zero() + } + + #[inline] + fn set_one(&mut self) { + self.re.set_one(); + self.im.set_zero(); + } +} + +macro_rules! write_complex { + ($f:ident, $t:expr, $prefix:expr, $re:expr, $im:expr, $T:ident) => {{ + let abs_re = if $re < Zero::zero() { + $T::zero() - $re.clone() + } else { + $re.clone() + }; + let abs_im = if $im < Zero::zero() { + $T::zero() - $im.clone() + } else { + $im.clone() + }; + + return if let Some(prec) = $f.precision() { + fmt_re_im( + $f, + $re < $T::zero(), + $im < $T::zero(), + format_args!(concat!("{:.1$", $t, "}"), abs_re, prec), + format_args!(concat!("{:.1$", $t, "}"), abs_im, prec), + ) + } else { + fmt_re_im( + $f, + $re < $T::zero(), + $im < $T::zero(), + format_args!(concat!("{:", $t, "}"), abs_re), + format_args!(concat!("{:", $t, "}"), abs_im), + ) + }; + + fn fmt_re_im( + f: &mut fmt::Formatter<'_>, + re_neg: bool, + im_neg: bool, + real: fmt::Arguments<'_>, + imag: fmt::Arguments<'_>, + ) -> fmt::Result { + let prefix = if f.alternate() { $prefix } else { "" }; + let sign = if re_neg { + "-" + } else if f.sign_plus() { + "+" + } else { + "" + }; + + if im_neg { + fmt_complex( + f, + format_args!( + "{}{pre}{re}-{pre}{im}i", + sign, + re = real, + im = imag, + pre = prefix + ), + ) + } else { + fmt_complex( + f, + format_args!( + "{}{pre}{re}+{pre}{im}i", + sign, + re = real, + im = imag, + pre = prefix + ), + ) + } + } + + #[cfg(feature = "std")] + // Currently, we can only apply width using an intermediate `String` (and thus `std`) + fn fmt_complex(f: &mut fmt::Formatter<'_>, complex: fmt::Arguments<'_>) -> fmt::Result { + use std::string::ToString; + if let Some(width) = f.width() { + write!(f, "{0: >1$}", complex.to_string(), width) + } else { + write!(f, "{}", complex) + } + } + + #[cfg(not(feature = "std"))] + fn fmt_complex(f: &mut fmt::Formatter<'_>, complex: fmt::Arguments<'_>) -> fmt::Result { + write!(f, "{}", complex) + } + }}; +} + +// string conversions +impl fmt::Display for Complex +where + T: fmt::Display + Num + PartialOrd + Clone, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write_complex!(f, "", "", self.re, self.im, T) + } +} + +impl fmt::LowerExp for Complex +where + T: fmt::LowerExp + Num + PartialOrd + Clone, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write_complex!(f, "e", "", self.re, self.im, T) + } +} + +impl fmt::UpperExp for Complex +where + T: fmt::UpperExp + Num + PartialOrd + Clone, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write_complex!(f, "E", "", self.re, self.im, T) + } +} + +impl fmt::LowerHex for Complex +where + T: fmt::LowerHex + Num + PartialOrd + Clone, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write_complex!(f, "x", "0x", self.re, self.im, T) + } +} + +impl fmt::UpperHex for Complex +where + T: fmt::UpperHex + Num + PartialOrd + Clone, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write_complex!(f, "X", "0x", self.re, self.im, T) + } +} + +impl fmt::Octal for Complex +where + T: fmt::Octal + Num + PartialOrd + Clone, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write_complex!(f, "o", "0o", self.re, self.im, T) + } +} + +impl fmt::Binary for Complex +where + T: fmt::Binary + Num + PartialOrd + Clone, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write_complex!(f, "b", "0b", self.re, self.im, T) + } +} + +#[allow(deprecated)] // `trim_left_matches` and `trim_right_matches` since 1.33 +fn from_str_generic(s: &str, from: F) -> Result, ParseComplexError> +where + F: Fn(&str) -> Result, + T: Clone + Num, +{ + let imag = match s.rfind('j') { + None => 'i', + _ => 'j', + }; + + let mut neg_b = false; + let mut a = s; + let mut b = ""; + + for (i, w) in s.as_bytes().windows(2).enumerate() { + let p = w[0]; + let c = w[1]; + + // ignore '+'/'-' if part of an exponent + if (c == b'+' || c == b'-') && !(p == b'e' || p == b'E') { + // trim whitespace around the separator + a = &s[..=i].trim_right_matches(char::is_whitespace); + b = &s[i + 2..].trim_left_matches(char::is_whitespace); + neg_b = c == b'-'; + + if b.is_empty() || (neg_b && b.starts_with('-')) { + return Err(ParseComplexError::expr_error()); + } + break; + } + } + + // split off real and imaginary parts + if b.is_empty() { + // input was either pure real or pure imaginary + b = if a.ends_with(imag) { "0" } else { "0i" }; + } + + let re; + let neg_re; + let im; + let neg_im; + if a.ends_with(imag) { + im = a; + neg_im = false; + re = b; + neg_re = neg_b; + } else if b.ends_with(imag) { + re = a; + neg_re = false; + im = b; + neg_im = neg_b; + } else { + return Err(ParseComplexError::expr_error()); + } + + // parse re + let re = from(re).map_err(ParseComplexError::from_error)?; + let re = if neg_re { T::zero() - re } else { re }; + + // pop imaginary unit off + let mut im = &im[..im.len() - 1]; + // handle im == "i" or im == "-i" + if im.is_empty() || im == "+" { + im = "1"; + } else if im == "-" { + im = "-1"; + } + + // parse im + let im = from(im).map_err(ParseComplexError::from_error)?; + let im = if neg_im { T::zero() - im } else { im }; + + Ok(Complex::new(re, im)) +} + +impl FromStr for Complex +where + T: FromStr + Num + Clone, +{ + type Err = ParseComplexError; + + /// Parses `a +/- bi`; `ai +/- b`; `a`; or `bi` where `a` and `b` are of type `T` + fn from_str(s: &str) -> Result { + from_str_generic(s, T::from_str) + } +} + +impl Num for Complex { + type FromStrRadixErr = ParseComplexError; + + /// Parses `a +/- bi`; `ai +/- b`; `a`; or `bi` where `a` and `b` are of type `T` + /// + /// `radix` must be <= 18; larger radix would include *i* and *j* as digits, + /// which cannot be supported. + /// + /// The conversion returns an error if 18 <= radix <= 36; it panics if radix > 36. + /// + /// The elements of `T` are parsed using `Num::from_str_radix` too, and errors + /// (or panics) from that are reflected here as well. + fn from_str_radix(s: &str, radix: u32) -> Result { + assert!( + radix <= 36, + "from_str_radix: radix is too high (maximum 36)" + ); + + // larger radix would include 'i' and 'j' as digits, which cannot be supported + if radix > 18 { + return Err(ParseComplexError::unsupported_radix()); + } + + from_str_generic(s, |x| -> Result { + T::from_str_radix(x, radix) + }) + } +} + +impl Sum for Complex { + fn sum(iter: I) -> Self + where + I: Iterator, + { + iter.fold(Self::zero(), |acc, c| acc + c) + } +} + +impl<'a, T: 'a + Num + Clone> Sum<&'a Complex> for Complex { + fn sum(iter: I) -> Self + where + I: Iterator>, + { + iter.fold(Self::zero(), |acc, c| acc + c) + } +} + +impl Product for Complex { + fn product(iter: I) -> Self + where + I: Iterator, + { + iter.fold(Self::one(), |acc, c| acc * c) + } +} + +impl<'a, T: 'a + Num + Clone> Product<&'a Complex> for Complex { + fn product(iter: I) -> Self + where + I: Iterator>, + { + iter.fold(Self::one(), |acc, c| acc * c) + } +} + +#[cfg(feature = "serde")] +impl serde::Serialize for Complex +where + T: serde::Serialize, +{ + fn serialize(&self, serializer: S) -> Result + where + S: serde::Serializer, + { + (&self.re, &self.im).serialize(serializer) + } +} + +#[cfg(feature = "serde")] +impl<'de, T> serde::Deserialize<'de> for Complex +where + T: serde::Deserialize<'de> + Num + Clone, +{ + fn deserialize(deserializer: D) -> Result + where + D: serde::Deserializer<'de>, + { + let (re, im) = serde::Deserialize::deserialize(deserializer)?; + Ok(Self::new(re, im)) + } +} + +#[derive(Debug, PartialEq)] +pub struct ParseComplexError { + kind: ComplexErrorKind, +} + +#[derive(Debug, PartialEq)] +enum ComplexErrorKind { + ParseError(E), + ExprError, + UnsupportedRadix, +} + +impl ParseComplexError { + fn expr_error() -> Self { + ParseComplexError { + kind: ComplexErrorKind::ExprError, + } + } + + fn unsupported_radix() -> Self { + ParseComplexError { + kind: ComplexErrorKind::UnsupportedRadix, + } + } + + fn from_error(error: E) -> Self { + ParseComplexError { + kind: ComplexErrorKind::ParseError(error), + } + } +} + +#[cfg(feature = "std")] +impl Error for ParseComplexError { + #[allow(deprecated)] + fn description(&self) -> &str { + match self.kind { + ComplexErrorKind::ParseError(ref e) => e.description(), + ComplexErrorKind::ExprError => "invalid or unsupported complex expression", + ComplexErrorKind::UnsupportedRadix => "unsupported radix for conversion", + } + } +} + +impl fmt::Display for ParseComplexError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self.kind { + ComplexErrorKind::ParseError(ref e) => e.fmt(f), + ComplexErrorKind::ExprError => "invalid or unsupported complex expression".fmt(f), + ComplexErrorKind::UnsupportedRadix => "unsupported radix for conversion".fmt(f), + } + } +} + +#[cfg(test)] +fn hash(x: &T) -> u64 { + use std::collections::hash_map::RandomState; + use std::hash::{BuildHasher, Hasher}; + let mut hasher = ::Hasher::new(); + x.hash(&mut hasher); + hasher.finish() +} + +#[cfg(test)] +pub(crate) mod test { + #![allow(non_upper_case_globals)] + + use super::{Complex, Complex64}; + use super::{ComplexErrorKind, ParseComplexError}; + use core::f64; + use core::str::FromStr; + + use std::string::{String, ToString}; + + use num_traits::{Num, One, Zero}; + + pub const _0_0i: Complex64 = Complex::new(0.0, 0.0); + pub const _1_0i: Complex64 = Complex::new(1.0, 0.0); + pub const _1_1i: Complex64 = Complex::new(1.0, 1.0); + pub const _0_1i: Complex64 = Complex::new(0.0, 1.0); + pub const _neg1_1i: Complex64 = Complex::new(-1.0, 1.0); + pub const _05_05i: Complex64 = Complex::new(0.5, 0.5); + pub const all_consts: [Complex64; 5] = [_0_0i, _1_0i, _1_1i, _neg1_1i, _05_05i]; + pub const _4_2i: Complex64 = Complex::new(4.0, 2.0); + pub const _1_infi: Complex64 = Complex::new(1.0, f64::INFINITY); + pub const _neg1_infi: Complex64 = Complex::new(-1.0, f64::INFINITY); + pub const _1_nani: Complex64 = Complex::new(1.0, f64::NAN); + pub const _neg1_nani: Complex64 = Complex::new(-1.0, f64::NAN); + pub const _inf_0i: Complex64 = Complex::new(f64::INFINITY, 0.0); + pub const _neginf_1i: Complex64 = Complex::new(f64::NEG_INFINITY, 1.0); + pub const _neginf_neg1i: Complex64 = Complex::new(f64::NEG_INFINITY, -1.0); + pub const _inf_1i: Complex64 = Complex::new(f64::INFINITY, 1.0); + pub const _inf_neg1i: Complex64 = Complex::new(f64::INFINITY, -1.0); + pub const _neginf_infi: Complex64 = Complex::new(f64::NEG_INFINITY, f64::INFINITY); + pub const _inf_infi: Complex64 = Complex::new(f64::INFINITY, f64::INFINITY); + pub const _neginf_nani: Complex64 = Complex::new(f64::NEG_INFINITY, f64::NAN); + pub const _inf_nani: Complex64 = Complex::new(f64::INFINITY, f64::NAN); + pub const _nan_0i: Complex64 = Complex::new(f64::NAN, 0.0); + pub const _nan_1i: Complex64 = Complex::new(f64::NAN, 1.0); + pub const _nan_neg1i: Complex64 = Complex::new(f64::NAN, -1.0); + pub const _nan_nani: Complex64 = Complex::new(f64::NAN, f64::NAN); + + #[test] + fn test_consts() { + // check our constants are what Complex::new creates + fn test(c: Complex64, r: f64, i: f64) { + assert_eq!(c, Complex::new(r, i)); + } + test(_0_0i, 0.0, 0.0); + test(_1_0i, 1.0, 0.0); + test(_1_1i, 1.0, 1.0); + test(_neg1_1i, -1.0, 1.0); + test(_05_05i, 0.5, 0.5); + + assert_eq!(_0_0i, Zero::zero()); + assert_eq!(_1_0i, One::one()); + } + + #[test] + fn test_scale_unscale() { + assert_eq!(_05_05i.scale(2.0), _1_1i); + assert_eq!(_1_1i.unscale(2.0), _05_05i); + for &c in all_consts.iter() { + assert_eq!(c.scale(2.0).unscale(2.0), c); + } + } + + #[test] + fn test_conj() { + for &c in all_consts.iter() { + assert_eq!(c.conj(), Complex::new(c.re, -c.im)); + assert_eq!(c.conj().conj(), c); + } + } + + #[test] + fn test_inv() { + assert_eq!(_1_1i.inv(), _05_05i.conj()); + assert_eq!(_1_0i.inv(), _1_0i.inv()); + } + + #[test] + #[should_panic] + fn test_divide_by_zero_natural() { + let n = Complex::new(2, 3); + let d = Complex::new(0, 0); + let _x = n / d; + } + + #[test] + fn test_inv_zero() { + // FIXME #20: should this really fail, or just NaN? + assert!(_0_0i.inv().is_nan()); + } + + #[test] + #[allow(clippy::float_cmp)] + fn test_l1_norm() { + assert_eq!(_0_0i.l1_norm(), 0.0); + assert_eq!(_1_0i.l1_norm(), 1.0); + assert_eq!(_1_1i.l1_norm(), 2.0); + assert_eq!(_0_1i.l1_norm(), 1.0); + assert_eq!(_neg1_1i.l1_norm(), 2.0); + assert_eq!(_05_05i.l1_norm(), 1.0); + assert_eq!(_4_2i.l1_norm(), 6.0); + } + + #[test] + fn test_pow() { + for c in all_consts.iter() { + assert_eq!(c.powi(0), _1_0i); + let mut pos = _1_0i; + let mut neg = _1_0i; + for i in 1i32..20 { + pos *= c; + assert_eq!(pos, c.powi(i)); + if c.is_zero() { + assert!(c.powi(-i).is_nan()); + } else { + neg /= c; + assert_eq!(neg, c.powi(-i)); + } + } + } + } + + #[cfg(any(feature = "std", feature = "libm"))] + pub(crate) mod float { + use super::*; + use num_traits::{Float, Pow}; + + #[test] + fn test_cis() { + assert!(close(Complex::cis(0.0 * f64::consts::PI), _1_0i)); + assert!(close(Complex::cis(0.5 * f64::consts::PI), _0_1i)); + assert!(close(Complex::cis(1.0 * f64::consts::PI), -_1_0i)); + assert!(close(Complex::cis(1.5 * f64::consts::PI), -_0_1i)); + assert!(close(Complex::cis(2.0 * f64::consts::PI), _1_0i)); + } + + #[test] + #[cfg_attr(target_arch = "x86", ignore)] + // FIXME #7158: (maybe?) currently failing on x86. + #[allow(clippy::float_cmp)] + fn test_norm() { + fn test(c: Complex64, ns: f64) { + assert_eq!(c.norm_sqr(), ns); + assert_eq!(c.norm(), ns.sqrt()) + } + test(_0_0i, 0.0); + test(_1_0i, 1.0); + test(_1_1i, 2.0); + test(_neg1_1i, 2.0); + test(_05_05i, 0.5); + } + + #[test] + fn test_arg() { + fn test(c: Complex64, arg: f64) { + assert!((c.arg() - arg).abs() < 1.0e-6) + } + test(_1_0i, 0.0); + test(_1_1i, 0.25 * f64::consts::PI); + test(_neg1_1i, 0.75 * f64::consts::PI); + test(_05_05i, 0.25 * f64::consts::PI); + } + + #[test] + fn test_polar_conv() { + fn test(c: Complex64) { + let (r, theta) = c.to_polar(); + assert!((c - Complex::from_polar(r, theta)).norm() < 1e-6); + } + for &c in all_consts.iter() { + test(c); + } + } + + pub(crate) fn close(a: Complex64, b: Complex64) -> bool { + close_to_tol(a, b, 1e-10) + } + + fn close_to_tol(a: Complex64, b: Complex64, tol: f64) -> bool { + // returns true if a and b are reasonably close + let close = (a == b) || (a - b).norm() < tol; + if !close { + println!("{:?} != {:?}", a, b); + } + close + } + + // Version that also works if re or im are +inf, -inf, or nan + fn close_naninf(a: Complex64, b: Complex64) -> bool { + close_naninf_to_tol(a, b, 1.0e-10) + } + + fn close_naninf_to_tol(a: Complex64, b: Complex64, tol: f64) -> bool { + let mut close = true; + + // Compare the real parts + if a.re.is_finite() { + if b.re.is_finite() { + close = (a.re == b.re) || (a.re - b.re).abs() < tol; + } else { + close = false; + } + } else if (a.re.is_nan() && !b.re.is_nan()) + || (a.re.is_infinite() + && a.re.is_sign_positive() + && !(b.re.is_infinite() && b.re.is_sign_positive())) + || (a.re.is_infinite() + && a.re.is_sign_negative() + && !(b.re.is_infinite() && b.re.is_sign_negative())) + { + close = false; + } + + // Compare the imaginary parts + if a.im.is_finite() { + if b.im.is_finite() { + close &= (a.im == b.im) || (a.im - b.im).abs() < tol; + } else { + close = false; + } + } else if (a.im.is_nan() && !b.im.is_nan()) + || (a.im.is_infinite() + && a.im.is_sign_positive() + && !(b.im.is_infinite() && b.im.is_sign_positive())) + || (a.im.is_infinite() + && a.im.is_sign_negative() + && !(b.im.is_infinite() && b.im.is_sign_negative())) + { + close = false; + } + + if close == false { + println!("{:?} != {:?}", a, b); + } + close + } + + #[test] + fn test_exp2() { + assert!(close(_0_0i.exp2(), _1_0i)); + } + + #[test] + fn test_exp() { + assert!(close(_1_0i.exp(), _1_0i.scale(f64::consts::E))); + assert!(close(_0_0i.exp(), _1_0i)); + assert!(close(_0_1i.exp(), Complex::new(1.0.cos(), 1.0.sin()))); + assert!(close(_05_05i.exp() * _05_05i.exp(), _1_1i.exp())); + assert!(close( + _0_1i.scale(-f64::consts::PI).exp(), + _1_0i.scale(-1.0) + )); + for &c in all_consts.iter() { + // e^conj(z) = conj(e^z) + assert!(close(c.conj().exp(), c.exp().conj())); + // e^(z + 2 pi i) = e^z + assert!(close( + c.exp(), + (c + _0_1i.scale(f64::consts::PI * 2.0)).exp() + )); + } + + // The test values below were taken from https://en.cppreference.com/w/cpp/numeric/complex/exp + assert!(close_naninf(_1_infi.exp(), _nan_nani)); + assert!(close_naninf(_neg1_infi.exp(), _nan_nani)); + assert!(close_naninf(_1_nani.exp(), _nan_nani)); + assert!(close_naninf(_neg1_nani.exp(), _nan_nani)); + assert!(close_naninf(_inf_0i.exp(), _inf_0i)); + assert!(close_naninf(_neginf_1i.exp(), 0.0 * Complex::cis(1.0))); + assert!(close_naninf(_neginf_neg1i.exp(), 0.0 * Complex::cis(-1.0))); + assert!(close_naninf( + _inf_1i.exp(), + f64::INFINITY * Complex::cis(1.0) + )); + assert!(close_naninf( + _inf_neg1i.exp(), + f64::INFINITY * Complex::cis(-1.0) + )); + assert!(close_naninf(_neginf_infi.exp(), _0_0i)); // Note: ±0±0i: signs of zeros are unspecified + assert!(close_naninf(_inf_infi.exp(), _inf_nani)); // Note: ±∞+NaN*i: sign of the real part is unspecified + assert!(close_naninf(_neginf_nani.exp(), _0_0i)); // Note: ±0±0i: signs of zeros are unspecified + assert!(close_naninf(_inf_nani.exp(), _inf_nani)); // Note: ±∞+NaN*i: sign of the real part is unspecified + assert!(close_naninf(_nan_0i.exp(), _nan_0i)); + assert!(close_naninf(_nan_1i.exp(), _nan_nani)); + assert!(close_naninf(_nan_neg1i.exp(), _nan_nani)); + assert!(close_naninf(_nan_nani.exp(), _nan_nani)); + } + + #[test] + fn test_ln() { + assert!(close(_1_0i.ln(), _0_0i)); + assert!(close(_0_1i.ln(), _0_1i.scale(f64::consts::PI / 2.0))); + assert!(close(_0_0i.ln(), Complex::new(f64::neg_infinity(), 0.0))); + assert!(close( + (_neg1_1i * _05_05i).ln(), + _neg1_1i.ln() + _05_05i.ln() + )); + for &c in all_consts.iter() { + // ln(conj(z() = conj(ln(z)) + assert!(close(c.conj().ln(), c.ln().conj())); + // for this branch, -pi <= arg(ln(z)) <= pi + assert!(-f64::consts::PI <= c.ln().arg() && c.ln().arg() <= f64::consts::PI); + } + } + + #[test] + fn test_powc() { + let a = Complex::new(2.0, -3.0); + let b = Complex::new(3.0, 0.0); + assert!(close(a.powc(b), a.powf(b.re))); + assert!(close(b.powc(a), a.expf(b.re))); + let c = Complex::new(1.0 / 3.0, 0.1); + assert!(close_to_tol( + a.powc(c), + Complex::new(1.65826, -0.33502), + 1e-5 + )); + } + + #[test] + fn test_powf() { + let c = Complex64::new(2.0, -1.0); + let expected = Complex64::new(-0.8684746, -16.695934); + assert!(close_to_tol(c.powf(3.5), expected, 1e-5)); + assert!(close_to_tol(Pow::pow(c, 3.5_f64), expected, 1e-5)); + assert!(close_to_tol(Pow::pow(c, 3.5_f32), expected, 1e-5)); + } + + #[test] + fn test_log() { + let c = Complex::new(2.0, -1.0); + let r = c.log(10.0); + assert!(close_to_tol(r, Complex::new(0.349485, -0.20135958), 1e-5)); + } + + #[test] + fn test_log2() { + assert!(close(_1_0i.log2(), _0_0i)); + } + + #[test] + fn test_log10() { + assert!(close(_1_0i.log10(), _0_0i)); + } + + #[test] + fn test_some_expf_cases() { + let c = Complex::new(2.0, -1.0); + let r = c.expf(10.0); + assert!(close_to_tol(r, Complex::new(-66.82015, -74.39803), 1e-5)); + + let c = Complex::new(5.0, -2.0); + let r = c.expf(3.4); + assert!(close_to_tol(r, Complex::new(-349.25, -290.63), 1e-2)); + + let c = Complex::new(-1.5, 2.0 / 3.0); + let r = c.expf(1.0 / 3.0); + assert!(close_to_tol(r, Complex::new(3.8637, -3.4745), 1e-2)); + } + + #[test] + fn test_sqrt() { + assert!(close(_0_0i.sqrt(), _0_0i)); + assert!(close(_1_0i.sqrt(), _1_0i)); + assert!(close(Complex::new(-1.0, 0.0).sqrt(), _0_1i)); + assert!(close(Complex::new(-1.0, -0.0).sqrt(), _0_1i.scale(-1.0))); + assert!(close(_0_1i.sqrt(), _05_05i.scale(2.0.sqrt()))); + for &c in all_consts.iter() { + // sqrt(conj(z() = conj(sqrt(z)) + assert!(close(c.conj().sqrt(), c.sqrt().conj())); + // for this branch, -pi/2 <= arg(sqrt(z)) <= pi/2 + assert!( + -f64::consts::FRAC_PI_2 <= c.sqrt().arg() + && c.sqrt().arg() <= f64::consts::FRAC_PI_2 + ); + // sqrt(z) * sqrt(z) = z + assert!(close(c.sqrt() * c.sqrt(), c)); + } + } + + #[test] + fn test_sqrt_real() { + for n in (0..100).map(f64::from) { + // √(n² + 0i) = n + 0i + let n2 = n * n; + assert_eq!(Complex64::new(n2, 0.0).sqrt(), Complex64::new(n, 0.0)); + // √(-n² + 0i) = 0 + ni + assert_eq!(Complex64::new(-n2, 0.0).sqrt(), Complex64::new(0.0, n)); + // √(-n² - 0i) = 0 - ni + assert_eq!(Complex64::new(-n2, -0.0).sqrt(), Complex64::new(0.0, -n)); + } + } + + #[test] + fn test_sqrt_imag() { + for n in (0..100).map(f64::from) { + // √(0 + n²i) = n e^(iπ/4) + let n2 = n * n; + assert!(close( + Complex64::new(0.0, n2).sqrt(), + Complex64::from_polar(n, f64::consts::FRAC_PI_4) + )); + // √(0 - n²i) = n e^(-iπ/4) + assert!(close( + Complex64::new(0.0, -n2).sqrt(), + Complex64::from_polar(n, -f64::consts::FRAC_PI_4) + )); + } + } + + #[test] + fn test_cbrt() { + assert!(close(_0_0i.cbrt(), _0_0i)); + assert!(close(_1_0i.cbrt(), _1_0i)); + assert!(close( + Complex::new(-1.0, 0.0).cbrt(), + Complex::new(0.5, 0.75.sqrt()) + )); + assert!(close( + Complex::new(-1.0, -0.0).cbrt(), + Complex::new(0.5, -(0.75.sqrt())) + )); + assert!(close(_0_1i.cbrt(), Complex::new(0.75.sqrt(), 0.5))); + assert!(close(_0_1i.conj().cbrt(), Complex::new(0.75.sqrt(), -0.5))); + for &c in all_consts.iter() { + // cbrt(conj(z() = conj(cbrt(z)) + assert!(close(c.conj().cbrt(), c.cbrt().conj())); + // for this branch, -pi/3 <= arg(cbrt(z)) <= pi/3 + assert!( + -f64::consts::FRAC_PI_3 <= c.cbrt().arg() + && c.cbrt().arg() <= f64::consts::FRAC_PI_3 + ); + // cbrt(z) * cbrt(z) cbrt(z) = z + assert!(close(c.cbrt() * c.cbrt() * c.cbrt(), c)); + } + } + + #[test] + fn test_cbrt_real() { + for n in (0..100).map(f64::from) { + // ∛(n³ + 0i) = n + 0i + let n3 = n * n * n; + assert!(close( + Complex64::new(n3, 0.0).cbrt(), + Complex64::new(n, 0.0) + )); + // ∛(-n³ + 0i) = n e^(iπ/3) + assert!(close( + Complex64::new(-n3, 0.0).cbrt(), + Complex64::from_polar(n, f64::consts::FRAC_PI_3) + )); + // ∛(-n³ - 0i) = n e^(-iπ/3) + assert!(close( + Complex64::new(-n3, -0.0).cbrt(), + Complex64::from_polar(n, -f64::consts::FRAC_PI_3) + )); + } + } + + #[test] + fn test_cbrt_imag() { + for n in (0..100).map(f64::from) { + // ∛(0 + n³i) = n e^(iπ/6) + let n3 = n * n * n; + assert!(close( + Complex64::new(0.0, n3).cbrt(), + Complex64::from_polar(n, f64::consts::FRAC_PI_6) + )); + // ∛(0 - n³i) = n e^(-iπ/6) + assert!(close( + Complex64::new(0.0, -n3).cbrt(), + Complex64::from_polar(n, -f64::consts::FRAC_PI_6) + )); + } + } + + #[test] + fn test_sin() { + assert!(close(_0_0i.sin(), _0_0i)); + assert!(close(_1_0i.scale(f64::consts::PI * 2.0).sin(), _0_0i)); + assert!(close(_0_1i.sin(), _0_1i.scale(1.0.sinh()))); + for &c in all_consts.iter() { + // sin(conj(z)) = conj(sin(z)) + assert!(close(c.conj().sin(), c.sin().conj())); + // sin(-z) = -sin(z) + assert!(close(c.scale(-1.0).sin(), c.sin().scale(-1.0))); + } + } + + #[test] + fn test_cos() { + assert!(close(_0_0i.cos(), _1_0i)); + assert!(close(_1_0i.scale(f64::consts::PI * 2.0).cos(), _1_0i)); + assert!(close(_0_1i.cos(), _1_0i.scale(1.0.cosh()))); + for &c in all_consts.iter() { + // cos(conj(z)) = conj(cos(z)) + assert!(close(c.conj().cos(), c.cos().conj())); + // cos(-z) = cos(z) + assert!(close(c.scale(-1.0).cos(), c.cos())); + } + } + + #[test] + fn test_tan() { + assert!(close(_0_0i.tan(), _0_0i)); + assert!(close(_1_0i.scale(f64::consts::PI / 4.0).tan(), _1_0i)); + assert!(close(_1_0i.scale(f64::consts::PI).tan(), _0_0i)); + for &c in all_consts.iter() { + // tan(conj(z)) = conj(tan(z)) + assert!(close(c.conj().tan(), c.tan().conj())); + // tan(-z) = -tan(z) + assert!(close(c.scale(-1.0).tan(), c.tan().scale(-1.0))); + } + } + + #[test] + fn test_asin() { + assert!(close(_0_0i.asin(), _0_0i)); + assert!(close(_1_0i.asin(), _1_0i.scale(f64::consts::PI / 2.0))); + assert!(close( + _1_0i.scale(-1.0).asin(), + _1_0i.scale(-f64::consts::PI / 2.0) + )); + assert!(close(_0_1i.asin(), _0_1i.scale((1.0 + 2.0.sqrt()).ln()))); + for &c in all_consts.iter() { + // asin(conj(z)) = conj(asin(z)) + assert!(close(c.conj().asin(), c.asin().conj())); + // asin(-z) = -asin(z) + assert!(close(c.scale(-1.0).asin(), c.asin().scale(-1.0))); + // for this branch, -pi/2 <= asin(z).re <= pi/2 + assert!( + -f64::consts::PI / 2.0 <= c.asin().re && c.asin().re <= f64::consts::PI / 2.0 + ); + } + } + + #[test] + fn test_acos() { + assert!(close(_0_0i.acos(), _1_0i.scale(f64::consts::PI / 2.0))); + assert!(close(_1_0i.acos(), _0_0i)); + assert!(close( + _1_0i.scale(-1.0).acos(), + _1_0i.scale(f64::consts::PI) + )); + assert!(close( + _0_1i.acos(), + Complex::new(f64::consts::PI / 2.0, (2.0.sqrt() - 1.0).ln()) + )); + for &c in all_consts.iter() { + // acos(conj(z)) = conj(acos(z)) + assert!(close(c.conj().acos(), c.acos().conj())); + // for this branch, 0 <= acos(z).re <= pi + assert!(0.0 <= c.acos().re && c.acos().re <= f64::consts::PI); + } + } + + #[test] + fn test_atan() { + assert!(close(_0_0i.atan(), _0_0i)); + assert!(close(_1_0i.atan(), _1_0i.scale(f64::consts::PI / 4.0))); + assert!(close( + _1_0i.scale(-1.0).atan(), + _1_0i.scale(-f64::consts::PI / 4.0) + )); + assert!(close(_0_1i.atan(), Complex::new(0.0, f64::infinity()))); + for &c in all_consts.iter() { + // atan(conj(z)) = conj(atan(z)) + assert!(close(c.conj().atan(), c.atan().conj())); + // atan(-z) = -atan(z) + assert!(close(c.scale(-1.0).atan(), c.atan().scale(-1.0))); + // for this branch, -pi/2 <= atan(z).re <= pi/2 + assert!( + -f64::consts::PI / 2.0 <= c.atan().re && c.atan().re <= f64::consts::PI / 2.0 + ); + } + } + + #[test] + fn test_sinh() { + assert!(close(_0_0i.sinh(), _0_0i)); + assert!(close( + _1_0i.sinh(), + _1_0i.scale((f64::consts::E - 1.0 / f64::consts::E) / 2.0) + )); + assert!(close(_0_1i.sinh(), _0_1i.scale(1.0.sin()))); + for &c in all_consts.iter() { + // sinh(conj(z)) = conj(sinh(z)) + assert!(close(c.conj().sinh(), c.sinh().conj())); + // sinh(-z) = -sinh(z) + assert!(close(c.scale(-1.0).sinh(), c.sinh().scale(-1.0))); + } + } + + #[test] + fn test_cosh() { + assert!(close(_0_0i.cosh(), _1_0i)); + assert!(close( + _1_0i.cosh(), + _1_0i.scale((f64::consts::E + 1.0 / f64::consts::E) / 2.0) + )); + assert!(close(_0_1i.cosh(), _1_0i.scale(1.0.cos()))); + for &c in all_consts.iter() { + // cosh(conj(z)) = conj(cosh(z)) + assert!(close(c.conj().cosh(), c.cosh().conj())); + // cosh(-z) = cosh(z) + assert!(close(c.scale(-1.0).cosh(), c.cosh())); + } + } + + #[test] + fn test_tanh() { + assert!(close(_0_0i.tanh(), _0_0i)); + assert!(close( + _1_0i.tanh(), + _1_0i.scale((f64::consts::E.powi(2) - 1.0) / (f64::consts::E.powi(2) + 1.0)) + )); + assert!(close(_0_1i.tanh(), _0_1i.scale(1.0.tan()))); + for &c in all_consts.iter() { + // tanh(conj(z)) = conj(tanh(z)) + assert!(close(c.conj().tanh(), c.conj().tanh())); + // tanh(-z) = -tanh(z) + assert!(close(c.scale(-1.0).tanh(), c.tanh().scale(-1.0))); + } + } + + #[test] + fn test_asinh() { + assert!(close(_0_0i.asinh(), _0_0i)); + assert!(close(_1_0i.asinh(), _1_0i.scale(1.0 + 2.0.sqrt()).ln())); + assert!(close(_0_1i.asinh(), _0_1i.scale(f64::consts::PI / 2.0))); + assert!(close( + _0_1i.asinh().scale(-1.0), + _0_1i.scale(-f64::consts::PI / 2.0) + )); + for &c in all_consts.iter() { + // asinh(conj(z)) = conj(asinh(z)) + assert!(close(c.conj().asinh(), c.conj().asinh())); + // asinh(-z) = -asinh(z) + assert!(close(c.scale(-1.0).asinh(), c.asinh().scale(-1.0))); + // for this branch, -pi/2 <= asinh(z).im <= pi/2 + assert!( + -f64::consts::PI / 2.0 <= c.asinh().im && c.asinh().im <= f64::consts::PI / 2.0 + ); + } + } + + #[test] + fn test_acosh() { + assert!(close(_0_0i.acosh(), _0_1i.scale(f64::consts::PI / 2.0))); + assert!(close(_1_0i.acosh(), _0_0i)); + assert!(close( + _1_0i.scale(-1.0).acosh(), + _0_1i.scale(f64::consts::PI) + )); + for &c in all_consts.iter() { + // acosh(conj(z)) = conj(acosh(z)) + assert!(close(c.conj().acosh(), c.conj().acosh())); + // for this branch, -pi <= acosh(z).im <= pi and 0 <= acosh(z).re + assert!( + -f64::consts::PI <= c.acosh().im + && c.acosh().im <= f64::consts::PI + && 0.0 <= c.cosh().re + ); + } + } + + #[test] + fn test_atanh() { + assert!(close(_0_0i.atanh(), _0_0i)); + assert!(close(_0_1i.atanh(), _0_1i.scale(f64::consts::PI / 4.0))); + assert!(close(_1_0i.atanh(), Complex::new(f64::infinity(), 0.0))); + for &c in all_consts.iter() { + // atanh(conj(z)) = conj(atanh(z)) + assert!(close(c.conj().atanh(), c.conj().atanh())); + // atanh(-z) = -atanh(z) + assert!(close(c.scale(-1.0).atanh(), c.atanh().scale(-1.0))); + // for this branch, -pi/2 <= atanh(z).im <= pi/2 + assert!( + -f64::consts::PI / 2.0 <= c.atanh().im && c.atanh().im <= f64::consts::PI / 2.0 + ); + } + } + + #[test] + fn test_exp_ln() { + for &c in all_consts.iter() { + // e^ln(z) = z + assert!(close(c.ln().exp(), c)); + } + } + + #[test] + fn test_exp2_log() { + for &c in all_consts.iter() { + // 2^log2(z) = z + assert!(close(c.log2().exp2(), c)); + } + } + + #[test] + fn test_trig_to_hyperbolic() { + for &c in all_consts.iter() { + // sin(iz) = i sinh(z) + assert!(close((_0_1i * c).sin(), _0_1i * c.sinh())); + // cos(iz) = cosh(z) + assert!(close((_0_1i * c).cos(), c.cosh())); + // tan(iz) = i tanh(z) + assert!(close((_0_1i * c).tan(), _0_1i * c.tanh())); + } + } + + #[test] + fn test_trig_identities() { + for &c in all_consts.iter() { + // tan(z) = sin(z)/cos(z) + assert!(close(c.tan(), c.sin() / c.cos())); + // sin(z)^2 + cos(z)^2 = 1 + assert!(close(c.sin() * c.sin() + c.cos() * c.cos(), _1_0i)); + + // sin(asin(z)) = z + assert!(close(c.asin().sin(), c)); + // cos(acos(z)) = z + assert!(close(c.acos().cos(), c)); + // tan(atan(z)) = z + // i and -i are branch points + if c != _0_1i && c != _0_1i.scale(-1.0) { + assert!(close(c.atan().tan(), c)); + } + + // sin(z) = (e^(iz) - e^(-iz))/(2i) + assert!(close( + ((_0_1i * c).exp() - (_0_1i * c).exp().inv()) / _0_1i.scale(2.0), + c.sin() + )); + // cos(z) = (e^(iz) + e^(-iz))/2 + assert!(close( + ((_0_1i * c).exp() + (_0_1i * c).exp().inv()).unscale(2.0), + c.cos() + )); + // tan(z) = i (1 - e^(2iz))/(1 + e^(2iz)) + assert!(close( + _0_1i * (_1_0i - (_0_1i * c).scale(2.0).exp()) + / (_1_0i + (_0_1i * c).scale(2.0).exp()), + c.tan() + )); + } + } + + #[test] + fn test_hyperbolic_identites() { + for &c in all_consts.iter() { + // tanh(z) = sinh(z)/cosh(z) + assert!(close(c.tanh(), c.sinh() / c.cosh())); + // cosh(z)^2 - sinh(z)^2 = 1 + assert!(close(c.cosh() * c.cosh() - c.sinh() * c.sinh(), _1_0i)); + + // sinh(asinh(z)) = z + assert!(close(c.asinh().sinh(), c)); + // cosh(acosh(z)) = z + assert!(close(c.acosh().cosh(), c)); + // tanh(atanh(z)) = z + // 1 and -1 are branch points + if c != _1_0i && c != _1_0i.scale(-1.0) { + assert!(close(c.atanh().tanh(), c)); + } + + // sinh(z) = (e^z - e^(-z))/2 + assert!(close((c.exp() - c.exp().inv()).unscale(2.0), c.sinh())); + // cosh(z) = (e^z + e^(-z))/2 + assert!(close((c.exp() + c.exp().inv()).unscale(2.0), c.cosh())); + // tanh(z) = ( e^(2z) - 1)/(e^(2z) + 1) + assert!(close( + (c.scale(2.0).exp() - _1_0i) / (c.scale(2.0).exp() + _1_0i), + c.tanh() + )); + } + } + } + + // Test both a + b and a += b + macro_rules! test_a_op_b { + ($a:ident + $b:expr, $answer:expr) => { + assert_eq!($a + $b, $answer); + assert_eq!( + { + let mut x = $a; + x += $b; + x + }, + $answer + ); + }; + ($a:ident - $b:expr, $answer:expr) => { + assert_eq!($a - $b, $answer); + assert_eq!( + { + let mut x = $a; + x -= $b; + x + }, + $answer + ); + }; + ($a:ident * $b:expr, $answer:expr) => { + assert_eq!($a * $b, $answer); + assert_eq!( + { + let mut x = $a; + x *= $b; + x + }, + $answer + ); + }; + ($a:ident / $b:expr, $answer:expr) => { + assert_eq!($a / $b, $answer); + assert_eq!( + { + let mut x = $a; + x /= $b; + x + }, + $answer + ); + }; + ($a:ident % $b:expr, $answer:expr) => { + assert_eq!($a % $b, $answer); + assert_eq!( + { + let mut x = $a; + x %= $b; + x + }, + $answer + ); + }; + } + + // Test both a + b and a + &b + macro_rules! test_op { + ($a:ident $op:tt $b:expr, $answer:expr) => { + test_a_op_b!($a $op $b, $answer); + test_a_op_b!($a $op &$b, $answer); + }; + } + + mod complex_arithmetic { + use super::{_05_05i, _0_0i, _0_1i, _1_0i, _1_1i, _4_2i, _neg1_1i, all_consts}; + use num_traits::{MulAdd, MulAddAssign, Zero}; + + #[test] + fn test_add() { + test_op!(_05_05i + _05_05i, _1_1i); + test_op!(_0_1i + _1_0i, _1_1i); + test_op!(_1_0i + _neg1_1i, _0_1i); + + for &c in all_consts.iter() { + test_op!(_0_0i + c, c); + test_op!(c + _0_0i, c); + } + } + + #[test] + fn test_sub() { + test_op!(_05_05i - _05_05i, _0_0i); + test_op!(_0_1i - _1_0i, _neg1_1i); + test_op!(_0_1i - _neg1_1i, _1_0i); + + for &c in all_consts.iter() { + test_op!(c - _0_0i, c); + test_op!(c - c, _0_0i); + } + } + + #[test] + fn test_mul() { + test_op!(_05_05i * _05_05i, _0_1i.unscale(2.0)); + test_op!(_1_1i * _0_1i, _neg1_1i); + + // i^2 & i^4 + test_op!(_0_1i * _0_1i, -_1_0i); + assert_eq!(_0_1i * _0_1i * _0_1i * _0_1i, _1_0i); + + for &c in all_consts.iter() { + test_op!(c * _1_0i, c); + test_op!(_1_0i * c, c); + } + } + + #[test] + #[cfg(any(feature = "std", feature = "libm"))] + fn test_mul_add_float() { + assert_eq!(_05_05i.mul_add(_05_05i, _0_0i), _05_05i * _05_05i + _0_0i); + assert_eq!(_05_05i * _05_05i + _0_0i, _05_05i.mul_add(_05_05i, _0_0i)); + assert_eq!(_0_1i.mul_add(_0_1i, _0_1i), _neg1_1i); + assert_eq!(_1_0i.mul_add(_1_0i, _1_0i), _1_0i * _1_0i + _1_0i); + assert_eq!(_1_0i * _1_0i + _1_0i, _1_0i.mul_add(_1_0i, _1_0i)); + + let mut x = _1_0i; + x.mul_add_assign(_1_0i, _1_0i); + assert_eq!(x, _1_0i * _1_0i + _1_0i); + + for &a in &all_consts { + for &b in &all_consts { + for &c in &all_consts { + let abc = a * b + c; + assert_eq!(a.mul_add(b, c), abc); + let mut x = a; + x.mul_add_assign(b, c); + assert_eq!(x, abc); + } + } + } + } + + #[test] + fn test_mul_add() { + use super::Complex; + const _0_0i: Complex = Complex { re: 0, im: 0 }; + const _1_0i: Complex = Complex { re: 1, im: 0 }; + const _1_1i: Complex = Complex { re: 1, im: 1 }; + const _0_1i: Complex = Complex { re: 0, im: 1 }; + const _neg1_1i: Complex = Complex { re: -1, im: 1 }; + const all_consts: [Complex; 5] = [_0_0i, _1_0i, _1_1i, _0_1i, _neg1_1i]; + + assert_eq!(_1_0i.mul_add(_1_0i, _0_0i), _1_0i * _1_0i + _0_0i); + assert_eq!(_1_0i * _1_0i + _0_0i, _1_0i.mul_add(_1_0i, _0_0i)); + assert_eq!(_0_1i.mul_add(_0_1i, _0_1i), _neg1_1i); + assert_eq!(_1_0i.mul_add(_1_0i, _1_0i), _1_0i * _1_0i + _1_0i); + assert_eq!(_1_0i * _1_0i + _1_0i, _1_0i.mul_add(_1_0i, _1_0i)); + + let mut x = _1_0i; + x.mul_add_assign(_1_0i, _1_0i); + assert_eq!(x, _1_0i * _1_0i + _1_0i); + + for &a in &all_consts { + for &b in &all_consts { + for &c in &all_consts { + let abc = a * b + c; + assert_eq!(a.mul_add(b, c), abc); + let mut x = a; + x.mul_add_assign(b, c); + assert_eq!(x, abc); + } + } + } + } + + #[test] + fn test_div() { + test_op!(_neg1_1i / _0_1i, _1_1i); + for &c in all_consts.iter() { + if c != Zero::zero() { + test_op!(c / c, _1_0i); + } + } + } + + #[test] + fn test_rem() { + test_op!(_neg1_1i % _0_1i, _0_0i); + test_op!(_4_2i % _0_1i, _0_0i); + test_op!(_05_05i % _0_1i, _05_05i); + test_op!(_05_05i % _1_1i, _05_05i); + assert_eq!((_4_2i + _05_05i) % _0_1i, _05_05i); + assert_eq!((_4_2i + _05_05i) % _1_1i, _05_05i); + } + + #[test] + fn test_neg() { + assert_eq!(-_1_0i + _0_1i, _neg1_1i); + assert_eq!((-_0_1i) * _0_1i, _1_0i); + for &c in all_consts.iter() { + assert_eq!(-(-c), c); + } + } + } + + mod real_arithmetic { + use super::super::Complex; + use super::{_4_2i, _neg1_1i}; + + #[test] + fn test_add() { + test_op!(_4_2i + 0.5, Complex::new(4.5, 2.0)); + assert_eq!(0.5 + _4_2i, Complex::new(4.5, 2.0)); + } + + #[test] + fn test_sub() { + test_op!(_4_2i - 0.5, Complex::new(3.5, 2.0)); + assert_eq!(0.5 - _4_2i, Complex::new(-3.5, -2.0)); + } + + #[test] + fn test_mul() { + assert_eq!(_4_2i * 0.5, Complex::new(2.0, 1.0)); + assert_eq!(0.5 * _4_2i, Complex::new(2.0, 1.0)); + } + + #[test] + fn test_div() { + assert_eq!(_4_2i / 0.5, Complex::new(8.0, 4.0)); + assert_eq!(0.5 / _4_2i, Complex::new(0.1, -0.05)); + } + + #[test] + fn test_rem() { + assert_eq!(_4_2i % 2.0, Complex::new(0.0, 0.0)); + assert_eq!(_4_2i % 3.0, Complex::new(1.0, 2.0)); + assert_eq!(3.0 % _4_2i, Complex::new(3.0, 0.0)); + assert_eq!(_neg1_1i % 2.0, _neg1_1i); + assert_eq!(-_4_2i % 3.0, Complex::new(-1.0, -2.0)); + } + + #[test] + fn test_div_rem_gaussian() { + // These would overflow with `norm_sqr` division. + let max = Complex::new(255u8, 255u8); + assert_eq!(max / 200, Complex::new(1, 1)); + assert_eq!(max % 200, Complex::new(55, 55)); + } + } + + #[test] + fn test_to_string() { + fn test(c: Complex64, s: String) { + assert_eq!(c.to_string(), s); + } + test(_0_0i, "0+0i".to_string()); + test(_1_0i, "1+0i".to_string()); + test(_0_1i, "0+1i".to_string()); + test(_1_1i, "1+1i".to_string()); + test(_neg1_1i, "-1+1i".to_string()); + test(-_neg1_1i, "1-1i".to_string()); + test(_05_05i, "0.5+0.5i".to_string()); + } + + #[test] + fn test_string_formatting() { + let a = Complex::new(1.23456, 123.456); + assert_eq!(format!("{}", a), "1.23456+123.456i"); + assert_eq!(format!("{:.2}", a), "1.23+123.46i"); + assert_eq!(format!("{:.2e}", a), "1.23e0+1.23e2i"); + assert_eq!(format!("{:+.2E}", a), "+1.23E0+1.23E2i"); + #[cfg(feature = "std")] + assert_eq!(format!("{:+20.2E}", a), " +1.23E0+1.23E2i"); + + let b = Complex::new(0x80, 0xff); + assert_eq!(format!("{:X}", b), "80+FFi"); + assert_eq!(format!("{:#x}", b), "0x80+0xffi"); + assert_eq!(format!("{:+#b}", b), "+0b10000000+0b11111111i"); + assert_eq!(format!("{:+#o}", b), "+0o200+0o377i"); + #[cfg(feature = "std")] + assert_eq!(format!("{:+#16o}", b), " +0o200+0o377i"); + + let c = Complex::new(-10, -10000); + assert_eq!(format!("{}", c), "-10-10000i"); + #[cfg(feature = "std")] + assert_eq!(format!("{:16}", c), " -10-10000i"); + } + + #[test] + fn test_hash() { + let a = Complex::new(0i32, 0i32); + let b = Complex::new(1i32, 0i32); + let c = Complex::new(0i32, 1i32); + assert!(crate::hash(&a) != crate::hash(&b)); + assert!(crate::hash(&b) != crate::hash(&c)); + assert!(crate::hash(&c) != crate::hash(&a)); + } + + #[test] + fn test_hashset() { + use std::collections::HashSet; + let a = Complex::new(0i32, 0i32); + let b = Complex::new(1i32, 0i32); + let c = Complex::new(0i32, 1i32); + + let set: HashSet<_> = [a, b, c].iter().cloned().collect(); + assert!(set.contains(&a)); + assert!(set.contains(&b)); + assert!(set.contains(&c)); + assert!(!set.contains(&(a + b + c))); + } + + #[test] + fn test_is_nan() { + assert!(!_1_1i.is_nan()); + let a = Complex::new(f64::NAN, f64::NAN); + assert!(a.is_nan()); + } + + #[test] + fn test_is_nan_special_cases() { + let a = Complex::new(0f64, f64::NAN); + let b = Complex::new(f64::NAN, 0f64); + assert!(a.is_nan()); + assert!(b.is_nan()); + } + + #[test] + fn test_is_infinite() { + let a = Complex::new(2f64, f64::INFINITY); + assert!(a.is_infinite()); + } + + #[test] + fn test_is_finite() { + assert!(_1_1i.is_finite()) + } + + #[test] + fn test_is_normal() { + let a = Complex::new(0f64, f64::NAN); + let b = Complex::new(2f64, f64::INFINITY); + assert!(!a.is_normal()); + assert!(!b.is_normal()); + assert!(_1_1i.is_normal()); + } + + #[test] + fn test_from_str() { + fn test(z: Complex64, s: &str) { + assert_eq!(FromStr::from_str(s), Ok(z)); + } + test(_0_0i, "0 + 0i"); + test(_0_0i, "0+0j"); + test(_0_0i, "0 - 0j"); + test(_0_0i, "0-0i"); + test(_0_0i, "0i + 0"); + test(_0_0i, "0"); + test(_0_0i, "-0"); + test(_0_0i, "0i"); + test(_0_0i, "0j"); + test(_0_0i, "+0j"); + test(_0_0i, "-0i"); + + test(_1_0i, "1 + 0i"); + test(_1_0i, "1+0j"); + test(_1_0i, "1 - 0j"); + test(_1_0i, "+1-0i"); + test(_1_0i, "-0j+1"); + test(_1_0i, "1"); + + test(_1_1i, "1 + i"); + test(_1_1i, "1+j"); + test(_1_1i, "1 + 1j"); + test(_1_1i, "1+1i"); + test(_1_1i, "i + 1"); + test(_1_1i, "1i+1"); + test(_1_1i, "+j+1"); + + test(_0_1i, "0 + i"); + test(_0_1i, "0+j"); + test(_0_1i, "-0 + j"); + test(_0_1i, "-0+i"); + test(_0_1i, "0 + 1i"); + test(_0_1i, "0+1j"); + test(_0_1i, "-0 + 1j"); + test(_0_1i, "-0+1i"); + test(_0_1i, "j + 0"); + test(_0_1i, "i"); + test(_0_1i, "j"); + test(_0_1i, "1j"); + + test(_neg1_1i, "-1 + i"); + test(_neg1_1i, "-1+j"); + test(_neg1_1i, "-1 + 1j"); + test(_neg1_1i, "-1+1i"); + test(_neg1_1i, "1i-1"); + test(_neg1_1i, "j + -1"); + + test(_05_05i, "0.5 + 0.5i"); + test(_05_05i, "0.5+0.5j"); + test(_05_05i, "5e-1+0.5j"); + test(_05_05i, "5E-1 + 0.5j"); + test(_05_05i, "5E-1i + 0.5"); + test(_05_05i, "0.05e+1j + 50E-2"); + } + + #[test] + fn test_from_str_radix() { + fn test(z: Complex64, s: &str, radix: u32) { + let res: Result::FromStrRadixErr> = + Num::from_str_radix(s, radix); + assert_eq!(res.unwrap(), z) + } + test(_4_2i, "4+2i", 10); + test(Complex::new(15.0, 32.0), "F+20i", 16); + test(Complex::new(15.0, 32.0), "1111+100000i", 2); + test(Complex::new(-15.0, -32.0), "-F-20i", 16); + test(Complex::new(-15.0, -32.0), "-1111-100000i", 2); + + fn test_error(s: &str, radix: u32) -> ParseComplexError<::FromStrRadixErr> { + let res = Complex64::from_str_radix(s, radix); + + res.expect_err(&format!("Expected failure on input {:?}", s)) + } + + let err = test_error("1ii", 19); + if let ComplexErrorKind::UnsupportedRadix = err.kind { + /* pass */ + } else { + panic!("Expected failure on invalid radix, got {:?}", err); + } + + let err = test_error("1 + 0", 16); + if let ComplexErrorKind::ExprError = err.kind { + /* pass */ + } else { + panic!("Expected failure on expr error, got {:?}", err); + } + } + + #[test] + #[should_panic(expected = "radix is too high")] + fn test_from_str_radix_fail() { + // ensure we preserve the underlying panic on radix > 36 + let _complex = Complex64::from_str_radix("1", 37); + } + + #[test] + fn test_from_str_fail() { + fn test(s: &str) { + let complex: Result = FromStr::from_str(s); + assert!( + complex.is_err(), + "complex {:?} -> {:?} should be an error", + s, + complex + ); + } + test("foo"); + test("6E"); + test("0 + 2.718"); + test("1 - -2i"); + test("314e-2ij"); + test("4.3j - i"); + test("1i - 2i"); + test("+ 1 - 3.0i"); + } + + #[test] + fn test_sum() { + let v = vec![_0_1i, _1_0i]; + assert_eq!(v.iter().sum::(), _1_1i); + assert_eq!(v.into_iter().sum::(), _1_1i); + } + + #[test] + fn test_prod() { + let v = vec![_0_1i, _1_0i]; + assert_eq!(v.iter().product::(), _0_1i); + assert_eq!(v.into_iter().product::(), _0_1i); + } + + #[test] + fn test_zero() { + let zero = Complex64::zero(); + assert!(zero.is_zero()); + + let mut c = Complex::new(1.23, 4.56); + assert!(!c.is_zero()); + assert_eq!(c + zero, c); + + c.set_zero(); + assert!(c.is_zero()); + } + + #[test] + fn test_one() { + let one = Complex64::one(); + assert!(one.is_one()); + + let mut c = Complex::new(1.23, 4.56); + assert!(!c.is_one()); + assert_eq!(c * one, c); + + c.set_one(); + assert!(c.is_one()); + } + + #[test] + #[allow(clippy::float_cmp)] + fn test_const() { + const R: f64 = 12.3; + const I: f64 = -4.5; + const C: Complex64 = Complex::new(R, I); + + assert_eq!(C.re, 12.3); + assert_eq!(C.im, -4.5); + } +} diff --git a/src/pow.rs b/src/pow.rs new file mode 100644 index 0000000..569f01d --- /dev/null +++ b/src/pow.rs @@ -0,0 +1,186 @@ +use super::Complex; + +use core::ops::Neg; +#[cfg(any(feature = "std", feature = "libm"))] +use num_traits::Float; +use num_traits::{Num, One, Pow}; + +macro_rules! pow_impl { + ($U:ty, $S:ty) => { + impl<'a, T: Clone + Num> Pow<$U> for &'a Complex { + type Output = Complex; + + #[inline] + fn pow(self, mut exp: $U) -> Self::Output { + if exp == 0 { + return Complex::one(); + } + let mut base = self.clone(); + + while exp & 1 == 0 { + base = base.clone() * base; + exp >>= 1; + } + + if exp == 1 { + return base; + } + + let mut acc = base.clone(); + while exp > 1 { + exp >>= 1; + base = base.clone() * base; + if exp & 1 == 1 { + acc = acc * base.clone(); + } + } + acc + } + } + + impl<'a, 'b, T: Clone + Num> Pow<&'b $U> for &'a Complex { + type Output = Complex; + + #[inline] + fn pow(self, exp: &$U) -> Self::Output { + self.pow(*exp) + } + } + + impl<'a, T: Clone + Num + Neg> Pow<$S> for &'a Complex { + type Output = Complex; + + #[inline] + fn pow(self, exp: $S) -> Self::Output { + if exp < 0 { + Pow::pow(&self.inv(), exp.wrapping_neg() as $U) + } else { + Pow::pow(self, exp as $U) + } + } + } + + impl<'a, 'b, T: Clone + Num + Neg> Pow<&'b $S> for &'a Complex { + type Output = Complex; + + #[inline] + fn pow(self, exp: &$S) -> Self::Output { + self.pow(*exp) + } + } + }; +} + +pow_impl!(u8, i8); +pow_impl!(u16, i16); +pow_impl!(u32, i32); +pow_impl!(u64, i64); +pow_impl!(usize, isize); +pow_impl!(u128, i128); + +// Note: we can't add `impl Pow for Complex` because new blanket impls are a +// breaking change. Someone could already have their own `F` and `impl Pow for Complex` +// which would conflict. We can't even do this in a new semantic version, because we have to +// gate it on the "std" feature, and features can't add breaking changes either. + +macro_rules! powf_impl { + ($F:ty) => { + #[cfg(any(feature = "std", feature = "libm"))] + impl<'a, T: Float> Pow<$F> for &'a Complex + where + $F: Into, + { + type Output = Complex; + + #[inline] + fn pow(self, exp: $F) -> Self::Output { + self.powf(exp.into()) + } + } + + #[cfg(any(feature = "std", feature = "libm"))] + impl<'a, 'b, T: Float> Pow<&'b $F> for &'a Complex + where + $F: Into, + { + type Output = Complex; + + #[inline] + fn pow(self, &exp: &$F) -> Self::Output { + self.powf(exp.into()) + } + } + + #[cfg(any(feature = "std", feature = "libm"))] + impl Pow<$F> for Complex + where + $F: Into, + { + type Output = Complex; + + #[inline] + fn pow(self, exp: $F) -> Self::Output { + self.powf(exp.into()) + } + } + + #[cfg(any(feature = "std", feature = "libm"))] + impl<'b, T: Float> Pow<&'b $F> for Complex + where + $F: Into, + { + type Output = Complex; + + #[inline] + fn pow(self, &exp: &$F) -> Self::Output { + self.powf(exp.into()) + } + } + }; +} + +powf_impl!(f32); +powf_impl!(f64); + +// These blanket impls are OK, because both the target type and the trait parameter would be +// foreign to anyone else trying to implement something that would overlap, raising E0117. + +#[cfg(any(feature = "std", feature = "libm"))] +impl<'a, T: Float> Pow> for &'a Complex { + type Output = Complex; + + #[inline] + fn pow(self, exp: Complex) -> Self::Output { + self.powc(exp) + } +} + +#[cfg(any(feature = "std", feature = "libm"))] +impl<'a, 'b, T: Float> Pow<&'b Complex> for &'a Complex { + type Output = Complex; + + #[inline] + fn pow(self, &exp: &'b Complex) -> Self::Output { + self.powc(exp) + } +} + +#[cfg(any(feature = "std", feature = "libm"))] +impl Pow> for Complex { + type Output = Complex; + + #[inline] + fn pow(self, exp: Complex) -> Self::Output { + self.powc(exp) + } +} + +#[cfg(any(feature = "std", feature = "libm"))] +impl<'b, T: Float> Pow<&'b Complex> for Complex { + type Output = Complex; + + #[inline] + fn pow(self, &exp: &'b Complex) -> Self::Output { + self.powc(exp) + } +} -- 2.7.4