#
source:
trunk/third/gmp/mpz/fib_ui.c
@
18191

Revision 18191, 4.4 KB checked in by ghudson, 22 years ago (diff) |
---|

Line | |
---|---|

1 | /* mpz_fib_ui -- calculate Fibonacci numbers. |

2 | |

3 | Copyright 2000, 2001 Free Software Foundation, Inc. |

4 | |

5 | This file is part of the GNU MP Library. |

6 | |

7 | The GNU MP Library is free software; you can redistribute it and/or modify |

8 | it under the terms of the GNU Lesser General Public License as published by |

9 | the Free Software Foundation; either version 2.1 of the License, or (at your |

10 | option) any later version. |

11 | |

12 | The GNU MP Library is distributed in the hope that it will be useful, but |

13 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |

14 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |

15 | License for more details. |

16 | |

17 | You should have received a copy of the GNU Lesser General Public License |

18 | along with the GNU MP Library; see the file COPYING.LIB. If not, write to |

19 | the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |

20 | MA 02111-1307, USA. */ |

21 | |

22 | #include <stdio.h> |

23 | #include "gmp.h" |

24 | #include "gmp-impl.h" |

25 | #include "longlong.h" |

26 | |

27 | |

28 | /* change to "#define TRACE(x) x" to get some traces */ |

29 | #define TRACE(x) |

30 | |

31 | |

32 | /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low |

33 | limb because the result F[2k+1] is an F[4m+3] and such numbers are always |

34 | == 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7. (This is |

35 | the same as in mpn_fib2_ui.) |

36 | |

37 | In the F[2k+1] for k even, the +2 won't give a carry out of the low limb |

38 | in normal circumstances. This is an F[4m+1] and we claim that F[3*2^b+1] |

39 | == 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence |

40 | if n < 2^BITS_PER_MP_LIMB then F[n] cannot have a low limb of 0 or 1. No |

41 | proof for this claim, but it's been verified up to b==32 and has such a |

42 | nice pattern it must be true :-). Of interest is that F[3*2^b] == 0 mod |

43 | 2^(b+1) seems to hold too. |

44 | |

45 | When n >= 2^BITS_PER_MP_LIMB, which can arise in test setups with a small |

46 | limb, then the low limb of F[4m+1] can certainly be 1, and an mpn_add_1 |

47 | must be used. */ |

48 | |

49 | void |

50 | mpz_fib_ui (mpz_ptr fn, unsigned long n) |

51 | { |

52 | mp_ptr fp, xp, yp; |

53 | mp_size_t size, xalloc; |

54 | unsigned long n2; |

55 | mp_limb_t c, c2; |

56 | TMP_DECL (marker); |

57 | |

58 | if (n <= FIB_TABLE_LIMIT) |

59 | { |

60 | PTR(fn)[0] = FIB_TABLE (n); |

61 | SIZ(fn) = (n != 0); /* F[0]==0, others are !=0 */ |

62 | return; |

63 | } |

64 | |

65 | n2 = n/2; |

66 | xalloc = MPN_FIB2_SIZE (n2) + 1; |

67 | MPZ_REALLOC (fn, 2*xalloc+1); |

68 | fp = PTR (fn); |

69 | |

70 | TMP_MARK (marker); |

71 | TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc); |

72 | size = mpn_fib2_ui (xp, yp, n2); |

73 | |

74 | TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n", |

75 | n >> 1, size, n&1); |

76 | mpn_trace ("xp", xp, size); |

77 | mpn_trace ("yp", yp, size)); |

78 | |

79 | if (n & 1) |

80 | { |

81 | /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k */ |

82 | mp_size_t xsize, ysize; |

83 | |

84 | #if HAVE_NATIVE_mpn_addsub_n |

85 | xp[size] = mpn_lshift (xp, xp, size, 1); |

86 | yp[size] = 0; |

87 | ASSERT_NOCARRY (mpn_addsub_n (xp, yp, xp, yp, size+1)); |

88 | xsize = size + (xp[size] != 0); |

89 | ysize = size + (yp[size] != 0); |

90 | #else |

91 | c2 = mpn_lshift (fp, xp, size, 1); |

92 | c = c2 + mpn_add_n (xp, fp, yp, size); |

93 | xp[size] = c; |

94 | xsize = size + (c != 0); |

95 | c2 -= mpn_sub_n (yp, fp, yp, size); |

96 | yp[size] = c2; |

97 | ASSERT (c2 <= 1); |

98 | ysize = size + c2; |

99 | #endif |

100 | |

101 | size = xsize + ysize; |

102 | c = mpn_mul (fp, xp, xsize, yp, ysize); |

103 | |

104 | #if BITS_PER_MP_LIMB >= BITS_PER_ULONG |

105 | /* no overflow, see comments above */ |

106 | ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= MP_LIMB_T_MAX-2); |

107 | fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2)); |

108 | #else |

109 | /* this code only for testing with small limbs, limb<ulong is unusual */ |

110 | if (n & 2) |

111 | { |

112 | ASSERT (fp[0] >= 2); |

113 | fp[0] -= 2; |

114 | } |

115 | else |

116 | { |

117 | ASSERT (c != MP_LIMB_T_MAX); /* because it's the high of a mul */ |

118 | c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2)); |

119 | fp[size-1] = c; |

120 | } |

121 | #endif |

122 | } |

123 | else |

124 | { |

125 | /* F[2k] = F[k]*(F[k]+2F[k-1]) */ |

126 | |

127 | mp_size_t xsize, ysize; |

128 | c = mpn_lshift (yp, yp, size, 1); |

129 | c += mpn_add_n (yp, yp, xp, size); |

130 | yp[size] = c; |

131 | xsize = size; |

132 | ysize = size + (c != 0); |

133 | size += ysize; |

134 | c = mpn_mul (fp, yp, ysize, xp, xsize); |

135 | } |

136 | |

137 | /* one or two high zeros */ |

138 | size -= (c == 0); |

139 | size -= (fp[size-1] == 0); |

140 | SIZ(fn) = size; |

141 | |

142 | TRACE (printf ("done special, size=%ld\n", size); |

143 | mpn_trace ("fp ", fp, size)); |

144 | |

145 | TMP_FREE (marker); |

146 | } |

**Note:**See TracBrowser for help on using the repository browser.