| 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 -------------------------------------------------*/ |