#include <math.h>
    #include <errno.h>
    int findRoot(double a, double b, double c, double *root1, double *root2) {
        double discriminant, q;
        if (a == 0.0) {
            return EDOM;
        }
        discriminant = b*b-4*a*c;
        if (discriminant < 0) {
            return EDOM;
        }
        q = -(b+copysign(sqrt(discriminant), b)) / 2;
        *root1 = q/a;
        *root2 = c/q;
        return 0;
    }

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

  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));
    }
  }
  

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?

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.

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.

"Sometimes, the elegant implementation is just a function. Not a method. Not a class. Not a framework. Just a function."

https://twitter.com/#!/ID_AA_Carmack/status/5351230045120102...

I always found it absurd that some languages force you to write "Math.cos(...)".

I disagree. Namespaces are a good thing. The "Foo dot bar" syntax is just commonly associated with class member access. I think it's a lot better to have to explicitly refer to where a function comes from (as long as it doesn't devolve into retarded Java-style verbosity, ie "Foo.bar.baz.spam.eggs..."). Though it's nice to have a way to omit the namespace, ie C++'s "using namespace Foo" or Python's "from Foo import bar". I just think that the sane default is to use explicit namespaces.

But in this case "Math" isn't even a namespace! It's a class of which cos() is a static member. http://docs.oracle.com/javase/6/docs/api/java/lang/Math.html...

EDIT: I see below that Java has added "import static", but at that point you're in the situation where you have two things that are almost the same but not quite, which is not a particularly clean design in my opinion.

And Math is a particularly special case of this. Ponder this: have you ever had to say

    import java.lang.Math

?

http://docs.oracle.com/javase/1.5.0/docs/guide/language/stat...

From the linked docs: So when should you use static import? Very sparingly!

To me, it sort of confirms the point that these inelegant features are necessary because the language doesn't support fundamental primitives like plain functions and named constants without a superfluous object definition.

static class Math { static int Pow(int a,int b) {...} }

    constexpr tuple<double, double> Roots(double a, double b, double c)
    {
        return (a == 0.0) ? throw invalid_argument("a has value 0.0!") :
            tuple<double, double>{
                (-b + (sqrt((b * b) - 4 * a * c))) / 2*a,
                (-b - (sqrt((b * b) - 4 * a * c))) / 2*a
            };
    }

  return /* ... */ throw invalid_argument(""); 

[1] http://akrzemi1.wordpress.com/2011/05/06/compile-time-comput...

  return a == ALPHA ? x :
         a == BETA ? y :
         a == GAMMA ? z :
         throw ...;

  #include <iostream>
  #include <string>
  #include <cmath>
  
  using std::string;
  using std::cout;
  using std::endl;
  
  static bool findRoots(
     double a,
     double b,
     double c,
     double &rootOne,
     double &rootTwo,
     string &errMsg
  )
  {
     if(a == 0.0)
     {
        errMsg = "Parameter 'a' cannot be zero";
        return(false);
     }
  
     double discriminant = b * b - 4 * a * c;
  
     if(discriminant < 0.0)
     {
        errMsg = "Discriminant cannot be less than zero";
        return(false);
     }
  
     double dSqrt = sqrt(discriminant);
     double twoA  = 2 * a;
  
     rootOne = (-b + dSqrt) / twoA;
     rootTwo = (-b - dSqrt) / twoA;
  
     return(true);
  }
  
  int main()
  {
     double a[] = { 0.0, 2.0, 2.0 };
     double b[] = { 1.0, 3.0, 1.0 };
     double c[] = { 2.0, 2.0, 0.0 };
  
     double rootOne;
     double rootTwo;
     string errMsg;
  
     for(int n = sizeof(a) / sizeof(a[0]), i = 0; i < n; i ++)
     {
        cout << "a: " << a[i] << ", b: " << b[i] << ", c: " << c[i] << endl;
  
        if(findRoots(a[i], b[i], c[i], rootOne, rootTwo, errMsg))
        {
           cout << rootOne << " " << rootTwo << endl;
        }
        else
        {
           cout << errMsg << endl;
        }
  
        cout << endl;
     }
  
     return(0);
  }

     def findRoot(p:(Double,Double,Double)) = {
       assert(p._1 != 0)
       val n = sqrt(p._2*p._2-4*p._1*p._3)
       val d = 2*p._1
       ( (-p._2+n)/d, (-p._2-n)/d )
     }
      ----
      scala> findRoot((6,-6,-36))
      res0: (Double, Double) = (3.0,-2.0)

why should they?

In a larger context, you might want to actually have an object/struct that represents a quadratic equation. The weird thing though was having an object for a quadratic equation and not having another object for its root (or more likely, either have a binomial class or just one generic polynomial class). I'm trying to think of a context where that makes any sense at all... and not finding it.

That's usually a result of people coming from Java doing C++ work.