[UOJ 48][UR #3]核聚变反应堆
先去考虑两个变量的sgcd怎么求……我们先求出他们的gcd，然后在除掉一个最小的质因子即可。
那么我们先把\(a_1\)分解质因数，然后所有数和它求一个gcd，然后去去掉一个尽可能小的质因子。注意到一个数\(p\)的质因子数量是\(O(\log p)\)级别的，所以每一个数找到要除的那个最小的质因子只需要\(O(\log a_1)\)的时间。
代码：
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <vector>
#include <cmath>
const int maxn = 100005;
using ll = long long;
ll V[maxn]; int vcnt = 0;
void des(ll x) {
ll m = (ll)sqrt(x + 0.5);
for(ll i = 2; i <= m; i ++) {
if(x % i == 0) {
V[vcnt ++] = i;
while(x % i == 0) x /= i;
}
}
if(x > 1LL) V[vcnt ++] = x;
}
ll gcd(ll a, ll b) {
if(!b) return a;
else return gcd(b, a % b);
}
ll A[maxn];
int main() {
int n; scanf("%d", &n);
for(int i = 1; i <= n; i ++) {
scanf("%lld", &A[i]);
}
des(A[1]);
for(int i = 1; i <= n; i ++) {
if(i > 1) putchar(' ');
ll bs = gcd(A[i], A[1]);
if(bs == 1) {
printf("-1"); continue;
}
for(int j = 0; j < vcnt; j ++) {
if(A[i] % V[j] == 0) {
bs /= V[j];
break;
}
}
printf("%lld", bs);
}
puts("");
return 0;
}
[BZOJ 4816][SDOI2017]数字表格
当年在考场上一点反演都不会啊啊啊啊啊啊
这题虽然牵扯到了对积的反演，但其实也不是很难。先让我们大力推导一波（逃
考虑枚举柿子中的\(f[(i, j)]\)中的最大公约数（设为\(k\)），然后式子变成了（这里不妨偷个懒，钦定\(n\leq m\)）：
\[\prod_{k = 1}^n f[k]^{\sum_{i = 1}^n \sum_{j = 1}^m\:[(i, j) = k]}\]
然后发现上面那个指数就是Problem b的式子，显然可化为\(\sum_{d = 1}^n \mu(d) \lfloor\frac{n}{kd}\rfloor \lfloor\frac{m}{kd}\rfloor\)。然后似乎化简没有余地了，但是根据直觉，我们可以钦定\(T = kd\)，然后枚举\(T\)。
然后柿子变成了：
\[\prod_{T = 1}^n (\prod_{k\mid T} f[k]^{\mu(\tfrac{T}{k})})^{\lfloor\tfrac{n}{T}\rfloor \lfloor\tfrac{m}{T}\rfloor}\]
柿子中的\(\prod_{k\mid T} f[k]^{\mu(\tfrac{T}{k})}\)是一个类似狄利克雷卷积的东西（取完对数就是了），可以枚举约数预处理。并且这个东西很不错的一点就是上面的指数是莫比乌斯函数，可以取到的值只有那三种，不需要快速幂。
然后剩下的部分就和通常的反演题一般无二了……
代码（卡线过的蒟蒻……求轻喷……但我在LOJ上跑的还挺快的）：
/**************************************************************
Problem: 4816
User: danihao123
Language: C++
Result: Accepted
Time:50840 ms
Memory:33048 kb
****************************************************************/
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <cctype>
#include <cassert>
const int N = 1000000;
const int maxn = N + 5;
typedef long long ll;
const ll ha = 1000000007LL;
bool p_flag = false;
ll pow_mod(ll a, ll b) {
if(!b) return 1LL;
ll ans = pow_mod(a, b / 2LL);
ans = (ans * ans) % ha;
if(1LL & b) ans = (ans * a) % ha;
#ifdef LOCAL
if(p_flag) printf("%lld^%lld : %lld\n", a, b, ans);
if(b == ha - 2LL) p_flag = false;
#endif
return ans;
}
ll inv(ll x) {
// if(x == 1LL) p_flag = true;
return pow_mod(x, ha - 2LL);
}
int miu[maxn];
void sieve() {
static int prm[maxn]; int cnt = 0;
static bool vis[maxn];
miu[1] = 1;
for(int i = 2; i <= N; i ++) {
if(!vis[i]) {
miu[i] = -1;
prm[cnt ++] = i;
}
for(int j = 0; j < cnt && prm[j] * i <= N; j ++) {
int v = prm[j] * i;
vis[v] = true;
if(i % prm[j] == 0) {
miu[v] = 0;
break;
} else {
miu[v] = miu[i] * -1;
}
}
}
}
ll f[maxn], inv_d[maxn], F[maxn];
void process() {
sieve();
f[1] = 1LL; inv_d[1] = 1LL;
for(int i = 2; i <= N; i ++) {
f[i] = (f[i - 1] + f[i - 2]) % ha;
inv_d[i] = inv(f[i]);
assert(f[i] != 0LL); assert(inv_d[i] != 0LL);
}
for(int i = 1; i <= N; i ++) F[i] = 1LL;
for(int i = 1; i <= N; i ++) {
for(int j = i; j <= N; j += i) {
int res = j / i;
if(miu[res] == 1) {
F[j] = (F[j] * f[i]) % ha;
#ifdef LOCAL
if(j <= 3) printf("f[%d](%lld) -> F[%d]\n", i, f[i], j);
#endif
continue;
}
if(miu[res] == -1) {
F[j] = (F[j] * inv_d[i]) % ha;
#ifdef LOCAL
if(j <= 3) printf("inv_f[%d](%lld) -> F[%d]\n", i, inv_d[i], j);
#endif
}
}
}
F[0] = 1LL;
bool flag = true;
for(int i = 1; i <= N; i ++) {
F[i] = (F[i] * F[i - 1]) % ha;
if(flag && F[i] == 0LL) {
printf("SF[%d] has been 0.\n", i);
flag = false;
}
}
}
ll calc(ll n, ll m) {
if(n > m) std::swap(n, m);
ll ans = 1LL;
for(ll i = 1; i <= n;) {
ll nx = std::min(n / (n / i), m / (m / i));
#ifdef LOCAL
// printf("nx of %lld : %lld\n", i, nx);
#endif
ll mulv = pow_mod((F[nx] * inv(F[i - 1])) % ha, (n / i) * (m / i));
ans = (ans * mulv) % ha;
i = nx + 1LL;
}
return ans;
}
int main() {
process();
int T; scanf("%d", &T);
while(T --) {
int n, m; scanf("%d%d", &n, &m);
printf("%lld\n", calc(n, m));
}
return 0;
}
[BZOJ 3992][SDOI2015]序列统计
终于调过去了……（然而不就是道NTT+生成函数水题吗至于调半天）
首先积非常的恶心，考虑转成和。这个事需要求指标的……求指标的话枚举原根的若干次幂即可（恰好$m$是素数一定有原根……），判断原根用比较大力的手段即可（我搞了一个$O(n\sqrt{n}logn)$的……求轻喷）。
然后这题还算比较简单吧……用生成函数表示原来的集合，然后$n$次幂就可以了。
注意那事个循环卷积……所以要开两倍然后每次乘完了再把右半部分搞过去。
代码：
/**************************************************************
Problem: 3992
User: danihao123
Language: C++
Result: Accepted
Time:6652 ms
Memory:1864 kb
****************************************************************/
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <queue>
#include <utility>
#include <vector>
#include <cmath>
typedef long long ll;
const int maxn = 32005;
const ll ha = 1004535809LL;
const ll bs = 3LL;
ll n, m, gl; int sz;
ll pow_mod(ll a, ll b, ll p) {
if(!b) return 1LL;
ll ans = pow_mod(a, b >> 1, p);
ans = (ans * ans) % p;
if(b & 1LL) ans = (ans * a) % p;
return ans;
}
ll inv(ll x, ll p) {
return pow_mod(x, p - 2LL, p);
}
void break_up(ll x, std::vector<ll> &V) {
int bd = sqrt(x + 0.5);
for(ll i = 2; i <= bd; i ++) {
if(x % i == 0) {
V.push_back(i);
while(x % i == 0) x /= i;
}
if(x == 1LL) break;
}
if(x > 1LL) V.push_back(x);
}
ll get_phi() {
if(m == 2LL) return 1LL;
ll mi = m - 1LL;
std::vector<ll> V;
break_up(mi, V);
for(ll i = 2; i <= mi; i ++) {
bool ok = true;
for(int j = 0; j < V.size(); j ++) {
ll b = mi / V[j];
if(pow_mod(i, b, m) == 1LL) {
ok = false;
#ifdef LOCAL
printf("%lld not passed!\n", i);
#endif
break;
}
}
if(ok) {
return i;
}
}
}
int bi, len;
int flip(int x) {
int ans = 0;
for(int i = 0; i < bi; i ++) {
if((1 << i) & x) {
ans += (1 << (bi - i - 1));
}
}
return ans;
}
void ntt(ll *A, bool flag = false) {
for(int i = 0; i < len; i ++) {
int v = flip(i);
if(v < i) std::swap(A[i], A[v]);
}
for(int L = 1; L < len; L <<= 1) {
ll xi_n = pow_mod(3LL, (ha - 1LL) / (ll(L << 1)), ha);
if(flag) xi_n = inv(xi_n, ha);
for(int i = 0; i < len; i += (L << 1)) {
ll w = 1;
for(int j = i; j < i + L; j ++) {
ll p1 = A[j], p2 = A[j + L];
A[j] = (p1 + (p2 * w) % ha) % ha;
A[j + L] = (p1 - (p2 * w) % ha + ha) % ha;
w = (w * xi_n) % ha;
}
}
}
}
void poly_mul(ll *A, ll *B) {
static ll C[maxn]; memset(C, 0, sizeof(C));
#ifdef LOCAL
printf("A :");
for(int i = 0; i < len; i ++) printf(" %lld", A[i]);
puts("");
printf("B :");
for(int i = 0; i < len; i ++) printf(" %lld", B[i]);
puts("");
#endif
ntt(A); ntt(B);
for(int i = 0; i < len; i ++) C[i] = (A[i] * B[i]) % ha;
ntt(C, true);
ll inv_n = inv(len, ha);
for(int i = 0; i < len; i ++) {
C[i] = (C[i] * inv_n) % ha;
}
#ifdef LOCAL
printf("C (not processed) :");
for(int i = 0; i < len; i ++) printf(" %lld", C[i]);
puts("");
#endif
for(int i = 0; i < len; i ++) {
int v = i % (m - 1LL);
if(v != i) {
C[v] = (C[v] + C[i]) % ha;
C[i] = 0LL;
}
}
#ifdef LOCAL
printf("C :");
for(int i = 0; i < len; i ++) printf(" %lld", C[i]);
puts("");
#endif
std::copy(C, C + len, A);
}
void poly_pow(ll *A, ll b) {
static ll B[maxn];
static ll res[maxn];
std::copy(A, A + len, res);
std::fill(A, A + len, 0); A[0] = 1LL;
#ifdef LOCAL
printf("A :");
for(int i = 0; i < len; i ++) printf(" %lld", A[i]);
puts("");
printf("res : ");
for(int i = 0; i < len; i ++) printf("%lld ", res[i]);
puts("");
#endif
while(b) {
if(1LL & b) {
std::copy(res, res + len, B);
poly_mul(A, B);
std::fill(B, B + len, 0);
}
std::copy(res, res + len, B);
poly_mul(res, B);
std::fill(B, B + len, 0);
b >>= 1LL;
}
}
int main() {
static bool vis[maxn];
static ll A[maxn];
scanf("%lld%lld%lld%d", &n, &m, &gl, &sz);
ll phi = get_phi();
for(int i = 1; i <= sz; i ++) {
int v; scanf("%d", &v);
if(v == 0) continue;
vis[v] = true;
}
int logx = 0;
bi = 0; len = 1;
while(len <= (2 * m - 2)) {
bi ++; len <<= 1;
}
for(int i = 0; i < (m - 1); i ++) {
int v = pow_mod(phi, i, m);
if(vis[v]) {
A[i] ++;
#ifdef LOCAL
printf("log(%d) : %d\n", v, i);
#endif
}
if(v == gl) {
logx = i;
}
}
poly_pow(A, n);
printf("%lld\n", A[logx]);
return 0;
}
[坑]狄利克雷卷积和反演
开个坑，记录一些和反演以及狄利克雷卷积的东西。
首先积性函数、狄利克雷卷积等基本概念不写了，就讲讲性质吧。
有几个一定要记住的东西：
\[\mu \ast 1 = e\]
\[\varphi \ast 1 = id\]
\[\mu \ast id = \varphi\]
这几个在推式子的过程中都有很大的作用，务必要记住。
所谓莫比乌斯反演，其实就是：
\[F = f \ast 1\Leftrightarrow f = F \ast \mu\]
（谜之音：其实很多所谓“反演题”都没用到这俩性质啊……）
关于莫比乌斯函数本身，还有一个好康的性质：
\[(\mu\ast 1)(k) = \sum_{i = 0}^k (-1)^i C_k^i\]
[BZOJ 2186]沙拉公主的困惑
这个题啊……亦可赛艇！
答案是[tex]\varphi(m!)*n!/m![/tex]。原因很简单，把[tex][1,n!][/tex]分成长度为[tex]m![/tex]的若干段，除去第一段外每一段中与[tex]m![/tex]互质的数[tex]k[/tex]肯定满足[tex](k\bmod m!,m!)=1[/tex]（否则，[tex]k[/tex]和[tex]m![/tex]就会有大于一的公因子了)。所以说每一段内与[tex]m![/tex]互质的数都有[tex]\varphi(m!)[/tex]个。
麻烦好像就在于求一个阶乘的欧拉函数。考虑一个新乘上的数能给答案带来的贡献——如果这个数是合数，它的所有因子在前面都有了，只能把他自己贡献出来；如果这个数是质数（假设为[tex]p[/tex]），出了贡献自己以外还会贡献一个[tex](1-1/p)[/tex]，最后乘起来就是贡献了[tex]p-1[/tex]。筛一遍素数再递推一下就好辣~
并且……[tex]n-m[/tex]可能非常大，所以说除去[tex]m![/tex]那块要用逆元做。
（顺便说下阶乘也要递推）
代码：
/**************************************************************
Problem: 2186
User: danihao123
Language: C++
Result: Accepted
Time:9408 ms
Memory:166836 kb
****************************************************************/
#include <cstdio>
#include <cmath>
typedef unsigned long long ull;
const int maxn=10000000;
ull R;
bool vis[maxn+5];
inline void sievePrime(){
register int i,j,m=sqrt(maxn+0.5);
for(i=2;i<=m;i++)
if(!vis[i])
for(j=i*i;j<=maxn;j+=i)
vis[j]=true;
}
ull fac[maxn+5];
inline void sieveFac(){
register int i;
fac[0]=1%R;
for(i=1;i<=maxn;i++)
fac[i]=(fac[i-1]*(i%R))%R;
}
ull phifac[maxn+5];
inline void sievePhiFac(){
register int i;
phifac[1]=1%R;
for(i=2;i<=maxn;i++){
if(vis[i])
phifac[i]=(phifac[i-1]*(i%R))%R;
else
phifac[i]=(phifac[i-1]*((i%R-1%R+R)%R))%R;
}
}
void exgcd(ull a,ull b,ull& d,ull& x,ull& y){
if(!b){
d=a;
x=1;
y=0;
}else{
exgcd(b,a%b,d,y,x);
y-=x*(a/b);
}
}
ull inv(ull a){
ull d,x,y;
exgcd(a,R,d,x,y);
return (x+R)%R;
}
int main(){
int T;
int n,m;
scanf("%d%llu",&T,&R);
sievePrime();
sieveFac();
sievePhiFac();
while(T--){
scanf("%d%d",&n,&m);
printf("%llu\n",(phifac[m]*((fac[n]*inv(fac[m]))%R))%R);
}
return 0;
}
[BZOJ 2118]墨墨的等式
论如何把数论乱搞和图论乱搞出在一起……
这个题由于要求[tex]x\ge 0[/tex]，所以不能gcd乱搞。我们可以先取[tex]\{a_n\}[/tex]的最小值[tex]p[/tex]（忽略为0的情况，为啥取最小值待会再说），对方程两边模[tex]p[/tex]。然后对于任何能使某个同余方程成立的[tex]\{x_n\}[/tex]，将其中所有[tex]x_i[/tex]同时加任意个[tex]p[/tex]，同余方程都成立。
取模后，[tex]B\in Z_p[/tex]，所以说只要对于[tex]Z_p[/tex]中的每个数找出最小的一组合法解即能推出其他解（所以说，剩余系越少效率越高，这也就要求取的[tex]a_i[/tex]要小）。不过这个最小的一组合法解怎么找？
我们先找出任意一个合法[tex]B[/tex]（比如说0吧），然后尝试加上[tex]a_i[/tex]，就可以推出其他[tex]B\in Z_p[/tex]的最小解。这个应用当然是需要最短路辣。
求出来的最短路，表示的是取最小解时的[tex]B[/tex]。这样的话就可以推出某个前缀区间中合法[tex]B[/tex]的个数（能加多少[tex]p[/tex]，就有多少个，注意不要忽略加0个的情况），并且答案符合区间减法，最后做差分即可。
注意忽略[tex]a_i=0[/tex]的情况（相当于不存在）。
代码：
/**************************************************************
Problem: 2118
User: danihao123
Language: C++
Result: Accepted
Time:1952 ms
Memory:5508 kb
****************************************************************/
#include <cstdio>
#include <cstring>
#include <queue>
#include <algorithm>
#ifdef DEBUG
#include <cassert>
#endif
using namespace std;
typedef long long ll;
const int maxn=15;
const ll INF=0x7f7f7f7f7f7f7f7f;
ll A[maxn];
bool inQueue[500005];
ll d[500005];
int n;
ll minv;
inline void SPFA(){
register int i,u,to;
queue<int> Q;
memset(d,0x7f,sizeof(d));
d[0]=0;
Q.push(0);
inQueue[0]=true;
#ifdef DEBUG
assert(d[1]==INF);
#endif
while(!Q.empty()){
u=Q.front();
Q.pop();
inQueue[u]=false;
for(i=1;i<=n;i++){
to=(u+A[i])%minv;
if(d[u]<INF && d[u]+A[i]<d[to]){
d[to]=d[u]+A[i];
if(!inQueue[to]){
Q.push(to);
inQueue[to]=true;
}
}
}
}
}
inline ll Calc(ll x){
register ll ans=0;
register int i=0;
for(i=0;i<minv;i++)
if(d[i]<=x)
ans+=(x-d[i])/minv+1;
return ans;
}
int main(){
ll l,r;
register int i;
scanf("%d%lld%lld",&n,&l,&r);
minv=0x7fffffff;
for(i=1;i<=n;i++){
scanf("%d",&A[i]);
if(!A[i]){
i--;
n--;
continue;
}
minv=min(minv,A[i]);
}
SPFA();
printf("%lld\n",Calc(r)-Calc(l-1));
return 0;
}
[BZOJ 2190]仪仗队
这个题啊……有规律可循。
我们可以发现，关于答案[tex]a[/tex]有如下规律：
[tex]a+1=2\mul (\Sigma_{i=2}^{n-1}\varphi(i)+2)[/tex]
然后筛法求欧拉函数即可（我听说神犇们都用了杜教筛哎）
代码：
/**************************************************************
Problem: 2190
User: danihao123
Language: C++
Result: Accepted
Time:28 ms
Memory:976 kb
****************************************************************/
#include <cstdio>
#include <algorithm>
using namespace std;
const int maxn=40005;
#define CUS_REP(i,a,b) for(i=(a);i<=(b);i++)
int phi[maxn];
int n;
void sieve(){
register int i,j;
phi[1]=1;
CUS_REP(i,2,n)
if(!phi[i])
for(j=i;j<=n;j+=i){
if(!phi[j])
phi[j]=j;
phi[j]=phi[j]/i*(i-1);
}
}
int main(){
register int i,ans=2;
scanf("%d",&n);
sieve();
/*
#ifdef DEBUG
CUS_REP(i,1,n)
printf("%d\n",phi[i]);
#endif
*/
CUS_REP(i,2,n-1)
ans+=phi[i];
printf("%d\n",ans*2-1);
return 0;
}
[BZOJ 1441]Min
这个问题看似无从下手。
我们可以先取[tex]n=2[/tex]，然后你就发现你似乎要解[tex]A_{1}X_{1}+A_{2}X_{2}>0[/tex]，而且还要求[tex]S[/tex]最小？你想到了什么？扩展欧几里得？对头！
由扩欧推广可得答案就是所有数的最大公约数。
代码：
/**************************************************************
Problem: 1441
User: danihao123
Language: C++
Result: Accepted
Time:0 ms
Memory:820 kb
****************************************************************/
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
int gcd(int a,int b){
return b==0?a:gcd(b,a%b);
}
int main(){
register int ans,i;
int n,temp;
scanf("%d",&n);
scanf("%d",&temp);
ans=abs(temp);
for(i=2;i<=n;i++){
scanf("%d",&temp);
ans=gcd(ans,abs(temp));
}
printf("%d\n",ans);
return 0;
}
[CF 371B]Fox Dividing Cheese
狐狸的三种操作的实质是除二，除三，除五。由此我们可以猜想在最优策略下，两块蛋糕最后的重量应该是[tex]gcd(a,b)[/tex]。然后试除即可。
（大胆猜想，不用证明）
代码：
#include <iostream>
#include <algorithm>
using namespace std;
int gcd(int a,int b){
if(!b)
return a;
else
return gcd(b,a%b);
}
static int P[3]={2,3,5};
inline int min_fj(int source,int direction){
register int i,ans=0;
source/=direction;
if(source==1)
return 0;
for(i=2;i>=0;i--){
while(source%P[i]==0){
source/=P[i];
ans++;
}
}
if(source==1)
return ans;
else
return -1;
}
int main(){
register int direction,t1,t2;
int a,b;
cin>>a>>b;
if(a==b){
cout<<0<<endl;
return 0;
}
direction=gcd(a,b);
t1=min_fj(a,direction);
if(t1==-1){
cout<<-1<<endl;
return 0;
}
t2=min_fj(b,direction);
if(t2==-1){
cout<<-1<<endl;
return 0;
}
cout<<t1+t2<<endl;
return 0;
}
[CodeVS 1012]最大公约数和最小公倍数问题
很经典的问题了吧……然而现在才A……
应注意[tex]P*Q=x*y[/tex]，然而[tex]P[/tex]和[tex]Q[/tex]都可表示为[tex]x[/tex]的乘积，问题就好思考多了。答案就是[tex]y/x[/tex]的质因子的子集数（同一种质因子不能同时分配到P与Q中，否则gcd(P,Q)会不等于x）。
注意有可能无解！
代码：
#include <iostream>
using namespace std;
inline int fj(int x){
register int i,ans=0;
for(i=2;i<=x && x>1;i++)
if(!(x%i)){
ans++;
while(!(x%i))
x/=i;
}
return ans;
}
int main(){
int a,b;
register int ans;
cin>>a>>b;
if(b%a)
ans=0;
else
ans=a==1?1:1<<fj(b/a);
cout<<ans<<endl;
return 0;
}
