Solve quadratic equation using only 3 variables and arithmetic 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 + - * / % &amp; | xor ~ &amp;&amp; || ! &lt; &gt; &gt;= &lt;= == != &lt;&lt; &gt;&gt; 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.
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-- &gt; 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-- &gt; 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(&amp;float a, &amp;float b, &amp;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 &lt; 4c: negative determinant, no solutions
            //   b^2 &lt; -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 &lt; c) { // 26
                a = 0;
            }
            else {
                t = a + c; // 27
                if (t &lt; 0) { // 28
                    c = a - c; // 29, c holds the determinant

                    // Carefully optimised case split:
                    if (c &lt; 0.05) // 30
                        a = 0.025;
                    else if (c &lt; 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;
        }
    }
}