求调质因数分解
  • 板块学术版
  • 楼主Sya_Resory
  • 当前回复3
  • 已保存回复3
  • 发布时间2022/9/19 23:37
  • 上次更新2023/10/27 10:33:43
查看原帖
求调质因数分解
114082
Sya_Resory楼主2022/9/19 23:37

RT 题目链接

已经跑不过暴力了(悲

做法

code:

#include <algorithm>
#include <bitset>
#include <cstdio>
#include <random>
#include <vector>
#include <cmath>
#include <map>

using i64 = long long;
using u64 = unsigned long long;
using i128 = __int128;
using u128 = __uint128_t;
using std::min; using std::swap;

struct u256 {
    u128 lo,hi;
    u256(u128 low = 0,u128 high = 0): lo(low),hi(high) {}

    static u256 Mul128(u128 x,u128 y) {
        u64 x_hi,x_lo,y_hi,y_lo,t_hi; u128 p1,p2,p3;
        x_hi = x >> 64,x_lo = u64(x);
        y_hi = y >> 64,y_lo = u64(y);
        p1 = u128(x_lo) * y_lo,p2 = u128(x_hi) * y_lo + u64(p1 >> 64);
        t_hi = p2 >> 64,p2 = u128(x_lo) * y_hi + u64(p2);
        p3 = u128(x_hi) * y_hi + u64(p2 >> 64) + t_hi;
        return u256(u64(p1) | (p2 << 64),p3);
    }
};

struct Mont {
    u128 mod,inv,r2;
    Mont(u128 n = 1): mod(n) {
        inv = n;
        for(int i = 0;i < 6;i ++) inv *= 2 - inv * n;
        r2 = -n % n;
		for(int i = 0;i < 4;i ++) if ((r2 <<= 1) >= mod) r2 -= mod;
		for(int i = 0;i < 5;i ++) r2 = Mul(r2,r2);
    }

    inline u128 redc(u256 x) {
        u128 y = x.hi - u256::Mul128(x.lo * inv,mod).hi;
        return i128(y) < 0 ? y + mod : y;
    }
    inline u128 redc(u128 x) { return redc(u256(x,0)); }
    inline u128 init(u128 n) { return redc(u256::Mul128(n,r2)); }
    inline u128 Mul(u128 a,u128 b) { return redc(u256::Mul128(a,b)); }
} mont;

const int Y = 2e5,P = 1.5e4;
using vec = std::bitset < P >;

i128 N;
int pcnt; bool isp[Y]; u128 pri[P];
std::mt19937 rnd(1919810);
vec v[P]; u128 a[P];

template < class T >
inline T read() {
#define gc c = getchar()
    T d = 0; int f = 0,gc;
    for(;c < 48 || c > 57;gc) f |= (c == '-');
    for(;c > 47 && c < 58;gc) d = d * 10 + (c ^ 48);
#undef gc
    return f ? -d : d;
}
template < class T >
inline void write(T x) {
    if(x < 0) putchar('-'),x = -x;
    if(x > 9) write(x / 10);
    putchar((x % 10) ^ 48);
}

inline int ctz(u128 x) {
    return u64(x) ? __builtin_ctzll(x) : __builtin_ctzll(x >> 64) + 64;
}

inline u128 gcd(u128 a,u128 b) {
    if(!a || !b) return a | b;
    int d = ctz(a | b);
    for(a >>= ctz(a),b >>= ctz(b);a;a -= b,a >>= ctz(a)) if(a < b) swap(a,b);
    return b << d;
}

inline u64 randint(u64 N = 1ull << 63) { return std::uniform_int_distribution < u64 > (1,N - 1)(rnd); }

inline u128 randint(u128 N) { return randint() * randint() % N; }

inline void Init(u128 N) {
    mont = Mont(N);
    int y = exp(sqrt(log(N) * log(log(N))));
    for(int i = 2;i <= y;i ++) {
        if(!isp[i]) pri[++ pcnt] = i;
        for(int t,j = 1;j <= pcnt && (t = i * pri[j]) <= y;j ++) {
            isp[t] =  true; if(!(i % pri[j])) break;
        }
    }
}

inline bool gauss() {
    for(int i = 1;i <= pcnt;i ++) {
        if(!v[i][i]) continue;
        for(int j = 1;j <= pcnt + 1;j ++) {
            if(i == j || !v[j][i]) continue;
            for(int k = 1;k <= pcnt;k ++) v[j][k] = v[i][k] ^ v[j][k];
        }
    }
    u128 al = mont.init(1),be = al,ga;
    for(int i = 1;i <= pcnt + 1;i ++) if(v[i].none()) al = mont.Mul(al,a[i]);
    ga = mont.redc(mont.Mul(al,al)),al = mont.redc(al);
    for(int i = 1;i <= pcnt;i ++)
        for(;!(ga % pri[i]);ga /= pri[i] * pri[i]) be = mont.Mul(be,mont.init(pri[i]));
    be = mont.redc(be); u128 ans = gcd((al + be) % N,N);
    if(ans != 1 && ans != N) {
        u128 p = ans,q = N / p;
        if(p > q) swap(p,q);
        write(p),putchar(' '),write(q);
        return true;
    } return false;
}

inline void Solve(u128 N) {
    for(;;) {
        int cnt = 0;
        for(;cnt <= pcnt;) {
            u128 x = mont.init(randint(N)),y = mont.redc(mont.Mul(x,x)); vec z;
            for(int i = 1;i <= pcnt;i ++) for(;!(y % pri[i]);y /= pri[i]) z[i] = !z[i];
            if(y == 1) a[++ cnt] = x,v[cnt] = z;
        }
        if(gauss()) return ;
    }
}

int main() {
    N = read < u128 > ();
    Init(N);
    Solve(N);
    return 0;
}

@JS_TZ_ZHR 试图诈骗(

2022/9/19 23:37
加载中...