Commit | Line | Data |
---|---|---|
ab916894 MW |
1 | /* -*-c-*- |
2 | * | |
3 | * Compute the %$n$%th Fibonacci number | |
4 | * | |
5 | * (c) 2013 Straylight/Edgeware | |
6 | */ | |
7 | ||
8 | /*----- Licensing notice --------------------------------------------------* | |
9 | * | |
10 | * This file is part of Catacomb. | |
11 | * | |
12 | * Catacomb is free software; you can redistribute it and/or modify | |
13 | * it under the terms of the GNU Library General Public License as | |
14 | * published by the Free Software Foundation; either version 2 of the | |
15 | * License, or (at your option) any later version. | |
16 | * | |
17 | * Catacomb is distributed in the hope that it will be useful, | |
18 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
19 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
20 | * GNU Library General Public License for more details. | |
21 | * | |
22 | * You should have received a copy of the GNU Library General Public | |
23 | * License along with Catacomb; if not, write to the Free | |
24 | * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, | |
25 | * MA 02111-1307, USA. | |
26 | */ | |
27 | ||
28 | /*----- Header files ------------------------------------------------------*/ | |
29 | ||
30 | #include "mp.h" | |
31 | #include "mpint.h" | |
32 | ||
33 | /*----- About the algorithm -----------------------------------------------* | |
34 | * | |
35 | * Define %$F_0 = 0$% and %$F_1 = 1$%, and %$F_{k+1} = F_k + F_{k-1}$% for | |
36 | * all %$k$%. (This defines %$F_k$% for negative %$k$% too, by | |
37 | * %$F_{k-1} = F_{k+1} - F_k$%; in particular, %$F_{-1} = 1$% and | |
38 | * %$F_{-2} = -1$%.) We say that %$F_k$% is the %$k$%th Fibonacci number. | |
39 | * | |
40 | * We work in the ring %$\ZZ[t]/(t^2 - t -1)$%. Every residue class in this | |
41 | * ring contains a unique representative with degree at most 1. I claim that | |
42 | * %$t^k = F_k t + F_{k-1}$% for all %$k$%. Certainly %$t = F_1 t + F_0$%. | |
43 | * Note that %$t (F_{-1} t + F_{-2}) = t (t - 1) = t^2 - t = 1$%, so the | |
44 | * claim holds for %$k = -1$%. Suppose, inductively, that it holds for | |
45 | * %$k$%; then %$t^{k+1} = t \cdot t^k = F_k t^2 + F_{k-1} t = {}$% | |
46 | * %$(F_k + F_{k-1}) t + F_k = F_{k+1} t + F_k$%; and %$t^{k-1} = {}$% | |
47 | * %$t^{-1} t^k = (t - 1) t^k = t^{k+1} - t^k = {}$% | |
48 | * %$(F_{k+1} - F_k) t + (F_k - F_{k-1}) = F_{k-1} t + F_{k-2}$%, proving the | |
49 | * claim. | |
50 | * | |
51 | * So we can compute Fibonacci numbers by exponentiating in this ring. | |
52 | * Squaring and multiplication work like this. | |
53 | * | |
54 | * * Square: %$(a t + b)^2 = a^2 t^2 + 2 a b t + b^2 = {}$% | |
55 | * %$(a^2 + 2 a b) t + (a^2 + b^2)$% | |
56 | * | |
57 | * * Multiply: %$(a t + b)(c t + d) = a c t^2 + (a d + b c) t + b d = {}$% | |
58 | * %$(a c + a d + b c) t + (a c + b d)$%. | |
59 | */ | |
60 | ||
61 | /*----- Exponentiation machinery ------------------------------------------*/ | |
62 | ||
63 | /* --- @struct fib@ --- * | |
64 | * | |
65 | * A simple structure tracking polynomial coefficients. | |
66 | */ | |
67 | ||
68 | struct fib { | |
69 | int n; /* Exponent for this entry */ | |
70 | mp *a, *b; /* Coefficients: %$a t + b$% */ | |
71 | }; | |
72 | ||
73 | #define MAX 100 /* Saturation bounds for exponent */ | |
74 | #define MIN -100 | |
75 | ||
76 | /* --- @CLAMP@ --- * | |
77 | * | |
78 | * Clamp @n@ within the upper and lower bounds. | |
79 | */ | |
80 | ||
81 | #define CLAMP(n) do { \ | |
82 | if (n > MAX) n = MAX; else if (n < MIN) n = MIN; \ | |
83 | } while (0) | |
84 | ||
85 | /* --- Basic structure maintenance functions --- */ | |
86 | ||
87 | static void fib_init(struct fib *f) | |
88 | { f->a = f->b = MP_NEW; } | |
89 | ||
90 | static void fib_drop(struct fib *f) | |
91 | { if (f->a) MP_DROP(f->a); if (f->b) MP_DROP(f->b); } | |
92 | ||
93 | static void fib_copy(struct fib *d, struct fib *x) | |
94 | { d->n = x->n; d->a = MP_COPY(x->a); d->b = MP_COPY(x->b); } | |
95 | ||
96 | /* --- @fib_sqr@ --- * | |
97 | * | |
98 | * Arguments: @struct fib *d@ = destination structure | |
99 | * @struct fib *x@ = operand | |
100 | * | |
101 | * Returns: --- | |
102 | * | |
103 | * Use: Set @d@ to the square of @x@. | |
104 | */ | |
105 | ||
106 | static void fib_sqr(struct fib *d, struct fib *x) | |
107 | { | |
108 | mp *aa, *t; | |
109 | ||
110 | /* --- Special case: if @x@ is the identity then just copy --- */ | |
111 | ||
112 | if (!x->n) { | |
113 | if (d != x) { fib_drop(d); fib_copy(d, x); } | |
114 | return; | |
115 | } | |
116 | ||
117 | /* --- Compute the result --- */ | |
118 | ||
119 | aa = mp_sqr(MP_NEW, x->a); /* %$a^2$% */ | |
120 | ||
121 | t = mp_mul(d->a, x->a, x->b); /* %$t = a b$% */ | |
122 | t = mp_lsl(t, t, 1); /* %$t = 2 a b$% */ | |
123 | d->a = mp_add(t, t, aa); /* %$a' = a^2 + 2 a b$% */ | |
124 | ||
125 | t = mp_sqr(d->b, x->b); /* %$t = b^2$% */ | |
126 | d->b = mp_add(t, t, aa); /* %$b' = a^2 + b^2$% */ | |
127 | ||
128 | /* --- Sort out the exponent on the result --- */ | |
129 | ||
130 | d->n = 2*x->n; CLAMP(d->n); | |
131 | ||
132 | /* --- Done --- */ | |
133 | ||
134 | MP_DROP(aa); | |
135 | } | |
136 | ||
137 | /* --- @fib_mul@ --- * | |
138 | * | |
139 | * Arguments: @struct fib *d@ = destination structure | |
140 | * @struct fib *x, *y@ = operands | |
141 | * | |
142 | * Returns: --- | |
143 | * | |
144 | * Use: Set @d@ to the product of @x@ and @y@. | |
145 | */ | |
146 | ||
147 | static void fib_mul(struct fib *d, struct fib *x, struct fib *y) | |
148 | { | |
149 | mp *t, *u, *bd; | |
150 | ||
151 | /* --- Lots of special cases for low exponents --- */ | |
152 | ||
153 | if (y->n == 0) { | |
154 | copy_x: | |
155 | if (d != x) { fib_drop(d); fib_copy(d, x); } | |
156 | return; | |
157 | } else if (x->n == 0) { x = y; goto copy_x; } | |
158 | else if (y->n == -1) { | |
159 | dec_x: | |
160 | t = mp_sub(d->a, x->a, x->b); | |
161 | d->a = MP_COPY(x->b); if (d->b) MP_DROP(d->b); d->b = t; | |
162 | d->n = x->n - 1; CLAMP(d->n); | |
163 | return; | |
164 | } else if (y->n == +1) { | |
165 | inc_x: | |
166 | t = mp_add(d->b, x->a, x->b); | |
167 | d->b = MP_COPY(x->a); if (d->a) MP_DROP(d->a); d->a = t; | |
168 | d->n = x->n + 1; CLAMP(d->n); | |
169 | return; | |
170 | } else if (x->n == -1) { x = y; goto dec_x; } | |
171 | else if (x->n == +1) { x = y; goto inc_x; } | |
172 | ||
173 | /* --- Compute the result --- */ | |
174 | ||
175 | bd = mp_mul(MP_NEW, x->b, y->b); /* %$b d$% */ | |
176 | t = mp_add(MP_NEW, x->a, x->b); /* %$t = a + b$% */ | |
177 | u = mp_add(MP_NEW, y->a, y->b); /* %$u = c + d$% */ | |
178 | t = mp_mul(t, t, u); /* %$t = (a + b)(c + d)$% */ | |
179 | u = mp_mul(u, x->a, y->a); /* %$u = a c$% */ | |
180 | ||
181 | d->a = mp_sub(d->a, t, bd); /* %$a' = a c + a d + b c$% */ | |
182 | d->b = mp_add(d->b, u, bd); /* %$b' = a c + b d$% */ | |
183 | ||
184 | /* --- Sort out the exponent on the result --- */ | |
185 | ||
186 | if (x->n == MIN || x->n == MAX) d->n = x->n; | |
187 | else if (y->n == MIN || y->n == MAX) d->n = y->n; | |
188 | else { d->n = x->n + y->n; CLAMP(d->n); } | |
189 | ||
190 | /* --- Done --- */ | |
191 | ||
192 | MP_DROP(bd); MP_DROP(t); MP_DROP(u); | |
193 | } | |
194 | ||
195 | /* --- Exponentiation --- */ | |
196 | ||
197 | #define EXP_TYPE struct fib | |
198 | #define EXP_COPY(d, x) fib_copy(&d, &x) | |
199 | #define EXP_DROP(x) fib_drop(&x) | |
200 | #define EXP_FIX(d) | |
201 | ||
202 | #define EXP_SQR(x) fib_sqr(&x, &x) | |
203 | #define EXP_MUL(x, y) fib_mul(&x, &x, &y) | |
204 | #define EXP_SETSQR(d, x) fib_init(&d); fib_sqr(&d, &x) | |
205 | #define EXP_SETMUL(d, x, y) fib_init(&d); fib_mul(&d, &x, &y) | |
206 | ||
207 | #include "exp.h" | |
208 | ||
209 | /*----- Main code ---------------------------------------------------------*/ | |
210 | ||
211 | /* --- @mp_fibonacci@ --- * | |
212 | * | |
213 | * Arguments: @long n@ = index desired (may be negative) | |
214 | * | |
215 | * Returns: The %$n$%th Fibonacci number. | |
216 | */ | |
217 | ||
218 | mp *mp_fibonacci(long n) | |
219 | { | |
220 | struct fib d, g; | |
221 | mp *nn; | |
222 | ||
223 | d.n = 0; d.a = MP_ZERO; d.b = MP_ONE; | |
224 | if (n >= 0) { g.n = 1; g.a = MP_ONE; g.b = MP_ZERO; } | |
225 | else { g.n = -1; g.a = MP_ONE; g.b = MP_MONE; n = -n; } | |
226 | nn = mp_fromlong(MP_NEW, n); | |
227 | ||
228 | EXP_WINDOW(d, g, nn); | |
229 | ||
230 | MP_DROP(nn); fib_drop(&g); MP_DROP(d.b); | |
231 | return (d.a); | |
232 | } | |
233 | ||
234 | /*----- Test rig ----------------------------------------------------------*/ | |
235 | ||
236 | #ifdef TEST_RIG | |
237 | ||
238 | #include <mLib/testrig.h> | |
239 | ||
240 | static int vfib(dstr *v) | |
241 | { | |
242 | long x = *(long *)v[0].buf; | |
243 | mp *fx = *(mp **)v[1].buf; | |
244 | mp *y = mp_fibonacci(x); | |
245 | int ok = 1; | |
246 | if (!MP_EQ(fx, y)) { | |
247 | fprintf(stderr, "fibonacci failed\n"); | |
248 | MP_FPRINTF(stderr, (stderr, "fib(%ld) = ", x), fx); | |
249 | MP_EPRINT("result", y); | |
250 | ok = 0; | |
251 | } | |
252 | mp_drop(fx); | |
253 | mp_drop(y); | |
254 | assert(mparena_count(MPARENA_GLOBAL) == 0); | |
255 | return (ok); | |
256 | } | |
257 | ||
258 | static test_chunk tests[] = { | |
259 | { "fibonacci", vfib, { &type_long, &type_mp, 0 } }, | |
260 | { 0, 0, { 0 } } | |
261 | }; | |
262 | ||
263 | int main(int argc, char *argv[]) | |
264 | { | |
0f00dc4c | 265 | test_run(argc, argv, tests, SRCDIR "/t/mp"); |
ab916894 MW |
266 | return (0); |
267 | } | |
268 | ||
269 | #endif | |
270 | ||
271 | /*----- That's all, folks -------------------------------------------------*/ |