3 * Compute the %$n$%th Fibonacci number
5 * (c) 2013 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
33 /*----- About the algorithm -----------------------------------------------*
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.
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
51 * So we can compute Fibonacci numbers by exponentiating in this ring.
52 * Squaring and multiplication work like this.
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)$%
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)$%.
61 /*----- Exponentiation machinery ------------------------------------------*/
63 /* --- @struct fib@ --- *
65 * A simple structure tracking polynomial coefficients.
69 int n
; /* Exponent for this entry */
70 mp
*a
, *b
; /* Coefficients: %$a t + b$% */
73 #define MAX 100 /* Saturation bounds for exponent */
78 * Clamp @n@ within the upper and lower bounds.
81 #define CLAMP(n) do { \
82 if (n > MAX) n = MAX; else if (n < MIN) n = MIN; \
85 /* --- Basic structure maintenance functions --- */
87 static void fib_init(struct fib
*f
)
88 { f
->a
= f
->b
= MP_NEW
; }
90 static void fib_drop(struct fib
*f
)
91 { if (f
->a
) MP_DROP(f
->a
); if (f
->b
) MP_DROP(f
->b
); }
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
); }
96 /* --- @fib_sqr@ --- *
98 * Arguments: @struct fib *d@ = destination structure
99 * @struct fib *x@ = operand
103 * Use: Set @d@ to the square of @x@.
106 static void fib_sqr(struct fib
*d
, struct fib
*x
)
110 /* --- Special case: if @x@ is the identity then just copy --- */
113 if (d
!= x
) { fib_drop(d
); fib_copy(d
, x
); }
117 /* --- Compute the result --- */
119 aa
= mp_sqr(MP_NEW
, x
->a
); /* %$a^2$% */
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$% */
125 t
= mp_sqr(d
->b
, x
->b
); /* %$t = b^2$% */
126 d
->b
= mp_add(t
, t
, aa
); /* %$b' = a^2 + b^2$% */
128 /* --- Sort out the exponent on the result --- */
130 d
->n
= 2*x
->n
; CLAMP(d
->n
);
137 /* --- @fib_mul@ --- *
139 * Arguments: @struct fib *d@ = destination structure
140 * @struct fib *x, *y@ = operands
144 * Use: Set @d@ to the product of @x@ and @y@.
147 static void fib_mul(struct fib
*d
, struct fib
*x
, struct fib
*y
)
151 /* --- Lots of special cases for low exponents --- */
155 if (d
!= x
) { fib_drop(d
); fib_copy(d
, x
); }
157 } else if (x
->n
== 0) { x
= y
; goto copy_x
; }
158 else if (y
->n
== -1) {
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
);
164 } else if (y
->n
== +1) {
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
);
170 } else if (x
->n
== -1) { x
= y
; goto dec_x
; }
171 else if (x
->n
== +1) { x
= y
; goto inc_x
; }
173 /* --- Compute the result --- */
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$% */
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$% */
184 /* --- Sort out the exponent on the result --- */
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
); }
192 MP_DROP(bd
); MP_DROP(t
); MP_DROP(u
);
195 /* --- Exponentiation --- */
197 #define EXP_TYPE struct fib
198 #define EXP_COPY(d, x) fib_copy(&d, &x)
199 #define EXP_DROP(x) fib_drop(&x)
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)
209 /*----- Main code ---------------------------------------------------------*/
211 /* --- @mp_fibonacci@ --- *
213 * Arguments: @long n@ = index desired (may be negative)
215 * Returns: The %$n$%th Fibonacci number.
218 mp
*mp_fibonacci(long n
)
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
);
228 EXP_WINDOW(d
, g
, nn
);
230 MP_DROP(nn
); fib_drop(&g
); MP_DROP(d
.b
);
234 /*----- Test rig ----------------------------------------------------------*/
238 #include <mLib/testrig.h>
240 static int vfib(dstr
*v
)
242 long x
= *(long *)v
[0].buf
;
243 mp
*fx
= *(mp
**)v
[1].buf
;
244 mp
*y
= mp_fibonacci(x
);
247 fprintf(stderr
, "fibonacci failed\n");
248 MP_FPRINTF(stderr
, (stderr
, "fib(%ld) = ", x
), fx
);
249 MP_EPRINT("result", y
);
254 assert(mparena_count(MPARENA_GLOBAL
) == 0);
258 static test_chunk tests
[] = {
259 { "fibonacci", vfib
, { &type_long
, &type_mp
, 0 } },
263 int main(int argc
, char *argv
[])
265 test_run(argc
, argv
, tests
, SRCDIR
"/t/mp");
271 /*----- That's all, folks -------------------------------------------------*/