题目78：考察硬币划分的方式数量

欧拉计划

Coin partitions

Let p(n) represent the number of different ways in which n coins can be separated into piles. For example, five coins can be separated into piles in exactly seven different ways, so p(5)=7.

OOOOO

OOOO   O

OOO   OO

OOO   O   O

OO   OO   O

OO   O   O   O

O   O   O   O   O

Find the least value of n for which p(n) is divisible by one million.

题目：

用 p(n) 表示 n 个硬币可以被划分成堆的方式。例如，5 个硬币可以被用 7 种方式划分成堆，所以 p(5)=7。

OOOOO

OOOO   O

OOO   OO

OOO   O   O

OO   OO   O

OO   O   O   O

O   O   O   O   O

找到 p(n) 可以被一百万整除的最小的 n 。

jerryxjr1220

这题是参考了维基百科上的计算公式：
According Wikipedia's information:
p(k) = p(k-(3*n*n-n)/2) + p(k-(3*n*n+n)/2) - p(k-(3*n*n+5*n+2)/2) - p(k-(3*n*n+7*n+4)/2) +... (n from 1 to ...) while p(0) = 1 and p(1) = 1
有了迭代公式，写程序相对还是比较容易的：

1. pdic = {0:1,1:1}
   
2. def p(k):
   
3.         global pdic
   
4.         s = 0
   
5.         if k == 0 or k == 1:
   
6.                 return 1
   
7.         for n in range(1,int(k**0.5+1),2):
   
8.                 a = k-(3*n*n-n)/2
   
9.                 b = k-(3*n*n+n)/2
   
10.                 c = k-(3*n*n+5*n+2)/2
    
11.                 d = k-(3*n*n+7*n+4)/2
    
12.                 if a >= 0:
    
13.                         s += pdic[a]
    
14.                 else:
    
15.                         break
    
16.                 if b >= 0:
    
17.                         s += pdic[b]
    
18.                 else:
    
19.                         break
    
20.                 if c >= 0:
    
21.                         s -= pdic[c]
    
22.                 else:
    
23.                         break
    
24.                 if d >= 0:
    
25.                         s -= pdic[d]
    
26.                 else:
    
27.                         break
    
28.         pdic[k] = s
    
29.         return s
    
30. 
    
31. for i in range(100000):
    
32.         if p(i) % 1000000 == 0:
    
33.                 print (i)
    
34.                 break

复制代码

大约8秒钟可以出结果：
55374

guosl

这题也可以通过组合数学中的母函数展开来做。只是时间比较长。但可以用GPU来编程计算。在GTX1060上2秒多点可以得出结果。

1. /*
   
2. 用母函数展开来做
   
3. 答案：55374
   
4. 耗时：2.07287秒（GTX1060显卡）
   
5. */
   
6. #include <iostream>
   
7. #include <algorithm>
   
8. #include <cuda_runtime.h>
   
9. #include <device_launch_parameters.h>
   
10. #include <helper_cuda.h>
    
11. using namespace std;
    
12. 
    
13. const int N = 64000;
    
14. const int nBlocksSize = 256;
    
15. 
    
16. int *a, *b, c[N];
    
17. bool *bFind;
    
18. 
    
19. __global__ void Polynomial_Mul(int *p1, int *p2, int k, bool *bFlag)
    
20. {
    
21.   int i = blockIdx.x * nBlocksSize + threadIdx.x;
    
22.   if (i >= N)
    
23.     return;
    
24.   p2[i] = 0;
    
25.   for (int j = 0; j <= i; j += k)
    
26.   {
    
27.     p2[i] += p1[i - j];
    
28.     p2[i] %= 1000000;
    
29.   }
    
30.   if (i == k)
    
31.   {
    
32.     if (p2[k] == 0)
    
33.       *bFlag = true;
    
34.   }
    
35. }
    
36. 
    
37. void find_zero(int *p1, int *p2, int k)
    
38. {
    
39.   Polynomial_Mul <<<N / nBlocksSize, nBlocksSize >>> (p1, p2, k, bFind);
    
40. }
    
41. 
    
42. int main(void)
    
43. {
    
44.   cudaDeviceProp deviceProp;
    
45.   int devID = gpuGetMaxGflopsDeviceId();
    
46.   cudaSetDevice(devID);//设置计算卡
    
47.   cudaGetDeviceProperties(&deviceProp, devID);
    
48.   printf("GPU Device %d: "%s" with compute capability %d.%d\n\n",
    
49.     devID, deviceProp.name, deviceProp.major, deviceProp.minor);
    
50.   //在显卡上分配内存
    
51.   cudaMalloc((void**)&a, N * sizeof(int));
    
52.   cudaMalloc((void**)&b, N * sizeof(int));
    
53.   cudaMalloc((void**)&bFind, sizeof(int));
    
54.   //设置计算的初始向量
    
55.   for (int i = 0; i < N; ++i)
    
56.     c[i] = 1;
    
57.   cudaMemcpy(a, c, N * sizeof(int), cudaMemcpyHostToDevice);
    
58.   int *p1 = a, *p2 = b;
    
59.   for (int k = 2; k < N; ++k) //开始逐个展开母函数
    
60.   {
    
61.     bool bFlag = false;//设置返回值的初值
    
62.     cudaMemcpy(bFind, &bFlag, sizeof(bool), cudaMemcpyHostToDevice);
    
63.     find_zero(p1, p2, k);//计算第k次展开
    
64.     swap(p1, p2);
    
65.     cudaGetLastError();
    
66.     cudaDeviceSynchronize();//等待计算完成
    
67.     cudaMemcpy(&bFlag, bFind, sizeof(bool), cudaMemcpyDeviceToHost);//获取本次计算的结果
    
68.     if (bFlag)//找到结果
    
69.     {
    
70.       cout << k << endl;
    
71.       break;
    
72.     }
    
73.   }
    
74.   cudaFree(a);
    
75.   cudaFree(b);
    
76.   cudaFree(bFind);
    
77.   return 0;
    
78. }

复制代码

fc1735

1. from numba import jit
   
2. 
   
3. @jit
   
4. def solve():
   
5.   mod = 10 ** 6
   
6.   bound = 10 ** 5
   
7.   dp = [0] * bound
   
8.   dp[0] = 1
   
9.   for i in range(1, bound):
   
10.     for j in range(i, bound):
    
11.       dp[j] = (dp[j] + dp[j-i]) % mod
    
12.     if dp[i] == 0:
    
13.       print(i)
    
14.       break
    
15. 
    
16. solve()

复制代码

https://github.com/devinizz/project_euler/blob/master/page02/78.py

永恒的蓝色梦想

和76好像有点关系