====== 소인수분해 (Prime Factorization) ======
  * Integer Factorization 이라고도 불린다 (이게 더 흔히 쓰이는 표현같기도..)

===== n을 소인수분해하기 =====
  * 참고 - [[https://cp-algorithms.com/algebra/factorization.html]] [[https://aruz.tistory.com/entry/factorization-1]]
  * 소인수분해는 암호학이나 여러 분야에서 매우 중요하게 사용되는 알고리즘이고, 그러다보니 다양한 알고리즘이 있는데.. PS 범위에서는 trial division과 Pollard's rho 두가지만 알면 충분하다
    * 실제로는 더 효율적인 알고리즘들이 있지만.. 2^64 정도의 범위에서는  Pollard's rho 방법으로도 충분히 잘 돌아가는것 같다.
    * [[https://cp-algorithms.com/algebra/factorization.html#fermats-factorization-method|페르마의 방법]] 이나 [[https://cp-algorithms.com/algebra/factorization.html#pollards-p-1-method|Pollard's p-1]] 등은 Pollard's rho보다 장점이 없기 때문에 굳이 알 필요가 없지만, 소인수분해에 대한 이해를 위해서라면 읽어봐도 좋다
    * [[https://koosaga.com/239]]은 더 효율적일수도 있지만, 복잡해서 PS에서는 못쓸것 같은 알고리즘이다. 한번 봐둬도 도움은 될듯.
  * 소인수분해 결과는 p_1^k_1 * p_2^k_2 * ... * p^n^k_n 의 형태로 나온다. 코드를 구현할때는 리턴값으로 p_i를 키로 k_i를 밸류로 갖는 dict를 리턴하도록 만들었다. n이 1일때는 빈 딕셔너리를 리턴한다

==== Trial division ==== 
  * n이 적당히 작은 경우에 효율적이다. 
  * 그냥 i=2부터 시작해서 i*i<n이 될때까지 전부 나눠보는 방식. 소수일 경우에 최대 sqrt(n)번의 루프를 돌아야 한다. 시간복잡도는 O(sqrt(n))
  * 소수 체크를 할때와는 다르게 루프를 2~sqrt(n)까지 도는게 아니라, 매번 i*i를 계산해서 n을 넘으면 중단하는 식으로 구현한다. 그 이유는 소인수를 발견하면 n을 그값으로 나누는 작업때문에 n이 점점 작아지기 때문에 이렇게 처리하는것이 보통 더 효율적이기 때문이다. 하지만, n이 소수인 경우에는 사실 이 과정이 불필요하고, 그냥 sqrt(n)을 먼저 구해놓고 거기까지 루프를 도는게 빠르긴 하다. 즉.. 평균적으로는 i*i 로 브레이크하는게 빠르지만, 최악경우에는 sqrt(n) 까지 돌리는게 빠르다는 의미.. 그치만 속도 차이는 별로 없다 
  * 역시 wheel optimization 도 가능하다. 2의 배수만 따로 전처리하면 빨라진다. 그 이상은 안해봄
  * 사실 이 로직을 최적화하는 것이 중요해질 정도로 n이 커질 경우에는 아래에서 설명할 Pollard's Rho를 적용할 것이기 때문에, 굳이 상수최적화에 진심을 다할 필요는 없긴 하다..
  * 코드
    * [[:ps:teflib:numtheory#prime_factorization_small|teflib.numtheory.prime_factorization_small]]
      * 2-wheel만 적용.
  * 관련 문제 
    * [[ps:problems:boj:11653]] : n<=10^7

==== Pollard's Rho ====
  * 기본적인 설명은 [[https://blog.naver.com/jinhan814/222141831551]] 과 [[https://gratus907.github.io/advanced-algorithms/Pollards-Rho/]]를 참고.
  * Brent가 발전시킨 알고리즘에 대해서는 [[https://cp-algorithms.com/algebra/factorization.html#brents-algorithm]] 을 참고
  * Brent의 알고리즘은 cycle를 찾는 방식을 개선한것이다. 원래는 [[https://en.wikipedia.org/wiki/Cycle_detection#Floyd's_tortoise_and_hare|Floyd's cycle-finding algorithm]]을 썼는데, 그 대신 Brent's cycle detection 방식을 사용한 것이다. 시간복잡도가 달라지지는 않지만 f의 계산횟수와 gcd 계산 횟수를 줄일수 있다.
    * 원 논문에서는 플로이드에 비해서 36%정도 빠르게 사이클을 찾고, 폴라드로 알고리즘 전체에서는 24% 정도의 속도 향상을 가져왔다고 한다 [[https://en.wikipedia.org/wiki/Cycle_detection#Brent's_algorithm]]

=== 시간복잡도 ===
  * Pollard's Rho 알고리즘의 시간복잡도는 평균 O(n^1/4)라고 부르지만, 모든 수에 대해서 성립하는지는 증명되지 않았다고 한다
    * O(n^1/4) 번의 루프안에서는 매번 gcd 연산이 들어가는데, 이것을 포함해서 O(n^1/4 * logn)이라고 쓰는 경우도 있기는 하다. 하지만 Brent 방식을 쓰면 매 루프마다 gcd를 계산하지 않아도 된다.
  * Pollard's Rho 알고리즘 자체는 인수 하나를 찾아주는 알고리즘이므로, 소인수분해를 하기 위해서는 Pollard's Rho 를 반복해서 돌려야 한다. 인수 f를 찾은 다음에 n/f 와 f에 대해서 다시 Pollard's Rho 알고리즘을 돌리는 식으로 동작하게 되는데, 이렇게 반복해서 돌리는것을 고려해도 전체 소인수분해의 시간복잡도는 O(n^1/4)이 된다. 

=== 구현 ===
  * 구현상에서 주의할 점이 좀 많다.
  * Pollard's rho 알고리즘은 기본적으로 인수 하나를 찾는 알고리즘이다. 
  * 이 함수는 수도랜덤생성함수 f와 (보통 f(x)=x*x+c 로 준다), 초기값 x0를 갖고서 작업을 수행하게 되는데, 기본적으로는 아래처럼 구현된다. 
    * <dkpr py>
def pollard_rho(n, x0, c):
    x, y = x0, (x0 * x0 + c) % n            
    while (g := math.gcd(x - y, n)) == 1:
        x, y = (x * x + c) % n, (y * y + c) % n
        y = (y * y + c) % n
    return g
</dkpr>
      * cpp 계열의 gcd 함수는 두 인자가 모두 양수여야 해서, abs(x-y)처럼 써줘야 하지만, python의 math.gcd는 abs를 안붙여도 제대로 된 값을 준다.
      * cpp에서는 x*x를 계산할때 오버플로우되지 않게 하기 위해서 추가 구현이 필요하지만, 파이썬에서는 신경쓸 필요가 없다.
      * x0와 c는 랜덤으로 생성된 값을 주거나, 그냥 고정된 값을 주는 경우도 있다. 고정된 값을 쓴다면 c=1, x0=2 를 쓰는것이 일반적이다
  * 하지만 이 함수를 그대로 쓰기에는 문제가 있다. 
  * 우선, 운나쁘게도 찾아낸 인수 g가 n과 같은 경우가 있다. 이 경우에는 f와 x0를 바꿔서 다시 돌리는 것이 정석이다. f(x)=x*x+c 의 형태는 유지한 채로, c값과 x0를 바꿔서 돌리면 된다. 랜덤을 다시 돌려서 찾거나, 이전값에서 1을 증가시키던가 해서 처리한다.
    * 이때 주의할점은, x0만 바꿔서는 돌려서는 안된다. 대표적인 예가 124376107291 을 인수분해 하는 경우인데, f=x*x+1 을 사용해서는 초기값을 아무리 바꿔도 인수를 찾아내지 못하고, 무한루프로 빠지게 된다. 차라리, x0는 안바꾸고 c만 바꾸는 것은 동작한다.
    * [[https://lpha-z.hatenablog.com/entry/2023/01/15/231500]] 를 참고하라
    * [[https://github.com/cheran-senthil/PyRival/blob/master/pyrival/algebra/factors.py|pyrival]] 에서도 f를 고정시키고 x0만 바꾸는 방식의 잘못된 구현을 사용하고 있다.
  * 다음으로는 n이 소수일 경우이다. 이때는 인수가 진짜로 없으므로 f와 x0를 아무리 바꿔서 다시 돌려도 당연히 인수를 찾지 못한다. (g가 n으로 나온다). 따라서 폴라드 로를 돌리기 전에 미리 n이 소수인지 [[ps:소수 판별]]을 하고 소수가 아닌 경우에만 폴라드로를 돌리도록 한다.
    * 그래서, 폴라드로 알고리즘의 난이도는 폴라드로 자체보다도 오히려 소수판별을 위해 쓰는 밀러라빈의 난이도가 더 높다는 의견도 있다.
    * 지금처럼 밀러라빈+폴라드로를 세트로 묶어서 사용하는 것이 일반화되기 이전에는, 소수판별을 미리 하지 않고 그냥 돌려서 g==n일 경우 f와 x0를 바꿔서 몇번 더 돌려보고 계속 g==n이 나올경우에는 소수라고 간주하는 방식으로 돌리기도 했다고 한다.
  * 이부분을 추가하면 다음과 같아진다
    * <dkpr py>
def find_factor(n):
    if is_prime(n):
        return n
    x0 = 2
    for c in itertools.count(1):
        g = pollard_rho(n, x0, c)
        if g != n:
            return g
</dkpr>  
  * 소인수분해를 하기 위해서는 폴라드로 알고리즘으로 한개의 인수를 찾는것을 반복적으로 수행해야 한다
  * n의 인수 f를 찾고 나면, f와 n/f에 대해서 다시 폴라드로를 돌려야 한다. 가장 간단한 것은 재귀로 구현하는 것이다. 아래는 pyrival의 구현체이다
    * <dkpr py>
def prime_factors(n):
    """returns a Counter of the prime factorization of n"""
    if n <= 1:
        return Counter()
    f = find_factor(n)
    return Counter([n]) if f == n else prime_factors(f) + prime_factors(n // f)
</dkpr>
  * pollard_rho(n, x0, c) 함수를 brent의 방법으로 바꿔서 속도를 더 빠르게 할수있다.
  * [[https://cp-algorithms.com/algebra/factorization.html#brents-algorithm]] 의 코드를 거의 그대로 python으로 옮기면 아래처럼 된다.
    * <dkpr py>
def brent_pollard(n,x0,c):
    m, l, x, g, q = 128, 1, x0, 1, 1
    while g == 1:
        y = x
        for _ in range(l - 1):
            x = (x * x + c) % n
        k = 0
        while k < l and g == 1:
            xs = x
            for _ in range(min(m, l - k)):
                x = (x * x + c) % n
                q = q * abs(y - x) % n
            g = math.gcd(q, n)
            k += m
        l += l
    if g != n:
        return g
    g = 1
    while g == 1:
        xs = (xs * xs + 1) % n
        g = math.gcd(xs - y, n)
    return g
</dkpr>
    * m 값에 대해서는 여기에서는 128을 줬는데, 1000을 사용하는 코드도 있다.
  * [[https://qiita.com/t_fuki/items/7cd50de54d3c5d063b4a]] 에서는 좀더 python에 맞춰 정리된 brent 알고리즘 코드가 있다. m을 n^(1/8)로 잡고 있다.
  * [[https://judge.yosupo.jp/submission/124476]]에서는 이 코드를 다시 깔끔하게 최적화한 코드를 볼수 있다.

=== 코드 ===
  * [[:ps:teflib:numtheory#prime_factorization|teflib.numtheory.prime_factorization]]
  * 폴라드로 + brent 구현부분은 [[https://judge.yosupo.jp/submission/124476]]를 많이 참고했다
  * 폴라드로를 반복적으로 돌리는 부분은 재귀 대신 반복으로 구현했다

=== 관련 문제 ===
  * ...

===== 여러개의 수를 소인수분해하기 =====
  * x_1, x_2, ..., x_k 를 모두 소인수분해하기. 모든 수는 n이하라고 생각하자.
  * 소인수분해해야 할 수의 갯수 k에 따라서 다음중 적당한 방법들을 선택하면 된다 (n은 trial division으로 소인수분해가 가능한 범위라고 새생각하자)
    * k가 작으면, 그냥 위의 trial division 방법을 모든 x_i에 대해서 적용하면 된다. O(k*n^1/2)
    * k가 적당히 크면, 먼저 sqrt(n)이하의 소수 목록을 모두 계산해놓고, 그것에 대해서만 trial division을 사용하는 게 좀더 빠를수 있다
      * 에라토스테네스의 체를 이용해서 소수 목록을 구하는 데에 O(sqrt(n)*loglogn)
      * 목록의 소수들만 이용해서 소인수분해를 하는데 O(sqrt(n)/ln(n))
      * k개를 처리하는 총 시간복잡도는 대략 O(sqrt(n)*(loglogn + k/logn)). 
      * 관련 문제 - [[ps:problems:boj:2312]]
    * k가 크면 (n에 가까울 정도로), 아래의 방법을 사용하자 O(n + klogn)

===== [1,n]까지의 모든 수를 소인수분해하기 =====
  * 이 경우는, 모든 수에 대해서 최소인수를 미리 전처리해놓는 방식으로 구현한다.
  * 모든 수에 대해서 최소인수를 구하는 것은 보통 [[ps:소수 목록 구하기#linear sieve]]를 이용하는 것으로 설명하는 자료가 많은데, [[ps:소수 목록 구하기#에라토스테네스의 체]]를 써서도 가능하긴 하다. 하지만 [[ps:소수 목록 구하기#linear sieve]]를 이용하는 것이 조금 더 빠르다. 시간복잡도는 O(n).
  * 모든수의 최소인수를 구하고 나면, 임의의 수 x의 소인수분해를 O(Ω(x))에 처리할 수 있다. Ω(x)([[wp>Prime omega function]])은 x 소인수의 개수를 의미하는데, 최악의 경우 O(logx), 평균적으로는 O(loglogx)이다.
  * 결과적으로 1~n까지의 모든 수를 소인수분해하는 것도 O(nloglogn) 에 가능하다 
  * 코드
    * [[:ps:teflib:numtheory#PrimeFactorizationCalculator|teflib.numtheory.PrimeFactorizationCalculator]]
  * 관련 문제
    * [[ps:problems:boj:16563]]
    * [[ps:problems:boj:3142]]