cpp/math/rho.cpp

• View this file on GitHub
• Last update: 2024-01-29 18:19:52+09:00

Depends on

• cpp/math/binary-power-method.cpp
• cpp/math/miller-rabin.cpp

Required by

• cpp/geometry/two-points-distance.cpp

Verified with

• cpp/z_test/yosupo-factorize.test.cpp

Code

#include "./miller-rabin.cpp"

using namespace std;
template <typename T>
struct Rho{
  mt19937 mt; //32 bit version
  T N;
  vector<T> factor;
  //std::mt19937_64 mt(rnd()); //64 bit version
  Rho(T n):N(n){
    random_device rnd;
    mt.seed(rnd());
  }


  vector<T> run(){
    factor = factors(N);
    sort(factor.begin(), factor.end());
    return factor;
  }

private:
  __int128 c;
  T f(__int128 x, T n){
    return (x*x + c)%n;
  }

  T find_factor(T n){
    c = mt()%n;
    T base = mt()%n;
    T d = 1;
    T x = base;
    T y = base;
    while(true){
      x = f(x, n);
      y = f(f(y,n),n);
      d = __gcd(abs(x-y), n);
      if(d == n){
        return -1;
      }else if(d != 1){
        return d;
      }
    }
  }

  vector<T> factors(T n){
    if(n == 1)return {};
    if(n == 4)return {2, 2};
    if(MillerRabinCheck(n)){
      return {n};
    }

    T factor = -1;
    while(factor == -1){
      factor = find_factor(n);
    }

    vector<T> f1 = factors(factor);
    vector<T> f2 = factors(n/factor);
    f1.insert(f1.end(), f2.begin(), f2.end());

    return f1;
  }
};

#line 2 "cpp/math/binary-power-method.cpp"

template <typename T>
T uPow(T z,T n, T mod){
  T ans = 1;
  while(n != 0){
    if(n%2){
      ans*=z;
      if(mod)ans%=mod;
    }
    n >>= 1;
    z*=z;
    if(mod)z%=mod;
  }
  return ans;
}

#line 2 "cpp/math/miller-rabin.cpp"
/*
true: 素数
false: 合成数
*/
template<typename T>
bool MillerRabinCheck(T n){
  if(n == 1)return false;
  if(n%2 == 0){
    return n == 2;
  }
  __int128 d = n-1;
  __int128 s = 0;
  while(d%2 == 0){
    d/=2;
    s++;
  }
  vector<__int128> base = {2,3,5,7,11,13,17,19,23,29,31,37};
  for(int i = 0; base.size() > i; i++){
    if(base[i] >= n)break;
    long long current = (long long)uPow(base[i],d,(__int128)n);
    if(current == 1 || current == n-1)continue;
    bool ok = false;
    for(int j = 0; s > j; j++){
      current = ((__int128)current * (__int128)current) % n;
      if(current == n-1){
        ok = true;
        break;
      }
    }
    if(!ok)return false;
  }
  return true;
}
#line 2 "cpp/math/rho.cpp"

using namespace std;
template <typename T>
struct Rho{
  mt19937 mt; //32 bit version
  T N;
  vector<T> factor;
  //std::mt19937_64 mt(rnd()); //64 bit version
  Rho(T n):N(n){
    random_device rnd;
    mt.seed(rnd());
  }


  vector<T> run(){
    factor = factors(N);
    sort(factor.begin(), factor.end());
    return factor;
  }

private:
  __int128 c;
  T f(__int128 x, T n){
    return (x*x + c)%n;
  }

  T find_factor(T n){
    c = mt()%n;
    T base = mt()%n;
    T d = 1;
    T x = base;
    T y = base;
    while(true){
      x = f(x, n);
      y = f(f(y,n),n);
      d = __gcd(abs(x-y), n);
      if(d == n){
        return -1;
      }else if(d != 1){
        return d;
      }
    }
  }

  vector<T> factors(T n){
    if(n == 1)return {};
    if(n == 4)return {2, 2};
    if(MillerRabinCheck(n)){
      return {n};
    }

    T factor = -1;
    while(factor == -1){
      factor = find_factor(n);
    }

    vector<T> f1 = factors(factor);
    vector<T> f2 = factors(n/factor);
    f1.insert(f1.end(), f2.begin(), f2.end());

    return f1;
  }
};

Back to top page