模板部分:
#define LL long long
#define FF(i,n) for(i=0;i<n;i++)
LL f[N];
LL ans[N];
LL init[N];
LL buf[N];
void matrixMul(LL a[N],LL b[N],LL n,LL mod)
{
int i,j;
FF(i, n) buf[i]=0;
FF(i, n) FF(j,n) if(a[(i - j + n) % n] && b[j])
buf[i] = (buf[i] + a[(i - j + n) % n] * b[j]);
FF(i, n) a[i] = buf[i] % mod; // 算完后再取模,节省时间
}
void matrixMul(int n,int m,int mod)
{
int i,j;
FF(i,n) ans[i] = (i==0);
while(m > 1)
{
if(m&1) matrixMul(ans, init, n, mod);
matrixMul(init, init, n, mod);
m >>= 1;
}
matrixMul(init, ans, n, mod);
}
UVALive 4727 Jump 约瑟夫变形问题
模板部分:
LL gcd(LL a, LL b)
{
return b == 0 ? a : gcd(b, a % b);
}
void E_gcd(LL a, LL b, LL &d, LL& x, LL& y)
{
if(!b)
{
d = a;
x = 1, y = 0;
}
else
{
E_gcd(b, a % b, d, y, x);
y -= x*(a / b);
}
}
hdu 4497 GCD and LCM 质分解
模板部分:
const int N = 111111;
vector<int> hav;
bool isp[N];
int p[N], cnt;
void getp()
{
int i, j;cnt = 0;
isp[0] = isp[1] = 1;
for(i = 2; i < N; i ++)
{
if(!isp[i])
{
p[cnt ++] = i;
if(i <= 1111)for(j = i * i; j < N; j += i) isp[j] = 1;
}
}
}
void get_hav(int h) //求h的所有质因数(算重复)
{
int i;
for(i = 0; i < cnt && h > 1; i ++)
{
if(h % p[i] == 0)
{
h /= p[i];
hav.push_back(p[i]);
i --;
}
}
if(h != 1) hav.push_back(h);
}
hdu 4418 Time travel 高斯消元
const int N = 111111;
vector<int> hav;
bool isp[N];
int p[N], cnt;
void getp()
{
int i, j;cnt = 0;
isp[0] = isp[1] = 1;
for(i = 2; i < N; i ++)
{
if(!isp[i])
{
p[cnt ++] = i;
if(i <= 1111)for(j = i * i; j < N; j += i) isp[j] = 1;
}
}
}
void get_hav(int h) //求h的所有质因数(算重复)
{
int i;
for(i = 0; i < cnt && h > 1; i ++)
{
if(h % p[i] == 0)
{
h /= p[i];
hav.push_back(p[i]);
i --;
}
}
if(h != 1) hav.push_back(h);
}
模板部分:
int gauss(int N, int M){
int i, j, r, c, pvt;
double maxp;
for (r = 0, c = 0; r < N && c < M; ++ r, ++ c) {
for (maxp = 0, i = r; i < N; ++ i)
if (fabs(a[i][c])>fabs(maxp)) maxp = a[pvt=i][c];
if (sgn(maxp) == 0) {
r--;
continue;
}
if (pvt != r)
for (j = r; j <= M; ++j) swap(a[r][j], a[pvt][j]);
for(int i = r + 1; i < N; i ++)
{
double tmp = a[i][c] / a[r][c];
for(int j = c; j <= M; j ++)
{
a[i][j] -= tmp * a[r][j];
}
}
}
for (i = r; i < N; ++i)
if (sgn(a[i][M])) return -1;
if (r < M) return M-r;
for (i = M-1; i >= 0; --i){
for (j = i+1; j < M; ++j)
a[i][M] -= a[j][M]*a[i][j]; a[i][M] /= a[i][i];
}
return 0;
}
ZOJ 2320 Cracking' RSA 其次布尔线性方程组,高斯消元。
int gauss(int N, int M){
int i, j, r, c, pvt;
double maxp;
for (r = 0, c = 0; r < N && c < M; ++ r, ++ c) {
for (maxp = 0, i = r; i < N; ++ i)
if (fabs(a[i][c])>fabs(maxp)) maxp = a[pvt=i][c];
if (sgn(maxp) == 0) {
r--;
continue;
}
if (pvt != r)
for (j = r; j <= M; ++j) swap(a[r][j], a[pvt][j]);
for(int i = r + 1; i < N; i ++)
{
double tmp = a[i][c] / a[r][c];
for(int j = c; j <= M; j ++)
{
a[i][j] -= tmp * a[r][j];
}
}
}
for (i = r; i < N; ++i)
if (sgn(a[i][M])) return -1;
if (r < M) return M-r;
for (i = M-1; i >= 0; --i){
for (j = i+1; j < M; ++j)
a[i][M] -= a[j][M]*a[i][j]; a[i][M] /= a[i][i];
}
return 0;
}
int gauss(int N, int M)
{
int r, c, pvt;
bool flag;
for (r = 0, c = 0; r < N && c < M; ++ r, ++ c) {
flag = false;
for (int i = r; i < N; ++ i)
if (a[i][c]) {
flag = a[pvt=i][c];
break;
}
if (!flag) {
r--; continue;
}
if (pvt != r)
for (int j = r; j <= M; ++j) swap(a[r][j], a[pvt][j]);
for (int i = r+1; i < N; ++i)
if(a[i][c])
{
a[i][c] = false;
for (int j = c+1; j <= M; ++j) {
if(a[r][j]) a[i][j] = !a[i][j];
}
}
}
return r; //return的是系数矩阵的秩
}
int gauss(int N, int M)
{
int r, c, pvt;
bool flag;
for (r = 0, c = 0; r < N && c < M; ++ r, ++ c) {
flag = false;
for (int i = r; i < N; ++ i)
if (a[i][c]) {
flag = a[pvt=i][c];
break;
}
if (!flag) {
r--; continue;
}
if (pvt != r)
for (int j = r; j <= M; ++j) swap(a[r][j], a[pvt][j]);
for (int i = r+1; i < N; ++i)
if(a[i][c])
{
a[i][c] = false;
for (int j = c+1; j <= M; ++j) {
if(a[r][j]) a[i][j] = !a[i][j];
}
}
}
return r; //return的是系数矩阵的秩
}
hdu 4746 Mophues 莫比乌斯反演
模板部分:
void Mobius()
{
CLR(isp, 0);CLR(cnt, 0);
int tot = 0;mob[1] = 1;
for(int i = 2; i < N; i ++)
{
if(!isp[i])
{
p[tot ++] = i;
mob[i] = -1;cnt[i] = 1;
}
for(int j = 0; j < tot && p[j] * i < N; j ++)
{
isp[p[j] * i] = true;
cnt[i * p[j]] = cnt[i] + 1;
if(i % p[j] == 0)
{
mob[p[j] * i] = 0;
break;
}
else
{
mob[p[j] * i] = -mob[i];
}
}
}
}
分块处理:
void Mobius()
{
CLR(isp, 0);CLR(cnt, 0);
int tot = 0;mob[1] = 1;
for(int i = 2; i < N; i ++)
{
if(!isp[i])
{
p[tot ++] = i;
mob[i] = -1;cnt[i] = 1;
}
for(int j = 0; j < tot && p[j] * i < N; j ++)
{
isp[p[j] * i] = true;
cnt[i * p[j]] = cnt[i] + 1;
if(i % p[j] == 0)
{
mob[p[j] * i] = 0;
break;
}
else
{
mob[p[j] * i] = -mob[i];
}
}
}
}
for(int i = 1, j = 1; i < n; i = j + 1)
{
j = min(n / (n / i), m / (m / i));
ans += (LL)(mbs[j][p] - mbs[i - 1][p]) * (n / i) * (m / i);
}
容斥原理:
for(int i = 1, j = 1; i < n; i = j + 1)
{
j = min(n / (n / i), m / (m / i));
ans += (LL)(mbs[j][p] - mbs[i - 1][p]) * (n / i) * (m / i);
}
void dfs(int st, LL now) ///num表示元素个数。
{
ans.insert(now);
if(st >= tot) return;
LL tmp = 1;
for(int i = 0; i <= num[st]; i ++)
{
dfs(st + 1, now * tmp);
tmp *= bef[st];
}
return ;
}
错排公式。(摘自百科)
void dfs(int st, LL now) ///num表示元素个数。
{
ans.insert(now);
if(st >= tot) return;
LL tmp = 1;
for(int i = 0; i <= num[st]; i ++)
{
dfs(st + 1, now * tmp);
tmp *= bef[st];
}
return ;
}
当n个编号元素放在n个编号位置,元素编号与位置编号各不对应的方法数用M(n)表示,那么M(n-1)就表示n-1个编号元素放在n-1个编号位置,各不对应的方法数,其它类推.
第一步,把第n个元素放在一个位置,比如位置k,一共有n-1种方法;
第二步,放编号为k的元素,这时有两种情况⑴把它放到位置n,那么,对于剩下的n-1个元素,由于第k个元素放到了位置n,剩下n-2个元素就有M(n-2)种方法;⑵第k个元素不把它放到位置n,这时,对于这n-1个元素,有M(n-1)种方法;
综上得到
M(n)=(n-1)[M(n-2)+M(n-1)]
特殊地,M⑴=0,M⑵=1
下面通过这个递推关系推导
通项公式:
为方便起见,设M(k)=k!N(k),(k=1,2,…,n)
则N⑴=0,N⑵=1/2
n>=3时,n!N(n)=(n-1)(n-1)!N(n-1)+(n-1)!N(n-2)
即 nN(n)=(n-1)N(n-1)+N(n-2)
于是有N(n)-N(n-1)=-[N(n-1)-N(n-2)]/n=(-1/n)[-1/(n-1)][-1/(n-2)]…(-1/3)[N⑵-N⑴]=(-1)^n/n!
因此
N(n-1)-N(n-2)=(-1)^(n-1)/(n-1)!
N⑵-N⑴=(-1)^2/2!
相加,可得
N(n)=(-1)^2/2!+…+(-1)^(n-1)/(n-1)!+(-1)^n/n!
因此
M(n)=n![(-1)^2/2!+…+(-1)^(n-1)/(n-1)!+(-1)^n/n!]
可以得到
错排公式为M(n)=n!(1/2!-1/3!+…..+(-1)^n/n!)
梅森素数 from 百科
序号
|
p
|
Mp=2^p-1
|
Mp的位数
|
发现时间
|
发现者
|
1
|
2
|
3
|
1
|
古代
|
古人
|
2
|
3
|
7
|
1
|
古代
|
古人
|
3
|
5
|
31
|
2
|
古代
|
古人
|
4
|
7
|
127
|
3
|
古代
|
古人
|
5
|
13
|
8191
|
4
|
1456年
|
无名氏
|
6
|
17
|
131071
|
6
|
1588年
|
Cataldi
|
7
|
19
|
524287
|
6
|
1588年
|
Cataldi
|
8
|
31
|
2147483647
|
10
|
1772年
| |
9
|
61
|
23058430092
|
19
|
1883年
|
Pervushin
|
10
|
89
|
618970019…449562111
|
27
|
1911年
| |
11
|
107
|
162259276…010288127
|
33
|
1914年
|
Powers
|
12
|
127
|
170141183…884105727
|
39
|
1876年
| |
13
|
521
|
686479766…115057151
|
157
|
1952年1月30日
| |
14
|
607
|
531137992…031728127
|
183
|
1952年1月30日
|
Robinson
|
15
|
1,279
|
104079321…168729087
|
386
|
1952年6月25日
|
Robinson
|
16
|
2,203
|
147597991…697771007
|
664
|
1952年10月7日
|
Robinson
|
17
|
2,281
|
446087557…132836351
|
687
|
1952年10月9日
|
Robinson
|
18
|
3,217
|
259117086…909315071
|
969
|
1957年9月8日
|
Riesel
|
19
|
4,253
|
190797007…350484991
|
1,281
|
1961年11月3日
|
Hurwitz
|
20
|
4,423
|
285542542…608580607
|
1,332
|
1961年11月3日
|
Hurwitz
|
21
|
9,689
|
478220278…225754111
|
2,917
|
1963年5月11日
|
Gillies
|
22
|
9,941
|
346088282…789463551
|
2,993
|
1963年5月16日
|
Gillies
|
23
|
11,213
|
281411201…696392191
|
3,376
|
1963年6月2日
|
Gillies
|
24
|
19,937
|
431542479…968041471
|
6,002
|
1971年3月4日
| |
25
|
21,701
|
448679166…511882751
|
6,533
|
1978年10月30日
|
Noll&
Nickel
|
26
|
23,209
|
402874115…779264511
|
6,987
|
1979年2月9日
|
Noll
|
27
|
44,497
|
854509824…011228671
|
13,395
|
1979年4月8日
|
Nelson& Slowinski
|
28
|
86,243
|
536927995…433438207
|
25,962
|
1982年9月25日
|
Slowinski
|
29
|
110,503
|
521928313…465515007
|
33,265
|
1988年1月28日
|
Colquitt& Welsh
|
30
|
132,049
|
512740276…730061311
|
39,751
|
1983年9月20日
|
Slowinski
|
31
|
216,091
|
746093103…815528447
|
65,050
|
1985年9月6日
|
Slowinski
|
32
|
756,839
|
174135906…544677887
|
227,832
|
1992年2月19日
|
Slowinski&
Gage
|
33
|
859,433
|
129498125…500142591
|
258,716
|
1994年1月10日
|
Slowinski& Gage
|
34
|
1,257,787
|
412245773…089366527
|
378,632
|
1996年9月3日
|
Slowinski& Gage
|
35
|
1,398,269
|
814717564…451315711
|
420,921
|
1996年11月13日
|
GIMPS/Joel Armengaud
|
36
|
2,976,221
|
623340076…729201151
|
895,932
|
1997年8月24日
|
GIMPS/Gordon Spence
|
37
|
3,021,377
|
127411683…024694271
|
909,526
|
1998年1月27日
|
GIMPS/Roland Clarkson
|
38
|
6,972,593
|
437075744…924193791
|
2,098,960
|
1999年6月1日
|
GIMPS/Nayan Hajratwala
|
39
|
13,466,917
|
924947738…256259071
|
4,053,946
|
2001年11月14日
|
GIMPS/Michael Cameron
|
40
|
20,996,011
|
125976895…855682047
|
6,320,430
|
2003年11月17日
|
GIMPS/Michael Shafer
|
41
|
24,036,583
|
299410429…733969407
|
7,235,733
|
2004年5月15日
|
GIMPS/Josh Findley
|
42
|
25,964,951
|
122164630…577077247
|
7,816,230
|
2005年2月18日
|
GIMPS/Martin Nowak
|
43*
|
30,402,457
|
315416475…652943871
|
9,152,052
|
2005年12月15日
|
GIMPS/Curtis Cooper及Steven Boone
|
44*
|
32,582,657
|
124575026…053967871
|
9,808,358
|
2006年9月4日
|
GIMPS/Curtis Cooper及Steven Boone
|
45*
|
37,156,667
|
202254406…308220927
|
11,185,272
|
2008年9月6日
|
GIMPS/Hans-Michael Elvenich
|
46*
|
42,643,801
|
169873516…562314751
|
12,837,064
|
2009年4月12日
|
GIMPS/Odd M. Strindmo
|
47*
|
43,112,609
|
316470269…697152511
|
12,978,189
|
2008年8月23日
|
GIMPS/Edson Smith
|
48* | 57,885,161 | 581887266…724285951 | 17,425,170 | 2013年1月25日 | GIMPS/Curtis Cooper[1] |
大步小步攻击算法(Baby step Giant Step)AekdyCoin大牛的模板,太牛了。。。poj 3243
#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
typedef long long LL;
const int maxn = 65535;
struct hash
{
int a,b,next;
} Hash[maxn << 1];
int flg[maxn + 66];
int top,idx;
void ins(int a,int b)
{
int k = b & maxn;
if(flg[k] != idx)
{
flg[k] = idx;
Hash[k].next = -1;
Hash[k].a = a;
Hash[k].b = b;
return ;
}
while(Hash[k].next != -1)
{
if(Hash[k].b == b) return ;
k = Hash[k].next;
}
Hash[k].next = ++ top;
Hash[top].next = -1;
Hash[top].a = a;
Hash[top].b = b;
}
int find(int b)
{
int k = b & maxn;
if(flg[k] != idx) return -1;
while(k != -1)
{
if(Hash[k].b == b) return Hash[k].a;
k = Hash[k].next;
}
return -1;
}
int gcd(int a,int b)
{
return b?gcd(b,a%b):a;
}
int ext_gcd(int a,int b,int& x,int& y)
{
int t,ret;
if (!b)
{
x=1,y=0;
return a;
}
ret=ext_gcd(b,a%b,x,y);
t=x,x=y,y=t-a/b*y;
return ret;
}
int Inval(int a,int b,int n)
{
int x,y,e;
ext_gcd(a,n,x,y);
e=(LL)x*b%n;
return e<0?e+n:e;
}
int pow_mod(LL a,int b,int c)
{
LL ret=1%c;
a%=c;
while(b)
{
if(b&1)ret=ret*a%c;
a=a*a%c;
b>>=1;
}
return ret;
}
int BabyStep(int A,int B,int C)///A^x % C = B.
{
top = maxn;
++ idx;
LL buf=1%C,D=buf,K;
int i,d=0,tmp;
for(i=0; i<=100; buf=buf*A%C,++i)if(buf==B)return i;
while((tmp=gcd(A,C))!=1)
{
if(B%tmp)return -1;
++d;
C/=tmp;
B/=tmp;
D=D*A/tmp%C;
}
int M=(int)ceil(sqrt((double)C));
for(buf=1%C,i=0; i<=M; buf=buf*A%C,++i)ins(i,buf);
for(i=0,K=pow_mod((LL)A,M,C); i<=M; D=D*K%C,++i)
{
tmp=Inval((int)D,B,C);
int w ;
if(tmp>=0&&(w = find(tmp)) != -1)return i*M+w+d;
}
return -1;
}
int main()
{
int A,B,C;
while(scanf("%d%d%d",&A,&C,&B)!=EOF,A || B || C)
{
B %= C;
int tmp=BabyStep(A,B,C);
if(tmp<0)puts("No Solution");
else printf("%d\n",tmp);
}
return 0;
}
单个欧拉函数求值
const int N = 100100;
bool isp[N];
vector<int> p;
void get_P()
{
CLR(isp, true);p.clear();
for(int i = 2; i < N; i ++)
{
if(isp[i])
{
p.push_back(i);
if(i < 1111) for(int j = i * i; j < N; j += i)
{
isp[j] = false;
}
}
}
}
int Euler_phi(int n)
{
int ret = n;
for(int i = 0; p[i] * p[i] <= n; i ++) if(n % p[i] == 0)
{
ret = ret / p[i] * (p[i] - 1);
while(n % p[i] == 0) n /= p[i];
}
if(n > 1) ret = ret / n * (n - 1);
return ret;
}
筛法求欧拉函数。
const int N = 100100;
int phi[N];
void phi_table()
{
for(int i = 2; i < N; i ++) phi[i] = 0;
phi[1] = 1;
for(int i = 2; i < N; i ++) if(!phi[i])
for(int j = i; j < N; j += i)
{
if(!phi[j]) phi[j] = j;
phi[j] = phi[j] / i * (i - 1);
}
}
void Ex_gcd(LL a, LL b, LL& d, LL& x, LL& y)
{
if(!b) {d = a; x = 1; y = 0;}
else {Ex_gcd(b, a % b, d, y, x); y -= x * (a / b); }
}
int main()
{
int n;
LL a1, r1, a2, r2, ans, a, b, c, d, x, y;
while(scanf("%d", &n) != EOF)
{
bool flag = true;
scanf("%I64d%I64d", &a1, &r1);
for(int i = 1; i < n; i ++)
{
scanf("%I64d%I64d", &a2, &r2);
if(!flag) continue;
a = a1, b = a2, c = r2 - r1;
Ex_gcd(a, b, d, x, y);
if(c % d != 0) flag = false;
int t = b / d;
x = (x * (c / d) % t + t) % t;
r1 = a1 * x + r1;
a1 = a1 * (a2 / d);
}
if(!flag) puts("-1");
else printf("%I64d\n", r1);
}
}
求逆元
void Ex_gcd(LL a, LL b, LL& d, LL& x, LL& y)
{
if(!b) {d = a; x = 1; y = 0;}
else {Ex_gcd(b, a % b, d, y, x); y -= x * (a / b); }
}
LL Inv(LL a, LL mod)///mod任意
{
LL d, x, y;
Ex_gcd(a, mod, d, x, y);
return d == 1? (x + mod) % mod : -1;
}
LL Inv(LL x, LL mod)///mod为质数
{
LL r, y;
for(r = 1, y = mod - 2; y; x = x * x % mod, y >>= 1)
(y & 1) && (r = r * x % mod);
return r;
}
///n 个方程:x = a[i](mod m[i]) (0 <= i < n)
LL China(int n, int *a, int *m)
{
LL M = 1, d, y, x = 0;
for(int i = 0; i < n; i ++) M *= m[i];
for(int i = 0; i < n; i ++)
{
LL w = M / m[i];
Ex_gcd(m[i], w, d, d, y);
x = (x + y * w * a[i]) % M;
}
return (x + M) % M;
}
毕达哥拉斯三原组个数,(x,y,x <= n)
x = m^2 - n^2;
y = 2*n*m;
z = m^2 + n^2;
代码实现:
int Pythagoras(int n)
{
int tmp, ret = 0;
tmp = sqrt(n + 0.0);
for(int i = 1; i * i <= tmp; i ++)
{
for(int j = i + 1; j <= tmp; j ++)
{
if(i * i + j * j > n) break;
if(i % 2 == j % 2) continue;
if(gcd(i, j) == 1)
{
ret ++;
}
}
}
return ret;
}
佩尔方程 x^2 - d*y^2 = 1;(d为非完全平方数)。
迭代式:
x(n) = x(n-1) * x(1) + d * y(n-1) * y (1).
y(n) = x(n-1) * y(1) + y(n-1) * x(1).
暴力求解x1, y1。如下:
int x, y, d;
void search()求
{
y = 1;
while(true)
{
x = sqrt((double)d * y * y + 1);
if(x * x - d * y * y == 1)
break;
y ++;
}
}
欧拉定理: 对任何两个互质的正整数a,m(m > 1)有 a^(Euler(m)) = 1 (mod m)。
求整数n的因子个数:
设n = p1^a1 * p2^a2 *...* pn^an; 即将n质分解。
则count(n) = (a1 + 1)*(a2 + 1) *...*(an + 1).
int count(int n)
{
int ret = 1;
for(int i = 2; i * i <= n; i ++) if(n % i == 0)///可以用素数表优化。
{
int tmp = 0;
while(n % i == 0)
{
n /= i;
tmp ++;
}
ret *= tmp + 1;
}
if(n > 1) ret *= 2;
return ret;
}
求1-n,和n的最小公倍数的和~!
LL a[maxn + 1000];
LL phi[maxn+1000];
void initphi()
{
LL i,j;
for(i = 1; i <= maxn; i++) phi[i] = i;
for(i = 2; i <= maxn; i += 2) phi[i] /= 2;
for(i = 3; i <= maxn; i += 2) if(phi[i] == i)
{
for(j = i; j <= maxn; j += i)
phi[j] = phi[j] / i * (i - 1);
}
}
LL get(LL now)
{
LL res=0;
for(LL i=1; i*i<=now; i++)
{
if(now%i==0)
{
if(i*i==now) res+=i*phi[i];
else res+=i*phi[i]+(now/i)*phi[now/i];
}
}
return res;
}
void init()
{
initphi();
for(int i=1; i<=maxn; i++)
{
a[i]=i*(1ll+get(i))/2ll;
}
}
求整数n的所有正因子和:
设n = p1^a1 * p2^a2 *...* pn^an; 即将n质分解。
设ji(i) = 1 + pi + pi^2 + ... + pi^ai = (pi^(ai + 1) - 1) / (pi - 1);
则sum(n) = ji(1) * ji(2) * ... ji(n).
int sum(int n)
{
int ret = 1;
for(int i = 2; i * i <= n; i ++) if(n % i == 0)
{
int tmp = 1;
while(n % i == 0)
{
n /= i;
tmp *= i;
}
ret *= (tmp * i - 1) / (i - 1);
}
if(n > 1) ret *= (1 + n);
return ret;
}
const int maxn = 1001000;
LL dp[maxn];
int isp[maxn];
void init()
{
CLR(isp, -1);
for(int i = 2; i < maxn; i ++) if(isp[i] == -1)
{
isp[i] = i;
if(i < 11111) for(int j = i * i; j < maxn; j += i)
{
isp[j] = i;
}
}
dp[1] = 1;
for(int i = 2; i < maxn; i ++)
{
LL cnt = isp[i];
int a = i;
while(a % isp[i] == 0)
{
cnt *= isp[i];
a /= isp[i];
}
dp[i] = dp[a] * (cnt - 1) / (isp[i] - 1);
}
}