/* (x-1.23456789e10)*(x-1.23456789e-10) = 0
 * 上記式を展開して得られる
 * x^2 -1.23456789e10x +1.52415788 = 0
 * を例にとる。
 */
#include <stdio.h>
#include <math.h>
int main(void) {
    double  a,b,c;
    double  x1,x2;
    a = +1.0;
    b = -1.23456789e10;
    c = +1.52415788;
    
    /* 解の公式 */
    x1 = (-b+sqrt(b*b-4.0*a*c))/(2.0*a);
    x2 = (-b-sqrt(b*b-4.0*a*c))/(2.0*a);
    printf("解の公式\n");
    printf(" x1 : %e\n",x1);
    printf(" x2 : %e\n",x2);
    
    /* 改良版 */
    if (b>0) x1 = (-b-sqrt(b*b-4.0*a*c))/(2.0*a);
    else     x1 = (-b+sqrt(b*b-4.0*a*c))/(2.0*a);
    /* ここで (x-x1)(x-x2) = 0 を考える。
     * これを展開すると
     *  x^2 -(x1+x2)x +x1x2 = 0
     * だから、加減算の入らない最終項の x1x2 を使用する。
     * このケースでは x1x2 = c(1.52415788)だから・・・
     */
    x2 = c / x1;
    printf("改良版\n");
    printf(" x1 : %e\n",x1);
    printf(" x2 : %e\n",x2);
    
    return 0;
}
