pow 구현

double y = pow(double x, double z);

y = x ^ z는
y = (x ^ (1 / b)) ^ a, z = a / b입니다.

계산을 더 조금만 시키기 위해 a와 b를 서로소로 만들어 줍니다.

long long gcp(long long a, long long b) {
    if (b) return gcp(b, a % b);
    return a;
}

double lld_pow(double, long long);
double inv_lld_pow(double, long long);

double pow(double x, double z) {
	if (x < 0 || z < 0) return -1.0; // 이 경우는 고려 안 하겠습니다.
    long long b = 1000000000, a = z * b, tmp = gcp(a, b);
    a /= tmp, b /= tmp;
    return lld_pow(inv_lld_pow(x, b), a);
}

정수 거듭제곱 함수 double lld_pow(double, long long)를 분할정복으로 구현합니다.

double lld_pow(double x, long long n) {
    double result = 1.0;
    while (n) {
        if (n & 1) result *= x;
        x *= x;
        n >>= 1;
    }
    return result;
}

x ^ (1 / b)를 구하는 함수 double inv_lld_pow(double, long long)를 구현합니다.

y = x0 ^ (1 / b)인 y를 구해야 합니다. 보통의 방정식을 풀 때, x를 미지수로 놓기 때문에 x와 y를 바꿔 생각하겠습니다.
x = y0 ^ (1 / b) -> x ^ b = y0이고, x의 값은
y = x ^ b - y0가 x축과 만날 때 x의 값입니다.
이를 만족하는 x의 근사값은 그림과 같이 Newton-Raphson Method를 사용하여 구할 수 있습니다.

임의의 x0 (x0 > 0, 식에서는 1.0으로 지정)값부터 시작하여, 접선과 x축이 만나는 점을 계속해서 구하다보면, 구하고자 하는 근사값을 구할 수 있습니다.

double inv_lld_pow(double y0, long long n) {
    const double ERR = 0.000001;
    double x = 1.0, y = lld_pow(x, n) - y0;
    while (y > ERR || y < -ERR) {
        x -= y / (n * lld_pow(x, n - 1));
        y = lld_pow(x, n) - y0;
    }
    return x;
}

double pow(double, double)

위의 내용들을 종합하여 완성된 모습입니다.

long long gcp(long long a, long long b) {
    if (b) return gcp(b, a % b);
    return a;
}
double lld_pow(double x, long long n) {
    double result = 1.0;
    while (n) {
        if (n & 1) result *= x;
        x *= x;
        n >>= 1;
    }
    return result;
}
double inv_lld_pow(double y0, long long n) {
    const double ERR = 0.000001;
    double x = 1.0, y = lld_pow(x, n) - y0;
    while (y > ERR || y < -ERR) {
        x -= y / (n * lld_pow(x, n - 1));
        y = lld_pow(x, n) - y0;
    }
    return x;
}

double pow(double x, double z) {
    if (x < 0 || z < 0) return -1.0;
    long long b = 1000000000, a = z * b, tmp = gcp(a, b);
    a /= tmp, b /= tmp;
    return lld_pow(inv_lld_pow(x, b), a);
}