+/* -*-c-*-
+ *
+ * Compute the %$n$%th Fibonacci number
+ *
+ * (c) 2013 Straylight/Edgeware
+ */
+
+/*----- Licensing notice --------------------------------------------------*
+ *
+ * This file is part of Catacomb.
+ *
+ * Catacomb is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Library General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * Catacomb is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU Library General Public
+ * License along with Catacomb; if not, write to the Free
+ * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+ * MA 02111-1307, USA.
+ */
+
+/*----- Header files ------------------------------------------------------*/
+
+#include "mp.h"
+#include "mpint.h"
+
+/*----- About the algorithm -----------------------------------------------*
+ *
+ * Define %$F_0 = 0$% and %$F_1 = 1$%, and %$F_{k+1} = F_k + F_{k-1}$% for
+ * all %$k$%. (This defines %$F_k$% for negative %$k$% too, by
+ * %$F_{k-1} = F_{k+1} - F_k$%; in particular, %$F_{-1} = 1$% and
+ * %$F_{-2} = -1$%.) We say that %$F_k$% is the %$k$%th Fibonacci number.
+ *
+ * We work in the ring %$\ZZ[t]/(t^2 - t -1)$%. Every residue class in this
+ * ring contains a unique representative with degree at most 1. I claim that
+ * %$t^k = F_k t + F_{k-1}$% for all %$k$%. Certainly %$t = F_1 t + F_0$%.
+ * Note that %$t (F_{-1} t + F_{-2}) = t (t - 1) = t^2 - t = 1$%, so the
+ * claim holds for %$k = -1$%. Suppose, inductively, that it holds for
+ * %$k$%; then %$t^{k+1} = t \cdot t^k = F_k t^2 + F_{k-1} t = {}$%
+ * %$(F_k + F_{k-1}) t + F_k = F_{k+1} t + F_k$%; and %$t^{k-1} = {}$%
+ * %$t^{-1} t^k = (t - 1) t^k = t^{k+1} - t^k = {}$%
+ * %$(F_{k+1} - F_k) t + (F_k - F_{k-1}) = F_{k-1} t + F_{k-2}$%, proving the
+ * claim.
+ *
+ * So we can compute Fibonacci numbers by exponentiating in this ring.
+ * Squaring and multiplication work like this.
+ *
+ * * Square: %$(a t + b)^2 = a^2 t^2 + 2 a b t + b^2 = {}$%
+ * %$(a^2 + 2 a b) t + (a^2 + b^2)$%
+ *
+ * * Multiply: %$(a t + b)(c t + d) = a c t^2 + (a d + b c) t + b d = {}$%
+ * %$(a c + a d + b c) t + (a c + b d)$%.
+ */
+
+/*----- Exponentiation machinery ------------------------------------------*/
+
+/* --- @struct fib@ --- *
+ *
+ * A simple structure tracking polynomial coefficients.
+ */
+
+struct fib {
+ int n; /* Exponent for this entry */
+ mp *a, *b; /* Coefficients: %$a t + b$% */
+};
+
+#define MAX 100 /* Saturation bounds for exponent */
+#define MIN -100
+
+/* --- @CLAMP@ --- *
+ *
+ * Clamp @n@ within the upper and lower bounds.
+ */
+
+#define CLAMP(n) do { \
+ if (n > MAX) n = MAX; else if (n < MIN) n = MIN; \
+} while (0)
+
+/* --- Basic structure maintenance functions --- */
+
+static void fib_init(struct fib *f)
+ { f->a = f->b = MP_NEW; }
+
+static void fib_drop(struct fib *f)
+ { if (f->a) MP_DROP(f->a); if (f->b) MP_DROP(f->b); }
+
+static void fib_copy(struct fib *d, struct fib *x)
+ { d->n = x->n; d->a = MP_COPY(x->a); d->b = MP_COPY(x->b); }
+
+/* --- @fib_sqr@ --- *
+ *
+ * Arguments: @struct fib *d@ = destination structure
+ * @struct fib *x@ = operand
+ *
+ * Returns: ---
+ *
+ * Use: Set @d@ to the square of @x@.
+ */
+
+static void fib_sqr(struct fib *d, struct fib *x)
+{
+ mp *aa, *t;
+
+ /* --- Special case: if @x@ is the identity then just copy --- */
+
+ if (!x->n) {
+ if (d != x) { fib_drop(d); fib_copy(d, x); }
+ return;
+ }
+
+ /* --- Compute the result --- */
+
+ aa = mp_sqr(MP_NEW, x->a); /* %$a^2$% */
+
+ t = mp_mul(d->a, x->a, x->b); /* %$t = a b$% */
+ t = mp_lsl(t, t, 1); /* %$t = 2 a b$% */
+ d->a = mp_add(t, t, aa); /* %$a' = a^2 + 2 a b$% */
+
+ t = mp_sqr(d->b, x->b); /* %$t = b^2$% */
+ d->b = mp_add(t, t, aa); /* %$b' = a^2 + b^2$% */
+
+ /* --- Sort out the exponent on the result --- */
+
+ d->n = 2*x->n; CLAMP(d->n);
+
+ /* --- Done --- */
+
+ MP_DROP(aa);
+}
+
+/* --- @fib_mul@ --- *
+ *
+ * Arguments: @struct fib *d@ = destination structure
+ * @struct fib *x, *y@ = operands
+ *
+ * Returns: ---
+ *
+ * Use: Set @d@ to the product of @x@ and @y@.
+ */
+
+static void fib_mul(struct fib *d, struct fib *x, struct fib *y)
+{
+ mp *t, *u, *bd;
+
+ /* --- Lots of special cases for low exponents --- */
+
+ if (y->n == 0) {
+ copy_x:
+ if (d != x) { fib_drop(d); fib_copy(d, x); }
+ return;
+ } else if (x->n == 0) { x = y; goto copy_x; }
+ else if (y->n == -1) {
+ dec_x:
+ t = mp_sub(d->a, x->a, x->b);
+ d->a = MP_COPY(x->b); if (d->b) MP_DROP(d->b); d->b = t;
+ d->n = x->n - 1; CLAMP(d->n);
+ return;
+ } else if (y->n == +1) {
+ inc_x:
+ t = mp_add(d->b, x->a, x->b);
+ d->b = MP_COPY(x->a); if (d->a) MP_DROP(d->a); d->a = t;
+ d->n = x->n + 1; CLAMP(d->n);
+ return;
+ } else if (x->n == -1) { x = y; goto dec_x; }
+ else if (x->n == +1) { x = y; goto inc_x; }
+
+ /* --- Compute the result --- */
+
+ bd = mp_mul(MP_NEW, x->b, y->b); /* %$b d$% */
+ t = mp_add(MP_NEW, x->a, x->b); /* %$t = a + b$% */
+ u = mp_add(MP_NEW, y->a, y->b); /* %$u = c + d$% */
+ t = mp_mul(t, t, u); /* %$t = (a + b)(c + d)$% */
+ u = mp_mul(u, x->a, y->a); /* %$u = a c$% */
+
+ d->a = mp_sub(d->a, t, bd); /* %$a' = a c + a d + b c$% */
+ d->b = mp_add(d->b, u, bd); /* %$b' = a c + b d$% */
+
+ /* --- Sort out the exponent on the result --- */
+
+ if (x->n == MIN || x->n == MAX) d->n = x->n;
+ else if (y->n == MIN || y->n == MAX) d->n = y->n;
+ else { d->n = x->n + y->n; CLAMP(d->n); }
+
+ /* --- Done --- */
+
+ MP_DROP(bd); MP_DROP(t); MP_DROP(u);
+}
+
+/* --- Exponentiation --- */
+
+#define EXP_TYPE struct fib
+#define EXP_COPY(d, x) fib_copy(&d, &x)
+#define EXP_DROP(x) fib_drop(&x)
+#define EXP_FIX(d)
+
+#define EXP_SQR(x) fib_sqr(&x, &x)
+#define EXP_MUL(x, y) fib_mul(&x, &x, &y)
+#define EXP_SETSQR(d, x) fib_init(&d); fib_sqr(&d, &x)
+#define EXP_SETMUL(d, x, y) fib_init(&d); fib_mul(&d, &x, &y)
+
+#include "exp.h"
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- @mp_fibonacci@ --- *
+ *
+ * Arguments: @long n@ = index desired (may be negative)
+ *
+ * Returns: The %$n$%th Fibonacci number.
+ */
+
+mp *mp_fibonacci(long n)
+{
+ struct fib d, g;
+ mp *nn;
+
+ d.n = 0; d.a = MP_ZERO; d.b = MP_ONE;
+ if (n >= 0) { g.n = 1; g.a = MP_ONE; g.b = MP_ZERO; }
+ else { g.n = -1; g.a = MP_ONE; g.b = MP_MONE; n = -n; }
+ nn = mp_fromlong(MP_NEW, n);
+
+ EXP_WINDOW(d, g, nn);
+
+ MP_DROP(nn); fib_drop(&g); MP_DROP(d.b);
+ return (d.a);
+}
+
+/*----- Test rig ----------------------------------------------------------*/
+
+#ifdef TEST_RIG
+
+#include <mLib/testrig.h>
+
+static int vfib(dstr *v)
+{
+ long x = *(long *)v[0].buf;
+ mp *fx = *(mp **)v[1].buf;
+ mp *y = mp_fibonacci(x);
+ int ok = 1;
+ if (!MP_EQ(fx, y)) {
+ fprintf(stderr, "fibonacci failed\n");
+ MP_FPRINTF(stderr, (stderr, "fib(%ld) = ", x), fx);
+ MP_EPRINT("result", y);
+ ok = 0;
+ }
+ mp_drop(fx);
+ mp_drop(y);
+ assert(mparena_count(MPARENA_GLOBAL) == 0);
+ return (ok);
+}
+
+static test_chunk tests[] = {
+ { "fibonacci", vfib, { &type_long, &type_mp, 0 } },
+ { 0, 0, { 0 } }
+};
+
+int main(int argc, char *argv[])
+{
+ test_run(argc, argv, tests, SRCDIR "/tests/mp");
+ return (0);
+}
+
+#endif
+
+/*----- That's all, folks -------------------------------------------------*/
projects / u / mdw / catacomb / commitdiff