백준 27875번: 파이파이

https://www.acmicpc.net/problem/27875

27875번: 파이파이

---

인터넷 검색으로 어떻게든 레퍼런스를 찾기만 하면 매우 쉬워지는 문제였습니다.

 

일반적으로 어떤 실수 $x$의 소수점 아래 $N$번째 수를 구하기 위해서는 $10^{N} \times x$을 관찰합니다.

 

이 문제에서 $\pi ^{2}$의 소수점 아래 $N$번째 수는, 16진법이므로, $16^{N} \times \pi ^{2}$를 얻어야 할 것입니다.

 

여기서 $\pi ^{2}$의 16진법 표현을 얻는 방법이 문제의 핵심이었습니다. $\pi$의 16진법 표현을 안다면 어쩌면 가능할 수도 있습니다만

 

소스 코드 길이 제한에 걸릴 가능성이 큽니다. 따라서, $\pi ^{2}$의 16진법 표현과 곧바로 연결되는 어떤 수식의 존재를 찾아야 합니다.

---

1997년 David H. Bailey, Peter B. Borwein, Simon Plouffe 3명의 저자가 출판한 On The Rapid Computation of Various Polylogarithmic Constants에서 한 방법을 제시합니다.

$$ \pi ^{2} = \sum_{k=0}^{\infty} \left [ \frac{1}{16^{k}} \left \{ \frac{16}{(8k+1)^{2}} + \frac{-16}{(8k+2)^{2}} + \frac{-8}{(8k+3)^{2}} + \frac{-16}{(8k+4)^{2}} + \frac{-4}{(8k+5)^{2}} + \frac{-4}{(8k+6)^{2}} + \frac{2}{(8k+7)^{2}} \right \} \right ] $$

 

이러한 방법론을 사용하는 공식을 저자들의 이름으로부터 BBP-type formulas라고 합니다. 일반적으로 PSLQ Integer Relation Algorithm같은 알고리즘을 통해 이런 공식을 찾는다고 합니다.

 

저 공식을 verify하려고 했는데, $\pi$ 버전과 다르게 $\pi^{2}$ 버전은 난이도가 좀 어려웠습니다. 나중에 추가해보겠습니다. 저자들에 따르면 아이디어는 똑같다고 하는데...

 

일단 27875번을 풀기 위해 저 식의 양변에 $16^{N}$을 곱해봅니다.

$$ 16^{N} \times \pi ^{2} = \sum_{k=0}^{\infty} \left [ 16^{N-k} \left \{ \frac{16}{(8k+1)^{2}} + \frac{-16}{(8k+2)^{2}} + \frac{-8}{(8k+3)^{2}} + \frac{-16}{(8k+4)^{2}} + \frac{-4}{(8k+5)^{2}} + \frac{-4}{(8k+6)^{2}} + \frac{2}{(8k+7)^{2}} \right \} \right ] $$

$$ = 16 \sum_{k=0}^{\infty} \frac{16^{N-k}}{(8k+1)^{2}} - 16 \sum_{k=0}^{\infty} \frac{16^{N-k}}{(8k+2)^{2}} - 8 \sum_{k=0}^{\infty} \frac{16^{N-k}}{(8k+3)^{2}} - 16 \sum_{k=0}^{\infty} \frac{16^{N-k}}{(8k+4)^{2}} - 4 \sum_{k=0}^{\infty} \frac{16^{N-k}}{(8k+5)^{2}} - 4 \sum_{k=0}^{\infty} \frac{16^{N-k}}{(8k+6)^{2}} + 2 \sum_{k=0}^{\infty} \frac{16^{N-k}}{(8k+7)^{2}} $$

여기서 많이 반복되는 부분인

$$ \sum_{i=0}^{\infty} \frac{16^{N-i}}{(8i+j)^{2}} = S(i,j) $$

을 정밀하게 구하면 문제가 해결되겠습니다.

---

그런데 저 값을 정말로 모두 구하는 게 아니라, 소수점 아래 부분만 구해도 됩니다. $x$의 소수점 아래 부분, 즉 fractional part를 $\left \{ x \right \}$라고 하면

$$ \left \{ 16^{0} \pi ^{2} \right \} = 0.DE \cdots $$

이므로 $N=1$부터의 답이 나오게 됩니다. 다시 말해, $ \left \{ 16^{k} \pi ^{2} \right \}$ 을 정밀하게 구하면 27875번의 인풋 $N=k+1$일 때의 값을 구한 것입니다.

 

또한 fractional part는 비유하면 실수에 modulo 1.0을 취한 값입니다. 따라서 $S(i,j)$가 다음과 같이 reduce됩니다.

$$ \left \{ S(i,j) \right \} = \left \{ \sum_{i=0}^{\infty} \frac{16^{N-i} \;\; \textrm{mod} \;\; (8i+j)^{2}}{(8i+j)^{2}} \right \}$$

여기에서 $i \leq N$이면 익숙한 binary exponentiation을 통해서 빠르게 구할 수 있습니다.

 

그리고 $i > N$이면, C/C++ 환경 기준 machine epsilon 보다 구하는 값의 변화량이 작아지는 것은 매우 빠릅니다. 거의 $O(1)$ 급이라고 저는 생각합니다.

 

다시 말해,

$$ \left \{ S(i,j) \right \} = \left \{ \left \{ \sum_{i=0}^{N} \frac{16^{N-i} \;\; \textrm{mod} \;\; (8i+j)^{2}}{(8i+j)^{2}} \right \} + \left \{ \sum_{i=N+1}^{\infty} \frac{16^{N-i}}{(8i+j)^{2}} \right \} \right \} $$

을 구하는 시간 복잡도가 $O(N \log N)$라는 것입니다. (Stirling's Approximation) (cf. Bianry Exponentiation을 구현할 때 modulo가 변수로 들어가게 되는데 상수가 조금 커집니다.)

 

한편 machine epsilon을 언급했는데, 이 문제에서 분수 구조체를 사용하거나 정수 연산만으로 원하는 답을 구하기는 약간 힘듭니다. 물론 가능은 하지만, $N$의 범위가 double만으로 풀리는 수준이므로 double을 쓰는게 정신적으로 편합니다.

 

그리고 double을 쓸 경우, mantissa가 53비트이므로 해당 53비트를 16진수 정수로 바꿔 출력하면 됩니다.

---

이로써 문제를 풀기 위한 준비가 모두 끝났습니다. 구현할 때, C/C++ 기준 실수에서 fractional part를 구하기 위해서는 fmod를 쓰게 되는데, 이는 음숫값을 리턴할 수 있음에 주의해야 합니다. 추가로, Binary Exponentiation에서 modulo의 범위에 주의합니다.

 

제 코드는 다음과 같습니다.

/*
 * Author : thdtjdals3
 * Date : 2023-08-28
 * Description : David H. Bailey, 2023.
 *               "A Compendium of BBP-Type Formulas for Mathematical Constants"
 *               available at: https://www.davidhbailey.com/dhbpapers/bbp-formulas.pdf
 *
 *               Wikipedia
 *               https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
 *
 *               Implementation Hint
 *               https://literateprograms.org/pi_with_the_bbp_formula__python_.html
*/

#include <stdio.h>
#include <math.h>

long long power_modular (long long X, long long Y, long long M);

double S2(int j, int n);
double PI_Squared_BBP (int n);

int main() {
   
    int N;
    scanf("%d", &N);
   
    long long Answer = PI_Squared_BBP(N) * ((long long)1 << 56);
    printf("%llX", Answer >> 52);
    return 0;
}

long long power_modular (long long X, long long Y, long long M) {
    long long result = 1;
    while (Y) {
        if (Y % 2) {
            result = (__int128)result * X % M;
        }
        X = (__int128)X * X % M;
        Y /= 2;
    }
    return result;
}

double S2(int j, int n) {
    
    /*
     * summation of 16^(n-k) / 8k+j for k = 0 to infty, only fractional part
     *
     * first, seperate the summation at k = n
     *
     * For k = 0 to n, use exponentiation_by_squaring and time complexity is O(nlgn)
     * (Stirling's approximation)
     *
     * For k > n, use floating point number, exit the loop when any changes wont occur
    */
    double result_left = 0.0;
    for (int k = 0; k <= n; k++) {
        long long R = 8 * k + j;
        result_left = fmod(result_left + power_modular(16, n-k, R * R) / (double)(R * R), 1.0);
    }
   
    double result_right = 0.0;
    int k = n + 1;
    while(1) {
        long long R = 8 * k + j;
        double temp = result_right + pow(16, n-k) / (R * R);
        if (temp == result_right) {
            break;
        }
       
        result_right = temp;
        k++;
    }
   
    return result_left + result_right;
}

double PI_Squared_BBP (int n) {
    /*
     * "A Compendium of BBP-Type Formulas for Mathematical Constants"
     * formula 29 : pi^2 = P(2, 16, 8, <16, -16, -8, -16, -4, -4, 2, 0>)
     *
     * Calculate the fractional part of (16^n * pi^2) to find nth digit of pi^2
     *
     * Decrement of n is a convention.
     *
     * +48.0 is a bias so that the result of fmod can always be positive
    */
    n--;
    double result = fmod(16*S2(1,n) - 16*S2(2,n) - 8*S2(3,n) - 16*S2(4,n) - 4*S2(5,n) - 4*S2(6,n) + 2*S2(7,n) + 48.0, 1.0);
    return result;
}​

(C11, 1560ms, 1304KB, 제출번호 65748132)

---

물론 전체 문제의 시간 복잡도가 결국 $S(i,j)$의 복잡도와 같게 되지만, 상수 차이가 나름대로 크게 작용할 수 있습니다.

 

무슨 말이냐면, $S(i,j)$의 호출 횟수를 줄이거나, 더 빠르게 원하는 자릿수까지 접근하는 것입니다. 다음 공식을 보겠습니다.

$$ \pi ^{2} = \frac{9}{8} \times \sum_{k=0}^{\infty} \left [ \frac{1}{64^{k}} \left \{ \frac{16}{(6k+1)^{2}} + \frac{-24}{(6k+2)^{2}} + \frac{-8}{(6k+3)^{2}} + \frac{-6}{(6k+4)^{2}} + \frac{1}{(6k+5)^{2}} \right \} \right ] $$

이 역시 마찬가지로 On The Rapid Computation of Various Polylogarithmic Constants에서 나온 공식입니다.

 

이를 통해 $\pi^{2}$의 64진법 표현을 구하면, 16진법 표현을 훨씬 적은 자릿수로도 구하게 됩니다. 또한, $S(i,j)$의 호출 횟수 자체도 줄어들 것을 알 수 있습니다.

 

다만, 64진법 표현을 통해 16진법 표현을 구한다는 부분을 구현하는 게 약간 자명하지는 않습니다. 그 원인은 $64^{s} = 16^{n}$이더라도, $s$를 함수의 인풋 값으로 넣을 때는 정수로 넣고 싶어져서 생기는 문제입니다.

 

그 부분만 보정해주면 아까보다 훨씬 빠른 속도로 답을 얻게 됩니다. 제 코드는 다음과 같습니다.

/*
 * Author : thdtjdals3
 * Date : 2023-08-28
 * Description : David H. Bailey, 2023.
 *               "A Compendium of BBP-Type Formulas for Mathematical Constants"
 *               available at: https://www.davidhbailey.com/dhbpapers/bbp-formulas.pdf
 *
 *               Wikipedia
 *               https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
 *
 *               Implementation Hint
 *               https://literateprograms.org/pi_with_the_bbp_formula__python_.html
*/

#include <stdio.h>
#include <math.h>

long long power_modular (__int128 X, __int128 Y, __int128 M);

double S2(int j, int n);
double PI_Squared_BBP (int n);

int main() {
   
    int N;
    scanf("%d", &N);
    
    if (N < 10) {
        printf("%c", "DE9E64DF22"[N-1]);
        return 0;
    }
    
    long long Answer;
    Answer = PI_Squared_BBP(N-1) * ((long long)1 << 56);
    printf("%llX\n", Answer >> 52);
    return 0;
}

long long power_modular (__int128 X, __int128 Y, __int128 M) {
    __int128 result = 1;
    while (Y) {
        if (Y % 2) {
            result = result * X % M;
        }
        X = X * X % M;
        Y /= 2;
    }
    return result;
}

double S2(int j, int n) {
    
    /*
     * summation of 64^(n-k) / (6k+j)^2 for k = 0 to infty, only fractional part
     *
     * first, seperate the summation at k = n
     *
     * For k = 0 to n, use exponentiation_by_squaring and time complexity is O(nlgn)
     * (Stirling's approximation)
     *
     * For k > n, use floating point number, exit the loop when any changes wont occur
    */
    double result_left = 0.0;
    for (int k = 0; k <= n; k++) {
        long long R = 6 * k + j;
        result_left = fmod(result_left + power_modular(64, n-k, R * R) / (double)(R * R), 1.0);
    }
   
    double result_right = 0.0;
    int k = n + 1;
    while(1) {
        long long R = 6 * k + j;
        double temp = result_right + pow(64, n-k) / (R * R);
        if (temp == result_right) {
            break;
        }
       
        result_right = temp;
        k++;
    }
   
    return result_left + result_right;
}

double PI_Squared_BBP (int n) {
    /*
     * "A Compendium of BBP-Type Formulas for Mathematical Constants"
     * formula 30 : pi^2 = (9/8) * P(2, 64, 6, <16, -24, -8, -6, 1, 0>)
     *
     * Calculate the fractional part of (16^n * pi^2) to find nth digit of pi^2
     * but using "64^s * pi^2" formula, let s = 2(n-1)/3
     *
     * if (n-1) % 3 == 1 :   {16 * 64^s * pi^2} => {16^(n-1) * pi^2}
     * elif (n-1) % 3 == 2 : { 4 * 64^s * pi^2} => {16^(n-1) * pi^2}
     * else :                {     64^s * pi^2} => {16^(n-1) * pi^2}
     *
     * Decrement of n is a convention.
     *
     * +38.0 is a bias so that the result of fmod can always be positive
    */
    n--;
    
    int s = n * 2 / 3;
    double result = fmod(16*S2(1,s) - 24*S2(2,s) - 8*S2(3,s) - 6*S2(4,s) + S2(5,s) + 38.0, 1.0);
    
    if (n % 3 == 1) {
        result = fmod(result * 288.0, 1.0);
    }
    else if (n % 3 == 2) {
        result = fmod(result * 72.0, 1.0);
    }
    else {
        result = fmod(result * 18.0, 1.0);
    }
    
    return result;
}

(C11, 716ms, 1304KB, 제출번호 65753824)

---

여기서 unsigned 타입을 써서 약간 더 커팅할 수 있습니다. (680ms, 제출번호 65754346, 단 math.h의 pow를 직접 구현하는 건 큰 도움이 되진 않았습니다.)

 

만약 더 정밀한 타입을 원하는 경우, long double을 써볼 수 있습니다. 그러나 문제에서 필수적으로 요구하는 사항도 아니라서 제출해보지는 않았습니다.

 

일단 주의할 점은, long double을 쓰면 math.h 특성상 powl, fmodl 같이 함수를 바꿔줘야 합니다. 또한 mantissa가 64비트까지 늘어나서, 이 부분도 구현 난이도를 증가시킵니다.

 

__float128 을 쓰는 것은 이보다 더 구현 난이도가 올라갈 것입니다. 저는 여기서부터는 시도도 해보지 못했습니다.

---

한편, 더 큰 실수 자료형을 쓰는 게 아니라 정수 연산을 써도 정밀도를 크게 올릴 수 있습니다. 이 사이트를 참고할 수 있을 것입니다.

---

언젠가 $\pi ^{2}$에 대한 BBP-type formula을 verify하고 정수 연산만을 활용한 버전의 코드를 구현해보고 싶습니다.