标签：线性筛 inline 同余方程 resid -- eps break euc 论文
/*** ****
数论模板
1.EXGCD
2.CRT (互质 与 非互质)
3.逆元(线性预处理 欧拉（费小）定理)
4.筛 （线性筛与数论函数 状压筛）
5.大素数计数 (Meisell-Lehmer O2/3 O3/4 )
6.SG函数 mex方法
7.Gauss消元
8.FFT NTT FWT
9.二次剩余 (二次域方法)
10.素性测试 (Miller_rabin)
11.LUCAS 大组合取模
/** ****/
//1.EXGCD
//
//用于解 ax + by = c mod n x与y的值
//前提: (a,b) = d | c
//先求 ax‘+by‘=(a,b) 得特解x‘ y‘
//其中一组解
//x0 = x‘·c/d
//y0 = y‘·c/d
//通解:
//x = (x‘+ k·b/d)c/d
//y = (y‘- k·a/d)c/d
//最后最小 x=(x‘ mod b/d) c/d
LL exgcd(LL a, LL b, LL &x, LL &y)
{
LL d = a;
if(b == 0)
x = 1, y = 0;
else
{
d = exgcd(b, a % b, y, x);
y -= (a / b) * x;
}
return d;
}
/********************************************************/
//2.CRT
//
//用于求同余方程的 x ≡ r[i] mod m[i] x的值
//可用于化简 模数较大但可分解成多个小模数的公式,通过列出同余方程得到最终的X值最后取模
//
//(1)模数互质
//设M = PI(i:1->n){ mod[i] } 即模数之积
//定义M[i] = M/mod[i] inv[i] = mod[i]^-1
//则X=(r[1]·M[1]·inv[1]) mod M 实则M除以mod[i]的平方.
//param1:同余方程数量
//param2:余数
//param3:模数
LL CRT(LL n, LL rem[], LL mod[])
{
LL M = 1, x, y;
for(int i = 0; i < n; i++)
M *= mod[i];
LL res = 0;
for(int i = 0; i < n; i++)
{
LL t = M / mod[i];
exgcd(t, mod[i], x, y);//EXGCD求mod[i]逆元
res = (res + t * rem[i] * x) % M;
}
return (res % M + M) % M;
}
//(2)模数非互质 (合并方程)
//利用EXGCD性质 m1·x + m2·y = r2 - r1 求出特解x‘ y‘
//再求出x最小解 x = x‘ * (c/g) mod (b/d)
//X ≡ (m1 * x + r1) mod (LCM(m1, m2))
//当前余数 模数 下个传入的余数 模数
LL crtCD(LL &rem, LL &mod, LL newRem, LL newMod)
{
LL c = newRem - rem;// c= r2 -r1
LL g = __gcd(mod, newMod);
if(c % g!=0)
return 0;
LL x, y;//特解x‘ y‘
exgcd(mod, newMod, x, y);
LL t = newMod / g;// b/d
LL minx = (c / g * x % t + t) % t;//x = x‘ * (c/g) mod (b/d) 最小x
rem = minx * mod + rem;// newrem = mod * x + rem,
mod = mod / g * newMod; //newmod = lcm(mod, newMod)
return 1;
}
/********************************************************/
//3.逆元
//
//大多线性预处理用于组合数
void init(){
inv[1] = inv[0] = 1;
for(int i = 1; i < N; i++)
inv[i] = (mod - mod / i) * inv[mod % i];
}
//单次求解
//等价方程 a/b % m = a mod (m·b)/b
LL Inv(LL a, LL b)
{
LL d = __gcd(a, b);
if(d != 1)
return -1;
LL x, y;
extend_Euclid(a, b, x, y);
return (x % b + b) % b;
}
//费小定理 略
/********************************************************/
//4.筛
//
//(1)线性筛
//常用于筛出数论函数
LL pri[N];
LL mu[N];
LL phi[N];
int c = 0;
bool vis[N];
void prime()
{
MMF(vis);
MMF(sum);
mu[1] = 1;
for(int i = 2; i < N; i++) {
if(!vis[i])
pri[c++] = i,
mu[i] = -1,
phi[i] = i - 1;
for(int j = 0; j < c && i * pri[j] < N; j++) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0){//倍数i=kp
phi[i* pri[j]] = 1;
mu[i * pri[j]] = 0;
break;
}
else
phi[i* pri[j]] = phi[i] * (pri[j] - 1);
mu[i * pri[j]] = -mu[i];
}
}
}
//(2)状压素数筛 不知道拿来干嘛用。
int pri[5800000];
unsigned int sum[5800000];
int vis[N/32+10];
int c = 0;
void prime()
{
pri[0] = 2;
sum[0] = 2;
c = 1;
for(int i = 3; i < 100000007; i+=2)
{
if(!(vis[i/32] & (1 << ((i/2)%16))) )
{
pri[c++] = i;
for(int j = 3*i; j < 100000007; j+=2*i)
vis[j/32] |= (1 << ((j/2)%16));
}
}
}
/********************************************************/
//5.大素数计数
//
//用于论文题
//(1) Meisell-Lehmer
const int N = 5e6 + 2;
bool np[N];
int prime[N], pi[N];
int getprime()
{
int cnt = 0;
np[0] = np[1] = true;
pi[0] = pi[1] = 0;
for(int i = 2; i < N; ++i)
{
if(!np[i]) prime[++cnt] = i;
pi[i] = cnt;
for(int j = 1; j <= cnt && i * prime[j] < N; ++j)
{
np[i * prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
return cnt;
}
const int M = 7;
const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;
int phi[PM + 1][M + 1], sz[M + 1];
void init()
{
getprime();
sz[0] = 1;
for(int i = 0; i <= PM; ++i) phi[i][0] = i;
for(int i = 1; i <= M; ++i)
{
sz[i] = prime[i] * sz[i - 1];
for(int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
}
}
int sqrt2(LL x)
{
LL r = (LL)sqrt(x - 0.1);
while(r * r <= x) ++r;
return int(r - 1);
}
int sqrt3(LL x)
{
LL r = (LL)cbrt(x - 0.1);
while(r * r * r <= x) ++r;
return int(r - 1);
}
LL getphi(LL x, int s)
{
if(s == 0) return x;
if(s <= M) return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
if(x <= prime[s]*prime[s]) return pi[x] - s + 1;
if(x <= prime[s]*prime[s]*prime[s] && x < N)
{
int s2x = pi[sqrt2(x)];
LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
for(int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];
return ans;
}
return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
}
LL getpi(LL x)
{
if(x < N) return pi[x];
LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x / prime[i]) - i + 1;
return ans;
}
LL lehmer_pi(LL x)//答案
{
if(x < N) return pi[x];
int a = (int)lehmer_pi(sqrt2(sqrt2(x)));int bit[64];
int a[N], cnt;
void gauss()
{
cnt = 1;
for(int k = 30; k >= 0; k--)
{
for(int i = cnt; i <= n; i++)
{
if(a[i] & (1 << k))
{
for(int j = 1; j <= n; j++)
{
if(j != i && (a[j] & (1 << k)))
a[j] ^= a[i];
}
swap(a[i], a[cnt]);//因为存的是从高到低 放到前面
cnt++;
break;
}
}
}
cnt--;
/*for(int i = 1; i <= n; i++)
printf("%d", a[i]);
cout << endl;
cout << "~" << cnt << endl;*/
}
int b = (int)lehmer_pi(sqrt2(x));
int c = (int)lehmer_pi(sqrt3(x));
LL sum = getphi(x, a) +(LL)(b + a - 2) * (b - a + 1) / 2;
for (int i = a + 1; i <= b; i++)
{
LL w = x / prime[i];
sum -= lehmer_pi(w);
if (i > c) continue;
LL lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1);
}
return sum;
}
//O 3/4 略
/********************************************************/
//6.SG
//
//用于组合游戏SG局面值求解
int sg[N];
int dfs(int x)
{
if(x < 0)
return 0;
if(sg[x] != -1)
return sg[x];
int vis[2010];//注意函数内声明
MMF(vis);//注意清空
for(int i = 1; i <= x; i++)
vis[dfs(i - 3) ^ dfs(x - i - 2)] = 1;
for(int i = 0; ;i++)//后继局面最小未出现
if(!vis[i])
return sg[x] = i;
return sg[x];
}
/********************************************************/
//7.Gauss
//
//线代用于求多元方程组解
//求线性基等
int bit[64];
int a[N], cnt;
void gauss()
{
cnt = 1;
for(int k = 30; k >= 0; k--)
{
for(int i = cnt; i <= n; i++)
{
if(a[i] & (1 << k))//消元部分
{
for(int j = 1; j <= n; j++)
{
if(j != i && (a[j] & (1 << k)))
a[j] ^= a[i];
}
swap(a[i], a[cnt]);//因为存的是从高到低 放到前面
cnt++;
break;
}
}
}
cnt--;
/*for(int i = 1; i <= n; i++)
printf("%d", a[i]);
cout << endl;*/
}
/********************************************************/
//8.快速傅里叶XX变化
//
//高精度 快速符号卷积运算
//(1) FFT
const double eps = 1e-8;
const double pi = acos(-1.0);
struct Complex
{
double a, b;
Complex(){}
Complex(double _a, double _b):a(_a),b(_b){}
Complex(double _a):a(_a),b(0.0){}
inline Complex operator +(const Complex &z)const{
return Complex(a + z.a, b + z.b);
}
inline Complex operator -(const Complex &z)const{
return Complex(a - z.a, b - z.b);
}
inline Complex operator *(const Complex &z)const{
return Complex(a * z.a - b * z.b, a * z.b + b * z.a);
}
inline Complex operator /(const Complex &z)const{
double m = z.a * z.a + z.b * z.b;
return Complex((a * z.a + b * z.b) / m, (z.a * b - z.b * a) / m);
}
}A[N], B[N];
int L;
int rev[N];
int ans[N];
char a[N], b[N];
void init()
{
MMF(A);
MMF(B);
MMF(ans);
L = 0;//?
}
int reverse(int x,int r) //蝴蝶操作
{
int ans=0;
for(int i=0; i<r; i++){
if(x&(1<<i)){
ans+=1<<(r-i-1);
}
}
return ans;
}
void rader(LL f[], int len)//雷德二进制反转
{
for(int i = 1, j = len >> 1; i < len - 1; i++)
{
if(i < j) swap(f[i], f[j]);
int k = len >> 1;
while(j >= k)
{
j -= k;
k >>= 1;
}
if(j < k) j += k;
}
}
void FFT(Complex c[], int nlen, int on) //主要部分 on为正逆开关
{
Complex wn, w, x, y;
for(int i = 0; i < nlen; i++)
if(i < rev[i]) swap(c[i], c[rev[i]]);
for(int i = 1; i < nlen; i <<= 1)
{
wn = Complex(cos(pi/i), sin(pi/i) * on);
for(int p = i << 1, j = 0; j < nlen; j+= p)
{
w = Complex(1, 0);
for(int k = 0; k < i; k++, w = w * wn)
{
x = c[j + k];
y = w * c[j + k + i];
c[j + k] = x + y;
c[j + k + i] = x - y;
}
}
}
if(on == -1)
for(int i = 0; i < nlen; i++)
c[i].a /= (double)nlen;
}
void work(Complex a[], Complex b[], int len)
{
FFT(a, len, 1); //
FFT(b, len, 1); //
for(int i = 0; i < len; i++)
a[i] = a[i] * b[i];
FFT(a, len, -1);
}
int main()
{
while(~scanf("%s%s", a, b))
{
init();
int alen = strlen(a);
int blen = strlen(b);
int len = alen + blen;
//getlen
int n;
for(n = 1; n < len - 1; n <<= 1, L++);
//rev
for(int i = 0; i < n; i++)
rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (L - 1));
//deal
for(int i = 0; i < alen; i++)
A[i] = a[alen - i - 1] - ‘0‘;
for(int i = 0; i < blen; i++)
B[i] = b[blen - i - 1] - ‘0‘;
work(A, B, n);
///
for(int i = 0; i <= len; i++)
ans[i] = (int)(A[i].a + 0.5);
for(int i = 0; i <= len; i++)
ans[i + 1] += ans[i] / 10, ans[i] %= 10;
while(ans[len] == 0 && len)
len--;
for(int i = len; ~i; i--)
printf("%c", ans[i] + ‘0‘);
printf("\n");
}
return 0;
}
//(2)NTT
//适用于有取模运算 模数为费马素数 p = c·2^k + 1
//HDU 6061
const LL N = 1e5 + 10;
const LL mod = 998244353;
const int g = 3;//原根注意
const int maxlen = 1 << 18;
LL wn[maxlen], fac[N], inv[N];
LL fpow(LL a, LL n)
{
LL ans = 1;
while(n)
{
if(n & 1)
ans = (ans * a % mod + mod) %mod;
a = (a * a + mod) % mod;
n >>= 1;
}
return ans;
}
void init()
{
wn[0] = 1;
wn[1] = fpow(g, ((mod - 1)>>18));
for(int i = 2; i < maxlen; i++)
wn[i] = wn[i - 1] * wn[1] % mod;
fac[1] = inv[1] = 1;
fac[0] = inv[0] = 1;
for(int i = 2; i < N; i++)
{
fac[i] = fac[i - 1] * i % mod;
inv[i] = (mod - mod / i) * inv[mod % i] % mod;
}
for(int i = 1; i < N; i++)
(inv[i] *= inv[i - 1]) %= mod;
}
void rader(LL f[], int len)
{
for(int i = 1, j = len >> 1; i < len - 1; i++)
{
if(i < j) swap(f[i], f[j]);
int k = len >> 1;
while(j >= k)
{
j -= k;
k >>= 1;
}
if(j < k) j += k;
}
}
void ntt(LL f[], int len, int on)
{
/*for(int i = 0, j = 0; i < len;i++)
{
if(i > j) swap(f[i], f[j]);
for(int l = len >> 1; (j^=l) < l; l>>=1);
}*/
rader(f, len);
for(int i = 1, d = 1; d < len; i++, d <<= 1)
{
//LL wnn = fpow(g, (mod-1)/(d<<1));
for(int j = 0; j < len; j += (d << 1))
{
//LL w = 1;
for(int k = 0; k < d; k++)
{
LL t = wn[(maxlen >> i) * k] * f[j + k + d] % mod;
//LL t = w*f[j+k+d]%mod;
//w = w*wnn % mod;
f[j + k + d] = ((f[j + k] - t) % mod + mod) % mod;
f[j + k] = ((f[j + k] + t) % mod + mod) % mod;
}
}
}
if(on == -1)
{
reverse(f + 1, f + len);
LL inv2 = fpow(len, mod - 2);
for(int i = 0; i < len; i++)
f[i] = f[i] % mod * inv2 % mod;
}
}
void work(LL a[], LL b[], int len)
{
ntt(a, len, 1);
ntt(b, len, 1);
for(int i = 0; i < len; i++)
a[i] = (a[i] * b[i] % mod + mod) % mod;
ntt(a, len, -1);
}
LL A[maxlen], B[maxlen], Suma, c[N];
int n, m;
int main()
{
init();
while(~scanf("%d", &n))
{
for(int i = 0; i <= n; i++) scanf("%lld", c + i);
scanf("%d", &m);
Suma = 0;
LL t;
for(int i = 0; i < m; i++)
scanf("%lld", &t), Suma -= t;
Suma = (Suma + mod) % mod;
while(Suma < 0)
Suma += mod;
if(Suma == 0)
{
for(int i = 0; i <= n; i++)
printf("%lld ", c[i]);
printf("\n");
continue;
}
//getLength
int len = 1;
while(len <= 2*n) len <<= 1;
//
LL ae = 1;
for(int i = 0; i < len; i++)
{
if(i <= n)
{
B[i] = ae * inv[i] % mod;
A[i] =(fac[n - i] * c[n - i]) % mod;
}
else A[i] = B[i] = 0;
ae = (ae * Suma) % mod;
//cout << A[i] << "~"<< B[i] << endl;
}
work(A, B, len);
for(int i = 0; i <= n; i++)
A[i] = A[i] * inv[n - i] % mod;
for(int i = 0; i <= n; i++)
printf("%lld ", A[n - i]);
printf("\n");
}
}
//(3)FWT 暂略
//用于快速计算位间操作符号的卷积
/********************************************************/
//9.二次剩余解
//
//常用于数论意义下的化简式子
//可用于复杂方程中替代有多重开根的值（一般为平方）
//因此需要求多重剩余的解x
//
//二次复数域方法求解
LL mod;
//二次域乘法????
//O1乘法取模黑科技
LL mul(LL x,LL y)
{
return (x * y-(LL)(x /(long double)mod * y + 1e-3) * mod + mod) % mod;
}
LL fpow(LL a, LL n)
{
LL ans = 1;
while(n)
{
if(n & 1)
//ans = (a * ans + mod)% mod;
ans = mul(ans, a);
n >>= 1;
//a = (a * a % mod + mod) % mod;
a = mul(a, a);
}
return ans;
}
LL QuadraticResidue()
{
LL rem = (-3 % mod + mod) % mod;
if(rem == 0)//特判mod==3
return 0;
if(mod == 2)//特判非奇素数
return 1;
if(fpow(rem, (mod - 1)/2) + 1 == mod)//欧拉判别条件 非剩余
return -1;
LL a = -1, b;
while(1)
{
a = (1 + rand() % (mod - 1));//注意模板是错的 原来a可能取到0
b = (a * a - rem) % mod;
while(b < 0)
b += mod;
if(fpow(a, (mod-1)/2) + 1 == mod)//找一个非剩余的b 推导上一阶剩余
break;
}
a = fpow(b, (mod + 1)/2);
return a;
}
LL QuadraticResidue2()//求解二次剩余根 x
{
LL rem = ((-3) % mod + mod) % mod;
if(rem == 0)//特判mod==3
return 0;
if(mod == 2)//特判非奇素数
return rem % mod;
if(fpow(rem, (mod - 1)/2) != 1)//欧拉判别条件 非剩余
return -1;
LL x = -1, b;
LL i, k;
if(mod % 4 == 3)
x = fpow(rem, (mod + 1)/4);
else
{
for(b = 1; fpow(b,(mod - 1)/2) == 1;
b = 1 + rand() % (mod - 1));
i = (mod - 1) / 2;
k = 0;
do
{
i /= 2,k /= 2;
if(( mul(fpow(rem, i),fpow(b, k)) + 1)% mod==0)
k += (mod - 1)/2;
}
while(i % 2 == 0);
x = mul(fpow(rem,(i+1)/2),fpow(b,k/2));
}
if(x * 2 > mod)
x = mod - x;
return x;
}
vectora;
int main()
{
int T;
cin >> T;
while(T--)
{
a.clear();
LL n, p;
scanf("%lld%lld", &n, &mod);
for(int i = 0; i < n; i++)
{
LL t;
scanf("%lld", &t);
if(t > 0)//注意有0的...
a.PB(t);
}
LL cnt = a.size();
sort(a.begin(), a.end());
///////////
LL ans = 0;
if(mod == 2)//特殊情况
ans = cnt * (cnt - 1) / 2LL;
else
{
LL t = QuadraticResidue2();
if(t == -1)
{
printf("0\n");
continue;
}
LL inv = (mod + 1) >> 1;
LL x = mul((mod + t - 1LL)%mod, inv);
LL y = mul((mod - t - 1LL)%mod, inv);
for(int i = 0; i < cnt; i++)
{
LL tmp = mul(x , a[i]);
ans += upper_bound(a.begin(), a.begin() + i, tmp)
- lower_bound(a.begin(), a.begin() + i, tmp);
}
if(x != y)//两个解
{
for(int i = 0; i < cnt; i++)
{
LL tmp = mul(y, a[i]);
ans += upper_bound(a.begin(), a.begin() + i, tmp)
- lower_bound(a.begin(), a.begin() + i, tmp);
}
}
}
printf("%lld\n", ans);
}
return 0;
}
/********************************************************/
//10.Miller-rabin 素性测试
//
//40个素数足以判别出64位整数
bool Miller_Rabbin(int n,int a)//米勒拉宾素数测试
{
int r = 0, s = n - 1;
if(!(n % a))
return false;
while(!(s & 1))
{
s>>=1;
r++;
}
LL k = fpow(a, s, n);
if(k == 1)
return true;
for(int j = 0; j < r; j++, k = k * k % n)
if(k == n-1)
return true;
return false;
}
bool IsPrime(int n)//判断是否是素数
{
int tab[14] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};
for(int i = 0; i < 13; i++)
{
if(n == tab[i])
return true;
if(!Miller_Rabbin(n, tab[i]))
return false;
}
return true;
}
/********************************************************/
//11.Lucas
//
//
LL C(LL n, LL m, LL mod)
{
if(m > n)
return 0;
LL ans = 0;
ans = ((fac[n] * inv[m] % mod)* inv[n - m]) % mod;
return ans;
}
LL lucas(LL n, LL m, LL mod)
{
if(m == 0)
return 1;
return C(n % mod, m % mod, mod) * lucas(n / mod, m / mod, mod) % mod;
}
标签：线性筛 inline 同余方程 resid -- eps break euc 论文
原文地址：http://www.cnblogs.com/Yumesenya/p/7747737.html
