chiark - git - mdw - catacomb/blob - math/mp-fibonacci.c

   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 {
 265   test_run(argc, argv, tests, SRCDIR "/t/mp");
 266   return (0);
 267 }
 268
 269 #endif
 270
 271 /*----- That's all, folks -------------------------------------------------*/