Project Euler 78

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进行分割的情况数。看到这样的题目首先脑海中浮现的就是递归法。

一、递归

    令f(sum,max)=对整数sum分割中最大值为max的情况，显然有p(n)=f(n,1)+f(n,2)+...+f(n,n)

    而根据定义，不难得到f(s,m)=f(s-m,1)+f(s-m,2)+...+f(s-m,s-m)

    再写出递归边界f(s,m)=1(s=m时),f(s,m)=0,（s<m），就可以算出p(n)

    这个方法好处是简洁易懂，好想。然而时间复杂度十分高，虽然我不太会计算递归的复杂度，不过起码是指数级数的时间了。用这种方法来虽然能得到p（n），却几乎不可能得到结果。

二、动态规划

    前面的递归算法效率实在太低，不得不对此进行优化。经过分析发现，我们重复计算了许多值，不妨以空间换时间，做一个三角形的数组，以储存f(sum,max)，于是我们每次就只需要将前面的结果相加，就可以得到新的f(sum,max)。

    这种算法的复杂度是o(n^3)，也十分不理想。（我在自家计算机上跑了20分钟才得到前10000的p（n）mod1000000的值）然而还没跑到答案。

三、母函数法

    这是我曾经在高中数学竞赛中接触的方法。。当年学了几乎没用过（捂脸）没想到在这种场合又有了用武之地。

    考虑函数G(x)=a0+a1*x+a2*x^2+...+an*x^n，则称G(x)是数列an的母函数。

    对数列an=p（n），设母函数为G1(x)=a0+a1*x+a2*x^2+...+an*x^n+...。

    考虑G2(x)=(1+x+x^2+...+x^n)(1+x^2+x^4+...+x^n+...)(1+x^3+...+x^n+..)....

    第一个括号表示 整数n的分割中1的个数，若取x^2，则1在n的分割中出现两次。

    分析G2每一项所表达的意义，可以发现G2(x)中x^k的系数恰恰就是p(k)。

    于是通过构造多项式相乘，我们就能够得到p（n）。时间复杂度同样是o（n^3），只不过可能相比第二种方法系数比较小，所以会快上一些，然而同样不能跑出结果。（跑出结果可能要八个小时左右？。。。不现实）

四、数论与五边形数定理

    在屡次尝试失败之后，查了一些资料。发现数学家Euler早已对这个问题给过比较完整的解答。Euler 的解答分三个部分：

    一、构造母函数（像上面那样）并化简成1/[(1-x)*(1-x^2)*(1-x^3)*...]（化简可参见无穷级数）

    二、分析(1-x)*(1-x^2)*(1-x^3)*...（这个函数也叫欧拉函数）的系数。（结论：系数的通项为 k*(3k±1）/2）（这个东西证明起来十分复杂，自行百度五边形数定理）

    三、p（n）与欧拉函数的关系，通过分析n次项系数可得p(n)-p(n-1)-p(n-2)+p(n-5)+p(n-7)-p(n-12)-p(n-15)+…-…=0

    于是我们可以用前几项的p（n）来得到新的p（n），这是效率相当高的方法。时间复杂度o(n^3/2)。答案：55374。贴上代码：

#include <stdio.h> 
#include <malloc.h>
#include <string.h>

const int wu_size=10000;
const int max=1000000;

int main() {
	int *guangwu,*ans,flag=1;
	ans=(int *)malloc(max*sizeof(int));
	memset(ans,0,max*sizeof(int));
	guangwu=(int *)malloc(wu_size*sizeof(int));
	guangwu[0]=0;
	for (int i=1,k=1;i<wu_size;i+=2) {
		guangwu[i]=k*(3*k-1)/2;
		k++;
	}
	for (int i=2,k=1;i<wu_size;i+=2) {
		guangwu[i]=k*(3*k+1)/2;
		k++;
	}
	
	ans[0] = 1;
	for (int i=1;i<max;i++) {
		for (int j=1;guangwu[j]<=i;j++) {
			if (j%4==1||j%4==2) flag=1;
			else flag=-1;
			ans[i] += flag*ans[i-guangwu[j]];
			if (ans[i]>=1000000&&ans[i]>0) ans[i] %= 1000000;
			if (ans[i]<0) ans[i]+=1000000;
		}
		//if (i<=10) printf("%d\n",ans[i]); //for test
		if (ans[i]%1000000==0) {
			printf("%d\n",i);
			break;
		}
	}
}