Hacker News Reader

Comment by ansgri | Hacker News Reader

Make it a template (parameterized with coefficient and solution type) and you'll get very generic and efficient solution across all mentioned 'dialects' of c(++). It has the added bonus that this routine can be converted to vector form almost without breaking the interface: just make a, b, c pointers and specify length. And you can specialize templates when you'll want to rewrite float32_t-version in assembly.

Numeric code is really inappropriate for demonstration of new C++ features, as good (fast!) numeric code is most often C with templates.

View on HN

A slight variation with neither error codes nor exceptions.

  template< typename T >
  auto findRoots(T a, T b, T c) -> pair<maybe<complex<T>>, maybe<complex<T>>> {
    typedef complex<T>       C;
    typedef maybe<C>         Root;
    typedef pair<Root, Root> Roots;

    if (a == 0) {
      if (b == 0) 
        return Roots(nothing,      nothing);
      else
        return Roots(Root(-c / b), nothing);
    } else {
      auto d     = C(b*b - a*c*4);
      auto two_a = a*2;
      return Roots(Root((-b + sqrt(d)) / two_a), 
                   Root((-b - sqrt(d)) / two_a));
    }
  }

maybe is simply a templatized std::pair wrapper that lets you check if the returned value is valid (either via maybe.valid() or if(maybe)). All imaginary results are handled automatically via std::complex. nothing decomposes into an invalid maybe value of the appropriate type.

The minor nit first: the original specification did not support complex numbers, so adding them here is even less relevant than my return codes.

More importantly, your implementation is numerically less correct than my version, which itself is less correct than it should be.

Consider when b is near sqrt(d). In that case, b-sqrt(d) can lose precision. That's why I used the copysign function, so that I'm always adding two numbers of equal sign and size.

Mine is incorrect if b2 is near the size of 4ac since there too b2-4ac loses precision. Ideally this intermediate result should be done in quad precision if the input is in double.

Question: How does C++ handle copysign() (vs. copysignf() for floats), and support templates which want to use an higher precision intermediate value?

My reply wasn't meant as a numerically superior solution, simply as a fairly simple alternative without error codes...

As for "original specification", I'm not sure why you even bring that up. complex is certainly in the C++11 specification, against which my snippet is compliant.

Oh, I see I didn't explain well enough. I'm curious about how one would write the "numerically superior solution." That is, does C++11 have a templated copysign function which takes different numeric types (at least float, double, and quad)? If not, then there's some unneeded type promotion (or downcast) going on.

And if the template uses a float, how do I get the type with double precision to use as my intermediate? (And the same where the template uses a double and I want the intermediate to use a quad.)

There must surely be a way to handle these, but I haven't done C++ programming for over a decade and I don't know the modern way of doing things.

As regards "original specification" - I mean the original article from feabhas.com, which has since disappeared. As I recall, it only supported real roots, and not imaginary ones.

BTW, _Complex is also in C99.

Reply on news.ycombinator.com