2차 방정식의 근을 구하는 c code에 관하여

개인적으로 네이버 지식인에서 코딩관련 답변을 즐겨? 다는 편인데 최근에 2차방정식의 근을 구하는 코드를 물어보는 질문이 올라왔다.

해당 방정식은 수직으로 던진 물체의 시간에 따른 높이에 관한 방정식인데 아래와 같다.

h = vt - gt^2 * 1/2

g라는 중력가속도 상황에서 v 라는 속도로 던져진 물체가 t 라는 시간이 지난뒤의 높이를 구하는 공식이다.

이때 h=100m, v=50m/s, g=9.8m/s^2 일때 t 를 구하는 코드를 만들어보라는 것이 문제이다.

해당 방정식을 t 에 대하여 푸는 것을 갓울프럼알파에게 물어보면(졸업한지 하도 오래되어서 저런거 손으로 푸는것도 이젠 어렵다 ㅠㅠ;;) 아래와 같다.

울프럼알파 짱...2차방정식이므로 해는 2개이다.

저렇게 나온 식에 값을 넣으면 아주 간단하겠지만...ㅎㅎㅎ 그게 문제가 아니라 값을 넣어가면서 근사값을 찾도록 알고리즘을 구성하라는 것이 문제이다.

일단 정답을 구해보면 아래와 같다.

그래프까지 아주 친절하게 그려주는 갓울프럼알파.

t 의 값은 t=2.7324, t=7.4612 2개가 정답이다.

이 값을 코드로 근사값을 구해보면 알고리즘을 아래와 구성하면 될듯 싶어서 해보았다.

h 값을 우변으로 옮기면 아래와 같은데...

0 = -h + vt - gt^2 * 1/2

이 공식에 주어진 값들을 넣어서 계산하고 t 값이 어디있는지 모르므로 일단 -100부터 100까지 돌면서 값을 구해본다.(초기 증감값을 1) 이때 각 결과값의 부호를 검사하여 부호가 바뀌는 순간의 직전값이 1자리에서의 정답이 되겠다. 그러나 우리고 구하고 싶은 값은 적당한 소수점 자리수까지의 값이므로 위에서 구한 2가지의 1자리 값에서 이후로는 0.1씩 증감값을 바꿔가면서 역시 부호가 바뀌는 값을 찾는다.

이후로는 계속 구하고 싶은 소수점 자리수까지 증감값을 1/pow(10, 구하고싶은자리수) 만큼으로 변경해가면서 loop를 돌면서 값을 구하면 된다.

이런 방식으로 소수점 6째 자리까지 구한 답은 아래와 같다.

c의 실수연산의 한계로 소수점 7자리 이상은 구하기 힘들다..ㅠㅠ

소수점 6째 자리까지 밖에 구하지 못했는데... 구한값이 t 는 시간값이므로.. 뭐 0.000001초 정도의 정밀도면... 충분하다고 보여진다 ㅋㅋ(더 작은 수 혹은 큰수를 다루는 수학 라이브러리나 파이썬의 수학 라이브러리 등을 이용하면 훨씬 작은 값으로도 구할 수 있을듯...)

아래는 사용한 코드이다.

int main(int argc, char* argv[])
{
    long double v0 = 50, h = 100, g = 9.81, t, possible_answer[2], offset, prev_sign, current_sign, increment_value, divide_value, one=1;
    int index = 0, decimal_point=0, i;

    std::setprecision(16);
    //t 값을 -100부터 100까지 돌면서 방정식에서 h를 우변으로 보냈을때 좌변이 0이 되므로 계산값이 부호가 바뀌는 시점을 2개 찾는다.
    for (t = -100; t < 100; t++) {
        offset = (h - (v0 * t - g * t * t * 0.5));
        if (offset > 0) current_sign = 1;
        else current_sign = -1;
        //printf("t %Lf  offset %Lf\n", t, offset);
        if (t > 0) {
            if (prev_sign != current_sign) {
                possible_answer[index++] = t - 1;
            }
        }
        prev_sign = current_sign;
    }
    //위에서 찾은 2개의 값에서 소수점 n번째 자리까지 돌면서 부호가 바뀌는 시점의 값을 계속 찾아나간다. 그 값이 해임.
    for (decimal_point = 1; decimal_point <= 6; decimal_point++) {
        divide_value = 1;
        for (i = 0; i < decimal_point; i++)
            divide_value *= 10;
        increment_value = one / divide_value;
        for (index = 0; index < 2; index++) {
            t = possible_answer[index];
            for (i = 1; i < 9; i++) {
                t = t+ increment_value;
                offset = (h - (v0 * t - g * t * t * 0.5));
                //printf("t %Lf  offset %Lf divide_value %Lf increment_value %Lf\n", t, offset, divide_value, increment_value);
                if (offset > 0) current_sign = 1;
                else current_sign = -1;
                if (i > 1) {
                    if (prev_sign != current_sign) {
                        possible_answer[index] = t - increment_value;
                    }
                }
                prev_sign = current_sign;
            }
        }
    }

    printf("Possible Answer 1 : %Lf\nPossible Answer 2 : %f", possible_answer[0], possible_answer[1]);

    return 0;
}