number.cpp

ACM经典算法ACM经典算法ACM经典算法

/*============================================================*\
 | 递推求phi(i)                                               |
\*============================================================*/
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);
}
/*============================================================*\
 | 单独求phi(x)                                               |
\*============================================================*/
unsigned euler(unsigned x)
{
	unsigned i, res=x;
	for (i = 2; i < (int)sqrt(x * 1.0) + 1; i++) if(x%i==0) {
		res = res / i * (i - 1);
		while (x % i == 0) x /= i;
	}
	if (x > 1) res = res / x * (x - 1);
	return res;
}
/*============================================================*\
 | GCD 最大公约数                                             |
\*============================================================*/
int gcd(int x, int y)
{
	if (!x || !y) return x > y ? x : y;
	for (int t; t = x % y; x = y, y = t);
	return y;
}
/*============================================================*\
 | 快速 GCD                                                   |
\*============================================================*/
int kgcd(int a, int b)
{
	if (a == 0) return b;
	if (b == 0) return a;
	if (!(a & 1) && !(b & 1)) return kgcd(a>>1, b>>1) << 1;
	else if (!(b & 1)) return kgcd(a, b>>1);
	else if (!(a & 1)) return kgcd(a>>1, b);
	else return kgcd(abs(a - b), min(a, b));
}
/*============================================================*\
 | 扩展 GCD                                                   |
 | 求x, y使得gcd(a, b) = a * x + b * y;                       |
\*============================================================*/
int extgcd(int a, int b, int & x, int & y)
{
   if (b == 0) { x=1; y=0; return a; }
   int d = extgcd(b, a % b, x, y);
   int t = x; x = y; y = t - a / b * y;
   return d;
}
/*============================================================*\
 | 模线性方程 a * x = b (% n)                                 |
\*============================================================*/
void modeq(int a, int b, int n)   // ! n > 0
{
	int e, i, d, x, y;
	d = extgcd(a, n, x, y);
	if (b % d > 0) printf("No answer!\n");
	else {
		e = (x * (b / d)) % n;
		for (i = 0; i < d; i++) // !!! here x maybe < 0
			printf("%d-th ans: %d\n", i+1, (e+i*(n/d))%n);
	}
}
/*============================================================*\
 | 模线性方程组                                               |
 | a=B[1](% W[1]); a=B[2](% W[2]); ......  a=B[k](mod W[k]);  |
 | 其中W, B已知，W[i]>0且W[i]与W[j]互质, 求a; (中国余数定理)  |
\*============================================================*/
int china(int b[], int w[], int k)
{
	int i, d, x, y, m, a = 0, n = 1;
	for (i = 0; i < k; i++) n *= w[i];   // ! 注意不能overflow
	for (i = 0; i < k; i++) {
		m = n / w[i];
		d = extgcd(w[i], m, x, y);
		a = (a + y * m * b[i]) % n;
	}
	if (a > 0) return a;
	else return (a + n);
}
/*============================================================*\
 | 筛素数 [1..n]                                              |
\*============================================================*/
char is[N]; int prm[N];
int getprm(int n)
{
	int i, j, k = 0;
	int s, e = (int)(sqrt(0.0 + n) + 1);
	memset(is, 1, sizeof(is));
	prm[k++] = 2; is[0] = is[1] = 0;
	for (i = 4; i < n; i += 2) is[i] = 0;
	for (i = 3; i < e; i += 2) if (is[i]) {
		prm[k++] = i;
		for (s = i * 2, j = i * i; j < N; j += s) is[j] = 0;
	}
	for ( ; i < n; i += 2) if (is[i]) prm[k++] = i;
	return k;               // 返回素数的个数
}
/*============================================================*\
 | 随机素数测试                                               |
 | CALL: bool res = miller(n);                                |
 | 快速测试n是否满足素数的'必要'条件, 出错概率很小;           |
 | 对于任意奇数 n>2 和正整数 s, 算法出错概率 <= 2^(-s);       |
\*============================================================*/
int witness(int a, int n)
{
	int x, d=1, i = ceil(log(n - 1.0) / log(2.0)) - 1;
	for ( ; i >= 0; i--) {
		x = d; d = (d * d) % n;
		if (d==1 && x!=1 && x!=n-1) return 1;
		if (((n-1) & (1<<i)) > 0) d = (d * a) % n;
	}
	return (d == 1 ? 0 : 1);
}
int miller(int n, int s = 50)
{
	if (n == 2) return 1;
	if ((n % 2) == 0) return 0;
	int j, a;
	for (j = 0; j < s; j++) {
		a = rand() * (n-2) / RAND_MAX + 1;
		if (witness(a, n)) return 0;
	}
	return 1;
}
/*============================================================*\
 | 组合数学相关                                               |
\*============================================================*/
1. {1, 2, ... n}的r组合a1, a2, ... ar出现在所有r组合中的字典
序位置编号, C(n, m)表示n中取m的组合数; index =
C(n, r) - C(n - a1, r) - C(n - a2, r-1) - ... - C(n - ar, 1)

2. k * C(n, k) = n * C(n-1, k-1);
   C(n, 0) + C(n, 2) + ... = C(n, 1) + C(n, 3) + ...
   1 * C(n, 1) + 2 * C(n, 2) + ... + n * C(n, n) = n * 2^(n-1)

3. Catalan数:   C_n = C(2*n, n) / (n+1)
                C_n = (4*n-2)/(n+1) * C_n-1
                C_1 = 1

4. 第二类Stirling数: S(p, k) = k * S(p-1, k) + S(p-1, k-1).
   S(p, 0) = 0, (p>=1);  S(p, p) = 1, (p>=0);
   且有 S(p, 1) = 1, (p>=1);
        S(p, 2) = 2^(p-1) - 1, (p>=2);
        S(p, p-1) = C(p, 2);
   含义: 将p个元素划分到k个同样的盒子, 每个盒子非空的方法数.
              1      k
   S(p, k) = --- * sigma((-1)^t * C(k, t) * (k-t)^p)
              k!    t=0

5. Bell数: B_p = S(p, 0) + S(p, 1) + ... + S(p, p)
   B_p = C(p-1, 0)*B_0 + C(p-1, 1)*B_1 + ... C(p-1, p-1)*B_(p-1)

6. 第一类stirling数:
   s(p, k)是将p个物体排成k个非空的循环排列的方法数.
   (或者: 把p个人排成k个非空圆圈的方法数)
   s(p, k) = (p-1) * s(p-1, k) + s(p-1, k-1);
/*============================================================*\
 | Polya计数                                                  |
 | c种颜色的珠子, 组成长为s的项链, 项链没有方向和起始位置;    |
\*============================================================*/
int gcd (int a, int b) { return b ? gcd(b,a%b) : a; }
int main (void)
{
	int c, s;
	while (scanf("%d%d", &c, &s)==2) {
		int k;
		long long p[64]; p[0] = 1; // power of c
		for (k=0 ; k<s ; k++) p[k+1] = p[k] * c;
		// reflection part
		long long count = s&1 ? s * p[s/2 + 1] :
		                  (s/2) * (p[s/2] + p[s/2 + 1]);
		// rotation part
		for (k=1 ; k<=s ; k++) count += p[gcd(k, s)];
		count /= 2 * s;
		cout << count << '\n';
	}
	return 0;
}