Unit Fraction Partition

さいしょに通分して，あとは総当りで求めました．
途中で continue してるのは，元の分数より単位分数の和のほうが大きくなったり，単位分数の数が n より多くなっちゃうときはもう見つかるわけないからです．

#include <stdio.h>

int main(void)
{
	//	ファイル出力準備
	FILE	*pfI;	//	入力ファイルポインタ
	FILE	*pfO;	//	出力ファイルポインタ
	pfI = fopen("c.in", "r");
	if(pfI == NULL)
	{
		printf("ERR: input file open\n");
		return -1;
	}
	pfO = fopen("c.out","w");
	if(pfO == NULL)
	{
		printf("ERR: output file open\n");
		fclose(pfI);
		return -2;
	}

	//
	//	入力
	//
	int p, q, a, n;
	while(1)
	{
		fscanf(pfI, "%d %d %d %d ", &p, &q, &a, &n);
		//printf("%d %d %d %d\n", p, q, a, n);
		if(p == 0 && q == 0 && a == 0 && n == 0)
			break;

		//
		//	処理
		//
		int c = 0;	//	単位分数の数

		//	分数が通分できたらする
		if(q % p == 0)
		{
			q /= p;
			p /= p;
		}
		//printf("%d, %d\n", p, q);

		for(int i1 = 1; i1 <= a; i1++)
		{
			if(p == 1 && q == i1)
			{
				c++;
				//printf("1/%d\n", i1);
			}
			if((double)p/(double)q < 1.0/(double)i1)
				continue;
			if(n < 2)
				continue;
			for(int i2 = i1; i1*i2 <= a; i2++)
			{
				if(p * (i1*i2) == q*i2 + q*i1)
				{
					c++;
					//printf("1/%d + 1/%d\n", i1, i2);
				}
				if((double)p/(double)q < 1.0/(double)i1 + 1.0/(double)i2)
					continue;
				if(n < 3)
					continue;
				for(int i3 = i2; i1*i2*i3 <= a; i3++)
				{
					if(p * (i1*i2*i3) == q*i2*i3 + q*i1*i3 + q*i1*i2)
					{
						c++;
						//printf("1/%d + 1/%d + 1/%d\n", i1, i2, i3);	
					}
					if(n < 4)
					if(n < 4)
						continue;
					for(int i4 = i3; i1*i2*i3*i4 <= a; i4++)
					{
						if(p * (i1*i2*i3*i4) == q*i2*i3*i4 + q*i1*i3*i4 + q*i1*i2*i4 + q*i1*i2*i3)
						{
							c++;
							//printf("1/%d + 1/%d + 1/%d + 1/%d\n", i1, i2, i3, i4);	
						}
						if((double)p/(double)q < 1.0/(double)i1 + 1.0/(double)i2 + 1.0/(double)i3 + 1.0/(double)i4)
							continue;
						if(n < 5)
							continue;
						for(int i5 = i4; i1*i2*i3*i4*i5 <= a; i5++)
						{
							if(p * (i1*i2*i3*i4*i5) == q*i2*i3*i4*i5 + q*i1*i3*i4*i5 + q*i1*i2*i4*i5 + q*i1*i2*i3*i5 + q*i1*i2*i3*i4)
							{
								c++;
								//printf("1/%d + 1/%d + 1/%d + 1/%d + 1/%d\n", i1, i2, i3, i4, i5);	
							}
							if((double)p/(double)q < 1.0/(double)i1 + 1.0/(double)i2 + 1.0/(double)i3 + 1.0/(double)i4 + 1.0/(double)i5)
								continue;
							if(n < 6)
								continue;
							for(int i6 = i5; i1*i2*i3*i4*i5*i6 <= a; i6++)
							{
								if(p * (i1*i2*i3*i4*i5*i6) == q*i2*i3*i4*i5*i6 + q*i1*i3*i4*i5*i6 + q*i1*i2*i4*i5*i6 + q*i1*i2*i3*i5*i6 + q*i1*i2*i3*i4*i6 + q*i1*i2*i3*i4*i5)
								{
									c++;
									//printf("1/%d + 1/%d + 1/%d + 1/%d + 1/%d + 1/%d\n", i1, i2, i3, i4, i5, i6);	
								}
								if((double)p/(double)q < 1.0/(double)i1 + 1.0/(double)i2 + 1.0/(double)i3 + 1.0/(double)i4 + 1.0/(double)i5 + 1.0/(double)i6)
									continue;
								if(n < 7)
									continue;
								for(int i7 = i6; i1*i2*i3*i4*i5*i6*i7 <= a; i7++)
								{
									if(p * (i1*i2*i3*i4*i5*i6*i7) == q*i2*i3*i4*i5*i6*i7 + q*i1*i3*i4*i5*i6*i7 + q*i1*i2*i4*i5*i6*i7 + q*i1*i2*i3*i5*i6*i7 + q*i1*i2*i3*i4*i6*i7 + q*i1*i2*i3*i4*i5*i7 + q*i1*i2*i3*i4*i5*i6)
									{
										c++;
										//printf("1/%d + 1/%d + 1/%d + 1/%d + 1/%d + 1/%d + 1/%d\n", i1, i2, i3, i4, i5, i6, i7);	
									}
									if((double)p/(double)q < 1.0/(double)i1 + 1.0/(double)i2 + 1.0/(double)i3 + 1.0/(double)i4 + 1.0/(double)i5 + 1.0/(double)i6 + 1.0/(double)i7)
										continue;
								}
							}
						}
					}
				}
			}
		}

		//
		//	出力
		//
		fprintf(pfO, "%d\n", c);
	}

	//ファイルポインタを閉じる
	fclose(pfI);
	fclose(pfO);

	return 0;
}