Annotation of OpenXM_contrib/gmp/mpz/fib_ui.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpz_fib_ui(result, n) -- Set RESULT to the Nth Fibonacci number.
2:
3: Copyright (C) 1998, 1999, 2000 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 "gmp.h"
23: #include "gmp-impl.h"
24:
25: /* This is fast, but could be made somewhat faster and neater.
26: The timing is somewhat fluctuating for even/odd sizes because
27: of the extra hair used to save variables and operations. Here
28: are a few things one might want to address:
29: 1. Avoid using 4 intermediate variables in mpz_fib_bigcase.
30: 2. Call mpn functions directly. Straightforward for these functions.
31: 3. Merge the three functions into one.
32:
33: Said by Kevin:
34: Consider using the Lucas numbers L[n] as an auxiliary sequence, making
35: it possible to do the "doubling" operation in mpz_fib_bigcase with two
36: squares rather than two multiplies. The formulas are a little more
37: complicated, something like the following (untested).
38:
39: F[2n] = ((F[n]+L[n])^2 - 6*F[n]^2 - 4*(-1)^n) / 2
40: L[2n] = 5*F[n]^2 + 2*(-1)^n
41:
42: F[2n+1] = (F[2n] + L[2n]) / 2
43: L[2n+1] = (5*F[2n] + L[2n]) / 2
44:
45: The Lucas number that comes for free here could even be returned.
46:
47: Maybe there's formulas with two squares using just F[n], but I don't
48: know of any.
49: */
50:
51: /* Determine the needed storage for Fib(n). */
52: #define FIB_SIZE(n) (((mp_size_t) ((n)*0.695)) / BITS_PER_MP_LIMB + 2)
53:
54: static void mpz_fib_bigcase _PROTO ((mpz_t, mpz_t, unsigned long int));
55: static void mpz_fib_basecase _PROTO ((mpz_t, mpz_t, unsigned long int));
56:
57:
58: #ifndef FIB_THRESHOLD
59: #define FIB_THRESHOLD 60
60: #endif
61:
62: void
63: #if __STDC__
64: mpz_fib_ui (mpz_t r, unsigned long int n)
65: #else
66: mpz_fib_ui (r, n)
67: mpz_t r;
68: unsigned long int n;
69: #endif
70: {
71: if (n == 0)
72: mpz_set_ui (r, 0);
73: else
74: {
75: mpz_t t1;
76: mpz_init (t1);
77: if (n < FIB_THRESHOLD)
78: mpz_fib_basecase (t1, r, n);
79: else
80: mpz_fib_bigcase (t1, r, n);
81: mpz_clear (t1);
82: }
83: }
84:
85: static void
86: #if __STDC__
87: mpz_fib_basecase (mpz_t t1, mpz_t t2, unsigned long int n)
88: #else
89: mpz_fib_basecase (t1, t2, n)
90: mpz_t t1;
91: mpz_t t2;
92: unsigned long int n;
93: #endif
94: {
95: unsigned long int m, i;
96:
97: mpz_set_ui (t1, 0);
98: mpz_set_ui (t2, 1);
99: m = n/2;
100: for (i = 0; i < m; i++)
101: {
102: mpz_add (t1, t1, t2);
103: mpz_add (t2, t1, t2);
104: }
105: if ((n & 1) == 0)
106: {
107: mpz_sub (t1, t2, t1);
108: mpz_sub (t2, t2, t1); /* trick: recover t1 value just overwritten */
109: }
110: }
111:
112: static void
113: #if __STDC__
114: mpz_fib_bigcase (mpz_t t1, mpz_t t2, unsigned long int n)
115: #else
116: mpz_fib_bigcase (t1, t2, n)
117: mpz_t t1;
118: mpz_t t2;
119: unsigned long int n;
120: #endif
121: {
122: unsigned long int n2;
123: int ni, i;
124: mpz_t x1, x2, u1, u2;
125:
126: ni = 0;
127: for (n2 = n; n2 >= FIB_THRESHOLD; n2 /= 2)
128: ni++;
129:
130: mpz_fib_basecase (t1, t2, n2);
131:
132: mpz_init (x1);
133: mpz_init (x2);
134: mpz_init (u1);
135: mpz_init (u2);
136:
137: for (i = ni - 1; i >= 0; i--)
138: {
139: mpz_mul_2exp (x1, t1, 1);
140: mpz_mul_2exp (x2, t2, 1);
141:
142: mpz_add (x1, x1, t2);
143: mpz_sub (x2, x2, t1);
144:
145: mpz_mul (u1, t2, x1);
146: mpz_mul (u2, t1, x2);
147:
148: if (((n >> i) & 1) == 0)
149: {
150: mpz_sub (t1, u1, u2);
151: mpz_set (t2, u1);
152: }
153: else
154: {
155: mpz_set (t1, u1);
156: mpz_mul_2exp (t2, u1, 1);
157: mpz_sub (t2, t2, u2);
158: }
159: }
160:
161: mpz_clear (x1);
162: mpz_clear (x2);
163: mpz_clear (u1);
164: mpz_clear (u2);
165: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>