Solve quadratic equation using only 3 variables and arithmetic

Question

Write function which solves equation of degree less than or equal to 2 (ax^2 + bx + c = 0).

Input

Function should get 3 parameters - floating point numbers which are coefficients of this equation. Coefficients may be zero (that means function must able to solve bx + c = 0 and c = 0 too).
Coefficients are between -1000 and 1000 with 3 significant decimal digits but maximum 3 digits after comma, for example, 999. or 3.14 or 0.233 but not 1001 or 0.0123.

Rules

• You may use only + - * / % & | xor ~ && || ! < > >= <= == != << >> operators
• No built-in functions (that means no square root or power or similar functions)
• No variables except 3 function parameters
• You can't use complicated expressions - only statements of the form a = b @ c or a = @ b are allowed (@ - any allowed operator and a, b, c - any variables, may be same or constants)
• No loops, recursion, goto's, custom types
• Casts and type conversion between integral data types are allowed
• Using features of dynamic typing (changing floating point number to array/object/something else that is not number) is not allowed
• If's are allowed with one or none operators in condition (if (a) {...} else {...} or if (a==b) {...} else {...}). Note that if itself and its condition scores separately
• You may use try{}catch{}

Output

• Output should be passed by reference to input variables, if your language doesn't have such feature, you anyway should put result to input variables and print their values
• In variable a - number of real roots - 0, 1, 2, 3 (value greater than 2 stands for infinitely many roots)
• In variables b and c - 1st and 2nd roots (if none or infinitely many roots than any values, if one root then c - any value)
• Note that 2 equal roots is not the same as 1 root
• Precision should be at least 3 significant decimal digits

Score

• Score is [number of operators] + [number of if's]. Assignment operators have score 0.
• Answer with lowest score wins.
• In case of same score answer with greater precision wins. If precision is same then first answer wins.
• if(a != 0) and if(a == 0) has score 1 (instead of 2) because if(a) or if(a)else{} can be used, however not all languages support it.
• For each catch{} +2 to score
• Use of 1 temporary variable of numeric type is +20 to score.

Answer 1

C, 60

I happen to use the same square root algorithm (Newton) as Peter. However, I think we need to use two different starting values in order to get away with even 8 iterations. The maximum value of b²/4-c (the normalised discriminant) is for b = -c = 1000. The minimum value is for b = 0.001; c = 0 (in fact, there are smaller values but these would be beyond the significant digits of b and c). If we split the range of values at 20, we can use 0.1 as the starting guess for the lower range and 100 as the starting guess for the upper range and always converge in 8 steps. If I could be bothered to keep searching I think I could get this down to 7, and with a few more subdivisions probably even below 6 while still saving operations (but I can't be bothered right now).

void qe(&float a, &float b, &float c)
{
    float t; // 20
    if (a)   // 21
    {
        b = b / a;  // 22
        c = c / a;  // 23
        b = b / -2; // 24
        t = b * b;  // 25
        c = t - c;  // 26
        if (c<0)    // 28
        {
            a = 0;
        }
        else
        {
            if (c<20)  // 30
                a = .1
            else
                a = 100

            t = c / a; // 31
            a = a + t; // 32
            a = a / 2; // 33

            t = c / a; // 34
            a = a + t;
            a = a / 2;

            t = c / a; // 37
            a = a + t;
            a = a / 2;

            t = c / a; // 40
            a = a + t;
            a = a / 2;

            t = c / a; // 43
            a = a + t;
            a = a / 2;

            t = c / a; // 46
            t = a + t;
            t = t / 2;

            t = c / a; // 49
            t = a + t;
            t = t / 2;

            t = c / a; // 52
            t = a + t;
            t = t / 2;

            a = 2;
            c = b + t; // 55
            b = b - t; // 56
        }           
    }
    else if (b) // 57
    {
        a = 1;
        b = b / c; // 58
        b = -b;    // 59
    }
    else if (c) // 60
    {
        a = 0;
    }
    else
    {
        a = 3;
    }
}

Answer 2

Hybrid of iterative approaches; score: 62

I originally posted an answer which purely used Newton-Raphson; however, in the cases of large b and c this requires a lot of iterations. Consider only the most interesting quadratics: those which have two solutions. If we normalise them to the form x^2 + bx + c = 0 then we can plot the number of iterations required for various approaches as a function of b and c. For example, we can modify the continued fraction route to get the direct approach:

c = -c;
a = b * 0.5;
while (loopCount-- > 0) { // This loop would be unrolled in the final solution
    a = a + b;
    a = c / a;
}
// Roots are a and -c/a

That gives a loop count diagram of

2236...................6322
2224n.................n4222
22236.................63222
22224n...............n42222
222236...............632222
222224m.............m422222
2222236.............6322222
2222224j...........j4222222
22222234P.........P43222222
22222223d.........d32222222
222222222t.......t222222222
2222222225Z.....Z5222222222
22222222219/////91222222222
/////////////./////////////
222222222*********222222222
222222222*********222222222
22222222***********22222222
22222223***********32222222
2222223*************3222222
2222224*************4222222
222223***************322222
222224***************422222
22223*****************32222
22224*****************42222
2223*******************3222
2224*******************4222
223*********************322

where c is negative in the first row and increases downwards, b increases from left to right, * means negative discriminant, / means only one root was within tolerance in 62 iterations, . means no root was within tolerance in 62 iterations, and numbers are in base-62 with digits 0-9a-zA-Z. This is pretty decent for a certain range, but fails completely for another range.

bmarks' original implementation is similar to (but not the same; this is a reimplementation from the continued fraction):

a = b;
while (loopCount-- > 0) {
    a = c / a;
    a = b - a;
}
// Roots are -a, -c/a (which is better than a-b)

with loop count diagram

1117...................7111
1112n.................n2111
11117.................71111
11112n...............n21111
111117...............711111
111112m.............m211111
1111117.............7111111
1111112j...........j2111111
11111115Q.........Q51111111
11111111d.........d11111111
111111113u.......u311111111
1111111116/...../6111111111
1111111111a/////a1111111111
1111111111111.1111111111111
111111111*********111111111
111111111*********111111111
11111111***********11111111
11111112***********21111111
1111111*************1111111
1111112*************2111111
111111***************111111
111113***************311111
11111*****************11111
11113*****************31111
1111*******************1111
1113*******************3111
111*********************111

My naïve N-R was quite bad, but the really interesting thing to note about it is where the worst areas are:

njfdcccccccccccccccccccdfjn
njgcbbbbbbbbbbbbbbbbbbbcgjn
njgc9999999999999999999cgjn
njgc9777777777777777779cgjn
njgc9666666666666666669cgjn
njgc9544444444444444459cgjn
njgc9521111111111111259cgjn
njgc9534444444444444359cgjn
njgc9535555555555555359cgjn
njgc9536777777777776359cgjn
njgc9536888888888886359cgjn
njgc9536899999999986359cgjn
njgc95368aaaaaaaaa86359cgjn
njgc95368aaaaaaaaa86359cgjn
njgc95368*********86359cgjn
njgc95369*********96359cgjn
njgc9536***********6359cgjn
njgc9536***********6359cgjn
njgc953*************359cgjn
njgc954*************459cgjn
njgc95***************59cgjn
njgc94***************49cgjn
njgc9*****************9cgjn
njgc8*****************8cgjn
njgc*******************cgjn
njgb*******************bgjn
njf*********************fjn

So the thing to do is to hybridise. For the areas where it's good, I use the continued fraction; for the rest, I use N-R, but I've done some optimisation: I've experimentally found two cutoffs and three initial values which converge in 5 iterations to a worst-case mixed error abs(value - correct_value) / (1 + abs(correct_value) of 0.000052. (Thanks to Martin Büttner for the idea of doing a case split). The final loop count (mixing the two cases into one diagram) is:

111555555555555555555555111
111233333333333333333332111
111133333333333333333331111
111125555555555555555521111
111115555555555555555511111
111112333333333333333211111
111111333333333333333111111
111111255555555555552111111
111111144444444444441111111
111111112111111111111111111
111111113222222222211111111
111111111333333333111111111
111111111444444444111111111
111111111111111111111111111
111111111*********111111111
111111111*********111111111
11111111***********11111111
11111112***********21111111
1111111*************1111111
1111112*************2111111
111111***************111111
111113***************311111
11111*****************11111
11113*****************31111
1111*******************1111
1113*******************3111
111*********************111

and the code is

solve_quad(&float a, &float b, &float c) {
    float t; // 20 points
    if (a != 0) { // 21
        if (c == 0) {
            // x(x+b) = 0
            b = -b; // 22
            a = 2;
        }
        else {
            // Normalise so that we're solving x^2 + bx + c = 0 and have `a` free to use
            // Hereon, mentions of b and c in comments refer to these normalised values.
            b = b / a; // 22
            c = c / a; // 23

            // Three cases:
            //   b^2 < 4c: negative determinant, no solutions
            //   b^2 < -4c: use Newton-Raphson
            //   Otherwise: use the continued fraction method that bmarks was the first to post

            a = b / 2; // 24
            a = a * a; // 25
            if (a < c) { // 26
                a = 0;
            }
            else {
                t = a + c; // 27
                if (t < 0) { // 28
                    c = a - c; // 29, c holds the determinant

                    // Carefully optimised case split:
                    if (c < 0.05) // 30
                        a = 0.025;
                    else if (c < 250) // 31
                        a = 2;
                    else a = 193;

                    t = c / a; a = a + t; a = a / 2; // 34
                    t = c / a; a = a + t; a = a / 2; // 37
                    t = c / a; a = a + t; a = a / 2; // 40
                    t = c / a; a = a + t; a = a / 2; // 43
                    t = c / a; a = a + t; a = a / 2; // 46

                    // At this point a ~= sqrt(determinant)
                    t = b / 2; // 47
                    a = a + t; // 48
                    t = t * t; // 49
                    c = t - c; // 50
                    // At this point a ~= b/2 + sqrt(determinant) and c holds its original value
                }
                else {
                    a = b;

                    a = c / a; a = b - a; // 52
                    a = c / a; a = b - a; // 54
                    a = c / a; a = b - a; // 56

                    // At this point a ~= b/2 + sqrt(determinant)
                }

                b = -a; // 57
                c = c / b; // 58
                a = 2;
            }
        }
    }
    else {
        // bx + c = 0
        if (b != 0) { // 59
            a = 1;
            c = -c; // 60
            b = c / b; // 61
        }
        else {
            // c = 0
            if (c != 0) a = 0; // 62
            else a = 3;
        }
    }
}

Answer 3

C# (46)

Used a different way of computing sqrt(x) (continued fractions http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Generalized_continued_fraction)

It fails for b=0 due to large rounding errors. I'm not sure how to fix that, but I hope I can help fellow contestants with a new way of doing the square roots.

Edit: Thanks to @Peter Taylor for the suggestion on how to handle b=0.

Edit 2: Thanks to @steveverrill for pointing out that there are only >2 solutions if c is also 0, thus requiring an extra if check.

       public static void doStuff(ref double a, ref double b, ref double c)
    {
        if (a == 0) //1
        {
            if (b == 0) //2
            {
                if (c == 0){ //3
                   a = 3;
                }
                else {
                   a = 0;
                }
            }
            else
            {
                b = -1 / b;//4
                b = c * b;//5
                a = 1;
            }
        }
        else
        {
            b = b / a;//6
            c = c / a;//7
            if (b == 0) //8
            {
                b = 1 + 1E-7;// still just one constant; since 0.001/1000>1E-7, b can never equal this value normally
                c = c + 0.25; //9
            }
            a = 2 * c;//10
            a = a / b;//11
            a = b - a;//12
            a = c / a; a = b - a;//14
            a = c / a; a = b - a;//16
            a = c / a; a = b - a;//18
            a = c / a; a = b - a;//20
            a = c / a; a = b - a;//22
            a = c / a; a = b - a;//24
            a = c / a; a = b - a;//26
            a = c / a; a = b - a;//28
            a = c / a; a = b - a;//30
            a = c / a; a = b - a;//32
            a = c / a; a = b - a;//34
            a = 2 / a;//35
            a = c * a;//36
            a = b - a;//37
            if (b == 1 + 1E-7) //38
            {
                b = 0;
                c = c - 0.25;//39
            }
            if (a == 0) //40
            {
                b = -1.0 / 2 * b;//41
                c = b;
                a = 2;
            }
            else
            {
                if (a < 0) //42
                {
                    a = 0;
                }
                else
                {
                    b = -1.0 / 2 * b;//43
                    a = a / 2;//44
                    c = b - a;//45
                    b = b + a;//46
                    a = 2;
                }
            }
        }
    }

Answer 4

C, 46 excluding handling a=0 cases, with "fast squareroot" (not "fast inverse squareroot")

I hope I've scored what I've done so far correctly. I don't know how much more I'll do on this because:

1. I'm not sure if casting's legal, or how to score it (the union used in another answer is more efficient, but I've not used them before.)

2. I'm a little unclear on the meaning of "2 equal roots are not the same as 1 root." I assume this means that a quadratic will always have either 0 o 2 real roots (which may or may not be equal) and a linear equation (a=0) will have 1 root. I would point out the only infinite case is a=b=c=0. If a=b=0 and c is nonzero, there are no roots.

Anyway, my main reason for posting is that I wanted to show that it's possible to define a "fast squareroot" (which is more suitable for the needs of this question than simply copying the famous "fast inverse squareroot" as per some other answers.) Instead of subtracting i>>1 from a constant, we add i>>1 to a constant.

The principle of the algorithm is as follows: Assume a float is stored in an imaginary little computer as 2^AAA*1.BBBB (the 2 and 1 are not stored. AAA is known as the exponent and BBBB as the fraction.)

Performing a rightshift on the number 2^100*1.0000=16 gives 2^010*1.0000=4 which is exactly the correct square root. Repeating the process gives 2^001*1.0000=2. So far we have seen perfectly accurate results.

Now we have an odd power of 2. When we try to calculate the square root of the number 2 by rightshifting, the odd digit of the exponent underflows into the fraction and we get 2^000*1.100 which in decimal is 1.5, an error of approximately +6% on the correct value of sqrt(2)=1.414.

Moving to practical implemenations the 32-bit IEEE float (and the 64-bit, 128 bit etc.) has a bias on the exponent, so we have to add a constant to the raw shifted value. Careful choice of this constant enables the error to be limited to a range of approximately +/-3% which is a pretty good starting guess for The Newton-Raphson method. I found the constant by trial and error (there is room for some improvement, on my test cases it is more likely to give a positive error than a negative error.)

Three iterations of Newton-Raphson cleans it up very well (and may be more than necessary.)

void quad(float &a, float &b, float &c) {
  printf("a=%f b=%f c=%f",a,b,c);

  int i;                             //score 20

  c=c/a;                              // c/a
  b=b/a;
  b=b/-2;                             // -b/2a
  a=b*b;
  a=a-c;                              // (b/2a)^2-c/a), score 25

  if (a<0) a=0; else{                 // score 27

    i  = * ( int * ) &a;                       
    i  = i >>1;
    i  = i+0x1FB90000;        
    c  = * ( float * ) &i;            //sqrt( (-b/2a)^2-c/a ) score 31


    i  = * ( long * ) &a;             //don't lose the value (b/2a)^2-c/a
    a=a/c;
    c=c+a;
    a=c/2;
    c  = * ( float * ) &i;            //second iteration
    c=c/a;
    a=a+c;
    a=a/2;
    c  = * ( float * ) &i;            //third iteration  
    c=c/a;
    a=a+c;
    a=a/2;                            //score 43 after 3 iteratations

    c=b-a;                            //-b/2a in each
    b=b+a;             
    a=2;                              //score 46
  }
}                     

Answer 5

C++ (29 45 44)

Note that each return can be removed, simply by using else branch.. it's just more readable. The same with if(!sth) can be simply changed using else and empty if branch....

edit: Now, with casts allowed, this answer should be correct.

void quadr(float &a, float &b, float &c) {  
    b = - b;                //1
    if(!(b && c)) {         //3
        if(!a) {            //4
            a = 3;          //
            return;
        }
        a = 1;              //
        return;
    }
    if(!a) {                //5
        a = 1;              //
        b = c / b;          //6
        return;
    }

    b = b / a;              //7
    c = c / a;              //8
    c = c * 4;              //9
    a = b * b;              //10
    a = a - c;              //11
    c = a / 2;              //12

    int t = *(int*)&a;      //32
    t = t >> 1;             //33
    t = 0x5f3759df - t;     //34
    a = *(float*)&t;        //

    c = c * a;              //35
    c = c * a;              //36
    c = 1.5f - c;           //37
    a = a * c;              //38

    a = 1 / a;              //39
    if(!a) {                //40
        a = 0;              //
        return;     
    }
    c = b - a;              //41
    b = b + a;              //42
    c = c / 2;              //43
    b = b / 2;              //44
    a = 2;                  //
}

I might have missed some case, feel free to prove me wrong :)

Answer 6

C - 50 48

I use the 32-bit float shortcut for getting a good initial seed for the square root. That allows me to get away with just a few iterations of Newton's method. I haven't found a set of coefficients yet that yield unacceptable results (outside the specified precision).

Also, I am not sure how the spec treats pointer operations. For all I know this is disallowed, but hey, I had fun!

void primitiveQuad(float *a, float *b, float *c) {
    if(0.0 == *a) { // 1
        if(0.0 != *b) { // 2
            *a = 1.0;
            *b = -1.0/(*b); // 3
            *b = (*b)*(*c); // 4
        }
        else if(*c) // 5
            *a = 0.0; 
        else
            *a = 3.0;
    }
    else {
        *b = (*b)/(*a); // 6
        *c = (*c)/(*a); // 7
        *c = 4.0*(*c); // 8
        *a = (*b)*(*b); // 9
        *a = (*a) - (*c); // 10
        *c = *a; // store discriminant here for later
        if(0.0 > *a) // 11
            *a = 0.0;
        else {
            int x = *((int*)(a)); // 31
            x = x - (0x800000); // 32
            x = x >> 1; // 33
            x = x + (0x20000000); // 34
            *a = *((float*)(&x));
            {
                // add precision here through iterations

                // iteration
                x = *((int*)(a));
                *a = (*c)/(*a); // 35
                *a = *a + *((float*)(&x)); // 36
                *a = 0.5*(*a); // 37

                // iteration
                x = *((int*)(a));
                *a = (*c)/(*a); // 38
                *a = *a + *((float*)(&x)); // 39
                *a = 0.5*(*a); // 40

                // iteration
                x = *((int*)(a));
                *a = (*c)/(*a); // 41
                *a = *a + *((float*)(&x)); // 42
                *a = 0.5*(*a); // 43

            }
            *b = 0.0 - *b; // 44
            *c = *b;
            *b = *b + *a; // 45
            *c = *c - *a; // 46
            *b = 0.5*(*b); // 47
            *c = 0.5*(*c); // 48
            *a = 2.0;
        }
    }
}