테일러 급수로 구현한 ln() 함수

2020/03/19 - 온도 센서의 출력 저항값으로 원래 온도값 계산 삽질기

 

이전 포스팅에서 로그와 간단한 식을 이용해서 센서 데이터를 처리하는 방식들을 얘기했다.

그런데, 새로운 문제가 제기되었다.

 

구현하는 환경에서 log() 함수가 지원되지 않으면 어떻게 하지?

 

쿼드코어 CPU가 발길에 채이는 세상이지만, 놀랍게도 이런 문제는 상존한다.

예컨데, Cortex 시리즈의 경우 M0는 물론 M3까지도 HW 부동소수점 기능이 제공되지 않는다.

 

 

1. 로그(ln(x))의 구현 방안

 

로그의 계산은 당연하게도(!) 테일러 급수에서 시작된다.

 

자연로그는 아래 두 가지 형태로 기술된다.

 

 

 

뭔가 좀 복잡해보이니 간단하게(?) 전개하면 아래와 같다.

 

 

 

그런데, 이 두 식은 조금만 생각해보면 홀수차 항만 남기도록 정리가 가능하다.

연산량을 가능한 줄이는 것이 좋다는 점을 고려하면 이 정리는 상당히 유용한 과정이다.

 

 

여기서, y를 아래와 같이 두면

 

 

x는 이렇게 정리가 가능하다.

 

 

주의해야 할 점은 -1 < x < 1 이므로 0 < y < 2 라는 점.

 

 

2. 실제 구현

 

여기까지 왔으면 이 내용을 그대로 코드에 담으면 되긴 하지만, 또 하나의 문제가 남는다.

과연 몇 자리까지 저 계산을 할 것인가 하는 점.

 

const static double E = 2.71828182845904523536;

bool LN_e_v1(double y, double limit, double& result, int& count) {
    result = 0;
    count = 0;

    if (y <= 0) {
        return false;
    }

    int N = 0;
    while (y >= 2) {
        ++N;
        y /= E;
    }

    double x = (y - 1) / (y + 1);
    double xx = x * x;

    double xxT1 = x;
    double xxT2;

    result = x;
    count = 0;
    while (true) {
        ++count;
        xxT1 *= xx;
        xxT2 = xxT1 / (count * 2 + 1);

        result += xxT2;

        if (fabs(xxT2) < limit) {
            break;
        }
    }

    result = result + result + N;
    return true;
}

 

이렇게 해서 항의 크기가 10^{-10}가 될 때까지 전개해보니 대략 10~12개 까지 전개하면 상당히 정확한 결과를 얻을 수 있었다.

하지만, 여긴 또 하나의 문제가 있는데, 입력값이 0에 가까워질수록 연산량이 심각하게 늘어난다는 점.

 

이런 점들을 고려해서 작성한 결과는 아래와 같다.

총 20회의 루프를 돌리는데 실제로는 10회만 돌려도 꽤 정확한 값이 나온다.

 

const static double E = 2.71828182845904523536;

double LN_e_v2(double y) {
    double result = 0;

    if (y <= 0) {
        return INFINITY;
    }

    int N = 0;
    while (y >= 2) { ++N; y /= E; }
    while (y <= 0.4) { --N; y *= E; }

    double x = (y - 1) / (y + 1);
    double xx = x * x;

    double xxT1 = x;
    double xxT2;

    result = x;
    int i2 = 3;
    for (int i = 0; i < 20; ++i, i2 += 2) {
        xxT1 *= xx;
        xxT2 = xxT1 / i2;

        result += xxT2;
    }

    result = result + result + N;
    
    return result;
}

 

같이 만들어본 sin(), pi()의 결과