+++ /dev/null
-/* -*-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 / blobdiff