【AHOI2013复仇】数论之神

Posted on 2013-03-15 19:24 Mato_No1 阅读(1185) 评论(2)  编辑 收藏 引用 所属分类: 数论
原题地址
这题真是太神犇了……可以让人完全搞懂数论同余部分的全部内容……

题解……由于虹猫大神已经在空间里写得很详细了,所以就不肿么写了囧……
主要说一下一些难想的和容易搞疵的地方:
(1)中国剩余定理的那个推论(多个同余方程的模数互质,则整个方程组在小于所有模数之积的范围内的解数等于各个方程解数之积)其实是很强大的,不光对线性同余方程有用,对这种非线性的同余方程也有用,只需要所有方程都满足:若模数为MOD,则a是解当且仅当(a+MOD)是解……本题显然满足,因此,只要在按质因数拆分后求出各个方程的解数,再相乘即可(本沙茶就是这里木有想起来,结果看了虹猫的题解囧……);
(2)对于余数不为0且和模数不互质的情况要特别注意(这个好像很多标程都疵了,比如虹猫给的标程,不过数据弱,让它们过了囧),首先必须是余数含p(p为该方程模数的质因数)因子的个数j是a的倍数(也就是余数是p^a的倍数)才能有解,然后,当有解时,转化为解必须是p^(j/a)的倍数以及x/(p^(j/a))满足一个模数指数为原来指数减j的方程,这里需要注意,这个新方程的解数乘以p^(j-j/a)才是原来方程的解数!!道理很简单,因为模数除以了p^j,而x只除以了p^(j/a)……可以用一组数据检验:3 330750 6643012,结果是135而不是15;
(3)原根只能暴力求(不过最小原根都很小,1000以内的所有质数最小原根最大只有19……),但在求的时候有一个小的优化:首先p的原根也是p的任意整数次方的原根,然后求p的原根时,将(p-1)的非自身因数(预先求出)递减排序,这样可以比较快地排除不合法解;
(4)求逆元时一定要注意,如果得到的逆元是负数,要转化为正数,另外要取模;
(5)BSGS的时候一定要注意去重,在保留重复元素的情况下即使使用另一种二分查找也会疵的;
(6)数组不要开小了;

代码:
#include <iostream>
#include 
<stdio.h>
#include 
<stdlib.h>
#include 
<string.h>
#include 
<math.h>
#include 
<algorithm>
using namespace std;
#define re(i, n) for (int i=0; i<n; i++)
#define re1(i, n) for (int i=1; i<=n; i++)
#define re2(i, l, r) for (int i=l; i<r; i++)
#define re3(i, l, r) for (int i=l; i<=r; i++)
#define rre(i, n) for (int i=n-1; i>=0; i--)
#define rre1(i, n) for (int i=n; i>0; i--)
#define rre2(i, r, l) for (int i=r-1; i>=l; i--)
#define rre3(i, r, l) for (int i=r; i>=l; i--)
#define ll long long
const int MAXN = 110, MAXP = 50010, INF = ~0U >> 2;
int P_LEN, _P[MAXP + 1], P[MAXP + 1];
int A, B, M, n, DS[MAXN], DK[MAXN], R[MAXN], KR[MAXP], res;
struct sss {
    
int v, No;
    
bool operator< (sss s0) const {return v < s0.v || v == s0.v && No < s0.No;}
} Z[MAXP];
void prepare0()
{
    P_LEN 
= 0int v0;
    
for (int i=2; i<=MAXP; i++) {
        
if (!_P[i]) P[P_LEN++= _P[i] = i; v0 = _P[i] <= MAXP / i ? _P[i] : MAXP / i;
        
for (int j=0; j<P_LEN && P[j]<=v0; j++) _P[i * P[j]] = P[j];
    }
}
void prepare()
{
    n 
= 0int M0 = M;
    re(i, P_LEN) 
if (!(M0 % P[i])) {
        DS[n] 
= P[i]; DK[n] = 1; M0 /= P[i]; while (!(M0 % P[i])) {DK[n]++; M0 /= P[i];} n++;
        
if (M0 == 1break;
    }
    
if (M0 > 1) {DS[n] = M0; DK[n++= 1;}
    
int x;
    re(i, n) {
        x 
= 1; re(j, DK[i]) x *= DS[i];
        R[i] 
= B % x;
    }
}
ll pow0(ll a, 
int b, ll MOD)
{
    
if (b) {ll _ = pow0(a, b >> 1, MOD); _ = _ * _ % MOD; if (b & 1) _ = _ * a % MOD; return _;} else return 1;
}
void exgcd(int a, int b, int &x, int &y)
{
    
if (b) {
        
int _x, _y; exgcd(b, a % b, _x, _y);
        x 
= _y; y = _x - a / b * _y;
    } 
else {x = 1; y = 0;}
}
int gcd(int a, int b)
{
    
int r = 0while (b) {r = a % b; a = b; b = r;} return a;
}
void solve()
{
    
int x, y; res = 1;
    re(i, n) 
if (!R[i]) {
        
if (DK[i] < A) x = 1else x = (DK[i] - 1/ A + 1;
        re2(j, x, DK[i]) res 
*= DS[i];
    } 
else if (!(R[i] % DS[i])) {
        x 
= 0while (!(R[i] % DS[i])) {R[i] /= DS[i]; x++;}
        
if (x % A) {res = 0return;} else {
            DK[i] 
-= x; y = x / A;
            re2(j, y, x) res 
*= DS[i];
        }
    }
    
int phi, m0, m1, KR_len, _r, v0, _left, _right, _mid, T; bool FF;
    re(i, n) 
if (R[i]) {
        x 
= DS[i] - 1; KR_len = 0;
        
for (int j=2; j*j<=x; j++if (!(x % j)) {
            KR[KR_len
++= j;
            
if (j * j < x) KR[KR_len++= x / j;
        }
        KR[KR_len
++= 1;
        re2(j, 
2, DS[i]) {
            FF 
= 1;
            rre(k, KR_len) {
                _r 
= (int) pow0(j, KR[k], DS[i]);
                
if (_r == 1) {FF = 0break;}
            }
            
if (FF) {x = j; break;}
        }
        phi 
= DS[i] - 1; re2(j, 1, DK[i]) phi *= DS[i]; v0 = phi / (DS[i] - 1* DS[i];
        m0 
= (int) ceil(sqrt(phi) - 1e-10);
        Z[
0].v = 1; Z[0].No = 0; re2(j, 1, m0) {Z[j].v = (ll) Z[j - 1].v * x % v0; Z[j].No = j;}
        _r 
= (ll) Z[m0 - 1].v * x % v0; sort(Z, Z + m0);
        m1 
= 1; re2(j, 1, m0) if (Z[j].v > Z[j - 1].v) Z[m1++= Z[j];
        exgcd(_r, v0, x, y); 
if (x < 0) x += v0; y = R[i];
        re(j, m0) {
            _left 
= 0; _right = m1 - 1; T = -1;
            
while (_left <= _right) {
                _mid 
= _left + _right >> 1;
                
if (y == Z[_mid].v) {T = j * m0 + Z[_mid].No; break;}
                
else if (y < Z[_mid].v) _right = _mid - 1else _left = _mid + 1;
            }
            
if (T >= 0breakelse y = (ll) y * x % v0;
        }
        x 
= gcd(A, phi); if (T % x) {res = 0break;} else res *= x;
    }
}
int main()
{
    
int tests;
    scanf(
"%d"&tests);
    prepare0();
    re(testno, tests) {
        scanf(
"%d%d%d"&A, &B, &M); M += M + 1; B %= M;
        
if (!A) {
            
if (B == 1) res = M; else res = 0;
        } 
else {
            prepare();
            solve();
        }
        printf(
"%d\n", res);
    }
    
return 0;
}

Feedback

# re: 【AHOI2013复仇】数论之神  回复  更多评论   

2013-04-08 10:30 by zrz1996
能发一下原题吗?zrz96@163.com 谢谢

# re: 【AHOI2013复仇】数论之神  回复  更多评论   

2013-06-14 17:02 by SillyJ
能发一下原题吗?sillyhook00@126.com 谢谢

只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理