Project_Euler-12 题解
题目
思路
我们可以从小到大枚举每一个三角形数,然后计算他们的约数个数,从而得到结果。
代码
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
//#include "./DBG_text.h"
#define MAX_N 100
typedef long long ll;
ll get_factory(ll n) {
ll cnt = 0;
for (ll i = 1, I = sqrt(n); i <= I; i++) {
if (n % i) continue;
cnt += 1;
cnt += (i != n / i);
}
return cnt;
}
ll triangle(ll n) {
return n * (n + 1) >> 1;
}
ll main() {
ll n = 0;
while (1) {
n++;
// 求第n个三角形数
ll val = triangle(n);
// 计算因子个数
ll cnt = get_factory(val);
//DBG("cnt = %lld\n", cnt);
if (cnt <= 500) continue;
printf("val = %lld, cnt = %lld\n", val, cnt);
break;
}
return 0;
}
优化思路
这里不得不提到 一个定理了:
约数个数定理
约数个数定理:如果一个大于1的正整数 n
可以分解为质因数的形式:n = p1^a1 * p2^a2 * ... * pk^ak
,其中 p1, p2, ..., pk
是 n
的不同质因数,a1, a2, ..., ak
是对应的幂次,那么 n
的约数个数就是 (a1+1) * (a2+1) * ... * (ak+1)
。
举个例子,假设 n = 2^3 * 3^2 * 5^1
,那么 n
的约数个数就是 (3+1) * (2+1) * (1+1) = 24
。
有了这个定理可以得到什么呢?当一个数字n
被确定时,我们通过素数快速计算出它的约数个数。
但是前提,我们需要先知道它的素因子是什么,这就需要我们使用线性筛来求解。
不管先不着急,因为我们要求的是三角形数字的因数个数,而不是n
的,因此我们需要进一步处理:
具体是这样的做的:
我们设 F ( n ) F(n) F(n) 表示数字 n n n 的因数个数, f ( n ) f(n) f(n) 表示第 n n n 个三角形数字,那么我们可以得到下面的关系:
F ( f ( n ) ) = F ( n ∗ ( n + 1 ) 2 ) F(f(n)) = F(\frac {n * (n + 1)} {2}) F(f(n))=F(2n∗(n+1))
此时我们发现 n n n 和 n + 1 n + 1 n+1 一定是一奇一偶的关系,并且他们一定互素。
两个互素的数相乘得到的数的因数 ( F ( n ∗ ( n + 1 ) ) F(n * (n + 1)) F(n∗(n+1))) 的个数等于什么?
F ( n ∗ ( n + 1 ) ) = F ( n ) ∗ F ( n + 1 ) F(n*(n + 1)) = F(n) * F(n + 1) F(n∗(n+1))=F(n)∗F(n+1)
为什么呢?因为两个数的最大公因数为1,所以他们没有相同的因数,因此他们两个相乘得到的因数的个数就是他们自己拆解后的因数不断排列组合得到的因数的个数。
举个例子: 5 5 5 和 6 6 6 。 5 5 5 的因数个数为 2 2 2, 6 6 6 的因数个数为 4 4 4 ,那么他们的乘积就是 30 30 30 , 30 30 30 的因数个数就是 2 ∗ 4 = 8 2 * 4 = 8 2∗4=8 个。
那么回到上面的式子:
F ( f ( n ) ) = F ( n ∗ ( n + 1 ) 2 ) F(f(n)) = F(\frac {n * (n + 1)} {2}) F(f(n))=F(2n∗(n+1))
就可以化简为:
F ( f ( n ) ) = F ( n ) ∗ F ( n + 1 ) 2 F(f(n)) = \frac{F(n) * F(n + 1)}{2} F(f(n))=2F(n)∗F(n+1)
对n分奇偶情况讨论,可以将2化简掉:
当 n 为偶数时: F ( f ( n ) ) = F ( n 2 ) ∗ F ( n + 1 ) 当 n 为奇数时: F ( f ( n ) ) = F ( n ) ∗ F ( n + 1 2 ) 当n为偶数时: F(f(n)) = F(\frac {n}{2}) * F(n + 1) \\ 当n为奇数时: F(f(n)) = F(n) * F(\frac {n + 1} {2}) 当n为偶数时:F(f(n))=F(2n)∗F(n+1)当n为奇数时:F(f(n))=F(n)∗F(2n+1)
得到这个公式,我们就可以通过 n n n 来表达第 n n n 个三角形数的因数个数了。
程序的主框架就完成了:
int main() {
// ... init();
ll n = 0;
while (1) {
n++;
ll val;
val = (n & 1) ? f[n] * f[(n + 1) >> 1] : f[n >> 1] * f[n + 1];
if (val <= 500) continue;
printf("n = %lld\n", n * (n + 1) >> 1);
break;
}
return 0;
}
接下来就是怎么确定 F ( n ) F(n) F(n) ,我们必须保证,当我们拿出一个 n n n 的时候,要得到 n n n 对应的因数的个数。
我们可以将 n n n 对应的因数个数存放在数组里。
我们需要找到 n n n 的最小素因子(prime[j]),这个可以通过遍历素数表来实现,但是对于 n n n 本身,我们需要在线性筛中实现( n = i ∗ p r i m e [ j ] n = i * prime[j] n=i∗prime[j]), i i i 是线性筛的遍历值,在筛中有以下几种情况讨论:
- 当 n n n为素数时:其因数个数只有1和它本身,因此就是2
- 当
i
i
i与某个素数
j
j
j互素的时候,
f[prime[j] * i] = f[prime[j]] * f[i];
- 当
i
i
i与某个素数
j
j
j非互素的时候
f[prime[j] * i] = (f[i] / (cnt[i] + 1)) * (cnt[i] + 2);
优化代码
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include "./DBG_text.h"
#define MAX_N 100000
typedef long long ll;
int prime[MAX_N + 5];
int f[MAX_N + 5];
int cnt[MAX_N + 5]; // 表示的是数字i中最小素因子的幂次
void init() {
for (int i = 2; i <= MAX_N; i++) {
if (!prime[i]) {
prime[++prime[0]] = i;
f[i] = 2;
cnt[i] = 1;
}
for (int j = 1; j <= prime[0]; j++) {
if (i * prime[j] > MAX_N) break;
prime[i * prime[j]] = 1;
// == 0 说明不互素并且prime[j]是i最小的素因数,否则说明互素
if (i % prime[j] == 0) {
// 这里的prime[j]就是n的最小素因数,cnt[i]中记录的就是它的幂数
// prime[j] * i说明n这个数相对于i这个数的最小素因数增加了1幂,此时需要先让f[i]先除掉原来的最小幂
// 然后再乘以新的幂数也就是+2
// 这个2中有一个1表示prime[j] * i之后,增加的1幂
// 另一个1表示根据约数个数定理而必须+1
f[prime[j] * i] = (f[i] / (cnt[i] + 1)) * (cnt[i] + 2); // 这句话较难理解
cnt[prime[j] * i] = cnt[i] + 1;
break;
} else {
f[prime[j] * i] = f[prime[j]] * f[i];
cnt[prime[j] * i] = 1;
}
}
}
return ;
}
int main() {
init();
ll n = 0;
while (1) {
n++;
ll val;
val = (n & 1) ? f[n] * f[(n + 1) >> 1] : f[n >> 1] * f[n + 1];
if (val <= 500) continue;
printf("n = %lld\n", n * (n + 1) >> 1);
break;
}
return 0;
}
补充
对于这段代码:
f[prime[j] * i] = (f[i] / (cnt[i] + 1)) * (cnt[i] + 2);
首先要清楚, p r i m e [ i ] prime[i] prime[i] 一定是 i i i 的最小素因数
之后可以用下面的例子来进行理解:
假设 i = p a i = p^a i=pa,那么 f [ i ] = a + 1 f[i] = a + 1 f[i]=a+1, c n t [ i ] = a cnt[i] = a cnt[i]=a。当 p r i m e [ j ] = p prime[j] = p prime[j]=p 时, p r i m e [ j ] ∗ i = p ( a + 1 ) prime[j] * i = p^{(a+1)} prime[j]∗i=p(a+1),其约数个数 f [ p r i m e [ j ] ∗ i ] f[prime[j] * i] f[prime[j]∗i] 应该是 a + 2 a + 2 a+2。这正好等于 f [ i ] / ( c n t [ i ] + 1 ) ∗ ( c n t [ i ] + 2 ) f[i] / (cnt[i] + 1) * (cnt[i] + 2) f[i]/(cnt[i]+1)∗(cnt[i]+2)。