QOJ.ac

QOJ

IDProblemSubmitterResultTimeMemoryLanguageFile sizeSubmit timeJudge time
#841394#9923. Ma Meilleure Ennemieucup-team3586#WA 884ms6552kbC++237.4kb2025-01-03 18:16:092025-01-03 18:16:10

Judging History

你现在查看的是最新测评结果

  • [2025-01-03 18:16:10]
  • 评测
  • 测评结果:WA
  • 用时:884ms
  • 内存:6552kb
  • [2025-01-03 18:16:09]
  • 提交

answer

#include<bits/stdc++.h>
#pragma GCC optimize("Ofast")
#pragma GCC optimize("unroll-loops")
#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2,tune=native")
uint64_t gcd_stein_impl( uint64_t x, uint64_t y ) {
    if( x == y ) { return x; }
    const uint64_t a = y - x;
    const uint64_t b = x - y;
    const int n = __builtin_ctzll( b );
    const uint64_t s = x < y ? a : b;
    const uint64_t t = x < y ? x : y;
    return gcd_stein_impl( s >> n, t );
}

uint64_t gcd_stein( uint64_t x, uint64_t y ) {
    if( x == 0 ) { return y; }
    if( y == 0 ) { return x; }
    const int n = __builtin_ctzll( x );
    const int m = __builtin_ctzll( y );
    return gcd_stein_impl( x >> n, y >> m ) << ( n < m ? n : m );
}

// ---- is_prime ----

uint64_t mod_pow( uint64_t x, uint64_t y, uint64_t mod ) {
    uint64_t ret = 1;
    uint64_t acc = x;
    for( ; y; y >>= 1 ) {
        if( y & 1 ) {
            ret = __uint128_t(ret) * acc % mod;
        }
        acc = __uint128_t(acc) * acc % mod;
    }
    return ret;
}

bool miller_rabin( uint64_t n, const std::initializer_list<uint64_t>& as ) {
    return std::all_of( as.begin(), as.end(), [n]( uint64_t a ) {
        if( n <= a ) { return true; }

        int e = __builtin_ctzll( n - 1 );
        uint64_t z = mod_pow( a, ( n - 1 ) >> e, n );
        if( z == 1 || z == n - 1 ) { return true; }

        while( --e ) {
            z = __uint128_t(z) * z % n;
            if( z == 1 ) { return false; }
            if( z == n - 1 ) { return true; }
        }

        return false;
    });
}

bool is_prime( uint64_t n ) {
    if( n == 2 ) { return true; }
    if( n % 2 == 0 ) { return false; }
    if( n < 4759123141 ) { return miller_rabin( n, { 2, 7, 61 } ); }
    return miller_rabin( n, { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 } );
}

// ---- Montgomery ----

class Montgomery {
    uint64_t mod;
    uint64_t R;
public:
    Montgomery( uint64_t n ) : mod(n), R(n) {
       for( size_t i = 0; i < 5; ++i ) {
          R *= 2 - mod * R;
       }
    }

    uint64_t fma( uint64_t a, uint64_t b, uint64_t c ) const {
        const __uint128_t d = __uint128_t(a) * b;
        const uint64_t    e = c + mod + ( d >> 64 );
        const uint64_t    f = uint64_t(d) * R;
        const uint64_t    g = ( __uint128_t(f) * mod ) >> 64;
        return e - g;
    }

    uint64_t mul( uint64_t a, uint64_t b ) const {
        return fma( a, b, 0 );
    }
};

// ---- Pollard's rho algorithm ----

uint64_t pollard_rho( uint64_t n ) {
    if( n % 2 == 0 ) { return 2; }
    const Montgomery m( n );

    constexpr uint64_t C1 = 1;
    constexpr uint64_t C2 = 2;
    constexpr uint64_t M = 512;

    uint64_t Z1 = 1;
    uint64_t Z2 = 2;
retry:
    uint64_t z1 = Z1;
    uint64_t z2 = Z2;
    for( size_t k = M; ; k *= 2 ) {
        const uint64_t x1 = z1 + n;
        const uint64_t x2 = z2 + n;
        for( size_t j = 0; j < k; j += M ) {
            const uint64_t y1 = z1;
            const uint64_t y2 = z2;

            uint64_t q1 = 1;
            uint64_t q2 = 2;
            z1 = m.fma( z1, z1, C1 );
            z2 = m.fma( z2, z2, C2 );
            for( size_t i = 0; i < M; ++i ) {
                const uint64_t t1 = x1 - z1;
                const uint64_t t2 = x2 - z2;
                z1 = m.fma( z1, z1, C1 );
                z2 = m.fma( z2, z2, C2 );
                q1 = m.mul( q1, t1 );
                q2 = m.mul( q2, t2 );
            }
            q1 = m.mul( q1, x1 - z1 );
            q2 = m.mul( q2, x2 - z2 );

            const uint64_t q3 = m.mul( q1, q2 );
            const uint64_t g3 = gcd_stein( n, q3 );
            if( g3 == 1 ) { continue; }
            if( g3 != n ) { return g3; }

            const uint64_t g1 = gcd_stein( n, q1 );
            const uint64_t g2 = gcd_stein( n, q2 );

            const uint64_t C = g1 != 1 ? C1 : C2;
            const uint64_t x = g1 != 1 ? x1 : x2;
            uint64_t       z = g1 != 1 ? y1 : y2;
            uint64_t       g = g1 != 1 ? g1 : g2;

            if( g == n ) {
                do {
                    z = m.fma( z, z, C );
                    g = gcd_stein( n, x - z );
                } while( g == 1 );
            }
            if( g != n ) {
                return g;
            }

            Z1 += 2;
            Z2 += 2;
            goto retry;
        }
    }
}

void factorize_impl( uint64_t n, std::vector<uint64_t>& ret ) {
    if( n <= 1 ) { return; }
    if( is_prime( n ) ) { ret.push_back( n ); return; }

    const uint64_t p = pollard_rho( n );

    factorize_impl( p, ret );
    factorize_impl( n / p, ret );
}

std::vector<uint64_t> factorize( uint64_t n ) {
    std::vector<uint64_t> ret;
    factorize_impl( n, ret );
    std::sort( ret.begin(), ret.end() );
    return ret;
}
using namespace std;
#define int long long
inline int read(){
   int s=0,w=1;
   char ch=getchar();
   while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
   while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
   return s*w;
}
const int p=998244353;
int qp(int x,int y)
{
	int res=1;
	for(int t=x; y; y>>=1,t=1ll*t*t%p)
		if(y&1) res=1ll*res*t%p;
	return res;
}
int pr[103],k[103],idx[103],sz;
int n=read(),m=read();
struct custom_hash {
    static uint64_t splitmix64(uint64_t x) {
        // http://xorshift.di.unimi.it/splitmix64.c
        x += 0x9e3779b97f4a7c15;
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
        x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
        return x ^ (x >> 31);
    }

    size_t operator()(uint64_t x) const {
        static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
        return splitmix64(x + FIXED_RANDOM);
    }
};
// unordered_map<int,int,custom_hash> d,f,g;
// unordered_set<int,custom_hash> st;
int f[1<<20],d[1<<20];
void dfs(int x,int y,int z,int px)
{
	if(y>sz){d[px]=z,f[px]=m/x;return;}
	for(int i=0; i<=k[y]; ++i,x*=pr[y])
		dfs(x,y+1,z*(i+1)%p,px+idx[y+1]*i);
}
void build()
{
	auto vec=factorize(n);
	int len=vec.size();
	reverse(vec.begin(),vec.end());
	for(int i=0; i<len; ++i)
	{
		int j=1;
		while(i+1<len&&vec[i]==vec[i+1]) ++i,++j;
		pr[++sz]=vec[i],k[sz]=j;
	}
	idx[sz+1]=1;
	for(int i=sz; i>=1; --i)
		idx[i]=idx[i+1]*(k[i]+1);
	dfs(1,1,1,0);
}
void calc(int v,int x,int y)
{
	if(x>sz)
	{
		f[v]=(f[v]+p-f[v+idx[y+1]])%p;
		return ;
	}
	for(int i=0; i<=k[x]-(x==y); ++i) calc(v+idx[x+1]*i,x+1,y);
}
void solve(int x,int y,int z,int px,int py)
{
	if(z>sz)
	{
		if(x==y) return ;
		f[px]=(f[px]+qp(f[py],x/y)*(p+3-d[px-py]))%p;
		return ;
	}
	for(int i=0; i<=k[z]; ++i,x*=pr[z])
		for(int j=0,ty=y; j<=i; ++j,ty*=pr[z])
			solve(x,ty,z+1,px+i*idx[z+1],py+j*idx[z+1]);
}
int cur[103];
void solve2(int val,int x,int y,int z)
{
	if(z>sz)
	{
		if(x==y) return ;
		f[x]=(f[x]+val*(p+3-d[x-y]))%p;
		return ;
	}
	for(int i=cur[z]; i<=k[z]; ++i)
	{
		if(i>cur[z]) val=qp(val,pr[z]);
		solve2(val,x+(i-cur[z])*idx[z+1],y,z+1);
	}
}
void solve1(int x,int y,int px)
{
	if(y>sz) return solve2(f[px],px,px,1);
	for(int i=0; i<=k[y]; ++i,x*=pr[y])
		cur[y]=i,solve1(x,y+1,px+i*idx[y+1]);
}
signed main()
{
	if(n==1) printf("%lld\n",m%p),exit(0);
	build();
	for(int i=1; i<=sz; ++i) calc(0,1,i);
	reverse(f,f+idx[1]);
	// solve(1,1,1,0,0);
	solve1(1,1,0);
	printf("%lld\n",f[idx[1]-1]);
	return 0;
}

Details

Tip: Click on the bar to expand more detailed information

Test #1:

score: 100
Accepted
time: 1ms
memory: 5844kb

input:

4 4

output:

6

result:

ok 1 number(s): "6"

Test #2:

score: 0
Accepted
time: 1ms
memory: 6116kb

input:

2338 1470

output:

18530141

result:

ok 1 number(s): "18530141"

Test #3:

score: -100
Wrong Answer
time: 884ms
memory: 6552kb

input:

941958815880242160 945059392259579928

output:

628117270

result:

wrong answer 1st numbers differ - expected: '57894579', found: '628117270'