欧拉计划 71~80
71
考虑小于 \(\frac{3}{7}\) 且分母小于等于 \(10^6\) 以内的最大分数是什么，注意到 \(7 \mid 999999\)，故我们通分 \(\frac{3}{7}\) 得到 \(\frac{428571}{999999}\)，由不等关系可知，合法的数字为 \(\frac{428570}{999999}\)。
# python
print(3 * (999999 // 7) - 1)
72
求分母小于 \(10^6\) 的最简真分数个数，我们转化为求对于所有分母，比它小且与它互质的数的个数，于是筛出欧拉函数求和即可。
// cpp
#include<bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10;
int primes[N], st[N], cnt;
int phi[N];
int main() {
phi[1] = 1;
for (int i = 2; i < N; i ++ ) {
if (!st[i]) primes[cnt ++ ] = i, phi[i] = i - 1;
for (int j = 0; i * primes[j] < N; j ++ ) {
st[i * primes[j]] = 1;
if (i % primes[j] == 0) {
phi[i * primes[j]] = phi[i] * primes[j];
break;
}
phi[i * primes[j]] = phi[i] * (primes[j] - 1);
}
}
long long ans = 0;
for (int i = 2; i <= 1e6; i ++ ) ans += phi[i];
cout << ans;
return 0;
}
73
Solution 1
在看题解前，你需要了解 Farey_sequence。
对 \(n\) 阶 Farey 数列，相邻两项 \(\frac{a}{b} < \frac{c}{d}\) 满足 \(bc - ad = 1\)，因此对于数列的下一项我们可以求出一个最小的 \(\frac{c^*}{d^*}\) 满足 \(bc^* - ad^* = 1\)，那么其解系为
于是求出下一项之后，根据相邻两项 \(\frac{a}{b}, \, \frac{c}{d}\) 均在 Stern-Brocot 树上，我们有下一项 \(\frac{p}{q}\) 满足 \(\frac{c}{d} = \frac{a + p}{b + q}\)，也就是 \((p + a)d = (q + b)c\)，因此存在 \(k\) 使得
为了使 \(p, \, q\) 精度高，那么 \(k\) 的值应该尽可能大，所以对上界 \(n\) 来说，\(p, \, q\) 由下式定义：
显然由真分数定义知 \(q > p\)，所以最大 \(k = \left\lfloor\frac{n + b}{d}\right\rfloor\)，因此
根据 \(n\) 阶 Farey 序列的个数的估计为 \(\frac{3}{\pi^2}n^2\)，本题取 \(n = 12000\) 且具有 \(\frac{1}{6}\) 常数，统计出来有 \(7295372\) 个。
# python
n = 12000
ans = 0
a, b = 1, 3
c, d = 1, 2
k = min((n - c) // a, (n - d) // b)
c, d = c + k * a, d + k * b
while c * 2 < d:
k = (n + b) // d
a, b, c, d = c, d, k * c - a, k * d - b
ans += 1
print(ans)
# takes 3.355603100033477s
solution 2
因为 \(\frac{1}{2}\) 和 \(\frac{1}{3}\) 是 \(2\) 层 Stern-Brocot 树上的相邻节点，因此我们直接考虑构建出两者中的 Stern-Brocot 树，统计其答案即可。
// cpp
#include<bits/stdc++.h>
using namespace std;
int build(int a, int b, int c, int d, int n) {
if (b + d > n) return 0;
return 1 + build(a, b, a + c, b + d, n) + build(a + c, b + d, c, d, n);
}
int main() {
int ans = build(1, 3, 1, 2, 12000);
cout << ans << "\n";
return 0;
}
// CPU time used: 0.028 s
Solution 3
\(O(n\log{n})\) 与 \(O(n^{\frac{2}{3}} + \sqrt{n}\log{n})\)，详见Farey 数列。
74
Solution 1
很暴力的做法，直接爆搜。
// cpp
#include<bits/stdc++.h>
using namespace std;
int fact[10];
int get_sum(int x) {
string s = to_string(x);
int t = 0;
for (auto c : s) t += fact[c - '0'];
return t;
}
int main() {
fact[0] = 1;
for (int i = 1; i < 10; i ++ ) fact[i] = fact[i - 1] * i;
int ans = 0;
for (int i = 1; i <= 1000000; i ++ ) {
unordered_map<int, int> mp;
int len = 0, x = i;
while (!mp.count(x) && len <= 60) mp[x] = 1, len ++ , x = get_sum(x);
if (len == 60) ans ++ ;
}
cout << ans << "\n";
return 0;
}
Solution 2
考虑到 \(999,999\) 的下一个数字会是 \(2,177,280\)，我们处理出 \(2,200,000\) 以内的所有数的下一个数，通过拓扑排序留下所有环，接下来我们处理出所有的环的长度，那么对于一个数字，我们走到环上就可以立即停止了。
// cpp
#include<bits/stdc++.h>
using namespace std;
int fact[10], ans;
int ring[2200000], nxt[2200000];
int d[2200000], st[2200000];
int get_sum(int x) {
string s = to_string(x);
int t = 0;
for (auto c : s) t += fact[c - '0'];
return t;
}
int get_len(int x, int dep) {
st[x] = 1;
if (st[nxt[x]]) return dep;
return get_len(nxt[x], dep + 1);
}
void dfs(int x, int dep) {
if (ring[x]) {
if (dep + ring[x] == 60) ans ++ ;
return;
}
dfs(nxt[x], dep + 1);
}
void solve() {
fact[0] = 1;
for (int i = 1; i < 10; i ++ ) fact[i] = fact[i - 1] * i;
for (int i = 1; i < 2200000; i ++ ) nxt[i] = get_sum(i), d[nxt[i]] ++ ;
queue<int> q;
for (int i = 1; i < 2200000; i ++ ) if (!d[i]) q.push(i);
while (q.size()) {
int t = q.front(); q.pop();
d[nxt[t]] -- ;
if (!d[nxt[t]]) q.push(nxt[t]);
}
for (int i = 1; i < 2200000; i ++ ) {
if (d[i] && !st[i]) {
int len = get_len(i, 1), x = nxt[i];
ring[i] = len;
while (x != i) ring[x] = len, x = nxt[x];
}
}
for (int i = 1; i <= 1000000; i ++ ) dfs(i, 0);
cout << ans << "\n";
}
int main() {
clock_t start = clock();
solve();
clock_t end = clock();
double cpu_time_used = static_cast<double>(end - start) / CLOCKS_PER_SEC;
cout << "CPU time used: " << cpu_time_used << " s\n";
return 0;
}
75
考虑本原毕达哥拉斯三元数，利用 #39 中的生成矩阵生成本原三角形，枚举其倍数覆盖即可。
# python
from numpy import *
A = matrix([[1, -2, 2],
[2, -1, 2],
[2, -2, 3]])
B = matrix([[1, 2, 2],
[2, 1, 2],
[2, 2, 3]])
C = matrix([[-1, 2, 2],
[-2, 1, 2],
[-2, 2, 3]])
O = matrix([[3], [4], [5]])
triangles = [O]
maxn = 1500000
L = []
for ori in triangles:
L.append(sum(ori))
nxtA, nxtB, nxtC = A * ori, B * ori, C * ori
if sum(nxtA) <= maxn: triangles.append(nxtA)
if sum(nxtB) <= maxn: triangles.append(nxtB)
if sum(nxtC) <= maxn: triangles.append(nxtC)
cnt = [0 for _ in range(maxn + 1)]
for ay in L:
for i in range(ay, maxn + 1, ay):
cnt[i] += 1
ans = 0
for i in range(maxn + 1):
if cnt[i] == 1:
ans += 1
print(ans)
发现 python 有点太卡了，且 sympy 比 numpy 慢三倍以上，我重新用 c++ 写了一份。
// cpp
#include<bits/stdc++.h>
using namespace std;
vector<vector<int>> A = {{1, -2, 2},
{2, -1, 2},
{2, -2, 3}};
vector<vector<int>> B = {{1, 2, 2},
{2, 1, 2},
{2, 2, 3}};
vector<vector<int>> C = {{-1, 2, 2},
{-2, 1, 2},
{-2, 2, 3}};
vector<int> O = {3, 4, 5};
int sum(vector<int> a) {
return accumulate(a.begin(), a.end(), 0);
}
vector<int> mul(vector<int> a, vector<vector<int>> b) {
vector<int> tmp(3);
for (int i = 0; i < 3; i ++ ) {
for (int j = 0; j < 3; j ++ ) {
tmp[i] += a[j] * b[i][j];
}
}
return tmp;
}
int maxn = 1500000;
queue<vector<int>> triangles;
unordered_map<int, int> st, cnt;
void solve() {
triangles.push(O);
while (!triangles.empty()) {
vector<int> ori = triangles.front(); triangles.pop();
st[sum(ori)] ++ ;
vector<int> nxtA = mul(ori, A), nxtB = mul(ori, B), nxtC = mul(ori, C);
if (sum(nxtA) <= maxn) triangles.push(nxtA);
if (sum(nxtB) <= maxn) triangles.push(nxtB);
if (sum(nxtC) <= maxn) triangles.push(nxtC);
}
int ans = 0;
for (auto [x, y] : st) {
for (int d = x; d <= maxn; d += x) {
cnt[d] += y;
}
}
for (auto [x, y] : cnt) if (y == 1) ans ++ ;
cout << ans << "\n";
}
int main() {
clock_t start = clock();
solve();
clock_t end = clock();
double cpu_time_used = static_cast<double>(end - start) / CLOCKS_PER_SEC;
cout << "CPU time used: " << cpu_time_used << " s\n";
return 0;
}
// CPU time used: 0.515 s
这个代码可以通过 HackerRank 上 5 \time 10^6 的超大数据。
// cpp
#include<bits/stdc++.h>
using namespace std;
vector<vector<int>> A = {{1, -2, 2},
{2, -1, 2},
{2, -2, 3}};
vector<vector<int>> B = {{1, 2, 2},
{2, 1, 2},
{2, 2, 3}};
vector<vector<int>> C = {{-1, 2, 2},
{-2, 1, 2},
{-2, 2, 3}};
vector<int> O = {3, 4, 5};
int sum(vector<int> &a) {
return a[0] + a[1] + a[2];
}
vector<int> mul(vector<int> &a, vector<vector<int>> &b) {
vector<int> tmp(3);
for (int i = 0; i < 3; i ++ ) {
for (int j = 0; j < 3; j ++ ) {
tmp[i] += a[j] * b[i][j];
}
}
return tmp;
}
#define PII pair<int, int>
int maxn = 5e6;
vector<int> ans(maxn + 1);
queue<vector<int>> triangles;
unordered_map<int, int> st, cnt;
priority_queue<PII, vector<PII>, greater<PII>> heap;
void init() {
triangles.push(O);
while (!triangles.empty()) {
vector<int> ori = triangles.front(); triangles.pop();
st[sum(ori)] ++ ;
vector<int> nxtA = mul(ori, A), nxtB = mul(ori, B), nxtC = mul(ori, C);
if (sum(nxtA) <= maxn) triangles.push(nxtA);
if (sum(nxtB) <= maxn) triangles.push(nxtB);
if (sum(nxtC) <= maxn) triangles.push(nxtC);
}
int res = 0;
for (auto [x, y] : st) heap.push({x, x});
for (int i = 12; i <= maxn; i ++ ) {
while (heap.size() && heap.top().first <= i) {
auto t = heap.top(); heap.pop();
int now = t.first, add = t.second, num = st[t.second];
if (!cnt.count(now) && num == 1) res ++ ;
else if (cnt[now] == 1) res -- ;
cnt[now] += num;
if (now + add <= maxn) heap.push({now + add, add});
}
ans[i] = res;
}
}
void solve() {
int n; cin >> n;
cout << ans[n] << "\n";
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
init();
int T = 1;
cin >> T;
while (T -- ) solve();
return 0;
}
// Used: 2249 ms, 86084 KB
76
简单完全背包
// cpp
#include<bits/stdc++.h>
using namespace std;
int dp[110];
int main() {
dp[0] = 1;
for (int i = 1; i < 100; i ++ ) {
for (int j = i; j <= 100; j ++ ) {
dp[j] += dp[j - i];
}
}
cout << dp[100];
return 0;
}
77
同样是完全背包问题
// cpp
#include<bits/stdc++.h>
using namespace std;
const int N = 1000;
int primes[N], st[N], cnt;
int dp[100];
int main() {
for (int i = 2; i < N; i ++ ) {
if (!st[i]) primes[cnt ++ ] = i;
for (int j = 0; i * primes[j] < N; j ++ ) {
st[i * primes[j]] = 1;
if (i % primes[j] == 0) break;
}
}
dp[0] = 1;
for (int i = 0; i < 100; i ++ ) {
for (int j = primes[i]; j < 100; j ++ ) {
dp[j] += dp[j - primes[i]];
}
}
for (int i = 1; i < N; i ++ ) {
if (dp[i] > 5000) {
cout << i << "\n";
return 0;
}
}
return 0;
}
78
普通完全背包干了 \(30000\) 以内没有搜到解，我们只能另辟奇径。
Attempt 1
考虑杨表，设 \(f[i][j]\) 表示考虑前 \(i\) 个数，和为 \(j\) 的选法，\(g[i][j]\) 表示考虑拆成了 \(i\) 个数，和为 \(j\) 的选法。前者就是一个简单的完全背包，后者则是利用了杨表的性质，我们显然可以给每个数增加 \(1\)，于是有 \(g[i][j - i]\) 可以转移到 \(g[i][j]\)；我们也可以新开一个数，使得 \(g[i - 1][j - 1]\) 转移到 \(g[i][j]\)。
考虑根号分治，对于不大于 \(m = \sqrt{n}\) 的数字，我们用第一种方法，对于大于 \(m\) 的数字，我们用杨表强行在开新数字的时候至少开 \(m + 1\) 即可，也就是说 \(g[i][j] = g[i][j - i] + g[i][j - m - 1]\)，那么最后统计答案只需要把两部分的答案并起来即可。
但是，很不幸的是，这个做法只能在 \(O(n^{1.5})\) 的时间复杂度内求出一个值，但是你可以注意到，对 \(g\) 的每一层求和后，答案会是一个卷积形式，我们可以，利用多项式手段来 \(O(n\log{n})\) 优化，但是用多项式显然不如我们直接去求生成函数了。
Solution 1
我们考虑欧拉五边形数定理：
设欧拉五边形函数 \(\phi(x) = \prod\limits_{n = 0}^{\infty}(1 - x^n)\)，它的展开形式为
其中我们可以得到 \(k\) 取奇数项时广义五边形数为 \(\frac{3k^2 + k}{2}\)，取偶数项时广义五边形数为 \(\frac{3k^2 - k}{2}\)，那么我们可以简单的推出 \(\phi(x)\) 的展开形式。
考虑分拆数 \(p(n)\) 的生成函数形式：
而 \(\frac{1}{\phi(x)} = \prod\limits_{n = 1}^{\infty}\frac{1}{1 - x^n}\)，容易看出，分拆数是欧拉五边形函数的倒数的第 \(n\) 次项系数，于是我们有
于是，我们整理出
这个东西的系数可以规约为欧拉五边形数 \(\frac{3k^2 + k}{2}\) 按顺序取 \(0, \, 1, \, -1, \, 2, \, -2, \, \cdots\) 的值。
按照模 \(4\) 分类，我们就可以通过递推的方式求出所有的 \(p(n)\)，因此，找到最小的一个模 \(10^6\) 余 \(0\) 的数即可。
// cpp
#include<bits/stdc++.h>
using namespace std;
const int N = 1e5 + 10;
int p[N], mod = 1e6;
int f(int x, int i) {
return x * (x * 3 + (i & 1 ? 1 : -1)) / 2;
}
int main() {
p[0] = 1;
for (int i = 1; i < N; i ++ ) {
for (int j = 0; f(j >> 1, j) <= i; j ++ ) {
if (j & 2) p[i] += p[i - f(j >> 1, j)];
else p[i] += mod - p[i - f(j >> 1, j)];
p[i] %= mod;
}
if (!p[i]) cout << i << "\n";
}
return 0;
}
Solution 2
这个完全背包可以考虑构建生成函数求，具体见 多项式例题。
79
这个是个 puzzle，直接瞪眼吧。
统计数字出现了 \(0, \, 1, \, 2, \, 3, \, 6, \, 7, \, 8, \, 9\)，注意到 \(762, \, 790, \, 718, \, 731\)，因此第一个数字为 \(7\)，接下来发现 \(319, \, 162\)，所以 \(0, \, 1, \, 2, \, 6, \, \, 8, \, 9\) 全在 \(3\) 后出现了，第二个数字为 \(3\)，然后再加上 \(316\) 这个数字，可以推出下一个数字为 \(1\)，库内还有 \(289, \, 162\)，所以后面的数字就确定为 \(62890\) 了。
答案为 \(73162890\)。
80
我该怎么说呢，python 王朝了。
# python
from decimal import *
getcontext().prec = 105
ans = 0
for i in range(2, 101):
if i ** 0.5 == int(i ** 0.5): continue
sqr = str(Decimal(i) ** Decimal(0.5))
sqr = sqr.replace('.', '')
ans += sum(int(j) for j in sqr[:100])
print(ans)
浙公网安备 33010602011771号