파이썬에서 팩토리얼을 구현한 방법 (파이썬 Factorial : binary splitting)

글을 읽기 전에 아래 글을 먼저 읽는 것이 도움됩니다

(팩토리얼에서 recursion 대신 divide and conquer를 쓰는 이유)

 

https://jhong92-pro.github.io/%EC%95%8C%EA%B3%A0%EB%A6%AC%EC%A6%98/%ED%8C%A9%ED%86%A0%EB%A6%AC%EC%96%BC-%EC%8B%9C%EA%B0%84%EB%B3%B5%EC%9E%A1%EB%8F%84-recursion/

팩토리얼 시간복잡도 - recursion

https://jhong92-pro.github.io/%EC%95%8C%EA%B3%A0%EB%A6%AC%EC%A6%98/%ED%8C%A9%ED%86%A0%EB%A6%AC%EC%96%BC-%EC%8B%9C%EA%B0%84%EB%B3%B5%EC%9E%A1%EB%8F%84-divide-and-conquer/

팩토리얼 시간복잡도 - divide and conquer

 

 

cpython에 대한 자세한 내용은 건너뛰고,

cpython에 팩토리얼이 라이브러리로 구현되어 있습니다.

 

현재까지 알려진 팩토리얼 구현 알고리즘 중 가장 빠른 알고리즘은 아니지만,

binary splitting이라는 방법으로 팩토리얼이 구현되어 있습니다

 

binary splitting은 두 가지 step으로 구성되어 있습니다

1. 2를 인수분해하여 비트 연산으로 처리

2. 인수분해 결과 2가 없는 것들은 divide and conquer로 곱셈 수행

예를 들어,

$4!=2^{3}*(1*1*3*1)$

이런 방식으로 처리합니다

 

관련 내용은

https://github.com/python/cpython/blob/main/Modules/mathmodule.c

여기서 확인할 수 있으며,

 

http://www.luschny.de/math/factorial/binarysplitfact.html

이 사이트를 참고했다고 설명되어 있습니다.

 

2를 인수분해 하기

1번부터 살펴봅시다

 

$factorial(n)에서$

2를 인수분해하여 $2^{k}$에서 n을 얻는데 $log(log(n))\approx1$ 만큼의 시간이 걸립니다

 

팩토리얼에서 k가 몇 번 등장하는지 살펴봅시다

18! 을 인수 분해한 것을 봅시다

 

2 인수분해

1~18까지의 행을 제외하고

첫 번째 줄을 먼저 보자면 2가 두 번 중 한 번 반복되고 있습니다

두 번째 줄은 2가 4번 중 한 번 반복되고 있습니다

세 번째 줄은 2가 8번 중 한 번 반복되고 있습니다

네 번째 줄은 2가 16번 중 한 번 반복되고 있습니다

 

위 규칙성을 조합해서 수식을 만들면

(2의 개수) = $\left [ \frac{n}{2}  \right ] + \left [ \frac{n}{4}  \right ] + \left [ \frac{n}{8}  \right ] + \left [ \frac{n}{16}  \right ] + ...$입니다

 

다행히 가우스 기호 $\left [ \frac{n}{2^{k}}  \right ]$  는 비트연산 $n>>k$와 같습니다

 

18을 이진수로 나타내면 10010입니다

비트 연산을 통해 2의 개수를 계산해봅시다

$\left [ \frac{18}{2}  \right ] = 1001$

$\left [ \frac{18}{4}  \right ] = 100$

$\left [ \frac{18}{8}  \right ] = 10$

$\left [ \frac{18}{16}  \right ] = 1$

 

다 더하면

 

   1 0 0 1

+    1 0 0

+       1 0

+          1

-----------

1 0 0 0 0

 

즉 16입니다

 

10010에서

첫 번째 자리를 추적해보면 첫 번째 자릿수 1에 의해 1111이 더해졌습니다

두 번째 자리는 0이지만 만약 1이었다면 두 번째 자릿수에 의해 111이 더해졌을 겁니다

따라서,

 

(2의개수) = 1*1111 + 0*111 + 0*11 + 1*1 + 0*0 입니다

또한, 1111 = ${2^{4}-1}$입니다

 

2의 개수를 n! 에 대해 일반화하여 표현해보겠습니다

n을 이진수로 표현하면

${n = A_0*2^{k}+A_1*2^{k-1}+A_2*2^{k-2}+...+A_k*2^{0}}$

입니다

 

n!에 대한 2의 개수는

${A_0*(2^{k}-1)+A_1*(2^{k-1}-1)+A_2*(2^{k-2}-1)+...+A_{k-1}*(2^{0}-1) \\
= (A_0*2^{k}+A_1*2^{k-1}+A_2*2^{k-2}+...+A_k*2^{0}) - (A_0+A_1+A_2+...+A_k)}$

이는 n - (n을 이진수로 표현했을 때 1의 개수) 와 같습니다

 

그러면 (n을 이진수로 표현했을 때 1의 개수)를 어떻게 구하는지 알아봅시다

 

자세한 내용은 Hamming weight 위키피디아를 참고하시면 됩니다

https://en.wikipedia.org/wiki/Hamming_weight

 

  int nminussumofbits(int v)
  {
      long w = (long)v;
      w -= (0xaaaaaaaa & w) >> 1;
      w = (w & 0x33333333) + ((w >> 2) & 0x33333333);
      w = w + (w >> 4) & 0x0f0f0f0f;
      w += w >> 8;
      w += w >> 16;
      return v - (int)(w & 0xff);
  }

비트 수는 32비트까지 한정한 코드입니다

 

      w -= (0xaaaaaaaa & w) >> 1;

 

이 부분만 이해하면 나머지는 쉽습니다

 

간단하게 비트 수가 32개가 아닌 2개라고 가정해봅시다

우리는 어떤 함수에 의해

00 -> 00

01 -> 01

10 -> 01

11 -> 10

이 되는 함수를 찾아야 합니다

 

첫 번째 비트를 두 번째 비트에 빼면 됩니다

 

w -= (0xaaaaaaaa & w) >> 1;

위 수식은 두 개씩 비트를 잘라서 1의 개수를 구한 것입니다

 

11 01 11 01 00 10 11 10 01 11 01 10 11 11 11 00 이라는 숫자가 있다면 위 수식에 의해

10 01 10 01 00 01 10 01 01 10 01 01 10 10 10 00 이 됩니다

 

두 번째는 4개씩 비트를 나눠서 1의 개수를 구하게 됩니다

(10 + 01)  (10 + 01) (00 + 01) (10 + 01) (01 + 10) (01 + 01) (10 + 10) (10 + 00)

0011 0011 0001 0011 0011 0010 0100 0010 이 됩니다

 

세 번째는 8개씩 비트를 나눠서 1의 개수를 구하게 됩니다

(0011 + 0011) (0001 + 0011) (0011 + 0010) (0100 + 0010)

00000110 00000100 00000101 00000110 이 됩니다

 

그다음은,

1번째 슬라이스와 2번째 슬라이스의 합을 2번째 슬라이스에 저장하고,

3번째 슬라이스와 4번째 슬라이스의 합을 4번째 슬라이스에 저장합니다

xxxxxxxx 00001010 xxxxxxxx 00001011

x친 부분은 연산으로 제거될 예정이므로 고려하지 않아도 됩니다

 

마지막으로는 2번째 슬라이스와 4번째 슬라이스를 더합니다

xxxxxxxx xxxxxxxx xxxxxxxx 00010101

w & 0xff 에 의해 x가 다 소거됩니다

 

나머지 곱셈 구하기

2를 인수분해 한 뒤, 나머지 곱셈을 구하는 방법을 알아보겠습니다

 

인수 분해하면 위와같은 모양이 됩니다

규칙성을 찾기 힘들어 변형을 해줍니다

 

 

초록색은 인수분해되어 나온 숫자입니다

 

세 번째 줄을 봅시다(4)

세번째 줄에는 18/4 보다 큰 홀수가 나올 수 없습니다

 

따라서, 다음과 같이 표현할 수 있습니다

n=18
result = 1
for i in range(n.bit_length()):
    result *= rec_product(1, (n>>i) + 1|1)
    # rec_product(1, 5) = 1*3 (exclude right border)

+1|1 은 n>>i 결과가 짝수면 +1, 홀수면 +2를 해줍니다

덧셈이 비트 연산보다 먼저 수행되니 참고해주세요

 

하지만 이렇게 되면,

초록색 1 줄과 초록색 2 줄을 곱할 때

1*3*5*7*9가 중복입니다

 

약간 최적화를 해주면,

def binsplit_factorial(n):
   """Factorial of nonnegative integer n, via binary split.
   """
   inner = outer = 1
   for i in range(n.bit_length(), -1, -1):
       inner *= rec_product((n >> i + 1) + 1 | 1, (n >> i) + 1 | 1)
       outer *= inner
   return outer << (n - num_of_set_bits(n))

이렇게 됩니다

 

또한, rec_product도 최적화를 위해 recursion이 아닌 divide and conquer 방식을 채택했습니다

 

질문 사항 있다면 댓글로 남겨주세요

감사합니다