Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Undue overflow or underflow on complex division #23

Closed
hydeik opened this issue May 17, 2018 · 6 comments
Closed

Undue overflow or underflow on complex division #23

hydeik opened this issue May 17, 2018 · 6 comments

Comments

@hydeik
Copy link

hydeik commented May 17, 2018

Complex division c1 / c2 might return an incorrect value when c2.norm_sqr() overflows or underflows.
For instance,

  extern crate num;
  use num::Complex;

  fn main() {
      let z1 = 1_f64 / Complex::new(1e300_f64, 1e300_f64);
      println!("z1 = {:e}", z1);
      // (expected output) z1 = 5e-301+5e-301i
      // (actual output)   z1 = 0e0+0e0i

      let z2 = 1_f64 / Complex::new(1e-301_f64, 1e-301_f64);
      println!("z2 = {:e}", z2);
      // (expected output) z2 = 5e300+5e300i
      // (actual output)   z2 = inf-infi
  }
@cuviper cuviper changed the title Undew overflow or underflow on complex division Undue overflow or underflow on complex division May 17, 2018
@cuviper
Copy link
Member

cuviper commented May 17, 2018

Yes, we use the naive algebraic formula, and it does have this downside. But the implementation of Div is in very generic code, for any T: Num + Clone, and this doesn't give us a lot to work with.

Also, the alternatives I see would fare poorly with non-float types or otherwise limited precision. For instance, with Gaussian integers like Complex<i32>, 25/(3+4i) is exactly 3-4i. But any division strategy that involves scaling will probably get a very different answer with floored integer division.

When Rust gets stable specialization, maybe we could implement a different Div just for T: Float. Or we could now implement an alternate method on impl<T: Float> Complex<T> -- div_slow or something -- it just wouldn't affect the / operator.

@cuviper
Copy link
Member

cuviper commented May 17, 2018

I know there are more sophisticated algorithms out there, but this is pretty accurate for your examples:

fn finv<T: Float>(c: Complex<T>) -> Complex<T> {
    let norm = c.norm();
    c.conj() / norm / norm
}

fn fdiv<T: Float>(a: Complex<T>, b: Complex<T>) -> Complex<T> {
    a * finv(b)
}

@hydeik
Copy link
Author

hydeik commented May 17, 2018

Thank you for detailed explanation.

I understand the reason why the current implementation uses the naive algebraic formula.
It would be nice if alternative method for T: Float is provided immediately in the library so that we could avoid this problem.
(fdiv seems to be good for method name.)

Of course, implementing specialized Div for T: Float would be more elegant solution after Rust stabilized specialization.

@benruijl
Copy link

Could you add the fdiv option for now? I am using this library for scientific computing and these overflows are causing issues with our numerical integrations.

@cuviper
Copy link
Member

cuviper commented Nov 15, 2018

@benruijl -- Have you tested whether the fdiv I gave above does work better for you? I can add them as methods, but they should work fine as the standalone functions given.

@benruijl
Copy link

@cupiver, I did. Your snippet works for me and I would like to see it included in the library.

cuviper added a commit to cuviper/num-complex that referenced this issue Feb 22, 2019
These use the floating point `norm()` to avoid cases where the more
generic `inv()` and `div()` might overflow in `norm_sqr()`.

Closes rust-num#23.
bors bot added a commit that referenced this issue Mar 28, 2019
41: Add `finv()` and `fdiv()` r=cuviper a=cuviper

These use the floating point `norm()` to avoid cases where the more
generic `inv()` and `div()` might overflow in `norm_sqr()`.

Closes #23.
cc @benruijl

Co-authored-by: Josh Stone <cuviper@gmail.com>
bors bot added a commit that referenced this issue Mar 28, 2019
41: Add `finv()` and `fdiv()` r=cuviper a=cuviper

These use the floating point `norm()` to avoid cases where the more
generic `inv()` and `div()` might overflow in `norm_sqr()`.

Closes #23.
cc @benruijl

Co-authored-by: Josh Stone <cuviper@gmail.com>
@bors bors bot closed this as completed in #41 Mar 28, 2019
adamnemecek pushed a commit to adamnemecek/num-complex that referenced this issue Mar 31, 2019
These use the floating point `norm()` to avoid cases where the more
generic `inv()` and `div()` might overflow in `norm_sqr()`.

Closes rust-num#23.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants