最优分解
源程序
下面就是求埃及分数最优分解的 C 语言源程序 EgyptianFraction.c：
1: #include
2: #include
3: #include
4:
5: const int SIZE = 64;
6: static int wp = -1, bp = -1;
7: static mpz_t z, w, work[SIZE], best[SIZE];
8: static mpq_t q1, q2;
9:
10: void fail(int code) { printf("***FAIL(%d)***\n", code); exit(code); }
11: void push(const mpz_t v) { if (wp > SIZE-2) fail(1); mpz_set(work[++wp], v); }
12: void pop() { if (wp < 0) fail(2); --wp; }
13: void copy() { bp = wp; for (int i = 0; i <= wp; i++) mpz_set(best[i], work[i]); }
14:
15: void sum(mpz_t sum, const mpz_t stack[], int sp)
16: {
17: mpz_set_si(sum, 0);
18: for (int i = 0; i <= sp; i++) mpz_add(sum, sum, stack[i]);
19: }
20:
21: void search(int depth, const mpq_t q, const mpz_t start)
22: {
23: if (depth == 0) fail(3);
24: mpz_t n; mpz_init(n); mpq_t q3; mpq_init(q3);
25: mpz_cdiv_q(n, mpq_denref(q), mpq_numref(q));
26: if (mpz_cmp(start, n) > 0) mpz_set(n, start);
27: while (1) {
28: mpq_set_den(q1, n);
29: int cmp = mpq_cmp(q, q1);
30: if (cmp > 0) {
31: mpq_set_si(q2, depth, 1); mpq_set_den(q2, n); mpq_canonicalize(q2);
32: if (mpq_cmp(q, q2) >= 0) break;
33: push(n);
34: mpz_add_ui(n, n, 1);
35: mpq_sub(q3, q, q1);
36: search(depth - 1, q3, n);
37: pop();
38: } else if (cmp == 0) {
39: push(n);
40: if (bp == -1) copy();
41: else {
42: cmp = mpz_cmp(best[bp], work[wp]);
43: if (cmp > 0) copy();
44: else if (cmp == 0) {
45: sum(z, best, bp); sum(w, work, wp);
46: if (mpz_cmp(z, w) > 0) copy();
47: }
48: }
49: pop();
50: break;
51: } else fail(4);
52: }
53: mpz_clear(n); mpq_clear(q3);
54: }
55:
56: int main(int argc, char *argv[])
57: {
58: if (argc != 3) {
59: puts("Usage: EgyptianFraction Numerator Denominator");
60: return 1;
61: }
62: mpz_t start; mpq_t q;
63: mpz_inits(start, z, w, 0);
64: mpz_set_str(z, argv[1], 10);
65: mpz_set_str(w, argv[2], 10);
66: mpq_inits(q, q1, q2, 0);
67: mpq_set_num(q, z);
68: mpq_set_den(q, w);
69: mpq_canonicalize(q);
70: if (mpq_sgn(q) <= 0) fail(5);
71: for (int i = 0; i < SIZE; i++) mpz_inits(work[i], best[i], 0);
72: gmp_printf("%Qd:", q); fflush(stdout);
73: mpq_set_si(q1, 1, 1);
74: mpz_set_si(start, 2);
75: for (int items = 1; bp == -1; items++) {
76: bp = -1;
77: search(items, q, start);
78: }
79: for (int i = 0; i <= bp; i++) gmp_printf(" %Zd", best[i]);
80: printf("\n");
81: for (int i = 0; i < SIZE; i++) mpz_clears(work[i], best[i], 0);
82: mpz_clears(start, z, w, 0);
83: mpq_clears(q, q1, q2, 0);
84: return 0;
85: }
编译和运行
编译：
$ clang -o EgyptianFraction EgyptianFraction.c -lgmp
运行：
$ ./EgyptianFraction
Usage: EgyptianFraction Numerator Denominator
$ ./EgyptianFraction 101 414
101/414: 6 18 46
$ time ./EgyptianFraction 732 733
732/733: 2 3 7 45 7330 20524 26388
real 7m8.216s
user 6m56.107s
sys 0m0.137s
源程序分析
这个源程序基本原理和上一篇“埃及分数(三)”中的 EgyptianFixItem.c 是一样的。 不同之处在于：
main 函数中，第 75 至 78 行的循环，从假设展开式只有一项开始，逐渐增加展开式的项数，调用 search 函数进行深度优先搜索。
search 函数中，第 38 至 50 行找到一个展开式后，将其和最优分解 best 比较，如果更优则替换之。
注意，不用比较 work 和 best 的项数是否相同，此时 wp 必定等于 bp。
