math/mpreduce.h: Missing include files.

[u/mdw/catacomb] / math / mp-fibonacci.c

/* -*-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 "/t/mp");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/