埃及分数c语言设计,埃及分数（四）

最优分解

源程序

下面就是求埃及分数最优分解的 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。