Quadratic equation in Ada

问题:

I just came around and decided to try some Ada. The downside is that the syntax and function strays quite away from C++. So I had to like cram various stuff to get this thing to work.

My question is if there are some better way to do this calculation that what I have done here

   IF(B < 0.0) THEN
      B := ABS(B);
      X1 := (B / 2.0) + Sqrt( (B / 2.0) ** 2.0 + ABS(C));
      X2 := (B / 2.0) - Sqrt( (B / 2.0) ** 2.0 + ABS(C));
   ELSE
      X1 := -(B / 2.0) + Sqrt( (B / 2.0) ** 2.0 - C);
      X2 := -(B / 2.0) - Sqrt( (B / 2.0) ** 2.0 - C);
   END IF;

I had some problem with negative numbers, that's why I did a IF statement and used ABS() to turn those into positive. But the weird thing is that it works perfectly for the other case, which is strange...

回答1:

Solving quadratic equations is not as simple as most people think.

The standard formula for solving a x^2 + b x + c = 0 is

delta = b^2 - 4 a c
x1 = (-b + sqrt(delta)) / (2 a)   (*)
x2 = (-b - sqrt(delta)) / (2 a)

but when 4 a c << b^2, computing x1 involves subtracting close numbers, and makes you lose accuracy, so you use the following instead

delta as above
x1 = 2 c / (-b - sqrt(delta))     (**)
x2 = 2 c / (-b + sqrt(delta))

which yields a better x1, but whose x2 has the same problem as x1 had above.

The correct way to compute the roots is therefore

q = -0.5 (b + sign(b) sqrt(delta))

and use x1 = q / a and x2 = c / q, which I find very efficient. If you want to handle the case when delta is negative, or complex coefficients, then you must use complex arithmetic (which is quite tricky to get right too).

Edit: With Ada code:

DELTA := B * B - 4.0 * A * C;

IF(B > 0.0) THEN
    Q := -0.5 * (B + SQRT(DELTA));
ELSE
    Q := -0.5 * (B - SQRT(DELTA));
END IF;

X1 := Q / A;
X2 := C / Q;

回答2:

Given ax² + bx + c = 0 the quadradic formula gives solutions for x = (-b +/- sqrt(b²-4ac) ) / 2a. The discriminant d = b²-4ac will be positive for real valued roots, negative for roots with a non-zero imaginary component (i.e. a non-real complex number) and will be 0 when the root is a double root.

So, the Ada code for this would be:

D := B ** 2.0 - 4.0 * A * C;
IF D >= 0.0 THEN
  X1 := (-B + Sqrt(D)) / (2.0 * A);
  X2 := (-B - Sqrt(D)) / (2.0 * A);
ELSE
  -- Deal with the fact that the result is a non-real complex number.
END IF;

Note: I'm a bit rusty in Ada, but this should be close to the proper syntax.

回答3:

The quadratic formula is x = ( -b +/- sqrt ( b ** 2 - 4*a*c ) ) / ( 2 * a )

I'm guessing a is 1.

so x = -( b/2 ) +/- sqrt ( ( ( b ** 2 ) / 4 ) - c )

Calculate d = ( b ** 2 ) * 0.25 - c then check its sign.

If the sign of d is negative you have complex roots; handle them as you wish.

Replacing - c with + abs ( c ) if b happens to be negative will give you rubbish.

Usually multiplying by 0.5 or 0.25 is better than dividing by 2.0 or 4.0.

回答4:

Though I don't know Ada, I see following things that can be optimized:

1. In the first branch of the IF instructiuon you already know that B is negative. So you can say B := -B instead of B := ABS(B). Or better: just use -B where you used B in the first branch.
2. You are using the subexpression B/2.0 four times. It could be more efficient (and also more clear) to assign B/2.0 to an auxilary variable B_2 (or assign it to B again if you don't want to spend another variable) and use it instead.
   Also the sqrt is calculated twice. Assigning it to an auxilariy variable saves runtime (and makes it explicite for the reader that exactly the same subexpresion is used twice).
3. Than it would probably be faster to use B_2*B_2 instead of **2.0; even better would be to use a dedicated square function, if there is one in Ada.

回答5:

To me, the question is more related to numeric algorithm than to Ada language. As always with numeric computing, one must often (if not -always-) refer to reference/academic papers.

These kinda questions always reminds me of this : https://en.wikipedia.org/wiki/Fast_inverse_square_root

You can only find the following tricks if you "do the maths" or find some paper that adresses your issue.

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

PS: as the wikiepdia article points out, this implementation is probably obsolete for most platforms now