标签：
题意：给出n(1<=n<=10^8)，求小于n的，求所有与n互质的数字的四次幂的加和是多少。
分析：容斥原理
首先要知道四次幂求和公式，1^4+2^4+...+n^4=n*(n+1)*(2n+1)*(3n^2+3n-1)/30
先求所有小于等于n的数字的四次幂和，然后减去那些不互质的即可。
这个减去的过程用到了容斥原理。
先对n分解质因子，每个不同的质因子只保留一个。
然后分别枚举这些质因子的组合情况，由奇数个因子组成的数要减去，由偶数个因子组成的数要加上。
对于一个因子组合的乘积a，我们需要一次性计算a^4+(2a)^4 + (3a)^4+...
将其转化为a^4 * (1^4+2^4+...)即可。
这道题还有一个难点，就是公式中有除法（除以30），却还要进行模运算。
除法是不支持模运算的，因此我们要将除法转化为乘法，除以30变为乘以30的逆元。
求逆元可以用扩展欧几里德gcd(30,MOD,x,y)，把x/gcd(30,MOD)整理到0～MOD-1范围内即为30的逆元。
具体原因查阅扩展欧几里德算法，和模运算中逆元的概念。
#include <cstdio> using namespace std; #define D(x) const int MOD = (int)(1e9) + 7; const int MAX_FACTOR = 40; int n; int factor_num; long long factor[MAX_FACTOR]; long long inverse; //n(n+1)(2n+1)(3n^2+3n-1)/30 long long to_forth(long long value) { long long ret = value; ret = ret * ret % MOD; ret = ret * ret % MOD; return ret; } long long cal(long long value) { long long num = n / value; long long ret = 1; ret = ret * num % MOD * (num + 1) % MOD; ret = ret * (2 * num + 1) % MOD; ret = ret * ((num * num % MOD * 3 % MOD + 3 * num % MOD - 1) % MOD) % MOD; if (ret / 30 != ret * inverse % MOD) { D(printf("#%lld %lld\n", ret / 30, ret * inverse % MOD)); }else { D(printf("**\n")); } ret = ret * inverse % MOD; ret = ret * to_forth(value) % MOD; return ret; } void get_factors() { factor_num = 0; int m = n; for (int i = 2; i * i <= m; i++) { if (m % i == 0) factor[factor_num++] = i; while (m % i == 0) { m /= i; } } if (m != 1) { factor[factor_num++] = m; } } long long work() { long long ans = 0; for (int i = 1; i < (1 << factor_num); i++) { int num = 0; long long temp = 1; int index = 0; for (int mask = 1; mask <= i; mask <<= 1, index++) { if ((mask & i) == 0) { continue; } num++; temp *= factor[index]; } D(printf("temp=%lld\n", temp)); if (num & 1) ans += cal(temp); else ans -= cal(temp); ans = (ans % MOD + MOD) % MOD; } ans = ((cal(1) - ans) % MOD + MOD) % MOD; return ans; } void gcd_extend(long long a,long long b,long long &g,long long &x,long long &y) { if (!b) { g = a; x = 1; y = 0; return; } gcd_extend(b, a % b, g, y, x); y -= a / b * x; } int main() { long long x, y, g; gcd_extend(30, MOD, g, x, y); D(printf("%lld %lld %lld\n", x, y, g)); x = (x % MOD + MOD) % MOD; inverse = x / g; D(printf("%lld\n", inverse)); D(printf("%lld\n", inverse * 30 % MOD)); int t; scanf("%d", &t); while (t--) { scanf("%d", &n); if (n == 1) { puts("0"); continue; } get_factors(); int ans_int = work(); printf("%d\n", ans_int); } return 0; }
标签：
原文地址：http://www.cnblogs.com/rainydays/p/4376075.html
